Comparison of diffusion approximation and Monte Carlo based finite element models for simulating thermal responses to laser irradiation in discrete vessels

Both diffusion approximation (DA) and Monte Carlo (MC) models have been used to simulate light distribution in multilayered human skin with or without discrete blood vessels. However, no detailed comparison of the light distribution, heat generation and induced thermal damage between these two models has been done for discrete vessels. Three models were constructed: (1) MC-based finite element method (FEM) model, referred to as MC-FEM; (2) DA-based FEM with simple scaling factors according to chromophore concentrations (SFCC) in the epidermis and vessels, referred to as DA-FEM-SFCC; and (3) DA-FEM with improved scaling factors (ISF) obtained by equalizing the total light energy depositions that are solved from the DA and MC models in the epidermis and vessels, respectively, referred to as DA-FEM-ISF. The results show that DA-FEM-SFCC underestimates the light energy deposition in the epidermis and vessels when compared to MC-FEM. The difference is nonlinearly dependent on wavelength, dermal blood volume fraction, vessel size and depth, etc. Thus, the temperature and damage profiles are also dramatically different. DA-FEM-ISF achieves much better results in calculating heat generation and induced thermal damage when compared to MC-FEM, and has the advantages of both calculation speed and accuracy. The disadvantage is that a multidimensional ISF table is needed for DA-FEM-ISF to be a practical modelling tool.


