## Abstract

In this work, we present a new analytical approach to model continuous wave laser induced temperature in highly homogeneous turbid media. First, the diffusion equation is used to model light transport and a comprehensive solution is derived analytically by obtaining a special Greens’ function. Next, the time-dependent bio-heat equation is used to describe the induced heat increase and propagation within the medium. The bio-heat equation is solved analytically utilizing the separation of variables technique. Our theoretical model is successfully validated using numerical simulations and experimental studies with agarose phantoms and ex-vivo chicken breast samples. The encouraging results show that our method can be implemented as a simulation tool to determine important laser parameters that govern the magnitude of temperature rise within homogenous biological tissue or organs.

© 2015 Optical Society of America

## 1. Introduction

The use of lasers has become considerably attractive for a variety of basic biological research and biomedical applications such as hypothermia laser therapy, photodynamic therapy (PDT) and thermal ablation therapy. Hypothermia laser therapy or low-level laser therapy (LLLT) has been used for the treatment of various diseases and injuries from muscle and brain injuries to cardiac disorders [1–3]. LLLT has a great potential to provide highly conformal thermal therapy without any harmful impact to the patient. Visible or near infrared (NIR) lasers (600–1064 nm) at the intensity level of about 0.001–5 W/cm^{2} are typically used for LLLT with an irradiation time of a few seconds or longer [4, 5]. Unfortunately, the mechanism of laser-tissue interaction during LLLT is not well understood and also the optimal treatment parameters have not been systematically established yet. Therefore, modeling of laser-tissue interaction for LLLT can play an important role in determining these parameters. Indeed, some recent work revealed the importance of the modeling of laser-tissue interaction since it showed that the skin temperature may reach as high as 44 °C during LLLT therapy [6]. When the tissue temperature rises beyond 42 °C, conformational changes of molecules take place and they are accompanied by bond destruction and membrane alterations [4]. At that rate, it is possible that a significant thermal damage may occur in biological tissues during LLLT therapy. Local tissue temperature is also one of the most important factors affecting the efficacy of PDT, which uses stationary laser [7, 8]. In this technique, a stationary laser beam is used to excite the tissue. The efficacy of the PDT has been shown with clinical trials and emerging data support the benefits of PDT for several applications involving management of nonmelanoma skin cancer and port wine stain birthmarks by combining PDT with pulsed dye laser therapy [9–14]. Since laser illumination can be utilized to warm up the tissue, modeling of laser-tissue interaction would be important for these type of applications.

Another important technique using the modeling of light due to the laser illumination is laser induced thermal ablation therapy [15–18]. This treatment method is a minimally invasive and alternative technique to the conventional surgery and can provide faster recovery times as well as less complication rates [19–21]. The increased temperature due to light absorption also promotes the tumor destruction. Two important factors should be considered in this technique, namely the magnitude and the duration of temperature increase in tissue. Generally, these two factors are monitored in real-time, e.g. using real-time magnetic resonance thermometry (MRT) monitoring. Nevertheless, simulation studies are performed prior to the laser irradiation in order to adjust the important laser parameters such as wavelength, power, and the irradiation time that govern the magnitude of temperature rise in tissue. Performing these simulation studies consists in the estimation of spatiotemporal distribution of heat deposited by the laser in the tissue. These simulations are performed by modeling the light transport within the tissue and the subsequent time-dependent temperature increase and propagation.

Biological tissue is usually considered a turbid medium characterized with given optical absorption and scattering. Light propagation is this type of media can be modeled using the radiative transfer equation [22, 23]. However, the diffusion approximation has been the preferred method since it is substantially easier to solve. The diffusion equation is generally solved using numerical methods such as the finite element method (FEM) [24–26]. Nevertheless, many researchers proposed analytical solutions to the diffusion equation [27–36]. Although structural inhomogeneities in tissue and its irregular geometry make it complex to obtain analytical solutions, these analytical methods can still be used on some regular geometries such as cylindrical, spherical and infinite or semi-infinite slabs [24, 36–41]. These analytical methods are usually based on the Greens’ function and series expansion methods [30, 41–43]. Meanwhile, the spatiotemporal temperature distribution in the medium can be described by the Pennes’ bio-heat equation. Like the diffusion equation, the bio-heat equation has been solved analytically, numerically and using Monte Carlo methods [44–57]. Technically, the diffusion and bio-heat equations are combined in order to model the laser induced temperature. To be able to achieve this combination, the source term of the bio-heat equation is expressed as a function of the solution of the diffusion equation and local absorption of the medium [58–66]. In these previous works, the combined diffusion and bio-heat equation system is solved analytically by using simple photon density expressions to describe the source term of the bio-heat equation.

In this study, we propose a new approach that not only provides a very comprehensive photon density expression but presents a very detailed analytical temperature expression for laser-induced temperature as well. For the first part, a detailed analytical solution is obtained for the diffusion equation based on an integral approach using the Robin boundary condition, which is more adequate for tissue modeling. For the second part, the heterogeneous heat equation is solved using the separation of variables technique considering that the temporal and spatial parts of the heat equation can be separable. Here, the heat convection boundary condition is used to derive the solution. In addition, our approach allows the implementation of any boundary condition making it a very flexible method. A simplified result for the zero (Dirichlet) boundary condition can also be obtained by just setting the refractive index mismatched constant in the Robin boundary condition to zero.

