Development of a Predictive Monte Carlo Radiative Transfer Model for Ablative Fractional Skin Lasers

It is possible to enhance topical drug delivery by pretreatment of the skin with ablative fractional lasers (AFLs). However, the parameters to use for a given AFL to achieve the desired depth of ablation or the desired therapeutic or cosmetic outcome are hard to predict. This leaves open the real possibility of overapplication or underapplication of laser energy to the skin. In this study, we developed a numerical model consisting of a Monte Carlo radiative transfer (MCRT) code coupled to a heat transfer and tissue damage algorithm. The simulation is designed to predict the depth effects of AFL on the skin, verified with in vitro experiments in porcine skin via optical coherence tomography (OCT) imaging. Ex vivo porcine skin is irradiated with increasing energies (50–400 mJ/pixel) from a CO2 AFL. The depth of microscopic treatment zones is measured and compared with our numerical model. The data from the OCT images and MCRT model complement each other well. Nonablative thermal effects on surrounding tissue are also discussed. This model, therefore, provides an initial step toward a predictive determination of the effects of AFL on the skin. Lasers Surg. Med. © 2020 The Authors. Lasers in Surgery and Medicine published by Wiley Periodicals LLC


INTRODUCTION
Ablative fractional lasers (AFLs), commonly Er:YAG (Erbium Yttrium Aluminum Garnet-2940 nm) or CO 2 (carbon dioxide−10600 nm), are widely used in both clinical and cosmetic applications. These lasers can vaporize tissue and create an array of microscopic treatment zones (MTZ) in the skin [1][2][3]. The formation of these MTZs presents a number of opportunities that may be taken advantage of to improve treatment outcomes.
Topical drug delivery, a cornerstone in dermatological medicine, is dependent on cutaneous drug penetration, which influences uptake and accumulation of the drug at the target site. However, drug penetration is typically reduced by the presence of the stratum corneum, the protective barrier of the skin [4]. Several studies have shown that pre-treatment of skin with AFL prior to topical drug application can increase cutaneous drug uptake [5][6][7]. This is particularly useful not only for increasing uptake and depth of penetration but also for allowing penetration of larger and more hydrophilic molecules that typically cannot pass through the stratum corneum.
Pre-treatment with AFLs has been used in applications of photodynamic therapy (PDT), where greater uptake of the light-activated drug and enhanced treatment efficacy have been shown [8][9][10].
Additionally, treatment with AFLs is widely used in aesthetic dermatology practices as this can result in improved cosmetic outcomes. Thermal stimulation of the dermis has been shown to enhance the production of collagen and restructure connective tissue [1]. Traditionally, laser skin resurfacing was carried out using continuous wave ablation of the top layer of the skin. However, CO 2 lasers, while effective, caused more severe side effects and prolonged downtime for the patient, and while the Er: YAG laser caused fewer side effects, this was less effective than the CO 2 laser [1,11,12]. AFLs have been developed to provide sufficient thermal damage to the dermis while sparing the majority of the epidermis, which provides an acceptable balance of efficacy and side effects [2]. The key to safe and effective use of AFLs is the understanding of the depth and severity of the effect on the tissue in relation to laser energy and heating.
Optical coherence tomography (OCT) is a noninvasive optical imaging method relying on the interference of two lasers and, when applied in dermatology, allows structural information of the skin to be reconstructed [13]. While histological analysis requires a biopsy sample, OCT may be carried out in vivo, and these high-resolution images have been shown to correlate well with histological data while investigating MTZs in the skin as a result of AFL [7,14,15]. Although OCT presents a noninvasive characterization modality for AFL, it is still inherently a postanalytical process.
Light-tissue interactions may be modeled by the use of Monte Carlo radiative transfer (MCRT) techniques [16][17][18]. In this paper, we present our initial numerical model to provide predictive capabilities for AFL skin interactions and we verify our code with OCT imaging of CO 2 -induced MTZs in ex vivo porcine skin.

Porcine Skin Irradiation
Samples are cut from fresh porcine belly skin sourced from a local butcher. The samples are cut to approximately 20 × 20 mm area and 5 mm thick. The porcine skin is then irradiated with increasing energies (50-400 mJ/pixel, 25 mJ/ pixel intervals) from a fractionated 70 W CO 2 (10.6 μm) laser (Pixel CO 2 ; Alma Lasers, Nürnberg, Germany). The Pixel CO 2 laser is a pulsed laser, which has a triangular temporal pulse shape and a beam width of 250 μm. The AFL outputs a 9 × 9 grid of microbeams (pixels), and as power is kept constant an increase in energy is obtained by lengthening the exposure duration (10-1000 milliseconds), thus a 50 mJ irradiation yields a pulselength ≈ 60 millisecond. The laser was operated in "SuperPulse" mode, which yields one laser pulse per operation. Laser irradiation is repeated on separate fresh samples. A smoke evacuator (ACU-EVAC-II; Schuco, Leavesden, UK) is used to remove laser-generated plumes from the air.

