Comprehensive analytical model for CW laser induced heat in turbid media

Hakan Erkol; Farouk Nouizi; Alex Luk; Mehmet Burcin Unlu; Gultekin Gulsen

doi:10.1364/OE.23.031069

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² 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]

- \nabla D \nabla Φ (r) + μ_{a} Φ (r) = S (r)

where Φ, D, μ_a, and S are the photon density, the diffusion coefficient, the absorption coefficient and the light source, respectively. The diffusion coefficient is

D (r) = \frac{1}{3 [μ_{a} (r) + μ_{s}^{'} (r)]}

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

ρ c \frac{\partial T (r, t)}{\partial t} = \nabla [k \nabla T (r, t)] + E (r)

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

- \nabla^{2} Φ (r) + \frac{μ_{a}}{D} Φ (r) = \frac{γ}{D} δ (r, r^{'})

for a point light source term modeled by the Dirac delta function δ(r, r′) where 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

Φ (r, θ) = \sum_{m = - \infty}^{\infty} [A_{m} (β) J_{m} (β r) + B_{m} (β) Y_{m} (β r)] [C_{m} (β) cos (m θ) + D_{m} (β) sin (m θ)]

where

β = i \sqrt{\frac{μ_{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.

Fig. 1 The schematic representing the geometry of a homogeneous medium with a δ function source located at r = r_i.

Download Full Size | PDF

Thus, the constant D_m(β) in Eq. (4) is zero and hence only cos(mθ) term contributes to the solution leading to

Φ (r, θ) = \sum_{m = - \infty}^{\infty} [a_{m} (β) J_{m} (β r) + b_{m} (β) Y_{m} (β r)] cos (m θ)

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

Φ_{<} (r, θ) = \sum_{m = - \infty}^{\infty} a_{m} (β) J_{m} (β r) cos (m θ)

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

Φ_{>} (r, θ) = \sum_{m = - \infty}^{\infty} [b_{m} (β) J_{m} (β r) + c_{m} (β) Y_{m} (β r)] cos (m θ)

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:

Φ (R, θ) + 2 ξ D {\frac{\partial Φ (r, θ)}{\partial r} |}_{r = R} = 0

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

Φ_{<} {(r, θ) |}_{r = r_{i}} = Φ_{>} {(r, θ) |}_{r = r_{i}} .

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

\sum_{m = - \infty}^{\infty} {\frac{1}{r} \frac{d}{d r} [r \frac{d}{d r} G_{m} (β r)] cos (m θ) + \frac{1}{r^{2}} \frac{d^{2}}{d θ^{2}} [cos (m θ)] G_{m} (β r) - \frac{μ_{a}}{D} G_{m} (β r) cos (m θ)} = - \frac{γ}{D} δ (r, r_{i})

where 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

\begin{array}{r} \sum_{m = \infty}^{\infty} \int_{r_{i} - ε}^{r_{i} + ε} [\frac{d}{d r} (r \frac{d G_{m} (β r)}{d r})] d r \int_{0}^{2 π} cos (m θ) cos (n θ) d θ \\ = - \frac{1}{D} γ \int_{r_{i} - ε}^{r_{i} + ε} \int_{0}^{2 π} δ (r, r_{i}) cos (n θ) r d r d θ . \end{array}

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}) = \frac{δ (r - r_{i}) δ (θ - θ_{i})}{r}

and the orthogonality relation for cosine functions,

\int_{0}^{2 π} 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

{\frac{d G_{m} (β r)}{d r} |}_{r = r_{i} + ε} - {\frac{d G_{m} (β r)}{d r} |}_{r = r_{i} - ε} = - \frac{γ}{π D r_{i}} cos (m θ_{i})

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

b_{m} (β) J_{m} (β R) + c_{m} (β) Y_{m} (β R) = - 2 ξ D {[b_{m} (β) \frac{d J_{m} (β r)}{d r} + c_{m} (β) \frac{d Y_{m} (β r)}{d r}] |}_{r = R},

b_{m} (β) J_{m} (β r_{i}) + c_{m} (β) Y_{m} (β r_{i}) = a_{m} (β) J_{m} (β r_{i})

and

{[b_{m} (β) \frac{d J_{m} (β r)}{d r} + c_{m} (β) \frac{d Y_{m} (β r)}{d r} - a_{m} (β) \frac{d J_{m} (β r)}{d r}] |}_{r = r_{i}} = - \frac{γ}{π D r_{i}} cos (m θ_{i}) .

Solving Eqs. (13)–(15) for the differentiation constants a_m(β), b_m(β) and c_m(β) yields

\begin{array}{l} a_{m} (β) = & \frac{γ cos (m θ_{i})}{2 D [2 D β ξ R J_{m - 1} (β R) + (R - 2 D m ξ) J_{m} (β R)]} \\ \times {J_{m} (β r_{i}) [2 D β ξ R Y_{m - 1} (β R) + (R - 2 D m ξ) Y_{m} (β R)] \\ + Y_{m} (β r_{i}) [- 2 D β ξ R J_{m - 1} (β R) - (R - 2 D m ξ) J_{m} (β R)]}, \end{array}

\begin{array}{l} b_{m} (β) = & γ & R & \frac{D β ξ [Y_{m - 1} (β R) - Y_{m + 1} (β R)] + Y_{m} (β R)}{2 D [2 D β ξ R J_{m - 1} (β R) + (R - 2 D m ξ) J_{m} (β R)]} \\ \times J_{m} (β r_{i}) cos (m θ_{i}) \end{array}

and

c_{m} (β) = - \frac{γ}{2 D} J_{m} (β r_{i}) cos (m θ_{i}),

respectively for r_i ≠ 0. Therefore, the photon density reaches its final form with the constants a_m(β), b_m(β) and c_m(β)

Φ (r, θ) = \sum_{m = - \infty}^{\infty} {\begin{matrix} a_{m} (β) J_{m} (β r) cos (m θ) & if & r \leq r_{i} \\ (b_{m} (β) J_{m} (β r) + c_{m} (β) Y_{m} (β r)) cos (m θ) & if & r \geq r_{i} \end{matrix} .

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

- k \frac{\partial T (r)}{\partial t} = h (T_{s} - T (r))

where h and T_s are the heat transfer coefficient and the surrounding temperature, respectively. For calculational simplicity, we define a shifted temperature T̃, T̃ = T − T_s. First, we consider the homogeneous case for the bio-heat thermal diffusion equation

ρ c \frac{\partial \tilde{T} (r, t)}{\partial t} = \nabla [k \nabla \tilde{T} (r, t)] .

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

\tilde{T} (t, r, θ) = {exp (\frac{- k}{ρ c} λ^{2} t), J_{p} (λ r), cos (p θ)} .

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 temperature T̃

- k \frac{\partial \tilde{T} (R)}{\partial t} = h \tilde{T} (R)

leading to

{J_{p} (λ r) |}_{r = R} - \frac{k}{h} {\frac{d J_{p} (λ r)}{d r} |}_{r = R} = 0 .

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

\tilde{T} (r, θ, t) = \sum_{p = - \infty}^{\infty} \sum_{l = 0}^{\infty} T_{p, l} exp (\frac{- k}{ρ c} λ_{l}^{2} t) J_{p} (λ_{l} r) cos (p θ)

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}

or

\tilde{T} (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

\frac{μ_{a} Φ (r, θ)}{ρ c} = \sum_{p = - \infty}^{\infty} \sum_{l = 0}^{\infty} ω (t) J_{p} (λ_{l} r) cos (p θ) .

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

Ω (t) = exp (\frac{- k}{ρ c} λ_{l}^{2} t) \int_{0}^{t} exp (\frac{k}{ρ c} λ_{l}^{2} t^{'}) ω (t^{'}) d t^{'} .

In our case, the source term is time independent so that ω can be taken out of the integral. Thus, Ω(t) takes the following form

Ω (t) = \frac{ρ c}{k λ_{l}^{2}} ω [1 - exp (\frac{- k}{ρ c} λ_{l}^{2} t)] .

To find ω, we write the solution of the diffusion equation Φ(r, θ) into Eq. (28)

\frac{μ_{a}}{ρ c} \sum_{m = - \infty}^{\infty} Λ_{m} (β r) cos (m θ) = \sum_{p = - \infty}^{\infty} \sum_{l = 0}^{\infty} ω J_{p} (λ_{l} r) cos (p θ)

where Λ_m(βr) is the radial solution of the tissue diffusion equation, which is

Λ_{m} (β r) = {\begin{matrix} a_{m} (β) J_{m} (β r) & if & r \leq r_{i} \\ b_{m} (β) J_{m} (β r) + c_{m} (β) Y_{m} (β r) & if & R \geq r \geq r_{i} \end{matrix} .

Multiplying Eq. (31) with rJ_p′ (λ_l′r) cos(p′θ) and then integrating between the intervals r ∈ [0, R] and θ ∈ [0, 2π] gives

ω_{m, l} = \frac{μ_{a}}{ρ c} \frac{\sum_{m = - \infty}^{\infty} \int_{0}^{R} Λ_{m} (β r) J_{m} (λ_{l} r) r d r \int_{0}^{2 π} cos (m θ) cos (p θ) d θ}{\int_{0}^{R} J_{p}^{2} (λ_{l} r) r d r \int_{0}^{2 π} {cos}^{2} (p θ) d θ} .

Utilization of the orthogonality relation of cosine functions drops the summation in Eq. (33). Calculating the other integrals in Eq. (33) yields

\begin{array}{l} ω_{m, l} = & \frac{μ_{a}}{ρ c} {- a_{m} (β) r_{i} [β J_{m - 1} (β r_{i}) J_{m} (r_{i} λ_{l}) - λ_{l} J_{m} (β r_{i}) J_{m - 1} (r_{i} λ_{l})] \\ + b_{m} (β) r_{i} [β J_{m - 1} (β r_{i}) J_{m} (r_{i} λ_{l}) - λ_{l} J_{m} (β r_{i}) J_{m - 1} (r_{i} λ_{l})] \\ - β R J_{m - 1} (R β) J_{m} (R λ_{l}) + R λ_{l} J_{m} (R β) J_{m - 1} (R λ_{l}) \\ + c_{m} (β) r_{i} [β Y_{m - 1} (β r_{i}) J_{m} (r_{i} λ_{l}) - λ_{l} Y_{m} (β r_{i}) J_{m - 1} (r_{i} λ_{l})] \\ - β R Y_{m - 1} (R β) J_{m} (R λ_{l}) + R λ_{l} Y_{m} (R β) J_{m - 1} (R λ_{l})} \\ \times \frac{1}{\frac{1}{2} R (β^{2} - λ_{l}^{2}) {R [J_{m - 1} {(R λ_{l})}^{2} + J_{m} {(R λ_{l})}^{2}] - \frac{2 m J_{m - 1} (R λ_{l}) J_{m} (R λ_{l})}{λ_{l}}}} \end{array}

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

T (r, θ, t) = T_{s} + \sum_{m = - \infty}^{\infty} \sum_{l = 0}^{\infty} \frac{ρ c}{k λ_{l}^{2}} ω_{m, l} J_{m} (λ_{l} r) cos (m θ) [1 - exp (\frac{- k}{ρ c} λ_{l}^{2} t)] .

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.

Fig. 2 The schematic of the interface used in the experiment. The laser beam is collimated using a GRIN lens and delivered to the homogeneous phantom from a hole prepared at the side of the custom built RF coil. The total power is distributed to three point light sources and simulation results for individual sources are superimposed.

Download Full Size | PDF

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.

Table 1. Agar phantom parameters used in the experiment.

View Table

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 (μ′_s) and the refractive index (n) are set to 0.0132 mm⁻¹, 0.8 mm⁻¹ 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³ and 4200 J (kg °C)⁻¹, 5.8 10⁻⁴ W (mm °C)⁻¹ and 8 10⁻⁶ W (mm² °C)⁻¹, 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).

Fig. 3 The experimental, analytical, and FEM based temperature maps for the phantom. The first and second rows show the linear, T(°C), and logarithmic, ln(T), scaled temperature maps, respectively. For the simulated temperature, the light sources are positioned at 1/μ′_s mm under the illumination site as imposed by the diffusion approximation. The simulation results are obtained 18 seconds after the beginning of the laser heating while the experimental data is acquired between 12 s and 24 s. Some artifacts are observable in the experimental data particularly towards the edges.

Download Full Size | PDF

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.

Fig. 4 The temperature profiles (a) linear, (b) logarithmic carried out along the x-axis on the temperature maps presented in Fig. 3. For the simulated temperature, the light sources are positioned at 1/μ′_s mm under the illumination site as imposed by the diffusion approximation. The simulation results are obtained 18 seconds after the beginning of the laser heating while the experimental data is acquired between 12 s and 24 s. All three methods are in good agreement above the MRT noise level, 0.1 °C (green dotted line).

Download Full Size | PDF

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.

Fig. 5 The temporal temperature profiles obtained between 0 s and 400 s at (a) x = 5 mm and (b) x = 12 mm. The temperature near the boundary (x = 12 mm) reaches a plateau much quicker due to heat convection.

Download Full Size | PDF

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.

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.

Download Full Size | PDF

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.

Fig. 7 The temperature profiles for chicken breast tissue (a) linear, (b) logarithmic carried out along the x-axis on the temperature maps presented in Fig. 6. 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. All three methods are in good agreement above the MRT noise level, 0.15 °C (green dotted line).

Download Full Size | PDF

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₂ 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]

Comprehensive analytical model for CW laser induced heat in turbid media

Abstract

1. Introduction

2. Methods

2.1. Analytical solution to the combined diffusion and heat equation

2.2. MRT spatiotemporal temperature measurements

2.3. Numerical solution to the combined diffusion and heat equation

3. Results and discussions

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Tables (1)

Equations (35)

Optics Express

Abbreviation	Parameter	Value	Unit
c	specific heat	4200	J (kg °C)⁻¹
k	thermal conductivity	5.8 10⁻⁴	W (mm °C)⁻¹
h	heat transfer coefficient	8.0 10⁻⁶	W (mm² °C)⁻¹
T_s	surrounding temperature	14	°C
ρ	density	1000	kg/m³