The results obtained by our method are validated using both numerical simulations as well as experimental phantom and ex-vivo tissue sample studies. In fact, the time-dependent temperature increase in the medium obtained with our method is in very good agreement with the FEM solutions and the experimental MRT measurements. Therefore, we believe that our approach will provide a suitable fast simulation platform for laser induced temperature rise and propagation in homogeneous organs such as liver [67, 68].

## 2. Methods

#### 2.1. Analytical solution to the combined diffusion and heat equation

As mentioned above, estimating laser induced temperature requires solution of a multiphysics problem, namely a combined diffusion and bio-heat equation system. The first one is the light propagation in tissue and the second one is the heat transfer. The photon transport in tissue can be described by the continuous wave diffusion equation [30, 69, 70]

where Φ, D,*μ*, and

_{a}*S*are the photon density, the diffusion coefficient, the absorption coefficient and the light source, respectively. The diffusion coefficient is $D(\mathbf{r})=\frac{1}{3\left[{\mu}_{a}(\mathbf{r})+{\mu}_{s}^{\prime}(\mathbf{r})\right]}$ where

*μ′*is the reduced scattering coefficient. On the other hand, the temperature distribution in tissue can be described by the Pennes’ bio-heat thermal diffusion equation. In this paper, the validation of the approach is carried out using agar phantoms and ex-vivo tissue samples. For this reason, the blood perfusion is considered null so that the Pennes’ bio-heat thermal equation reduces to

_{s}*ρ*,

*c*,

*k*represent the density, specific heat and thermal conductivity for tissue, respectively. Here,

*E*(

*r*) =

*μ*is the thermal energy absorbed from the laser heating [71]. Now, we first solve the diffusion equation for a spatially invariant diffusion coefficient so that the diffusion equation becomes

_{a}ϕ*δ*(

**,**

*r***) where**

*r′***and**

*r′**γ*are the position and the strength of the light source, respectively.

The solution of the homogeneous diffusion equation in 2D cylindrical polar coordinates is

*J*and

_{m}*Y*are the first and second kind Bessel functions, respectively and

_{m}*i*is the complex number. Here,

*A*and

_{m}*B*are the constants resulting from the radial solutions while

_{m}*C*and

_{m}*D*are the constants resulting from the angular ones. If the source is placed in the x axis (

_{m}*θ*= 0), the photon density is symmetric with respect to this axis. It is important to note that our approach is valid provided that the source position is different from the origin (

*r*≠ 0). The schematic of the problem can be seen in Fig. 1.

_{i}Thus, the constant *D _{m}*(

*β*) in Eq. (4) is zero and hence only cos(

*mθ*) term contributes to the solution leading to

*a*(

_{m}*β*) and

*b*(

_{m}*β*) are defined as

*a*(

_{m}*β*) =

*A*(

_{m}*β*).

*C*(

_{m}*β*) and

*b*(

_{m}*β*) =

*B*(

_{m}*β*).

*C*(

_{m}*β*) for calculational simplicity, respectively. Next, we solve the non-homogeneous diffusion equation with a source term represented by the Dirac

*δ*function. We consider the tissue domain as two sub-domains according to the source position. For

*r*≤

*r*, the Bessel function of first kind

_{i}*J*(

_{m}*βr*) only contributes to the solution since the Bessel function of the second kind tends to infinity when

*r*→ 0. Therefore, the solution takes the following form

*r*≤

*r*. On the other hand, for

_{i}*R*≥

*r*≥

*r*, the photon density consists of the two radial solutions since it has a finite value in this region. Therefore, the photon density becomes

_{i}*R*≥

*r*≥

*r*. We have three unknown differential constants

_{i}*a*(

_{m}*β*),

*b*(

_{m}*β*) and

*c*(

_{m}*β*). To find these constants, we introduce three expressions. The first one is the utilization of the Robin boundary condition at the tissue boundary: where

*ξ*is the refractive index mismatched constant [72–74]. The other two expressions are obtained using the mathematical properties of the Dirac delta source. Firstly, the photon densities given by Eqs. (6) and (7) are equal at the source position The second expression can be obtained based on the fact that the derivative of the photon density has an abrupt change at the source position. Considering this sudden change, the diffusion equation for the whole domain becomes

*G*(

_{m}*βr*) = {

*J*(

_{m}*βr*),

*Y*(

_{m}*βr*)} is the radial solutions. By multiplying Eq. (10) with

*r*cos(

*nθ*) and integrating over

*r*and

*θ*between the intervals (

*r*=

*r*−

_{i}*ε*,

*r*=

*r*+

_{i}*ε*) and (

*θ*= 0,

*θ*= 2

*π*) as

*ε*→ 0 leads to

*n*, which is an integer between −∞ and ∞, is introduced for another set of cosine functions (or angular solutions, cos(

*nθ*)) to drop the summation in the left hand side of Eq. (11). In Eq. (10), only the first term on the left hand side survives and the others become zero since the photon densities for the sub-domains are equal to each other at the source position. We use the expression of the Dirac delta function in 2D polar coordinates, $\delta (\mathbf{r},{\mathbf{r}}_{i})=\frac{\delta (r-{r}_{i})\delta (\theta -{\theta}_{i})}{r}$ and the orthogonality relation for cosine functions, ${\int}_{0}^{2\pi}\text{cos}(m\theta )\text{cos}(n\theta )d\theta =\pi {\delta}_{m,n}$ if both

*m*and

*n*are not equal to zero (the integral is simply 2

*π*if

*m*= 0 and

*n*= 0). Substituting these relations into Eq. (11) and using the equality of the photon density yields