OCT
Phase-sensitive OCT (PhS-OCT) was applied to obtain high-resolution structural information from the irradiated samples [19,20]. A more detailed description of the bespoke system is described elsewhere [21,22]. Briefly, the system produces structural images with an axial resolution of approximately 9 μm and a lateral resolution of approximately 15 μm across a 3.5 mm × 3.5 mm scan range.

MCRT
MCRT is the "gold-standard" for simulating the transport of light through biological tissue [23,24]. This is due to its flexibility in modeling three-dimensional (3D) geometries, various light sources, and micro-physics, such as absorption, scattering, fluorescence, and polarization [25][26][27][28][29]. It uses interaction probabilities and random numbers to model the "random walk" that photons undergo in a turbid medium. We simulate the propagation of power packets, which represent photons with a given power, derived from the incident radiant source. MCRT has been used to model light-tissue interactions in many different medical and biophotonics applications [30][31][32]. MCRT is used here to calculate the energy deposited by the laser, which is then passed to the heat transport simulation. The original MCRT code was developed for astronomy applications [27,33] and has since been adapted and validated for use in the medical field in previous works [30,34].
The tissue medium for the MCRT and heat transport simulations is a 3D voxel model. This allows the variation of optical and thermal properties from voxel to voxel, making it the ideal type of grid for modeling tissue ablation. A voxel model also allows us to update the medium as the simulation progression, allowing the modeling of a moving ablation front. We use 160 × 160 × 160 voxels, representing a tissue sample size of 0.06 cm × 0.06 cm × 0.12 cm. We assume that the porcine skin is uniform, so initially, our voxel model is uniform, and the optical properties of porcine skin at the wavelength of interest (10.6 μm) is mainly that of water mixed with protein [35]. Thus, we adopt an absorption coefficient of 850 cm −1 . As the absorption coefficient is large and absorption dominates at this wavelength, we assume that no scattering takes place, thus the scattering coefficient is 0 cm −1 . As we assume the scattering coefficient is zero, we could use an analytical method such as Beer's law to model the propagation of photons through the tissue. However, we use MCRT to allow simulations of different laser wavelengths where scattering may dominate.
To calculate the energy absorbed in the porcine tissue via the laser, we use the pathlength counter method devised by Lucy [36]. The energy absorbed in voxel i, E i abs , is therefore calculated as: where P is the power This grid of absorbed energy per timestep is then passed to the heat transport portion of the simulation so that the heat diffusion in the porcine tissue can be calculated.

Heat Simulation
We solve the nonlinear heat equation using the finite difference method in Cartesian coordinates with a heat sourceq [37].
where ( ) T x y z t , , , is the temperature as a function of space and time, κ is the thermal conductivity, ρ is the density, c p is the heat capacity, andq is the source/sink term as a function of time and space. As the heat equation does not have an analytic solution in arbitrary 3D geometries, we employ a finite difference method (FDM) to solve it numerically. The specific FDM we use is the simple explicit method, due to its ease of implementation and parallelization. Applying this scheme to Equation (2) yields where, in the x direction (the other directions have the same form, but are omitted for brevity): The κ ρ ± ± , , and ± c p terms occur, due to the nonlinearity of the heat equation, as the medium is inhomogeneous. Thus, we solve the nonlinear heat equation by averaging these thermal properties across the discontinuities in the medium. Cooling of the tissue is not currently considered. As we employ the simple explicit finite difference method, the time step is constrained. This can lead to large computational times, as the time step has to be kept small to avoid loss of numerical precision. These large computational times are offset by parallelizing the code using the 'halo swap' technique [38].