Introduction
In the treatment of port wine stain (PWS) birthmarks, advances in biomedical technology for imaging of discrete blood vessels in human skin have created a demand for accurate models which predict treatment results based on images that are acquired on an individual patient basis. Two subsequent processes need to be included in the models: (1) light propagation inside the tissue; and (2) transfer of heat generated by light absorption.
Both Monte Carlo (MC) and diffusion approximation (DA) models have been used to simulate light distribution in multilayered human skin with or without discrete vessels (Flock et al 1989, Keijzer et al 1989, Jacques and Wang 1995, Barton et al 1998, Shafirstein et al 2004. The MC approach is a rigorous but computationally intensive method to estimate the exact solution to the radiative transport equation (RTE) by sampling the set of all possible paths of photons through the turbid medium. It is universally agreed that MC is the 'gold standard' for simulating light distribution in tissue. Verkruysse et al (1997) and Kienle and Hibst (1995) used the MC method to calculate the total photon absorption inside a blood vessel in human skin and the subsequent heat generation and induced thermal damage. Theory DA is a low-order truncated solution to the RTE. The criteria for the validity of DA are that (1) light absorption is negligible compared to scattering, i.e., µ a /µ s (1 − g) 1; and (2) the radiance is only slightly anisotropic (Star 1995). In the case of PWS laser treatment, these criteria are invalid because the high absorption coefficients of melanin and blood violate the first criterion, and collimated laser light violates the second criterion in all but the most superficial layer of the skin.
The finite element method (FEM) is a numerical model used to solve partial differential equations that underlie most science and engineering applications. Hendrata and Franklin (2001) compared the performance of the FEM with MC in modelling light distributions in turbid media and concluded that, for a reasonably low error level (<1%), the FEM is at least one order of magnitude faster than MC when using a standard single processor computer. Other features of the FEM and associated software packages (e.g., Femlab 3.0, Comsol, Burlington, MA) also make the model a very attractive tool to simulate laser-tissue interaction. For example, the FEM is flexible enough to construct any user-defined geometry. It is also capable of solving coupled multi-physics problems-for our purposes, light distribution and heat transfer.
On the assumption that incident collimated light becomes diffuse once it enters tissue, DA-based FEM models were constructed to simulate heat generation and induced vessel damage profiles during laser-tissue interaction (Verhey et al 2003, Shafirstein et al 2004. Shafirstein et al (2004) scaled the laser energy according to chromophore concentration in order to compensate for the overestimation made by DA (details are illustrated below).
The question arises about the nature of the difference between the MC-based FEM and DA-based FEM models. Spott and Svaasand (2000) demonstrated that the DA model underestimates the fluence at the tissue surface due to the lack of a suitable description of collimated light sources. However, no detailed study has been done to compare different energy depositions inside a vessel using the DA and MC models, neither the resulting heat generation nor has the induced thermal damage been investigated. The aim of the current work is to clarify the feasibility of applying the DA-based FEM to calculate light distribution, heat generation and induced thermal damage by comparing it with a more accurate MC-based FEM. Based on the comparison, a new model that uses a scaling factor table to match the DA-based FEM with the MC-based FEM will be presented as a constructive remedy to the shortcomings of the current DA-based FEM models.

Skin models
Skin models used in this study consisted of a melanin-filled epidermis and dermis with certain blood volume fraction (BVF). For damage calculations, a blood vessel was embedded within the dermis in the centre of the laser beam. Optical and thermal properties of the epidermis, dermis and blood are listed in table 1. Note that in this study a simple geometry is used for illustration purpose, however, an important feature of the FEM is that irregular geometries can be easily simulated.

Monte Carlo method
In this study, we used the MC code originally designed by Keijzer (1993) and later modified and extended by Lucassen et al (1996). One of the capabilities of this code is to calculate energy deposition in model human skin containing multiple discrete vessels in a multilayered structure.

Diffusion approximation
In DA, light transport is described by the following differential equation: where is the optical diffusion coefficient (m 2 s −1 ) of skin, φ is the fluence rate (W m −2 ) of the incident laser light as a function of coordinates x and z, c e,d,b is the speed of light in the epidermis, dermis or blood vessel, µ a is the absorption coefficient (m −1 ) of skin, µ s is the scattering coefficient (m −1 ) of skin and g is the optical anisotropy factor.
The boundary condition on the skin surface receiving laser irradiation is given by (Shafirstein et al 2004) ( where P(t) laser is the laser irradiance (W m −2 ), r is the ratio of reflected light to the laser power output and n is the norm.
To model an infinitely wide laser beam, which is a reasonable assumption given the fact that the currently used clinical spot size (e.g., 10 mm) is much larger than the reduced scattering path length in skin (Keijzer 1993), the fluence at the lateral interfaces of the modelled skin and surrounding tissue is continuous, i.e., Since the depth of the modelled tissue is much greater than the light penetration depth, the same boundary condition (equation (4)) is used for the bottom interface.

Heat transfer equation
The heat transfer equation through conduction can be expressed as where T is the temperature ( • C) as a function of coordinates x, z and time t, ρ is the density

through light absorption and C(T)
is the specific heat capacity of the tissue. Shafirstein et al (2004) used a simplified numerical solution to approximate the discontinuity of phase transition at T 0 = 100 • C: where L is latent heat (J kg −1 ) for a liquid-to-gas phase transition, assuming that water evaporation occurs over 1 • C, i.e., T = 1. The boundary condition at the skin surface is modelled as convective heat transfer: where h = 150 W • C −1 m −2 is the heat transfer coefficient estimated for air flow over a flat plate for free convection and T ∞ = 26 • C is the ambient temperature. An Arrhenius-type integral was used to calculate the induced thermal damage (Pearce and Thomsen 1995): where (t) is the thermal damage parameter at time t, C d (0) is the original concentration of native tissue, C d (t) is the undamaged native tissue at time t, A is a frequency factor (s −1 ) or the molecule collision rate, A = 1.8 × 10 51 s −1 for the epidermis and A = 7.6 × 10 66 s −1 for blood, E a is an activation energy barrier, E a = 3.27 × 10 5 J mol −1 for the epidermis and E a = 4.55 × 10 5 J mol −1 for blood (Weaver and Stoll 1969), R is the universal gas constant (8.3144 J mol −1 K −1 ) and T is the absolute temperature (K). As illustrated in equation (8), tissue thermal damage is a complex nonlinear function of both temperature and time. By assuming that the coagulation threshold for irreversible damage occurs when = 1 (63.2% of the tissue is damaged), the coagulation profiles within specific vessels can be calculated (Weaver and Stoll 1969, Takata 1974, Pearce et al 1998, Shafirstein et al 2004.

Scaling factors based on chromophore concentration (SFCC)
Shafirstein et al (2004) introduced a scheme of simple scaling factors based on chromophore concentration, i.e., epidermal melanin concentration E e f and blood haematocrit E b f . The absorbed energy calculated from DA was multiplied by these scaling factors, such that the heat transfer equation (5) was modified into To be comparable with the parameters used in our MC-FEM model (  (9). An ISF value can be determined by varying it in a certain range and selecting the one that minimizes the error between the total deposited energy calculated by MC and unscaled DA-FEM for a certain skin structure, e.g., the epidermis or the blood vessel. This process can be easily done by any data analysis software. We employed a nonlinear optimization algorithm which is commercially implemented by Frontline Systems Inc. (Incline Village, NV) and available as add-in ('solver') in MS Office Excel R . This algorithm uses the generalized reduced gradient (GRG) method as implemented in an enhanced version of Lasdon and Waren's GRG2 code (Lasdon et al 1974, Lasdon and Waren 1983, Fylstra et al 1998. Detailed values of ISFs will be given in sections 3 and 4 for different situations.

Software implementation
A commercial software package (Femlab R 3.0 with Matlab R 6.5, Comsol, Burlington, MA) was used to construct the MC-FEM and DA-FEM models. In the MC-FEM model, a function in Matlab was programmed to interpolate the light distribution calculated by MC on regular spaced x-z coordinates into a user-defined FEM mesh, and heat transfer was calculated by the FEM solver that was provided by Femlab. In the DA-FEM model, both light distribution and heat transfer were calculated simultaneously by the FEM solver. The users only need to input the coefficients of equations (1)- (7) and (9). Figure 1 shows the comparison of energy deposition versus depth in a cross section through a vessel (radius R = 100 µm and the distance between the vessel top and the skin surface is 300 µm) resulting from 1 W cm −2 irradiation at two different wavelengths, i.e., 585 nm (figures 1(a) and (c)) and 595 nm (figures 1(b) and (d)). Two different values of BVF in the surrounding tissue were also considered: 1% representing normal skin (figures 1(a) and (b)) and 5% representing PWS skin (figures 1(c) and (d)). Three methods were used to calculate the energy deposition: MC, DA-FEM-SFCC and DA-FEM-ISF. For the current skin model, the ISF for epidermis, E e f , is 0.55 for both wavelengths of 585 and 595 nm. In the case of dermal BVF of 1%, the ISFs for the blood vessel, E b f , are 0.56 and 0.49 for 585 and 595 nm wavelengths, respectively. In the case of dermal BVF of 5%, the ISFs for the blood vessel, E b f , are 0.62 and 0.52 for 585 and 595 nm wavelengths, respectively.

Results
We further calculated the heat generation using equations (5)-(7) and induced epidermal and vessel damage using equation (8). Using the scenario of 5% BVF representing PWS skin, figures 2 and 3 compare the three methods by time-dependent temperature and damage profiles at wavelengths of 585 and 595 nm, respectively. A laser dosage of 8 J cm −2 was used for the simulation.

Effect of wavelength
From figure 1, it can be concluded that a DA-FEM without any scaling factors (equations (1)-(5)) will greatly overestimate energy deposition and heat generation. Figure 1 also shows that the DA-FEM-SFCC method generally results in a lower energy deposition in both the epidermis and blood vessel when compared to MC-FEM, while the DA-FEM-ISF results match much better with MC-FEM. However, when comparing the energy deposition in a vessel, all three methods have the same trend: while the energy deposition at the top of the vessel is higher at 585 than at 595 nm, a higher energy deposition is observed in the centre of the vessel at 595 than at 585 nm. Similar results are observed when we compare the effect of wavelength regardless of the method used, i.e., a high blood absorption coefficient at 585 nm causes a faster temperature rise at the top of the vessel (as shown in figures 2(a)-(c)) while a lower blood absorption coefficient at 595 nm causes slower but more uniform heating inside the vessel (as shown in figures 3(a)-(c)). dermal BVF is lower. This phenomenon is more pronounced at 585 nm than at 595 nm because the surrounding blood is a competing absorber and the blood absorption coefficient is higher at 585 nm than at 595 nm. It is concluded that the dermal BVF is an important parameter to be included when calculating thermal response induced by laser irradiation in a specific target vessel. This observation is similar to that reported in a previous study which compared the wavelengths of 577 and 585 nm using a layered tissue model (Pickering andvan Gemert 1991, Verkruysse et al 1997) and a recent study using the MC and finite difference methods to calculate thermal damage inside discrete vessels (Dai et al 2004). Table 2 lists E b f for different vessel depths and radii at a wavelength of 585 nm. Vessel depth is defined as the distance between the skin surface and top of the vessel. It is shown that E b f is nonlinearly dependent on vessel geometry. In other words, the discrepancy between MC-FEM and DA-FEM is a nonlinear function of vessel size and depth. By applying E b f as listed in table 2 into the DA-FEM-ISF model, the resulting damage profile inside a vessel is within ±5%, compared to that obtained by the MC-FEM model. Figure 4 compares the percentage of damage in different sized vessels calculated using the MC-FEM model immediately ( figure 4(a)) and 100 ms after laser irradiation ( figure 4(b)),  assuming that all other conditions were the same: light dosage 12 J cm −2 , pulse duration 1.5 ms, top of the vessel is 300 µm below the skin surface, and ignoring the epidermal damage which can be compensated for by cryogen spray cooling. Figure 4(b) shows that the 595 nm laser wavelength is superior to 585 nm in treating larger vessels (Svaasand et al 1995. For example, for a vessel radius of 100 µm, the 595 nm wavelength induces 100% damage in the vessels while 585 nm only induces 63% damage in vessels. Figure 4 also shows that a relaxation period after laser irradiation cannot be ignored when calculating thermal damage. For example, figure 4(a) shows that for a vessel radius of 100 µm, 585 and 595 nm wavelengths induce 15% and 25% damage in vessels, respectively, which are much smaller values than those shown in figure 4(b).

Effect of epidermal melanin absorption
As shown in figure 1, the DA-FEM-SFCC model severely underestimates the energy deposition in the epidermis-approximately 13% of the amount calculated by MC-FEM. Because epidermal safety is a critical factor for selecting the maximal permissible light dosage for PWS laser therapy (Dai et al 2004), an accurate estimation of light energy deposition is essential. The DA-FEM-ISF model with E e f = 0.55 achieves a much closer representation of MC results than DA-FEM-SFCC. We also found that E e f = 0.55 can be applied for other   (Jacques 1998) where µ melanosome a = 1.7 × 10 11 λ −3.48 (mm −1 ) is the absorption coefficient of a melanosome at wavelength λ expressed in nm, µ T a = 7.84 × 10 7 λ −3.255 (mm −1 ) is the baseline absorption coefficient of human skin (Saidi 1992) and f m is the volume fraction of melanosomes in the epidermis. It should be noted that table 3 was obtained assuming a laser irradiation of 585 nm, dermal BVF of 5%, vessel depth of 300 µm and vessel radius of 100 µm.

Multidimensional scaling factor table
Sections 4.1-4.4 also show that there are at least five parameters, i.e., wavelength, dermal BVF, vessel size and depth, and epidermal melanin concentration, that need to be taken into account for a complete ISF table. Table 2 varies two while fixing the other three parameters. Table 3 varies one while fixing the other four parameters. These two tables present only a small portion of a five-dimensional ISF space. It is impractical to give a complete five-dimensional ISF table in the current work; however, the readers can create their own ISF tables according to their particular needs following the same principle. Table 4 lists a detailed comparison on the ratios of deposited energy (as shown in figure 1) calculated by the DA-FEM-SFCC and DA-FEM-ISF to that calculated by the MC method. The ratios are different at the top, centre and bottom of the vessel. These inconsistent ratios make it impractical to use a simple scaling factor to exactly match the DA-FEM with the MC-FEM in details. If a more rigorous simulation is required, neither SFCC nor ISF are as good as MC-FEM. However, the speed of the FEM to simulate irregular blood vessel geometry is an advantage over the MC method. Thus, the DA-FEM-ISF model has the advantages of both speed and relative accuracy when compared to other two methods. Therefore, if approximate temperature increase and damage estimation (within error of ±5%) are of interest, DA-FEM-ISF is a good choice.

Conclusions
This study compares the differences between the diffusion approximation based finite element model (DA-FEM) and the Monte Carlo based finite element model (MC-FEM) to simulate the laser-induced thermal response of human skin containing discrete vessels. The DA-FEM-SFCC model that uses chromophore concentration as a scaling factor underestimates the deposited energies and thermal damages in both the epidermis and blood vessels. The ratios between DA and MC are a nonlinear function of laser wavelength, dermal BVF, vessel geometries and epidermal melanin concentration, such that a simple scaling factor becomes impractical. The MC-FEM model is recommended where rigorous simulation is required. The DA-FEM-ISF model uses two improved scaling factors that equalize the total energy deposition in the epidermis and blood vessels. It achieves a much closer representation of the MC-FEM model as compared to the DA-FEM-SFCC model. Although a multidimensional scaling factor space exists, once the ISF table is constructed as needed, this method has the advantages of both speed and relative accuracy when compared to the other two methods.