Numerical simulation of nanoparticles assisted laser photothermal therapy: a comparison of the P1-approximation and discrete ordinate methods

Photothermal therapy (PTT) with combined use of laser radiation and photon absorber nanoparticles is a promising technique to treat cancer. Treatment planning and devising appropriate protocols for cancer photo thermal therapy require the computational simulation of coupled physical phenomena, such as radiation, conduction, and blood perfusion. The P1-approximation is a numerical method to solve radiation heat transfer which features the advantage of being computationally fast and, therefore, desirable for PTT simulations. However, the method is known to become inaccurate under certain conditions. In this study, the P1-approximation and the accurate discrete ordinate method were applied to solve a set of test problems idealized to portray conditions encountered in PTT. The test problems were one-dimensional, and the radiation scattering was assumed as isotropic. Tissues composed by layers with different properties were considered, including cases in which gold nanoparticles were embedded in the tissue to increase photon absorption. For the problems considered here, the P1-approximation and discrete ordinate method results presented quite good agreement for the time-dependent temperature distribution, which is the quantity of interest in PTT.


Introduction
Photothermal therapy (PTT) is a procedure in which thermal radiation is used to treat diseases. It has become quite promising for the treatment of various types of cancer with the recent development of nanotechnology. The combined use of nanoparticles and collimated laser light allows to selectively heat and induce hyperthermia of cancerous cells, while preserving the surrounding healthy tissues, thus making the treatment minimally invasive. Due to this reason, there is a growing interest on using PTT to treat cancer, as can be verified in the review by Bayazitoglu [1].
Near-infrared (NIR) lasers are of common use in medicine, and research works involving their application for photothermal cancer therapy can be found in the literature [2][3][4][5][6][7][8][9][10][11][12][13]. In the range of the NIR spectrum encompassing the wavelengths between 700 and 1100 nm, the extinction coefficient of most biological tissues is low [1]. Therefore, laser radiation within such a spectral range propagates into the human body without significant attenuation. On the other hand, there are nanoparticles that strongly absorb NIR radiation, which can be embedded in tissues to increase the absorption in specific regions. Most of these particles are constituted of noble metals and can be designed for a specific resonance frequency [14]. Also, noble metals present low levels of toxicity [2]. More recently, interest in the application of carbon-based nanoparticles to PTT has increased, such as carbon nanotubes [15]. Differently of the nanoparticles made of noble metals, the ones composed of carbon-based nanomaterials are not tunable for being highly absorbant in a specific wavelength; instead, they have a high extinction coefficient along a large range of the visible and NIR spectra.
The application of PTT for cancer therapy has attracted the attention of researchers from different areas. In vivo experiments are of primary importance, since they allow to assay the feasibility of treating cancers with different PTT strategies, which can vary mostly due to differences in the laser configuration, type of nanoparticles, and methods to deliver the particles, as can be verified from the reviews in references [1,16,17].
Another field of interest concerns the nanoparticles themselves, since the success of the therapy depends on their optical properties. Examples of works in this area include those by Tuersun and Han [14], An et al. [18], and Lui et al. [19]. The first shows the relation of the backscattering and absorption properties of nanoshells with the dimensions of the silica core radii and gold shells, the second presents a numerical investigation about the radiative properties of silver nanorod dimers assemblies, while the last one is concerned with a model to calculate resonance properties of nanoshells with metal cores.
Although the photothermal therapy strategies must be ultimately validated for medical applications by in vivo experiments, numerical methods are a powerful tool, which can be used to simulate the heat transfer processes occurring during a PTT procedure. Numerical simulations can partially replace experimental investigation, with the advantage that some parameters can easily be modified in a computer code, allowing simulating an enormous variety of nanoparticles, tissues, and lasers with minor increase of costs, time, and manual labor. Therefore, the development of reliable numerical models for simulating PTT procedures can be considerably useful, and some models were proposed. The majority of the models solve the Pennes energy equation [20] coupled to the radiative transfer equation [21,22]. In such models, the most complex task is the solution of the radiative heat transfer problem. Tjahjono and Bayazitoglu [7], and Vera and Bayazitoglu [8] considered an one-dimensional problem with isotropic radiation scattering. For the radiative transfer solution, they employed the spherical harmonics with the first-order approximation (P1-approximation). This method is computationally fast, which is highly desirable for applications to simulate the complex processes of heat transfer in PTT, specially with irregular three-dimensional geometries such as those which can occur in real tumors. However, the method can be inaccurate in certain cases. For example, the P1-approximation becomes less accurate in multidimensional geometries with thin media associated to strongly anisotropic radiation scattering [22].
In spite of the fast growing of computers' processing power, simulating coupled radiation, conduction, and blood perfusion, which can occur in real three-dimensional PTT treatments, is still considerably expensive in terms of computational time, especially because of the radiation solution. Thus, using a computationally cheap method as the P1-approximation is highly desirable, and because of that, the method has been applied in researches aimed at gaining knowledge about PTT simulations, including cases in which simple one-dimensional problems are approached, examples are [7, 8, 10, and 23]. Nonetheless, the P1-approximation can be inaccurate in some cases. As demonstrated by Mengüç and Viskanta [24], under conditions for which the P1-approximation becomes inaccurate, the discrete ordinates method can still provide accurate results. The present study compares the P1-approximation with the more accurate discrete ordinate method with eighth order of angular discretization (S 8 ). The approached problems simulate cases of skin and breast cancer. In both cases, the geometry is approximated as one-dimensional, and the medium scatters radiation isotropically. Tissues with and without nanoparticles embedded in are considered. Results are presented for the temperature field obtained by solving the coupled radiation, conduction, and blood perfusion, as well as for the divergence of the radiative heat flux, to evaluate the methods for pure radiation solutions.

