## 1.

## Introduction

Laser is an effective tool in clinical applications. A particular use of laser energy is the selective photothermal interaction in cancer treatment. It has been shown that the combination of an
$805\text{-}\mathrm{nm}$
laser and *in situ* indocyanine green (ICG), a dye with an absorption peak around
$800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, can result in selective heating when treating tumors buried in normal tissue.^{1, 2, 3} Recent studies using nanoparticles as *in situ* dye also showed the selective photothermal effect.^{4} Since tumor tissue is more sensitive to temperature elevation than normal tissue,^{5, 6} the control of temperature in the target region during laser irradiation is crucial to the outcomes of the treatment.

Selective photothermal interaction can destroy tumor cells directly by raising the temperature of the target tissue. It can also serve as a precursor for immune activation by exposing tumor antigens, when combined with an immunotherapy-oriented treatment.^{7, 8} Tumor tissue responds to thermal treatment according to the treatment conditions. High-level thermal irradiation with extra-high temperature elevation can lead to total tissue destruction with denatured cells, which do not have much effect on the host immune stimulations. Low-level thermal irradiation may not destroy enough tumor cells and/or may not expose enough tumor antigens. The optical outcome of laser photothermal interaction is to kill as much target tumor cells as possible, and at the same time preserve the tumor proteins for the immunological stimulations. Such outcome depends on the temperature distribution in the tumor tissue.

The desired thermal effect strongly depends on the laser parameters, and the tissue optical and thermal properties. In the case of dye-enhanced laser therapy, the dye administration also becomes crucial. The concentrations of the absorbing dye in, and the optical and thermal properties of, a specific target region are important factors. Improper dye administration causes undesired tissue damages or insufficient heating of target tissue. To choose the treatment parameters properly and to predict the outcome of the photothermal effect, a reliable tissue model and simulation method is needed. Since biological tissue is a turbid medium and the mean free path of the light involved in irradiation is far less than the dimension of the medium, light propagation in tissue quickly becomes a random process. The ideal mathematical tool for such a process is the Monte Carlo method.^{9, 10, 11} The Monte Carlo method can effectively simulate absorption and scattering of the light in tissue, hence providing the spatial distribution of photon absorption. Heat transfer in the tissue during laser irradiation also affects the temperature distribution and tissue responses in the target tissue and surrounding tissue. Due to the relatively large physical dimensions and inhomogeneous characteristics of biological tissues, the diffusion process also needs to be modeled during laser irradiation. Using the absorbed energy in tissue as the heat source, the thermal energy diffusion process can be determined, hence providing information on thermal responses of tissue. In this study, Monte Carlo simulation was used to acquire the information of laser absorption in tissue, and the finite difference method for thermal diffusion in the tissue was used to determine the temperature distribution during and after laser irradiation. The dynamic processes of tissue thermal responses to laser irradiation under different tissue configurations and laser parameters are reported.

## 2.

## Monte Carlo Method for Laser-Tissue Interaction

Biological tissue is considered as a linear isotropic homogeneous medium. When a photon enters the tissue, the Monte Carlo method samples random propagation variables from a well-defined probability distribution. Selection of the step size of the photon traveling inside the tissue between two interaction sites is based on the photon’s mean free path. For a specific type of tissue with an absorption coefficient
${\mu}_{a}$
and scattering coefficient
${\mu}_{s}$
, the step size is given as^{9, 10, 11}

When a photon moves from one site to the next, it will have a deflection from its original direction. We define the deflection with an angle $\theta $ and an azimuthal angle $\phi $ . The probability distribution of $\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta $ is described by the scattering function:

## 2

$$p\left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta \right)=\frac{1-{g}^{2}}{2{(1+{g}^{2}-2g\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta )}^{3\u22152}},$$^{12}The anisotropy factor $g$ is the average value of $\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta $ .

For each random step, the deflection angle is determined using the distribution in Eq. 2:

## 3

$$\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta =\{\begin{array}{cc}\frac{1}{2g}[1+{g}^{2}-{\left(\frac{1-{g}^{2}}{1-g+2g\zeta}\right)}^{2}\phantom{)}\hfill & \hfill \phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.3em}{0ex}}g>0\hfill \\ 2\zeta -1\hfill & \hfill \phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.3em}{0ex}}g<0\hfill \end{array}\phantom{)}.$$The azimuth angle is uniformly distributed throughout the interval from 0 to $2\pi $ .