Tissue Damage Simulation
The final portion of the simulation is the tissue damage model. We split the tissue damage model into two separate categories: physical and coagulation. The coagulation portion accounts for all the thermal damage accrued under 100°C, where the physical damage model accounts for all damage above 100°C. To model coagulation damage, we use the Arrhenius damage model, widely used in the literature to model thermal damage [39][40][41].
where Ω is the damage value [−]; A is the frequency factor [s −1 ]; ΔE the activation energy [L mol −1 ]; R is the universal gas constant [J mol −1 K −1 ]; T is the temperature in kelvins, and t i and t f are the initial and final times that the temperature is above the threshold damage temperature.
The Arrhenius equation gives a number, which can be related to the severity of the burn [42]. We adopt values of 3.1 × 10 98 for A, and 6.3 × 10 5 for ΔE [43,40,44].
There are several proposed ablation models that explain the underlying mechanism [45][46][47][48]. To model the ablation process, we adopt an ablation model based upon water boiling within the tissue. Though this model is simplistic, we find it gives good results when compared with our experimental evidence. The physical portion of the damage model changes the medium as the temperatures increases. This damage is split into two areas: water boiling and ablation. When the temperature in a voxel reaches 100°C, we track the energy deposited in each voxel. As this energy is accrued, we start to boil off the water in the voxel as a function of energy in the voxel. As the water boils off, this changes the thermal and optical properties of the voxel [49].
where P ro is the amount of protein, and W is the amount of water in a voxel, Q current is the energy accrued in a voxel, Q vapor is the energy required to boil off the water in a voxel, mass M. We set P ro as 0.25 and W as 0.75. The updated optical properties are used in the next time step in the MCRT portion of the simulation, to ensure that light entering the simulation is correctly modeled. The final stage in the physical damage model is the ablation stage. This occurs between 300°C and 500°C [50,51,43,47]. At the ablation temperature T a , we remove the remaining mass in the voxel and change the optical and thermal properties to that of the air.

Validation
We validated our numerical heat simulation against an analytical case. We assumed that the medium is uniform, with thermal properties all set to unity; this gives the medium thermal properties of a metal comparable with copper. The medium is initially at 37°C in an ambient medium of 0°C. We then allowed the simulation to evolve and compared our numerical model against the analytical case. This comparison is shown in Figure 1.

Two-dimensional and 3D
OCT scans of the MTZs as a result of CO 2 laser irradiation are shown in Figure 2. The holes are spaced out in a 9 × 9 grid. However, to ensure high resolution of the reconstructed images, the scan was  conducted over a limited field of view under rotation of the galvo-mirror. As a result, nine MTZs (three per sample) were analyzed for each laser energy. OCT was able to resolve structural information to a depth of approximately 500 μm, sufficient to characterize the MTZs. Figure 3 shows the increasing trend of MTZ depth with increasing laser energy (50-400 mJ/pixel). MTZ depth and surrounding crater lip were seen to increase visually with increasing energy. Visually, the MTZs take a less welldefined shape at very high energies (400 mJ/pixel). From the OCT images, four parameters are determined in order to characterize MTZs: depth of MTZ channel, width of MTZ channel, vertical expansion of crater lip, and horizontal expansion of crater lip (Fig. 2). In comparing the results obtained from OCT imaging and from our MCRT methods, the depth of MTZ channels is used as the comparative metric as crater lips are currently not modeled in the MCRT simulation.
As the porcine tissue was kept on ice prior to the experiment, we assumed that the temperature throughout the tissue was 5°C. Due to computational constraints, some approximations to the experimental setup have to be made. The porcine skin was a large thin slice of the topmost layers of the skin. However, the area of interest is around the MTZs. Therefore, in our numerical model, we only model the volume around one MTZ. To ensure this is a valid approach to take, i.e thermal damage does not overlap with other MTZs, one full simulation of the 81-pixel beams was carried out. Figure 4 shows the results of this simulation, and clearly shows that the MTZs do not interfere with one another such that only simulating one MTZ is a valid approach. We, therefore, modeled a volume of 0.06 cm × 0.06 cm × 0.18 cm, using 160 3 voxels. The depth of MTZ for each laser energy is shown in Figure 3 alongside depth values obtained from the numerical model. As the literature does not have a consensus on the exact temperature that the skin ablates at, we modeled over the range of temperatures found in the literature. As shown in Figure 3, the ablation temperature that best fits the model is approximately 500°C. We also simulated the thermal damage around the MTZs, with an example of the thermal damage around an MTZ, as shown in Figure 5. The thermal damage shown in Figure 5 falls within the previously found thermal damage [52].
As the numerical model is very flexible, we can also change parameters that may require different equipment or new experiments to be carried out. To illustrate this, we varied several laser and tissue parameters from the experiment. The first parameter we varied was the beam profile of the laser. As the laser has an unknown beam profile, we trialed a tophat and Gaussian beam profile. We found that the tophat beam profile best fits the experimental data. The results of this trial are shown in Figures 3 and 6. The laser beam profile is not likely to be a tophat beam, but some hybrid of Gaussian/tophat beam. The second parameter we varied was the temporal pulse profile of the laser. The Pixel CO 2 laser has a triangular pulse profile, so we trialed Gaussian and a tophat profile. The results are shown in Figure 7. The final parameter we varied was the initial temperature of the porcine tissue. Originally, we assumed that the tissue was at ≈ ∘ 5 C, so we evaluated the model for a range of temperatures from 0°C to 25°C with the results illustrated in Figure 8.