Fundamental equations and modeling
In the problems examined here, a tissue is irradiated by collimated laser light. While the temperature of the tissue rises due to radiative energy absorption, heat is removed by conduction and blood perfusion. The quantity of major interest for PTT is the time and spatially dependent tissue temperature, T, since the success of the treatment depends on the magnitudes of temperature in the tissues and the time that they are subjected to high temperatures. To obtain the temperature distribution, the following form of the energy equation needs to be solved: where ρ and c p are, respectively, the tissue density and specific heat, t is time, k is the thermal conductivity of the tissue, and the last term accounts for the removal of heat from the tissue due to blood perfusion, as proposed by Pennes [20]. With the exception of T, all quantities appearing in this last term are relative to blood, as indicated by the subscript b, being ρ b , c b , T b ,and v b ,respectively, the density, specific heat, blood temperature, and blood perfusion rate. The laser radiation absorption is accounted for the divergence of the radiative heat flux, in the second term of the right-hand side of Eq. (1).
The problems considered here are one-dimensional. In this case, Eq. (1) can be reduced and written for the x direction as Obtaining q r,x requires the solution of the radiative transfer equation (RTE), which gives the variation of the intensity, I , as a radiation ray propagates in a particular direction, s. For quasi-steady conditions, the RTE takes the following form: where I b, is the blackbody spectral intensity given by the Planck function, κ is the absorption coefficient, σ sc, is the scattering coefficient, and β is the extinction coefficient, which is the sum of the absorption and scattering coefficients. The subscript indicates dependence with respect to the wavelength, while the subscript sc stands for scattering. The last term of the right-hand side of the RTE accounts for the portion of the radiative energy coming from all possible directions that is scattered toward the considered direction of propagation.
The I ,i is the intensity of radiation in the incidence direction, as indicated by the subscript i; and i are the coordinates of the solid angles in the directions of propagation and incidence, respectively; and � (�, � i ) is the scattering phase function, which describes the probability that a radiation ray coming from the direction i be scattered into the direction . The specific conditions regarding the PTT problems examined here allow some simplifications of Eq. (3). The temperature necessary to induce hyperthermia is low, about 43 °C [2]. Thus, emission of radiation from the tissue is insignificant compared to the laser power, and the term containing the blackbody intensity, which accounts for emission, can be neglected. Also, the scattering was assumed as isotropic, leading to = 1. Also, as lasers are monochromatic, the dependence of properties with the wavelength is not required. With these simplifications, Eq. (3) can be rewritten as The intensity in Eq. (4) has two components: one collimated coming directly from the laser, and another diffuse that arises due to scattering. However, as presented by Modest [22], the problem can be reduced to one with solely diffuse radiation by placing the effects of the collimated radiation in a source term, S c : In the above equation, s denotes the depth below the surface where the irradiation is imposed, q o is the laser power, and R e is the reflectivity of the boundary surface irradiated by the laser, which is external to the tissue slab, as indicated by the subscript e. Considering the boundaries as diffuse and gray, and neglecting surface emission, the boundary condition for Eq.
where R is the reflectivity of the surface boundary, which is internal to the tissue slab, n is a vector outward normal to the surface, ŝ is a unit vector with the direction of I i , such that n ·ŝ = cos(θ i ), being θ i relative to the surface normal.
The radiative heat flux crossing an element of area normal to the direction of x is a result of intensities incident from all directions. Thus, it can be computed as, where θ is relative to x.
Two different numerical methods were applied for solving the radiative transfer equation, namely, discrete ordinate method with eighth order of approximation for the angular discretization (S 8 ) and spherical harmonics with the first order of approximation (P1-approximation). Both S 8 and P1-approximation were implemented in a Fortran90 code and coupled to the same finite volume code also written in Fortran90, which contains routines for the conduction and perfusion problems to provide the transient solution of the coupled radiation, conduction and perfusion heat transfer. The finite volume method was implemented following the procedures presented by Maliska [25], while the discrete ordinate and P1-approximation were implemented according to the procedures presented by Howell et al. [21], and Modest [22]. The codes are able to consider nonuniform meshes, since it is important to handle with the nonhomogeneities of biological tissues and accurately represent the temperature gradients which occur in PTT.
The same finite volume code was applied to the two solution methods examined for the radiative transfer equation. Therefore, the results presented below allow a direct comparison between the discrete ordinate and P1-approximation, as shown in the next section.