An incoming photon may hit the layer boundary when it propagates. When a photon reaches an interface, the remaining photon step is changed to
${s}_{0}-{s}_{1}$
, where
${s}_{0}$
is the step size at the original layer and
${s}_{1}$
is defined as follows:^{9, 11}

## 4

$${s}_{1}=\{\begin{array}{cc}(z-{z}_{0})\u2215(-{\mu}_{z})\hfill & \hfill \phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.3em}{0ex}}{\mu}_{z}<0\hfill \\ ({z}_{1}-z)\u2215{\mu}_{z}\hfill & \hfill \phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.3em}{0ex}}{\mu}_{z}>0\hfill \end{array}\phantom{\}},$$If the photon bounces back, the photon would travel in the same layer with the same property; if the photon goes into another layer, the step size will be changed to ${s}^{\prime}={s}_{0}{\mu}_{t0}\u2215{\mu}_{t1}$ , where ${\mu}_{t1}$ and ${\mu}_{t0}$ are the interaction coefficients for the new layer and the original layer, respectively.

Snell’s law and Fresnel reflectance at the boundary are used to determine if the photon is reflected back or injected into the next layer.

The simulation process is as follows. A photon is launched and reaches the top of the tissue. It propagates a certain step size $s$ , with a rotation angle $\phi $ , and deflection angle $\theta $ , to an interaction site and deposits a portion of its energy to the site, determined by the tissue optical properties. The deposition of energy is determined by the weight $W$ that a photon leaves behind after an interaction at a specific site. The change of weight is defined as $\mathrm{\Delta}W=W{\mu}_{a}\u2215{\mu}_{t}$ . When this value is below a desired threshold, the photon is terminated. For this simulation, we terminated a photon if one ten thousandth of its original weight was left after its interactions with the surrounding media. After termination, a new photon is launched and the process is repeated.

Once the Monte Carlo simulation is complete for all the photons, an absorption power density matrix is generated for a given tissue configuration. This matrix represents the amount of laser power absorbed by the region in the form of power density. In a cylindrical system, the absorption power density matrix is a function of the coordinates $r$ and $z$ , with a cylindrical symmetry. Multiplying the absorption photon probability density by the desired laser power, the total power absorbed in the region is obtained. The absorbed power at each location can be used as the heat source for the thermal diffusion process.

## 3.

## Finite Difference Method

The diffusion equation with constant thermal properties and steady heat generation can be used for laser tissue interactions,

## 5

$${\nabla}^{2}T(r,z,t)+\frac{\stackrel{\u0307}{q}(r,z)}{k}=\frac{\rho c}{k}\frac{\partial T(r,z,t)}{\partial t},$$## 6

$$\stackrel{\u0307}{q}=\mathrm{photon}\phantom{\rule{0.3em}{0ex}}\mathrm{absorption}\phantom{\rule{0.3em}{0ex}}\mathrm{probability}*\mathrm{laser}\phantom{\rule{0.3em}{0ex}}\mathrm{power},$$From Eq. 6 one can find the temperature increase as a function of position and time. Using the Laplace equation in cylindrical coordinates and assuming azimuthal symmetry, Eq. 5 becomes

## 7

$$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)+\frac{{\partial}^{2}T}{\partial {z}^{2}}+\frac{q}{k}=\frac{\rho c}{k}\frac{\partial T}{\partial t}.$$^{13, 14, 15}Specifically, Eq. 7 can be rewritten as,

## 8

$$\frac{{T}_{i-1,j}^{n}-2{T}_{i,j}^{n}+{T}_{i+1,j}^{n}}{{\left(\mathrm{\Delta}r\right)}^{2}}+\frac{1}{i\mathrm{\Delta}r}\frac{{T}_{i+1,j}^{n}-{T}_{i-1,j}^{n}}{2\mathrm{\Delta}r}+\frac{{T}_{i,j-1}^{n}-2{T}_{i,j}^{n}+{T}_{i,j+1}^{n}}{{\left(\mathrm{\Delta}z\right)}^{2}}+\frac{{q}_{i,j}^{n}}{k}=\frac{\rho c}{k}\frac{{T}_{i,j}^{n+1}-{T}_{i,j}^{n}}{\mathrm{\Delta}t},$$^{13, 14, 15}For the solution to converge and be stable, the time step is limited towhere $\alpha =k\u2215\rho c$ . This value corresponds to a mesh Fourier number ${F}_{0m}\equiv \alpha \mathrm{\Delta}t\u2215(\mathrm{\Delta}r{)}^{2}\le 1\u22154$ . Values above 1/4 will cause the average temperature at location $(i,j)$ to overshoot the average of its nearest neighbors.