*a*(

_{m}*β*),

*b*(

_{m}*β*) and

*c*(

_{m}*β*) yields

*r*≠ 0. Therefore, the photon density reaches its final form with the constants

_{i}*a*(

_{m}*β*),

*b*(

_{m}*β*) and

*c*(

_{m}*β*)

Now, we introduce the following heat convection boundary condition to solve the Pennes’ bio-heat thermal diffusion equation

where*h*and

*T*are the heat transfer coefficient and the surrounding temperature, respectively. For calculational simplicity, we define a shifted temperature

_{s}*T̃*,

*T̃*=

*T*−

*T*. First, we consider the homogeneous case for the bio-heat thermal diffusion equation

_{s}*k*is homogeneous. In our case, the solution is symmetric with respect to the x axis since the light source is in the x axis and the phantom is a disk shaped with a radius

*R*. For this reason, the solutions

*Y*(

_{n}*λr*) and sin(

*pθ*) are excluded. Therefore, the solution becomes

*λ*is a constant (to be determined from the boundary condition) and

*p*is an integer ranging from −∞ to ∞. In order to determine the constant

*λ*, we use the heat convention boundary condition given by Eq. (20) for the shifted temperature

*T̃*leading to There are infinite number of solutions for this equation. In other words, the solutions are the roots of Eq. (24), which are

*λ*=

*λ*where

_{l}*l*= 0, 1, 2....

Therefore, the solution for the homogeneous part of the bio-heat equation is

*T*is constant. In order to find

_{p,l}*T*, we use the following initial condition at

_{p,l}*t*= 0 or As a result, both

*T*and the solution to the homogeneous part of the bio-heat equation become zero.

_{p,l}Now we consider the non-homogeneous bio-heat equation and define the non-homogeneous part in terms of the solutions of the homogenous part

*ω*can be taken out of the integral. Thus, Ω(

*t*) takes the following form

*ω*, we write the solution of the diffusion equation Φ(

*r*,

*θ*) into Eq. (28)

*(*

_{m}*βr*) is the radial solution of the tissue diffusion equation, which is

*rJ*(

_{p′}*λ*) cos(

_{l′}r*p′θ*) and then integrating between the intervals

*r*∈ [0,

*R*] and

*θ*∈ [0, 2

*π*] gives

*a*(

_{m}*β*),

*b*(

_{m}*β*) and

*c*(

_{m}*β*) have been presented in the previous section. Therefore, the temperature distribution can be calculated by the following solution

#### 2.2. MRT spatiotemporal temperature measurements

To validate our approach experimentally, we first perform experiments with homogeneous phantoms mimicking biological tissue. These phantoms are mainly prepared using agarose (OmniPur Agarose, Lawrence, KS). Intralipid and Indian ink are added in order to set the optical properties of the phantom. Secondly, the performance of our method is validated using chicken breast tissue sample. Although chicken breast tissue is considered homogeneous, it contains thick fibers which makes it moderately heterogenous. The spatiotemporal temperature distribution inside the media under investigation is measured using magnetic resonance thermometry (MRT), while it is illuminated with laser light as shown in Fig. 2.