Results
To compare the P1-approximation with the discrete ordinate method, a set of test problems were solved for two tissue slabs composed of different layers, as follows.

Skin cancer
This case was based on that presented by Dombrovsky et al. [10]. The geometry consists of an one-dimensional slab with five different layers, corresponding to epidermis, cancer tumor, papillary dermis, reticular dermis, and fat, as represented by Fig. 1. The slab is irradiated from the left by collimated laser radiation with wavelength of 0.6328   μm. The emissive power of the laser is 20 kW/m 2 , and is imposed on the tissue periodically, according to the function shown in Fig. 2. It is assumed that all radiative energy incident on the slab penetrates the tissue, that is, the reflectivity of the external surface of the slab, appearing in Eqs. (6) and (7), is R e = 0. Internally, the reflectivity of the surface boundaries are R = 1, at the left, and R = 0.49, at the right. The conduction boundary conditions for Eq. (2) are prescribed temperature of 37 °C, on the left surface, and convective heat transfer, on the right, for simulating heat transfer to inward body tissues. For the inner boundary, the convective heat transfer coefficient is h = 50 W/(m 2 K) [10], while the initial temperature and the temperature of the surrounding medium exchanging heat with the surface on the right are both equal to 37 °C. The properties of the tissue layers without nanoparticles are shown in Table 1. Figure 3a presents the temperatures as a function of the depth in the tissue at four different instants of the periodic heating process: 10, 25, 35, and 50 s. The temperatures were computed with the P1-approximation and discrete ordinate method with eighth order of angular discretization (S 8 ). The results obtained with the two methods agreed well, the temperature presented similar profile, the maximum difference was about 0.3 °C for 35 s of heating, and the relative difference was not higher than 0.75 %. Furthermore, an interesting observation which can be made by inspecting Fig. 3a is that the peak of temperature occurs about 1.5 mm deep in the tissue, away from the tumor. In this case, the laser heating is more likely to kill healthy cells than the cells in the tumor region. Thus, to reach temperatures sufficiently high to kill cancerous cells, healthy tissues would be overheated to temperatures even higher, which would be harmful for the patient.
The absorption of heat in the tumor region can be increased using nanoparticles. Figure 3b shows the temperature distribution when the first and second layers, representing the epidermis and the tumor are embedded with nanoparticles, which are gold nanoshells. The optical properties of the tissue layers with the nanoshells were obtained in accordance with the procedure presented by Dombrovsky et al. [10], using the following equations to compute, respectively, the absorption and scattering coefficients: where Q a and Q s are the dimensionless efficiency factor and transport efficiency factor of scattering for single particles, κ t and σ sc,t are the absorption and scattering coefficients of the tissue without nanoparticles, r is the radius of the nanoparticles, and f v is the volume fraction of the nanoparticles relative to the embedded tissue. Dombrovsky et al. [10] considered that the radiation source is a helium--neon laser with wavelength of 0.6328 μm, and that the nanoparticles are silica-core gold nanoshells with radius r = 20 nm (which give the maximum absorption at the refered wavelength). The radius of the silica core was considered equal to 0.725r. For such nanoparticles, and assuming a refractive index n = 1.45 relative to the ambient human tissues, Dombrovsky et al. [10] applied the Mie theory to estimate the following values for the efficiency factors: Q a = 7.828 and Q s = 1.144. Assuming a volume fraction of the nanoparticles f v = 10 −5 , Eqs. (9) and (10) provide κ = 3115.6 m −1 and σ sc = 2789 m −1 , for the epidermis, and κ = 2985.5 m −1 and σ sc = 1029 m −1 , for the tumor. In layers without nanoparticles, f v = 0, and thus, the absorption and scattering coefficients reduce to κ = κ t and σ sc = σ sc,t .
As can be seen in Fig. 3a, b, the peak temperature was moved to the tumor region (0.3 mm < x < 0.4 mm) when nanoparticles were used. Similarly to the case without nanoparticles, the P1-approximation and S 8 results presented good agreement, with a difference lower than 1 % for the peak temperature.
The results for the temperature shown in Fig. 3a, b were obtained with the presence of conduction and blood perfusion, since both are relevant in PTT. A comparison of the S 8 and P1-approximation predictions in a pure radiation problem is presented in Fig. 4a, b. The divergence of the radiative heat flux is shown: (a) in the tissue without nanoparticles, and (b) in the tissue embedded with nanoparticles. In the case without nanoparticles (see Fig. 4a), the S 8 predicted higher radiation absorption than that obtained with the P1-approximation in all layers, which is consistent with the higher temperatures predicted using the S 8 method (see Fig. 3a). The main differences in the divergence of the radiative heat flux occurred in the first two tissues; the differences were about 10% near the boundary and higher in the region of the tumor. However, the absorption coefficient is relatively small in the locations where the higher errors in the divergence of the radiative flux do occur, and thus have small influence on the temperature.
In the case of the tissue embedded with nanoparticles, the divergence of the radiative heat flux, shown in Fig. 4b, is considerably higher than that for the case without nanoparticles, since the particles increase radiation absorption. The results obtained with the S 8 and P1-approximation presented good agreement: the maximum differences occurred in the first two tissue layers and were about 5 %. Differently than for the tissue without nanoparticles, the divergence of the radiative heat flux obtained with the S 8 was lower than that obtained with the P1-approximation in the region where the peak temperature occurs.
A consistent behaviour was verified in the results for the temperature, since now is the S 8 method, instead of the P1-approximation, which predicts a higher peak temperature.
The addition of nanoparticles in the two first layers of tissue makes the results for the divergence of the radiation heat flux obtained with the S 8 and P1 methods compare better, as can be seen in Fig. 4. This is consistent with the results for the temperature, shown in Fig. 3, which also presented better agreement in the two first layers. Although the third, fourth, and fifth deeper tissue layers are not embedded with nanoparticles and, thus, their optical properties remain unchanged, the temperature profile continues presenting better agreement as compared with the case in which no nanoparticles are embedded in the two first layers. It occurs because the nanoparticles strongly increase the absorption in the region of the epidermis and tumor, which makes the radiation absorption relatively irrelevant in the three deeper tissue layers, and therefore, radiation absorption in this region without nanoparticles has no significant effect on the heat transfer process, which is dominated by conduction.
In general, the results presented by Figs. 3 and 4 show that the temperatures predicted by both methods examined here exhibit small discrepancies. In fact, although some significant differences in the divergence of the radiative heat flux were verified, they occurred when the radiation absorption was relatively low and did not strongly affected the temperature field.
To address the possibility of further model reduction for the radiative transfer problem, the diffusive contribution of the radiation intensity was neglected for the computation of the divergence of the radiative flux, which was computed only with the collimated intensity given by Beer-Lambert's law. The temperature variation for the case of the tissues containing nanoparticles is presented by Fig. 5, at time t = 35 s, for the radiation solutions given by discrete ordinates, P1-approximation and Beer-Lambert's law. The discrepancy between the temperature that results from Beer-Lambert's law is 0.75 °C in the region of the tumor, but much larger in deeper regions. Moreover, the temperature profile obtained with the Beer-Lambert law differs from that obtained with the two other methods. Therefore, while the P1-approximation provides an accurate prediction of the temperature in the region, Beer-Lambert's law results in inappropriate predictions for the hyperthermia treatment in this skin problem.
For the results presented above, the layers representing the dermis, tumor, papillary dermis, reticular dermis, and fat were discretized with 20, 40, 16, 32, and 64 finite volumes. Such number of volumes were selected based on a grid convergence analysis. Also, the S 8 solution with discrete ordinates were compared with those for lower orders (S 6 and S 4 ) and no significant differences were verified. All presented results were obtained using an Intel ® Core ™ i7-4770S @ 3.10 GHz, in a 64-bit system, with 8 GiB