^{13, 14, 15}

The boundaries of the cylinder are subjected to convection, which follows a Neumann boundary condition. The convection condition at the top of the cylinder is given by,

where $h$ is the convection constant, and ${T}_{b}$ and ${T}_{\infty}$ are the temperatures of the boundary and the environment, respectively. Similarly, at the bottom of the cylinder the same condition applies, but there is a sign change on the derivative due to the orientation of the coordinate system. The outer surface of the cylinder is also subjected to convection, and the heat transfer can be treated similarly.The center of the cylinder must be treated with care, since the Laplacian is singular at the origin. Symmetry can be used to avoid the singularity. Specifically, as $r$ approaches zero following a radial path, the slope of the temperature field vanishes,

The approximated value of this boundary can be found byThis can be validated by the fact that at the origin, the temperature from $i=-1$ to $i=1$ can be approximated to be the same. Therefore, the slope of the temperature distribution across the center of the cylinder vanishes, allowing the second derivative to vanish also. This condition requires ${T}_{i+1,j}={T}_{i-1,j}$ . With boundary conditions, all the nodes can be obtained by iterating over all of the individual space steps dictated by the physical size of the region to be simulated. Also, one should notice that Eq. 8 can be used for all space steps, but at the boundaries the appropriate substitutions must be made as dictated by the boundary conditions.## 4.

## Simulations of Energy Distribution and Thermal Diffusion

The coordinate system of the simulation is a cylinder with user-defined radius and height. For simplicity, the radius and the height were chosen to be the same in our study. The cylinder configuration may consist of different layers with varying optical and thermal properties, or with an embedded spherical tumor, as shown in Fig. 1. Using the method described in Sec. 2, algorithms are constructed to determine the energy deposition using Monte Carlo simulation. Then, the finite difference method is used for the determination of heat and temperature distribution. Figure 2 shows the flow of the simulation process. The beam properties must be taken into consideration, since it is the source of the energy. A flat circular beam, which can be achieved experimentally by using a fiber fitted with a diffusing lens, is usually used in laser thermal treatment. In this simulation, the flat incident beam is constructed of multiple concentric rings. Each ring consists of an equally weighted number of photons to produce an isotropic incident beam. For a desired beam radius, the number of rings may be increased to approach a continuum of incident energy. A total of 100,000 photons was used in each of the simulations in this study.

Given in Tables
1, 2, 3, 4
are the tissue optical and thermal properties, laser beam parameters, and tissue configurations used in this simulation, based on experimentally measured values for different types of tissue,^{9, 11} and on desired target tissue absorption properties.

## Table 1

Tissue optical properties used in the simulation.

Symbol | Physical variable | Unit | Values used in this study |
---|---|---|---|

${\mu}_{a}$ | Absorption coefficient | $1\u2215\mathrm{cm}$ | 0.01 for normal tissue 1.0 for dye-enhanced tissue |

${\mu}_{s}$ | Scattering coefficient | $1\u2215\mathrm{cm}$ | 100 for normal tissue 100 for dye-enhanced tissue |

$n$ | Index of refraction | 0 for normal tissue 1 for dye-enhanced tissue | |

$g$ | Anisotropy | 0.9 for normal tissue 0.9 for dye-enhanced tissue |

## Table 2

Tissue thermal properties used in the simulation.

Symbol | Physical variable | Unit | Values used in this study |
---|---|---|---|

$\rho $ | Density | $\mathrm{g}\u2215{\mathrm{cm}}^{3}$ | 1.07 |

$k$ | Thermal conductivity | $\mathrm{W}\u2215\mathrm{cm}\xb0\mathrm{C}$ | 0.0056 |

$c$ | Specific heat | $\mathrm{J}\u2215\mathrm{g}\xb0\mathrm{C}$ | 3.4 |

$h$ | Heat transfer coefficient | $\mathrm{W}\u2215{\mathrm{cm}}^{2}\xb0\mathrm{C}$ | 0.000433 |

${T}_{0}$ | Initial layer temperature | $\xb0\mathrm{C}$ | 0.0 |

## Table 3

Laser beam parameters used in the simulation.

Symbol | Physical variable | Unit | Values used in this study |
---|---|---|---|

$P$ | Power | W | 1.7 to 5.0 |