The laser diode emitting at 808 nm and its driver are both placed in the MRI control room. Light is transferred to the MRI bore using 15 meter long 400 *μ*m diameter optical fiber. After being collimated by a GRIN lens, the laser beam (1.8 mm diameter) is delivered to the phantom from a hole prepared at the side of the custom built RF coil, Fig. 2. The power of the laser is adjusted to 4.5 W. The light output of the fiber is collimated and the spot size of the beam illuminating the surface of the phantom is 1.8 mm. When using diffusion approximation, an isotropic point source located a transport mean free path (1/*μ′ _{s}*) beneath the surface is generally utilized to model a point light source at the surface of the phantom. Since our solution is valid for a point source, a finite number of point sources can be utilized to model the collimated light beam of 1.8 mm diameter. Our simulation studies showed that three point sources located a transport mean free path beneath the surface are enough for accurate modeling of this collimated beam with a small spot size. Hence, the solution of the combined equation is calculated for each source separately and summed to get the effect of the collimated beam. The optical and thermal properties of the homogeneous agar phantom are listed in Table 1.

The MRT experiments are conducted inside a Philips 3T Achieva system and data acquisition is synchronized with the laser driver. MRT measurements can be performed using any phase sensitive technique such as Proton Resonance Frequency (PRF). In this work, PRF images are obtained with a gradient echo sequence using 60 and 12 milliseconds as repetition (TR) and echo time (TE), respectively. First, a T1 weighted low resolution MR pilot image is obtained to determine the axial position of the laser probe using fiducial markers. After the axial plane is located, a baseline high resolution phase image *ϕ* is obtained using a gradient echo scan. Next, the phantom is illuminated by the laser beam for 24 seconds while two consecutive MRT images are acquired, each lasts 12 seconds. The difference between the corresponding phase maps for these two time points and the baseline phase image provide high resolution temperature images. Please note that these difference images reveal the relative temperature increase induced by only the laser illumination. Our extensive tests with phantoms and tissue samples confirmed that the temperature induced by MRT in the measurements is negligible.

#### 2.3. Numerical solution to the combined diffusion and heat equation

Our solution is also compared with the results obtained using our FEM solver [45]. In these simulations, numerical phantoms having the same geometry as the experimental phantoms are utilized. The phantoms are discretized using meshes consisting of 7268 triangular elements connected at 3758 nodes. Similar to analytical solutions, the spatiotemporal temperature distribution is obtained in two steps. First, the photon distribution is obtained by solving the diffusion equation. Next, this photon density is used to express the source term of the heat equation [45]. It is worth noticing that the same boundary conditions as the ones used in the analytical method are applied.

## 3. Results and discussions

To validate our analytical approach, we calculate the temperature map analytically and compare it first with the map obtained using FEM then with an experimentally measured temperature map. The analytical temperature distribution is obtained using Eq. (35). The optical properties, the absorbtion coefficient (*μ _{a}*), the reduced scattering coefficient (

*μ′*) and the refractive index (n) are set to 0.0132 mm

_{s}^{−1}, 0.8 mm

^{−1}and 1.4, respectively. For thermal properties, the density (

*ρ*), the specific heat of the phantom (

*c*), the thermal conductivity (

*k*) and the heat transfer coefficient (

*h*) are chosen as 1000 kg/m

^{3}and 4200 J (kg °C)

^{−1}, 5.8 10

^{−4}W (mm °C)

^{−1}and 8 10

^{−6}W (mm

^{2}°C)

^{−1}, respectively [45]. For both analytical and numerical simulations, the source of light is placed at

*r*=

_{i}*R*− 1/

*μ′*on the x axis (

_{s}*θ*= 0). The distance 1/

*μ′*represents the reduced scattering pathlength imposed by the approximation of the diffusion equation while R represents the radius. Moreover, the power level in the simulations is adjusted to match experimental results.

_{s}Technically, the temperature provided by MRT at each time point is the average temperature value over the duration of the data acquisition time which is equal to 12 s. Unlike MRT, the synthetic temperature distributions are provided at a specific time point. During this validation, the experimental data that correspond to the second time point (12 s – 24 s) are used. Therefore, the measured MRT temperature is compared with the synthetic temperatures obtained at *t* = 18 s, the midpoint for the MRT acquisition. The temperature maps obtained by the three methods are very similar. Figure 3 shows the cross-sectional linear (top) and logarithmic (bottom) temperature distributions corresponding to the optical fiber plane obtained experimentally, analytically, and numerically (FEM).

Solving the combined diffusion and heat equation system using FEM requires preparation of a mesh and iteratively inverting the system matrix describing the problem, which are both very time consuming. This disadvantage can be overcome using analytical methods. In fact, using our analytical approach, the temperature map is obtained in approximately three seconds. Using FEM, the computation time necessary to generate a temperature map utilizing a CPU parallelized solver is slightly higher than 60 s. This represents more than a 95% reduction in the computation time.

Figure 3 shows the experimental, analytical, and FEM based temperature maps for the 25 mm diameter agarose phantom (*R* = 12.5 mm). As specified in the experimental method, the temperature map provided by MRT does not represent the absolute temperature but the relative increase in temperature induced by the laser, Fig. 3. Thus, the room temperature, 14 °C, is subtracted from the absolute synthetic temperature maps in order to compare with the experimental one, Fig. 3. One can notice some slight artifacts at the lower boundary of the phantom in the measured temperature map, Fig. 3. Nevertheless, this measured temperature map is in a very good agreement with the simulated temperature maps. As expected, the highest increase in temperature can be observed near the source of light positioned at *θ* = 0 and *r _{i}* =

*R*− 1/

*μ′*.

_{s}Figures 4(a) and 4(b) show the linear and logarithmic scaled temperature profiles carried out along x axis, respectively. As seen from the profiles, temperature reaches up to 4 °C at the light source position, 1.5 mm below the boundary, 18 s after turning on the laser. The temperature decreases with the depth and goes down to sub-Celsius level approximately 5 mm below the boundary. The logarithmic profiles show that both our method and FEM successively match the experimental ones when temperature is above the MRT noise level (0.1 °C), Figs. 4(a) and 4(b). Accordingly, the MRT measurements are dependable up to 9.8 mm (x = 2.7 mm) deep under the illumination site.

Even though the heating of the phantom is realized using a CW laser, the heat propagation within the medium is simulated using the time-dependent bio-heat equation. Therefore, the time dependance of our method needs to be validated for. Figures 5(a) and 5(b) show the analytical and FEM simulated temperatures obtained between 0 s and 400 s with a temporal resolution of 10 s for different positions. Moreover, in order to validate the effect of our boundary conditions on the time dependent solution, the temporal profiles are performed close to (*x* = 12 mm) and slightly further away (*x* = 5 mm) from the source.

As you can see in Figs. 5(a) and 5(b), our temporal temperature profiles show a very good match with those obtained with FEM. As expected, the profiles obtained near the source at *x* = 12 mm show a higher increase in temperature compared to the ones performed at *x* = 5 mm. For the point far away from the boundary (*x* = 5 mm), the temperature increases almost linearly. However, it reaches to a plateau for the point closer to the boundary (*x* = 12 mm) after 200 s. For both points, the maximum error in the temporal profiles is less than 1.5% which validates the temporal aspect of our method. This validation shows that our analytical method is not only fast but as accurate as the widely used FEM method.

To validate our method on real biological tissue, we also measure spatiotemporal temperature maps on a chicken breast sample 40 seconds after turning on the laser. Again, the measured temperature map is compared to the ones computed both analytically and using FEM. To do so, the cross-section of the chicken breast sample is approximated by a 28 mm diameter circle. During this experiment, the measured data that correspond to the fourth time point (36 s – 48 s) is used. Therefore, the measured MRT temperature is compared with the synthetic temperature maps obtained at *t* = 42 s, the midpoint for the MRT acquisition. Figure 6 shows the experimental, analytical, and FEM temperature maps. Some artifacts are visible in the measured temperature map. Nevertheless, the measured temperature maps show a good agreement with both analytical and FEM maps.

A better comparison of these results can be made on the linear and logarithmic temperature profiles presented in Figs. 7(a) and 7(b). After heating the chicken breast sample for 40 seconds, the analytical and FEM results almost overlap. In fact, the maximum error between the two simulated temperature maps is approximately 0.05 °C. When compared to the measured temperature, the profiles obtained by both simulation methods show a good match. Here, the maximum error observed between the simulated profiles and the experimental one is slightly less than 0.17 °C (x = 3.7 mm) where the MRT signal is below the noise level. Overall, our method provides very encouraging results compared to not only the widely used FEM but also the experimental MRT results.

## 4. Conclusion

In this study, we present an analytical approach for CW laser induced temperature in homogeneous turbid media. Similar to all analytical solution proposed in the literature, the homogeneity is assumed in all optical and thermal properties within the medium. We first solve the diffusion equation analytically with an integral method by obtaining a special Greens’ function for the Robin boundary condition. Secondly, we obtain a solution for the Pennes’ bio-heat equation by utilizing the separation of variables technique for the heat convection boundary condition. This approximation is made based on the fact that the spatial and temporal parts of the heat equation can be separable. Moreover, our approach allows implementation of any boundary condition. Although we use the Robin boundary condition in this work, a simplified result can also be obtained by just setting refractive index mismatched constant to zero, which corresponds to zero boundary condition. Also, it should be noted that our method is suitable for cases with low or negligible blood perfusion since it is not accounted for.

Our analytical method is first validated with experimental studies. The experimental spatiotemporal temperature distribution within the tissue simulating phantom and chicken breast sample is obtained using magnetic resonance thermometry (MRT). In addition, we also validated our approach using numerical results (FEM). The comparison shows that our analytical approach is not only fast but also yields accurate results as well. The encouraging results show that our new analytic approach can be used for the fast determination of key laser parameters for any application in which a CW laser is utilized to heat a homogeneous turbid medium. These applications may include low level laser therapy, thermally modulated PDT and CW laser induced thermal ablation therapy. Since it has been recently shown that the optical properties of the tissue depend on temperature, it is important to take this phenomena into account during the modeling of laser-tissue interaction [6]. Our method can be iteratively used to estimate a more accurate temperature distribution by updating the optical properties of the tissue periodically at each time point. Our next aim will be measuring the temperature dependence of the optical properties of tissue and incorporating this into our model. Furthermore, we will work on expanding our technique for pulsed laser beams since they are extensively used in laser induced thermal ablation therapy. Finally, we are planning to add blood convention to our analytic model and validate it in vivo.

## Acknowledgments

This research is supported by National Institutes of Health (NIH) grants R01EB008716, F31CA171745, R21CA191389, R21EB013387 and P30CA062203, TUBITAK Grant 2219, Bogazici University Research funding Grant No. BAP 7126 and TUBITAK Grant No. 112T253.

## References and links

**1. **L. Assis, A. I. Soares-Moretti, T. B. Abrahão, H. P. de Souza, M. R. Hamblin, and N. A. Parizotto, “Low-level laser therapy (808 nm) contributes to muscle regeneration and prevents fibrosis in rat tibialis anterior muscle after cryolesion,” Lasers Med. Sci. **28**, 947–955 (2013). [CrossRef]

**2. **S. Eiichi, K. Hiroyasu, F. Hirosuke, F. Motoki, K. Tadashi, O. Yasutaka, Y. Susumu, T. Ryosuke, M. Tsuyoshi, and S. Michiyasu, “Diverse effects of hypothermia therapy in patients with severe traumatic brain injury based on the computed tomography classification of the traumatic coma data bank,” J. Neurotrauma **32**(5), 353–361 (2015). [CrossRef]

**3. **G. Debaty, M. Maignan, S. Ruckly, and J-F Timsit, “Impact of intra-arrest therapeutic hypothermia on outcome of prehospital cardiac arrest: response to comment by Saigal and Sharma,” Intensive Care Med. **41**(1), 172–173 (2015). [CrossRef]

**4. **M. H. Niemz, *Laser-Tissue Interactions: Fundamentals and Applications* (Biol. Med. Phys. Biomed., 2007).

**5. **Y. Y. Huang, “Biphasic dose response in low level light therapy an update,” Dose-Response **9**(4), 602–618 (2011). [CrossRef]

**6. **S. Kim and S. Jeong, “Effects of temperature-dependent optical properties on the fluence rate and temperature of biological tissue during low-level laser therapy,” Lasers Med. Sci. **9**(4), 637–644 (2014). [CrossRef]

**7. **A. Y. Seteikin, I. V. Krasnikov, E. Drakaki, and M. Makropoulou, “Dynamic model of thermal reaction of biological tissues to laser-induced fluorescence and photodynamic therapy,” J. Biomed. Opt. **18**(7), 075002 (2013). [CrossRef] [PubMed]

**8. **V. V. Barun and A. P. Ivanov, “Temperature regime of biological tissue under photodynamic therapy,” Biofizika **57**(1), 120–129 (2012). [PubMed]

**9. **L. R. Braathen, C. A. Morton, N. Basset-Seguin, R. Bissonnette, M. J. Gerritsen, Y. Gilaberte, P. Calzavara-Pinton, A. Sidoroff, H. C. Wulf, and R. M. Szeimies, “Photodynamic therapy for skin field cancerization: an international consensus. International society for photodynamic therapy in dermatology,” J. Eur. Acad. Dermatol. Venereol. **26**(9), 1063–1069 (2012). [CrossRef] [PubMed]

**10. **B. Zhao and Y. Y. He, “Recent advances in the prevention and treatment of skin cancer using photodynamic therapy,” Expert. Rev. Anticancer Ther. **10**(11), 1797–1809 (2010). [CrossRef] [PubMed]

**11. **W. J. Moy, S. J. Patel, B. S. Lertsakdadet, R. P. Arora, K. M. Nielsen, K. M. Kelly, and B. Choi, “Preclinical in vivo evaluation of NPe6-mediated photodynamic therapy on normal vasculature,” Lasers Surg. Med. **44**(2), 158–162 (2012). [CrossRef] [PubMed]

**12. **K. M. Kelly, W. J. Moy, A. J. Moy, B. S. Lertsakdadet, J. J. Moy, E. Nguyen, A. Nguyen, K. E. Osann, and B. Choi, “Talaporfin sodium-mediated photodynamic therapy alone and in combination with pulsed dye laser on cutaneous vasculature,” J. Invest. Dermatol. **135**(1), 302–306 (2015). [CrossRef]

**13. **J. Yang, A. C. Chen, Q. Wu, S. Jiang, X. Liu, L. Xiong, and Y. Xia, “The influence of temperature on 5–aminolevulinic acid-based photodynamic reaction in keratinocytes in vitro,” Photodermatol. Photoimmunol. Photomed. **26**(2), 83–91. (2010). [CrossRef] [PubMed]

**14. **A. Willey, R. R. Anderson, and F. H. Sakamoto, “Temperature-modulated photodynamic therapy for the treatment of actinic keratosis on the extremities: a pilot study,” Dermatol. Surg. **40**(10), 1094–1102 (2014). [CrossRef] [PubMed]

**15. **G. L. LeCarpentier, M. Motamedi, L. P. McMath, S. Rastegar, and A. J. Welch, “Continuous wave laser ablation of tissue: analysis of thermal and mechanical events,” IEEE Trans. Biomed. Eng. **40**(2), 188–200 (1993). [CrossRef] [PubMed]

**16. **G. Wendt-Nordahl, S. Huckele, P. Honeck, P. Alken, T. Knoll, M. S. Michel, and A. Häcker, “Systematic evaluation of a recently introduced 2-μm continuou-wave thulium laser for vaporesection of the prostate,” J. Endourol. **22**(5), 1041–1045 (2008). [CrossRef] [PubMed]

**17. **Y. Xu, D. Sun, Z. Wei, B. Hong, and Y. Yang, “Clinical study on the application of a 2-μm continuous wave laser in transurethral vaporesection of the prostate,” Exp. Ther. Med. **5**, 1097–1100 (2013). [PubMed]

**18. **T. Bach, N. Huck, F. Wezel, A. Häcker, A. J. Gross, and M. S. Michel, “70 vs 120 W thulium:yttrium-aluminium-garnet 2-μm continuous-wave laser for the treatment of benign prostatic hyperplasia: a systematic ex-vivo evaluation,” BJU Int. **106**, 368–372 (2009). [CrossRef]

**19. **E. Bay, A. Douplik, and D. Razansk, “Optoacoustic monitoring of cutting efficiency and thermal damage during laser ablation,” Lasers. Med. Sci. **29**, 1029–1035 (2014). [CrossRef]

**20. **E. Posadzka, R. Jach, K. Pitynski, and M. J. Jablonski, “Treatment efficacy for pain complaints in women with endometriosis of the lesser pelvis after laparoscopic electroablation vs. CO_{2} laser ablation,” Med. Phys. **30**(1), 147–152 (2015).

**21. **S. Sinha, E. Hargreaves, N. V. Patel, and S. F. Danish, “Assessment of irrigation dynamics in magnetic-resonance guided laser induced thermal therapy (MRgLITT),” Lasers Surg. Med. **47**, 273–280 (2015). [CrossRef] [PubMed]

**22. **C. Darne, Y. Lu, and E. M. Sevick-Muraca, “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. **59**(1), R1–R64 (2014). [CrossRef]

**23. **A. Liemert and A. Kienle, “Novel analytical solution for the radiance in an anisotropically scattering medium,” Appl. Opt. **54**(8), 1963–1969 (2015). [CrossRef] [PubMed]

**24. **S. R. Arridge, “A finite element approach for modelling photon transport in tissue,” Med. Phys. **20**, 299–309 (1993). [CrossRef] [PubMed]

**25. **H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. **25**, 711–743 (2008). [CrossRef] [PubMed]

**26. **M. B. Unlu, O. Birgul, R. Shafiha, G. Gulsen, and O. Nalcioğlu, “Diffuse optical tomographic reconstruction using multifrequency data,” J. Biomed. Opt. **11**, 054008 (2006). [CrossRef]

**27. **M. S. Patterson and S. J. Madsen, “Diffusion equation representation of photon migration in tissue,” IEEE MTT-S Digest **2**, 905–908 (1991).

**28. **M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the finite element method for the forward and inverse model in optical tomography,” J. Mat. Imaging Vis. **3**, 263–283 (1993). [CrossRef]

**29. **J. M. Schmitt, G. X. Zhou, E. C. Walker, and R. T. Wall, “Multilayer model of photon diffusion in skin,” J. Opt. Soc. Am. A **7**, 2141–2153 (1990). [CrossRef] [PubMed]

**30. **S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. **37**, 1531–1560 (1992). [CrossRef] [PubMed]

**31. **T. J. Farrel, M. S. Patterson, and B. C. Wilson, “A diffusion theory model of specially resolved steady state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. **19**, 879–888 (1992). [CrossRef]

**32. **D. A. Boas, M. A. O’leary, B. Chances, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications,” Proc. Natl. Acad. Sci. **91**, 4887–4891 (1994). [CrossRef] [PubMed]

**33. **S. A. Walker, D. A. Boas, and E. Gratton, “Photon density waves scattered from cylindrical inhomogeneities: theory and experiments,” Appl. Opt. **37**, 1935–1944 (1998). [CrossRef]

**34. **B. W. Pogue and M. S. Patterson, “Frequency-domain optical absorption spectroscopy of finite tissue volumes using diffusion theory,” Phys. Med. Biol. **39**, 1157–1180 (1994). [CrossRef] [PubMed]

**35. **D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation I. Theory,” Appl. Opt. **36**, 4587–4599 (1997). [CrossRef] [PubMed]

**36. **A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A. **14**, 246–254 (1997). [CrossRef]

**37. **F. Martelli, A. Sassaroli, Y. Yamada, and G. Zaccanti, “Analytical approximate solutions of the timedomain diffusion equation in layered slabs,” J. Opt. Soc. Am. A. **19**, 71–80 (2002). [CrossRef]

**38. **A. Kienle, “Light diffusion through a turbid parallelepiped,” J. Opt. Soc. Am. A **22**, 1883–1888 (2005). [CrossRef]

**39. **F. Martelli, A. Sassaroli, S. D. Bianco, and G. Zaccanti, “Solution of the time-dependent diffusion equation for a three-layer medium: application to study photon migration through a simplified adult head model,” Phys. Med. Biol. **52**, 2827–2843 (2007). [CrossRef] [PubMed]

**40. **A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. I. Homogeneous Case,” Opt. Express **18**, 9456–9473 (2010). [CrossRef] [PubMed]

**41. **A. Zhang, D. Piao, C. F. Bunting, and B. W. Pogue, “Photon diffusion in a homogeneous medium bounded externally or internally by an infinetely long circular cylindrical applicator. I. Steady-state theory,” J. Opt. Soc. Am. A **27**, 648–662 (2010). [CrossRef]

**42. **J. Sikora, A. Zacharopoulos, A. Douiri, M. Schweiger, L. Horesh, S. R. Arridge, and J. Ripoll, “Diffuse photon propagation in multilayered geometries,” Phys. Med. Biol. **51**, 497–516 (2006). [CrossRef] [PubMed]

**43. **A. Liemert and A. Kienle, “Light diffusion in N-layered turbid media: frequency and time domains,” J Biomed. Opt. **15**, 025002 (2010). [CrossRef] [PubMed]

**44. **H. H. Pennes, “Analysis of tissue and arterial blood temperatures in resting forearm,” J. Appl. Psychol. **1**, 93–122 (1948).

**45. **Y. Lin, Ha. Gao, D. Thayer, A. L. Luk, and G. Gulsen, “Photo-magnetic imaging: resolving optical contrast at MRI resolution,” Phys. Med. Biol. **58**, 3551–3562 (2013). [CrossRef] [PubMed]

**46. **H. Arkin, L. X. Xu, and K. R. Holmes, “Recent developments in modeling heat transfer in blood perfused tissues,” Phys. Med. Biol. **41**(2), 97–107 (1984).

**47. **G. Brix, M. Seebass, G. Hellwiga, and J. Griebela, “Estimation of heat transfer and temperature rise in partial-body regions during MR procedures: an analytical approach with respect to safety considerations,” Magn. Reson. Imaging **20**(1), 65–76 (2002). [CrossRef] [PubMed]

**48. **Z. S. Deng and J. Liu, “Mathematical modeling of temperature mapping over skin surface and its implementation in thermal disease diagnostics,” Comput. Biol. Med. **34**(6), 495–521 (2004). [CrossRef] [PubMed]

**49. **L. Jiang, W. Zhan, and M. H. Loew, “Modeling static and dynamic thermography of the human breast under elastic deformation,” Phys. Med. Biol. **56**, 187–202 (2011). [CrossRef]

**50. **J. Liu, X. Chen, and L. X. Xu, “New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating,” IEEE Trans. Biomed. Eng. **46**(4), 420–428 (1999). [CrossRef] [PubMed]

**51. **E. Neufeld, N. Chavannes, T. Samaras, and N. Kuster, “Novel conformal technique to reduce staircasing artifacts at material boundaries for FDTD modeling of the bioheat equation,” IEEE Trans. Biomed. Eng. **52**(4), 4371–4381 (2007).

**52. **M. N. Ozisik, *Heat Conduction* (Wiley-Interscience, 1993).

**53. **G. Wang, H. Shen, W. Cong, S. Zhao, and G. Wei, “Temperature-modulated bioluminescence tomography,” Opt. Express **14**(17), 7852–7871 (2006). [CrossRef] [PubMed]

**54. **J. J. Zhao, J. Zhang, N. Kang, and F. Yang, “A two level finite difference scheme for one dimensional Pennes bioheat equation,” Appl. Math. Comput. **171**(1), 320–331 (2005). [CrossRef]

**55. **S. A. Mirnezami, M. R. Jafarabadi, and M. Abrishami, “Temperature distribution simulation of the human eye exposed to laser radiation,” J. Lasers Med. Sci. **4**(4), 175–181 (2013). [PubMed]

**56. **S. Sarkar, C. Fan, J. C. Hsiang, and R. M. Dickson, “Modulated fluorophore signal recovery buried within tissue mimicking phantoms,” J. Phys. Chem. A **117**, 9501–9509 (2013). [CrossRef] [PubMed]

**57. **M. E. Kowalski and J. M. Jin, “Determination of electromagnetic phased-array driving signals for hyperthermia based on a steady-state temperature criterion,” IEEE Trans. Microw. Theory Techn. **48**(11), 1864–1873 (2000). [CrossRef]

**58. **Y. Feng and D. Fuentes, “Model-based planning and real-time predictive control for laserinduced thermal therapy,” Int. J. Hyperthermia **27**(8), 751–761 (2011). [CrossRef] [PubMed]

**59. **A. M. Elliott, J. Schwartz, J. Wang, A. M. Shetty, C. Bourgoyne, D. P. ONeal, J. D. Hazle, and R. J. Stafford, “Quantitative comparison of delta P1 versus optical diffusion approximations for modeling near-infrared gold nanoshell heating,” Med. Phys. **36**, 1351–1358 (2009). [CrossRef] [PubMed]

**60. **S. K. Cheong, S. Krishnan, and S. H. Cho, “Modeling of plasmonic heating from individual gold nanoshells for near-infrared laser-induced thermal therapy,” Med. Phys. **36**, 4664–4671 (2009). [CrossRef] [PubMed]

**61. **D. R. Wyman and W. M. Whelan, “Basic optothermal diffusion theory for interstitial laser photocoagulation,” Med. Phys. **21**, 1651–1656 (1994). [CrossRef] [PubMed]

**62. **R. K. Banerjee, L. Zhu, P. Gopalakrishnan, and M. J. Kazmierczak, “Influence of laser parameters on selective retinal treatment using single-phase heat transfer analyses,” Med. Phys. **34**(5), 1828–1841 (2007). [CrossRef] [PubMed]

**63. **M. Jaunicha, S. Rajea, K. Kimb, K. Mitraa, and Z. Guob, “Bio-heat transfer analysis during short pulse laser irradiation of tissues,” Int. J. Heat Mass Transfer **51**, 5511–5521 (2008). [CrossRef]

**64. **R. Zhang, W. Verkruysse, G. Aguilar, and J. S. Nelson, “Comparison of diffusion approximation and Monte Carlo based finite element models for simulating thermal responses to laser irradiation in discrete vessels,” Phys. Med. Biol. **50**(5), 4075–4086 (2005). [CrossRef] [PubMed]

**65. **M. N. Rylander, Y. Feng, J. Bass, and K. R. Diller, “Heat shock protein expression and injury optimization for laser therapy design,” Laser Surg. Med. **39**, 731–746 (2007). [CrossRef]

**66. **G. Shafirstein, W. Bäumler, M. Lapidoth, S. Ferguson, P. E. North, and M. Waner, “A new mathematical approach to the diffusion approximation theory for selective photothermolysis modeling and its implication in laser treatment of port-wine stains,” Lasers Surg. Med. **34**, 335–347 (2004). [CrossRef] [PubMed]

**67. **B. S. Hijmansa, A. Grefhorstb, M. H. Oosterveera, and A. K. Groena, “Zonation of glucose and fatty acid metabolism in the liver: mechanism and metabolic consequences,” Biochimie **96**, 121–129 (2014). [CrossRef]

**68. **A. Miranda, A. P. P López-Cardona, R. Laguna-Barraza, A. Calle, I. López-Vidriero, B. Pintado, and A. Gutiérrez-Adán, “Transcriptome profiling of liver of non-genetic low birth weight and long term health consequences,” BMC Genomics **15**(327), 1–12 (2014). [CrossRef]

**69. **L. H. V. Wang and H. Wu, *Biomedical Optics: Principles and Imaging* (Wiley, 2007).

**70. **F. Nouizi, M. Torregrossa, R. Chabrier, and P. Poulet, “Improvement of absorption and scattering discrimination by selection of sensitive points on temporal profile in diffuse optical tomography,” Opt. Express **19**(13), 12843–12854 (2011). [CrossRef] [PubMed]

**71. **S. H. Diaz, G. Aguilar, E. J. Lavernia, and B. J. F. Wong, “Modeling the thermal response of porcine cartilage to laser irradiation,” IEEE J. Sel. Top. Quantum Electron. **7**, 944–951 (2001). [CrossRef]

**72. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems **15**, 41–93 (1999). [CrossRef]

**73. **H. Erkol and M. B. Unlu, “Virtual source method for diffuse optical imaging,” Appl. Opt. **52**, 4933–4970 (2013). [CrossRef] [PubMed]

**74. **H. Erkol, A. Demirkiran, N. Uluc, and M. B. Unlu, “Analytical reconstruction of the bioluminescent source with priors,” Opt. Express **22**(16), 19758–19773 (2014). [CrossRef] [PubMed]