## 2 Materials and Methods  
### Image Reconstruction Algorithm  
Two steps were required in this study. The first step was to perform DOT. In DOT, both the absorption and scattering coefficient images were recovered from an iterative finite element (FEM)-based algorithm, as previously described in detail by [PERSON] et al. [37] (pages 9-15). The photon diffusion equation which is coupled with Robin boundary conditions,  
\[\
abla\cdot D(r)\
abla\Phi(r)-\mu_{a}(r)\Phi(r)=-S(r) \tag{1}\]  
\[-D(r)\
abla\Phi_{n}=\alpha\Phi, \tag{2}\]  
is discretized and solved by using finite element method, where \(\Phi(r)\) is the photon density, \(S(r)\) is the source term, \(\mu_{a}(r)\) is the absorption coefficient, and \(D\) is the diffusion coefficient. The diffusion coefficient can be written as \(D(r)=1/3\Big{[}\mu_{a}(r)+\mu_{s}^{{}^{\prime}}(r)\Big{]}\), where \(\mu_{s}^{{}^{\prime}}(r)\) is the reduced scattering coefficient and \(\alpha\) is the boundary condition (BC) coefficient related to the internal reflection at the boundary. Then, a regularized [PERSON]'s method is used to update an initial optical property distribution iteratively in order to minimize a weighted sum of the squared difference between computed and measured optical data along its boundary.  
The core procedure in our reconstruction algorithms is to iteratively solve the following regularized matrix equation:  
\[\Big{(}J^{T}J+\lambda I\Big{)}\Delta q=J^{T}\Big{(}\Phi^{(m)}-\Phi^{(c)}\Big{)} \tag{3}\]  
where \(J\) is the Jacobian matrix that is formed by \(\partial\Phi/\partial D\) and \(\partial\Phi/\partial\mu_{a}\) at the boundary measurement sites, \(I\) is the identity matrix, and \(\lambda\) is the regularization parameter determined by combined Marquardt and Tikhonov regularization schemes. \(\Delta q=\left(\Delta D_{1},\Delta D_{2},\ldots,\Delta D_{N},\Delta\mu_{a,1},\Delta\mu_{a,2},\ldots,\Delta\mu_{a,N}\right)^{T}\) is the update vector for the optical property profiles, where \(N\) is the total number of nodes in the finite element mesh. \(\Phi^{(m)}\) and \(\Phi^{(c)}\) are the measured and computed photon density, respectively. A three-dimensional cylindrical FE mesh of 700 nodes was used for the image reconstruction. The reconstruction algorithm converged within nine iterations in a GPU (graphics processing unit) parallelized computer with 8 GB RAM (random-access memory).  
The second step is to establish the relationship of the reduced scattering coefficient obtained by DOT and the tissue temperature measured by a thermocouple. Since the first change in optical properties was detected at 45 \({}^{\circ}\)C [24], only the data when the temperature was higher than 45 \({}^{\circ}\)C was utilized in this study. In this step, the DOT data from five tissue samples were fit to a two-degree polynomial (\(y=ax^{2}+bx+c\)) or two-term exponential model (\(y=ae^{bx}+ce^{dx}\)), in which a, b, c and d are the coefficients. Due to the non-linear relationship between the temperature and \(\mu_{s}^{{}^{\prime}}\) as observed in the experimental data, the entire heating process was divided into two steps in the fitting with the polynomial model (1 st step: 45 to 54 \({}^{\circ}\)C, 2 nd step: 54 to 59 \({}^{\circ}\)C), whereas only one step was needed in the exponential model (45 to 59 \({}^{\circ}\)C). Once this relationship was established, the temperature of a tissue could be estimated using this relationship. To verify the performance of each method, the temperatures predicted by different fitting models were compared with those measured by thermometer using other five separate tissue samples. Additionally, the root mean square error (RMSE) between the measured and predicted temperature was calculated to evaluate the precision of each fitting model.