$br$ | Beam radius | cm | 0.5 to 3.0 |

$nring$ | Number of rings within beam | 10 to 100 | |

$nphot$ | Number of photons to be propagated | ${10}^{5}$ | |

$bt$ | Beam type (flat or Gaussian) |

## Table 4

Geometrical tissue configurations used in the simulation.

Symbol | Physical variable | Unit | Values used in this study |
---|---|---|---|

$dr$ | Radial space increment | cm | 0.03 |

$dz$ | Axial depth increment | cm | 0.03 |

$nr$ | Number of radial increments | 100 | |

$nz$ | Number of axial depth increments | 100 | |

$na$ | Number of angular subdivisions | 30 | |

$nl$ | Number of tissue layers | 1–3 | |

$hi$ | Top of layer | cm | |

$hf$ | Depth of layer | (cm) | |

$d$ | Depth of spherical layer | (cm) | 1.0 or 3.0 |

$r$ | Radius of spherical layer | (cm) | 1.0 or 3.0 |

Using the finite difference method, the temperature increase as a function of space and time are obtained by a simple explicit iterative process. The temperature simulation follows the photon diffusion, as described in Fig. 2. The generator of the thermal energy is the product of the absorption power density matrix obtained by the Monte Carlo function and the desired laser power. This heat generation is taken to be constant throughout the irradiation process. The energy transport is governed by the thermal properties of the region as given in Table 2, which are accepted values for biological tissue.^{13}

## 5.

## Results

A cylindrical tissue of $1\text{-}\mathrm{cm}$ radius and $1\text{-}\mathrm{cm}$ height with a dye-enhanced spherical tumor of $0.3\text{-}\mathrm{cm}$ radius was used to mimic a buried tumor in normal tissue. This system was irradiated by a laser beam of $0.2\text{-}\mathrm{cm}$ radius. The photon absorption power density matrix is shown in Fig. 3 for a tumor of $0.2\text{-}{\mathrm{cm}}^{-1}$ absorption coefficient [Fig. 3a] and for a tumor of $0.6\text{-}{\mathrm{cm}}^{-1}$ absorption coefficient [Fig. 3b]. With the absorption enhancement, the buried tumor absorbed more photons than the surrounding normal tissue.

The geometrical configuration is important in photon absorption. In cancer treatment, different tissue layers with different optical properties may exist in and around the target region. Figure 4 depicts the photon absorption in a configuration of normal tissue (a cylinder of radius $1\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ and height $1\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ ) with a flat layer of dye-enhanced tissue ( $0.1\text{-}\mathrm{cm}$ thickness) $0.3\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ below the surface. The selective photon absorption can be observed in this configuration.

A larger cylindrical system, $3\text{-}\mathrm{cm}$ radius and $3\text{-}\mathrm{cm}$ height with a flat dye-enhanced layer of $0.2\text{-}\mathrm{cm}$ thickness placed in the center of the cylinder, was used for the temperature simulation. The laser beam was $3\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ in radius with $5\text{-}\mathrm{W}$ power and an irradiation duration of $4000\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ . The temperature changes inside the tissue at three different locations (surface of the cylinder, the center of the dye-enhanced layer, and the bottom of the cylinder) are shown in Fig. 5 . As shown in Fig. 5, the dye-enhanced layer had the most temperature increase. Figure 6 shows the temperature changes under the same conditions as those of Fig. 5, except for an irradiation duration of $2000\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ , followed by $2000\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ cooling time.

Different beam sizes were used to irradiate a cylindrical tissue system with a $3\text{-}\mathrm{cm}$ radius and $3\text{-}\mathrm{cm}$ height with a dye-enhanced tumor of $0.5\text{-}\mathrm{cm}$ radius at the center. The laser power was $1.7\phantom{\rule{0.3em}{0ex}}\mathrm{W}$ . The photon absorption density matrix, the temperatures at selected locations, and the temperature distributions as function of time were obtained for a laser beam of $0.5\text{-}\mathrm{cm}$ radius (Fig. 7 ) and of $1.5\text{-}\mathrm{cm}$ radius (Fig. 8 ). With dye enhancement, the tumor strongly absorbs the laser energy, as shown in Figs. 7a and 8a. When a smaller laser beam was used, the target temperature increase was much more significant [Fig. 7b] than that of a larger laser beam [Fig. 8b], although the thermal diffusion process for both cases followed a similar pattern as shown in Figs. 7c and 8c.