Breast cancer
For this case, the geometry once again was one-dimensional, and the scattering was assumed as isotropic. However, the tissue properties, dimensions and laser power were selected to simulate the heat transfer processes occurring in a photothermal treatment of breast cancer. The geometry is represented in Fig. 6, while the layers properties are given in Table 2. The specific heat, density, and thermal conductivity were obtained from He et al. [26], the metabolic heat of the tumor was estimated based on the analysis presented by Dombrovsky et al. [10], the metabolic heat generated by the gland was taken from Vo-Dinh [27], and the optical properties of the tumor and breast gland were obtained from Fatini [28]. The laser irradiates the left surface and is periodically applied according to the function shown in Fig. 2, as in the previous case. However, the laser emissive power is now 28 kW/m 2 . The boundary conditions are exactly the same as that of the case involving the skin cancer.
Results were obtained for the tissue without nanoparticles and when the 2/3 final part of the layer corresponding to the tumor was embedded with gold nanoshells with the goal of bringing the peak temperature to the tumor region. The nanoparticles are shells with the same features as for the skin cancer case. The properties of the tissue embedded with these particles were again determined according to the procedure presented by Dombrovsky et al. [10], using Eqs. (9) and (10). However, at this time the volume fraction of the nanoshells relative to the embedded tissue was f v = 2 · 10 −7 , which resulted in the following values of absorption and scattering coefficients: κ = 67.11 m −1 and σ sc = 1508.58 m −1 . Figure 7a, b present the temperature distributions at different times, obtained for the tumor without and loaded with nanoparticles, respectively. These figures show that, with the addition of nanoparticles, the peak temperature increases from approximately 41 °C to almost 44.5 °C, while shifts 2 mm to right, getting closer to the center of  the tumor region. For this case, the maximum temperature difference obtained with the P1-approximation and S 8 discrete ordinates is less than 1 %, both with the tumor containing or not containing nanoparticles.
Results for the divergence of the radiative heat flux are presented for the tumor without and with nanoparticles, respectively by Fig. 8a, b. For the tumor without nanoparticles, the larger discrepancies occur at the irradiated boundary and at the region of maximum divergence. However, the discrepancies between the two solutions decrease for deeper regions farther from the external surface. For the case of the tumor containing nanoparticles, two peaks of the divergence of the heat flux are observed: one ahead of the tumor and another ahead of the region containing nanoparticles inside the tumor. Moreover, the results obtained with S 8 discrete ordinates and P1-approximation are very similar in the region containing nanoparticles, with larger discrepancies at the boundary, and at the first peak of the divergence of the radiative heat flux. Anyhow, these local large discrepancies do not significantly affect the temperature distribution in the region, as can be noticed in Fig. 7a, b. Figure 9 shows results for the temperature at 35 s for the case in which the radiation solution was obtained with the Beer-Lambert law, discrete ordinates, and P1-approximation. The results obtained with the simpler Beer-Lambert law were totally discrepant from those obtained with the two other methods; the peak temperature differed considerably in its location and value, which was about 14 °C higher than those obtained with the S 8 and P1-approximation methods. These results demonstrate that the chosen radiation model has an important impact on the temperature distribution and indicate that the simplified Beer-Lambert's model, which neglects the diffusive contribution for the radiation intensity, cannot be used for the prediction of the temperature distribution in this case dealing with breast cancer.
Similarly to the case involving skin cancer, nonuniform finite volumes were used. The first three layers representing gland, tumor, and gland were discretized using 80, 1440, and 140 volumes, respectively. Such number of volumes were based on a grid convergence analysis. Also, the S 8 discrete ordinates solution was converged with respect to lower orders of angular discretization, such as, S 6 and S 4 . The computational time required to obtain the P1-approximation results with and without nanoparticles was, respectively, 1838 and 1832 s. In the case of the S 8 method, the computational time was 2420 s, for the case without nanoparticles, and 2412 s when nanoparticles were added.