DISCUSSION
While clinical studies exist demonstrating and characterizing the effects of AFLs on the skin for enhanced drug delivery and improved therapeutic and cosmetic outcomes, a method of predicting laser response has until now has not been created. Our numerical model shows a good match to the experimental OCT data, based on MTZ depth when an ablation temperature of 500°C is adopted. This is within the parameters set out in the literature for human skin. However, as ex vivo porcine skin was used in this work, in vivo parameters for human skin may vary.
Previous work has shown that OCT is a suitable technique for imaging MTZs and has been validated against histological measurements [14]. The benefits of a noninvasive imaging modality are apparent, particularly in imaging MTZ formation in vivo. Dark shadows are observed below the MTZs in the OCT images in Figure 2. These are thought to be a result of the change in optical properties of the coagulated tissue [53] at the surfaces of the MTZ, perhaps reflecting or scattering more light away and resulting in a "blind spot" on the OCT image. However, it is shown in Figure 5 that nonablative levels of laser energy are deposited in these deeper tissues, resulting in heating, so while there are likely still biological effects taking place, OCT is perhaps unable to resolve any information in this area. MTZ depth has been shown to correlate with drug uptake [54,55]. Assuming that the MTZ is filled with the topical formulation of a chosen drug, it then affects a shorter drug-target distance [56]. Therefore, knowing where the target is located in the skin, it would be possible to predict the required energy density from a given AFL to create MTZs of sufficient depth, while avoiding excess laser exposure and hence excess side effects associated with AFLs. Increasing the MTZ density, by decreasing the spacing between MTZs, also increases drug uptake. While not the primary focus of this work, simulating all 81 MTZs from the laser output was shown in Figure 4, and therefore it would be perfectly feasible to further investigate effects associated with MTZ density using our approach.
It is possible to model the diffusion of drugs or other topical formulations through the skin [57]. Previous authors have also successfully combined drug diffusion models with MCRT simulations [58]. Thus, as a future avenue of work, we could combine a drug diffusion model with our predictive ablation model to quantify the effect that laser ablation has on drug diffusion.
From Figures 3 and 6 it is seen that using a Gaussian beam profile will result in a shallower MTZ depth than that of a tophat profile when the same energy densities are simulated. This is due to the distribution of power in the Gaussian beam compared with that of the circular beam. Additionally, our MCRT results suggest that a Gaussian temporal pulse profile may create shallower MTZs than a triangular or tophat profile, as shown in Figure 7. This can be explained by the fact that the Gaussian pulse delivers energy over a longer time period, thus letting more heat diffuse away from the MTZ, yielding shallower MTZs in comparison with a triangular or tophat temporal pulse. This is a useful information as it reinforces the need to look beyond simply what power settings a laser is used at, and demonstrates the power of our MCRT model in directly comparing specific variables.
Aside from simply determining the depth of the MTZs, this MCRT model offers the ability to investigate surrounding nonablative heat transfer in the tissue. This is evident in Figure 5. It is known that thermal stimulation of the dermis induces collagen production and hence the  appeal in the context of aesthetic dermatology [59,60]. Alongside understanding what temperatures the dermis may be subject to during irradiation with AFLs, specific temperature-dependent biological mechanisms may then be further probed.

CONCLUSIONS
We have demonstrated an MCRT model that is able to predict the depth of MTZ formation as a result of AFLs on ex vivo porcine skin, verified via OCT. The authors believe this is the first time the techniques used in this study have been used to predict laser tissue ablation depth. This paves the way for studies comparing the predictions made by the MCRT model and OCT imaging of in vivo human skin following AFL treatment. Beyond the depth of the MTZs with changing laser parameters, the model may also give information on the nondestructive heat distribution as a result of laser exposure. This predictive information is invaluable in terms of providing important information on temperaturedependent mechanisms in the surrounding tissues. Also, in the provision of information to help inform optimization of AFL treatment parameters to attain the desired clinical or aesthetic outcomes, while minimizing adverse effects.

DATA AVAILABILITY STATEMENT
All code and data to reproduce the results of this study are openly available. The code [61] and data [62] are published online.