Figure 9 shows the simulation results of a dye-enhanced flat layer placed at the center of a cylinder. A $0.2\text{-}\mathrm{cm}$ flat layer was placed in the center with an absorption coefficient of $1.0\text{-}{\mathrm{cm}}^{-1}$ . A $3.0\text{-}\mathrm{cm}$ radius beam with power of $5.0\phantom{\rule{0.3em}{0ex}}W$ was used. The photon absorption power density is given in Fig. 9a, the temperatures at different locations are given in Fig. 9b, and the temperature distribution at different time frames is given in Fig. 9c. Although the laser power is high, the temperature increase is moderate due to the large laser beam.

## 6.

## Discussion

Using the Monte Carlo method and a finite difference algorithm, we studied the absorption of laser light by tissue of different optical properties, specifically for dye-enhanced target tissue buried in normal tissue. To theoretically predict the temperature increase, we designed an algorithm that calculated the temperature increase of a target tissue that was selectively heated. It takes into consideration the photon diffusion from a laser source using the Monte Carlo method and the thermal diffusion using the explicit finite difference method. The effectiveness of these methods can be easily coupled to model the entire photon/heat diffusion process.

The absorption power density matrices and the temperature distributions in different tissue configurations were obtained. The objective in laser photothermal therapy is to achieve the maximum light absorption in the target region and minimum damage to the surrounding tissue. We demonstrated the feasibility of the thermal selectivity using the simulation of a spherical target embedded within normal tissue. The absorption power density matrices shown in Figs. 3 and 4 illustrate the importance of the proper choice of absorption coefficient in selectively heating the target tissue.

The temperatures at different locations (Fig. 5) not only showed thermal selectivity in the dye-enhanced target tissue, but also showed the time course of thermal equilibrium during laser irradiation. Under the given conditions, the thermal equilibrium was reached $3000\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ after the start of the laser irradiation. The thermal relaxation process was shown in Fig. 6.

Figures 7, 8, 9 vividly illustrated the dynamics of the thermal responses of tissue systems containing dye-enhanced targets in different configurations. The impact of beam size is demonstrated in Figs. 7 and 8. Under the same laser power, a small, more concentrated laser beam apparently was more effective in raising the temperature of the target tissue (Fig. 8). An optimal combination of dye concentration, tissue configuration, and laser parameters can lead to a desired photothermal cancer treatment regimen.

When treating tumors *in vivo*, perfusion in the tumor and its surround structure will cause significant heat loss. Therefore, further studies to simulate the dye-enhanced photothermal interaction should include the heat loss term due to perfusion. In such a case, the thermal equilibrium shown in Figs. 5 and 6 will be different. However, with appropriate dye and irradiation doses, thermal equilibrium at a desired temperature can be achieved. Such a tradeoff can also be achieved by choosing the appropriate beam size. A smaller beam spot can result in a higher temperature increase, particularly in the central region of the target, as shown in Fig. 7. A wider beam with the same power can result in a much lower temperature increase in the target region, as shown in Fig. 8. This is due to the fact that the majority of the photons from the wide beam would miss the dye-enhanced target, hence reducing the overall absorption by the target. Therefore, final treatment protocol will be a function of target configuration, beam size, dye dose, and irradiation dose.

In the treatment of animal tumors, indocyanine green (ICG) has been used as the light absorbing dye for enhancement of selective photothermal interaction.^{1, 2, 3} The pharmacokinetic properties of ICG indicate a rapid hepatic excretion, and the differential uptake of ICG by tumor cells is not as significant as certain photosensitizers. Therefore, a systemic administration of ICG is not realistic for the proposed selective photothermal interaction for tumor destruction through prolonged laser irradiation. That is the reason why intratumoral administration of ICG has been adopted for tumor treatment and the desired selective photothermal effects have been observed.^{1, 2, 3}

The main goal of laser cancer treatment is to selectively destroy the cancer cells without denaturing the proteins. This allows the body’s immune system to recognize the tumor as “foreign” and to amount a systemic attack on the tumors. This goal can be achieved by a carefully designed treatment protocol with optimal operating parameters. The results obtained in this study demonstrate the feasibility of such a treatment regimen and the method developed in this study can be used effectively to guide treatment procedures.

## Acknowledgment

This work was supported in part by a grant from the College of Graduate Studies and Research of the University of Central Oklahoma and by a grant from the National Institute of Health (P20 RR016478 from the Ideal Networks for Biomedical Reseach Excellence (INBRE) Program of the National Center for Research Resources).