Conclusion
The P1-approximation is a simple and computationally cheap method to solve radiation heat transfer problems. Thus, its application for time-consuming coupled Fig. 8 Divergence of the radiative heat flux for the breast cancer: a tissue without nanoparticles and b tissue embedded with gold nanoshells Fig. 9 Temperature distribution for the breast cancer at 35 s in the tissue embedded with gold nanoshells simulations such that relative to photothermal therapy (PTT) is highly desirable. However, the method is known to become inaccurate under certain conditions. Aiming at evaluating the accuracy of the P1-approximation for PTT simulations, the method was applied to tackle the radiation heat transfer in two test problems, and the results were compared with other results obtained with the accurate discrete ordinate method. The test problems required the one-dimensional coupled solution of conduction, blood perfusion, and radiation in isotropically scattering media. Two tissue slabs were considered: one representing a case of skin cancer and another representing a case of breast cancer. In both cases, simulations were performed for tissues with and without gold nanoshells embedded in. The obtained results were presented for the time-dependent temperature distribution, which is usually the quantity of interests in PTT, as well as for the divergence of the radiative heat flux, to compare the methods in a pure radiation problem.
Comparing the divergence of the radiative heat flux with the predicted temperature, it was demonstrated how variations in the radiation results are perceived on the temperature profile. The differences in the divergence of the radiative heat flux obtained with the S 8 and P1-approximation were higher in regions of the tissue where the absorption was relatively low so that having minor influence on the temperature distribution. In general, the profiles of the divergence of the radiative heat flux obtained with the S 8 and P1-approximation were similar, and the values presented good agreement where the absorption was higher. Thus, results for the temperature computed using both methods presented good agreement. Additionally, results for the temperature were obtained using a much simpler approach for the radiation solution, the Beer-Lambert law. The results were discrepant from that obtained with the S 8 and P1-approximation, thus showing that the radiation model can have considerable influence on the temperature field.
For future works, we suggest to extend the analysis to multidimensional geometries and include anisotropic scattering of radiation, since those are factors which may lead to a loss of accuracy in the P1-approximation results.