C omprehensive analytical model for CW laser induced heat in turbid media

: 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. 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


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][2][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][10][11][12][13][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][16][17][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][20][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][25][26]. Nevertheless, many researchers proposed analytical solutions to the diffusion equation [27][28][29][30][31][32][33][34][35][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][37][38][39][40][41]. These analytical methods are usually based on the Greens' function and series expansion methods [30,[41][42][43]. Meanwhile, the spatiotemporal temperature distribution in the medium can be described by the Pennes' bioheat equation. Like the diffusion equation, the bio-heat equation has been solved analytically, numerically and using Monte Carlo methods [44][45][46][47][48][49][50][51][52][53][54][55][56][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][59][60][61][62][63][64][65][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 laserinduced 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].

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] −∇D∇Φ(r r r) + μ a Φ(r r r) = S(r r r) (1) where Φ, D, μ a , and S are the photon density, the diffusion coefficient, the absorption coefficient and the light source, respectively. The diffusion coefficient is where μ s 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 where ρ, c, k represent the density, specific heat and thermal conductivity for tissue, respectively. Here, E(r) = μ a φ 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 for a point light source term modeled by the Dirac delta function δ (r r r, r r r ) where r r r and γ are the position and the strength of the light source, respectively. The solution of the homogeneous diffusion equation in 2D cylindrical polar coordinates is where β = i μ a D , J m and Y m are the first and second kind Bessel functions, respectively and i is the complex number. Here, A m and B m are the constants resulting from the radial solutions while C m and D m are the constants resulting from the angular ones. If the source is placed in the x axis (θ = 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 i = 0). The schematic of the problem can be seen in Fig. 1.
where 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 i , the Bessel function of first kind 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 for r ≤ r i . On the other hand, for R ≥ r ≥ r i , the photon density consists of the two radial solutions since it has a finite value in this region. Therefore, the photon density becomes for R ≥ r ≥ r i . We have three unknown differential constants 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][73][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 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 Here, 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, δ (r, r i ) = δ (r−r i )δ (θ −θ i ) r and the orthogonality relation for cosine functions, 2π 0 cos(mθ ) cos(nθ )dθ = πδ 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 for the sudden increase of the derivative of the photon density. The Robin boundary condition and the photon densities for each sub-domain given by Eqs. (9) and (12) can be simplified further to give and Solving Eqs. and respectively for r i = 0. Therefore, the photon density reaches its final form with the constants 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 s are the heat transfer coefficient and the surrounding temperature, respectively. For calculational simplicity, we define a shifted temperatureT ,T = T − T s . First, we consider the homogeneous case for the bio-heat thermal diffusion equation Now, we use the separation of variables method and assume that 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 Here λ 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 temperatureT leading to There are infinite number of solutions for this equation. In other words, the solutions are the roots of Eq. (24), which are λ = λ l where l = 0, 1, 2.... Therefore, the solution for the homogeneous part of the bio-heat equation is where T p, l is constant. In order to find T p, l , we use the following initial condition at t = 0 T (r, θ , 0) = T s (26) orT (r, θ , 0) = 0.
As a result, both T p, l and the solution to the homogeneous part of the bio-heat equation become zero. Now we consider the non-homogeneous bio-heat equation and define the non-homogeneous part in terms of the solutions of the homogenous part Writing Eq. (28) into the Pennes' bio-heat thermal diffusion equation given by Eq. (2) and using the radial and the angular equations yields the temporal solution In our case, the source term is time independent so that ω can be taken out of the integral. Thus, Ω(t) takes the following form To find ω, we write the solution of the diffusion equation Φ(r, θ ) into Eq. (28) where Λ m (β r) is the radial solution of the tissue diffusion equation, which is Multiplying Eq. (31) with rJ p (λ l r) cos(p θ ) and then integrating between the intervals r ∈ [0, R] and θ ∈ [0, 2π] gives Utilization of the orthogonality relation of cosine functions drops the summation in Eq. (33).
Calculating the other integrals in Eq. (33) yields where 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

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.

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.

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 (μ s ) and the refractive index (n) are set to 0.0132 mm −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/μ s on the x axis (θ = 0). The distance 1/μ s 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.
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 meas-ured temperature map. Nevertheless, the measured temperature maps show a good agreement with both analytical and FEM maps.

Experimental
Analytical FEM T Fig. 6. The temperature maps for the chicken breast tissue obtained experimentally, analytically, and numerically (FEM), respectively. The simulation results are obtained 42 seconds after the beginning of the laser heating while the experimental data is acquired between 36 s and 48 s. The cross-section of the chicken breast sample is approximated by a 28 mm diameter circle indicated by the green dash-line.
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.

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.