A Three-Dimensional Modular Adaptable Grid Numerical Model for Light Propagation During Laser Irradiation of Skin Tissue

— Information regarding energy deposition during laser irradiation of structurally complex biological tissue is needed to understand and improve the results of clinical procedures. A modular adaptive geometry numerical model capable of simulating the propagation of laser light in a wide variety of multiple component tissues has been developed and tested. A material grid array is generated by assigning a value representing a tissue type to each of a large number of small voxels. The grid array is used to indicate optical properties in an existing variable step size, weighted-photon Monte Carlo algorithm that has been modiﬁed to account for voxels-to-voxels changes in optical properties. To test the model, simple geometric shapes and optical low coherence reﬂectometry images of rat skin have been used to create material grids consisting of epidermis, dermis, and blood. The model assumes 1-J/cm 2 irradiation of the tissue samples with a 1.0-mm diameter uniform beam at 585 nm. Computed results show good quantitative and qualitative agreement with published data. Various effects due to shading and scattering, similar to those suggested in the literature, are noted. This model provides a way to achieve more realistic representation of anatomical geometry as compared to other models, and produces accurate results. laser tissue interaction, thermal mechanisms of pulsed laser ablation, and transmyocardial laser revascularization.

A N UMERICAL modeling of laser tissue interaction can be effective tools for understanding and improving laser treatments. Models enable the reseacher to visualize basic physical mechanisms-such as light absorption, temperature rise, and thermal damage-and to make computational "experiments" by varying laser parameters, tissue geometry, and thermal boundary conditions.
Modeling approaches to the general problems of light transport and heat transfer have been researched extensively. Techniques such as the Monte Carlo method [1] for optical modeling and finite difference or finite element methods [2] for thermal diffusion have been shown to be quite effective.
However, the accuracy of laser tissue interaction numerical models is often called into question because they do not take into account the contribution of factors such as complex tissue geometry, water loss, thermal and mechanical damage, and dynamic changes in optical and thermal properties.
The distribution of light in turbid media was first described using the radiative transport equation [3] which cannot be solved analytically for most tissue geometries. Researchers have used techniques such as the diffusion approximation [4]- [6] and Beer's law [7], which are valid for highly scattering and nonscattering materials.
The Monte Carlo technique, an excellent overview of which is presented by Jacques and Wang [1], is a statistical method for simulating light transport that is valid for a wider variety of situations. Using experimentally determined absorption and scattering coefficients, Henyey-Greenstein anisotropy factors, indexes of refraction, and several physical and probabilistic relations, the photon "random walk" can be simulated repeatedly until a sufficient quantity of photons has been absorbed to establish reliable trends of energy deposition.
Previous research has shown that the Monte Carlo method can be used in different ways to effectively simulate laser irradiation of skin using homogeneous [8] and layered geometries [9]- [11]. Several studies have used the Monte Carlo technique to simulate treatment of port wine stain (PWS), combining simple geometric shapes to approximate the components of skin [8], [12]- [14]. Lucassen et al. [14] used layers, long cylinders, and long "curved vessels" to simulate the microstructure of skin. They reproduced some of the optical phenomena seen in laser-tissue interaction, including shading caused by more superficialvessels and how diffuse light fills partially shaded regions. Using layers and various sizes of straight, cylindrical vessels in a Monte Carlo model, Smithies et al. [13] demonstrated that changes in diameter and spatial distribution of blood vessels can affect energy deposition in skin more than a change in wavelength from 577 to 585 nm.
A two-dimensional (2-D) optical-thermal laser tissue interaction model has been developed by Kim [15] and Rastegar et al. [16]. The model iterates between Monte Carlo optical component and a finite-element thermal component to account for local changes in absorption and scattering coefficients due to temperature rise and thermal damage.
We have developed a three-dimensional (3-D) modular adaptive grid numerical model (MAGNUM) for light prop-agation in geometrically complex biological tissue. Although a modular representation of nonhomogeneous tissue would be useful in the simulation of any number of laser procedures, one obvious application is PWS treatment. Both clinical [17] and numerical studies [8], [13], [14] have shown that the specific vessel geometry plays a major role in determining the outcome of PWS treatment. Recent research has indicated that appropriate laser parameters depend on the unique 3-D mapping of PWS microstructure [18], [19].
This study has two goals: to digitally represent structurally complex biological tissue in a way that can be easily adapted to different geometries and to augment the standard Monte Carlo algorithm to handle three-dimensional (3-D), heterogeneous geometry.

II. METHODS
A complete, yet concise description of algorithms added to a standard Monte Carlo model are covered in detail, while the basic Monte Carlo techniques are summarized.

A. Material Grid Generation
The key to the augmentation is a material grid which corresponds in size one-to-one with the grid that accumulates absorbed photon weights. In this grid, locations can be specified using Cartesian coordinates or the indexes of a rectangular voxel . If all voxels are specified as cubes, a simple routine can be used to convert from Cartesian location to voxel index. The dimensions of the voxels and and the number of elements in each direction and are used to establish the grid skeleton. The material grid is represented by an array in which each voxel is assigned an integer value corresponding to a material (in the case of skin: epidermis, dermis, or blood vessel).
Material grids were specified using two different techniques. The first set of grids was taken directly from geometry specified in the literature. These included a layered geometry (simulation 1) and a combination of layers and straight cylinders (simulations 2-5). The second set (simulations 6 and 7) was based on an in vivo blood vessel pair described by Barton et al. [20]. Subsurface high-resolution images of rat skin were acquired using optical low-coherence reflectometry (OLCR). The blood vessel outlines were manually traced from the images, and an epidermal layer of constant thickness was estimated. During the image-to-material grid conversion, irregularities in the curvature of the vessels were smoothed. In simulation 7, the vessels pair was duplicated, rotated approximately 90 , and placed below the original pair to illustrate asymmetric shading effects.

B. Photon Movement
An example of photon movement in two dimensions ( and axes) is presented in Fig. 1. In the following description, references to Fig. 1 are enclosed in brackets.
The algorithm begins with the standard weighted photon "launch," which accounts for the beam profile and specular reflection. Initial directional cosines and are specified to represent light normal to the surface of the tissue.

Path lengths
are calculated using the local absorption in cm and scattering coefficients in cm and a random number. Using the present location and directional cosines and the location within homogeneous tissue of the next potential attenuation event is found The most significant augmentation of the Monte Carlo routine involved the movement of photons through elements and across their boundaries. Essentially, this required extending the algorithm used for layered tissue into three dimensions. Movement within any voxel follows the standard algorithm for homogeneous tissue [segment AB]. When the size of the voxels is small with respect to penetration depth , where , it becomes unlikely that movement of distance will be confined within a voxel. In the case that the element of origin and potential end point voxel are different [as in segment BE], each point at which the photon exits, one voxel and enters another must be found. Based on the direction of the photon, three of the six planes defining a voxel can be eliminated as possible exit locations. The distances to each of the remaining planes and the directional cosines are used to calculate the partial path lengths that would have to be traveled to reach the planes. In the direction, for example, if if photon will not exit voxel along direction where is the center of voxel and is the original photon location. and are found by replacing with and , respectively, in the above equations. and are shown as dashed lines; corresponds to segment length BC; corresponds to segment length BD; if the photon is traveling only in the plane, and . The shortest of these partial path lengths [ , segment length BC] defines the location of photon exit from the voxel . If, for example, was the shortest The assignment of the location exactly on the voxel boundary ( above) is made instead of using the equation in order to avoid problems due to potential rounding errors. The originally computed path length [BE] is then reduced by the distance traveled within element [distance traveled in voxel 1, segment length BC] If the new path length [segment CE] indicates that another voxel boundary must be crossed, the photon is propagated from one voxel boundary to the next [segment CD]; this sequence is repeated until the photon can be propagated from a voxel boundary to the location of the attenuation event [segment DE, assuming and ] or reaches the boundary of the tissue.
Upon entering a voxel with a different attenuation coefficient is rescaled according to the ratio of the attenuation coefficients (The effect of rescaling is seen in the difference in length between segments DF and DG, assuming ). The  [14] initially computed optical distance (in dimensionless units of optical depth) is equal to the sum of the optical distances traveled in each of the voxels Changes in the index of refraction are accounted for differently at tissue-tissue and air-tissue interfaces. If a photon encounters an index of refraction change between tissue types (interface between voxels 2 and 3 at point D) and the angle of the photon with the incident surface is larger than the critical angle, the photon is reflected. At angles less than the critical, the Fresnel equation and a randomly generated number are used to determine if the photon is reflected (segment DG) or transmitted through the interface and is refracted (segment DF). If the surface is parallel to the plane, for example, refraction alters the directional cosines as follows: where is calculated from Snell's Law: , and subscripts and indicate the incident and transmitted sides of the interface. At the air-tissue interface, the critical angle determines total internal reflection, whereas the Fresnel equation is used to determine the fraction of the photon weight that escapes. The remaining weight is reflected back into the tissue.
If a photon travels the entire path length without escaping, the standard Monte Carlo absorption, Henyey-Greenstein scattering, and roulette termination procedures are used. Each absorption is recorded and used to calculate the local energy deposition (J/cm ).

C. Simulations of Skin with Blood Vessels
In all of the vessel simulations (2-7), several basic parameters remain the same. Cubic voxels are used. The size of the tissue samples being modeled is 1.  Table I. The index of refraction of air ( 1.0) was used at the boundaries. Simulations were performed using 10 000 000 photons.

III. RESULTS
The seven simulations detailed in this section have been divided into two parts. Simulations 1-5 use highly simplified geometries (one with layers only, four with straight vessels) similar to that found in the literature. Simulations 6 and 7 employ more complex, anatomically realistic OLCR imagederived geometry. An overview of the parameters used in the six simulations using vessel geometry (simulations 2-7) is given in Table II.

A. Simplified Geometry Simulations
Simulation 1 involves a semi-infinite (in plane) threelayer slab configuration specified by Chan et al. [21], consisting of 1.  Table III, along with corresponding data experimental and modeled from Chan et al. [21].
The geometries for simulations 2-5 use one or two straight, cylindrical blood vessels running parallel to the axis, directly below the laser beam. The goal of these simulations is to assess MAGNUM's accuracy by analyzing the effect of a change in grid resolution and comparing MAGNUM's results with simulations from the literature. Since symmetric geometries are used in simulations 2-5, only half the tissue slice is shown in the surface plots.
The geometry for simulation 2 (Fig. 2) involves a 120m diameter vessel located at a depth of 250 m. Since the size of the voxels (20 m on a side) that specify the material grid are on the same order as the cylinder diameter, there is noticeable distortion from the desired circular shape. The distribution of energy deposition in the plane at m is shown in Fig. 3. The energy deposited in the column of cubes directly below the center of the beam ( 700 m) is illustrated in Fig. 4. The material grid cross section shown in Fig. 5 shows the geometry for simulations 3-5. Epidermal and dermal layers and two vertically aligned cylindrical vessels of the same size as in simulation 2 (120m diameter) are represented using cubes measuring 10 m on a side. Smaller voxel size results    Due to increased resolution, energy deposition results from simulation 5 (Fig. 6) show a more realistic distribution in the blood vessel than seen in simulation 2 (Fig. 3). The statistical variation in the dermis has increased due to a decrease in the photon deposition per voxel. The surface plot of simulation 5 results (Fig. 6) shows no effect of direct shading on the deeper vessel, but less energy is deposited than in the superficial vessel. Energy deposition results for simulations 3-5 are shown as a function of tissue depth below the beam center in Fig. 7. To minimize the statistical variation of the graphs in Fig. 7, for each data set the two columns of voxels on either side of the laser beam center (see thick line in Fig. 5) were averaged. Fig. 7(a) uses a linear scale to provide a direct comparison with data from the literature and to show the effect of shading on energy deposition in blood vessels. Fig. 7(b) uses a logarithmic scale to illustrate the effects of shading on the dermis and, indirectly, the scattering-induced recovery of the fluence.

B. Simulations Using OLCR Image-Based Geometry
The geometries for the final two simulations were adapted from OLCR images of blood vessels in rat skin. The material grid top view in Fig. 8(a) and cross section in Fig. 8(b) show the four vessels used in simulations 6 and 7. The two sets of roughly perpendicular vessel pairs were used in simulation 7, whereas only the superficial pair was used in simulation 6. Each vessel pair is comprised of a 90m diameter vessel and a 190-m diameter vessel (the diameters of the superficial and deep vessels are specified in the and planes, respectively). The superficial and deeper vessel pairs are centered at depths of 365 mand 615 m, respectively.
The results for simulation 6 were very similar to that of the superficial vessel pair in simulation 7 and are, therefore, not presented. This data can be inferred from Figs. 9 and 10.
Simulation 7 used the most complex geometry, involving two layers of OLCR-specified vessel pairs at different depths, with the primary intent of investigating optical mechanisms with more realistic, complex vessel geometry. Fig. 9 shows a

IV. DISCUSSION
There are three main issues to be addressed in evaluating MAGNUM: the effect of modularization of the material grid on the representation of vessel geometry and the results, the quality of the data insofar as the data compares to that found in the literature, and computational considerations.

A. Reproduction of Tissue Geometry
A main goal of this research is to maximize accuracy when developing representations of 3-D tissue geometries. Since MAGNUM relies on cubic voxels to build tissue structures, it cannot specify rounded edges as well as can be done using geometric shapes such as cylinders. To adequately represent the curvature of a region, voxel size must be much smaller than the region or the true geometry will be lost. Reducing the voxel size brought about a significantly better representation of the curvature of the vessel. Higher resolution also led to a more accurate, continuous distribution of absorbed energy within the vessel, although statistical variation increased (most noticeably in the dermis) due to a factor of eight decrease in the average number of photons absorbed in each voxel. In choosing a voxel yvsize, both the resolution needed to properly represent structures and the statistical variation due to small voxel size must be considered in order to find an appropriate balance. If computation time and memory size limits the minimum voxel size to 10 m, then the cross section of a 10-m vessel would be represented by a single voxel. Note that the smallest vessels are 10 m or larger. Choice of grid parameters is a basic issue of modeling which will be discussed further in the section on computational considerations.

B. Comparison with Modeled Results from the Literature
Simulations 1-5 give good evidence of MAGNUM's ability to simulate the basic optical mechanisms of absorption, scattering, reflection, and refraction. The transmission and diffuse reflection results produced by MAGNUM for the simplified case in simulation 1 were within 1% of data calculated using a well-established Monte Carlo model by Wang and Jacques [9] and within 2% of experimental values found by Chan et al. [21].
The results from simulations 2-5 can be compared to studies that used similar geometries with simple straight cylindrical vessels [13], [14], [22]. Due to differences in methodology, laser parameters, and inability to read specific data points from contour plots, it is difficult to make direct comparisons in most cases. Therefore, trends in energy deposition in the epidermis, dermis, and blood vessels will be compared.
Results from MAGNUM and those from two other studies [13], [14] exhibit numerous common features that can be attributed to absorption, scattering, and/or shading effects. One example is the energy deposition in the epidermis, which decreases steadily by a few percent over its 60m thickness, with its magnitude comparable to the peak in deeper vessels in simulations 4 and 5, whereas in superficial vessels from simulations 2, 3, and 5, the peak is about 3 to 4 times greater than the level in the epidermis. Additionally, the epidermal layer tends to be affected only minimally by adding vessels in the dermis.
In these models, distribution of energy absorbed in the vessels is seen to be affected by several mechanisms. The primary influence is the large of blood at 585 nm, which is an order of magnitude greater than the of epidermis. Because of the short penetration depth ( 15 m) in blood at 585 nm, a large quantity of energy is absorbed at the outer, top edge of the vessels and a radial gradient (sharp, exponential decrease toward the center) in absorbed energy is seen. The light directionality also plays an important role. If a vessel is directly below the beam center, the gradients are strongest in the direction (depth), whereas for a vessel located to one side of the beam, the gradients are directed at an angle (see Figs. 9 and 10). Smaller vessels exhibit a more constant distribution of absorbed energy over their cross section (see Fig. 10) because the number of photons that reache the middle of the vessel is proportional to the distance they must travel through highly absorbing blood.
Scattering also leads to features observed in all of the simulations. Because of scattering, collimated laser light becomes diffuse. Thus, photons also impact a vessel from the bottom and sides, causing a peak along the peripheral walls of the vessels. Below a highly absorbing vessel, scattering in the dermis quickly diminishes the effect of direct shading [ Fig. 7(b)]. Although this leads to a distributed decrease in the magnitude of energy deposition in the tissue, gradients remain unaffected except in a small region directly behind the vessel. The only instance of vessel shading is seen in Figs. 9 and 10, in which the larger of the deeper vessels shows a significant decrease in energy deposition over a region directly below the larger superficial vessel.
Probably the most interesting comparison can be made with Fig. 7(a) from this paper and Fig. 4(b) from Lucassen et al. [14]. The graphs show agreement within a few percent for the distribution of energy deposition in the epidermis for all three cases and in the vessels for the single superficial and single deep vessel cases. However, according to MAGNUM the energy deposition in the superficial blood vessel for the two vessel case is only slightly reduced from the single vessel case. Data published by Lucassen et al. [14] indicates that the addition of a deeper vessel significantly affects energy deposition in the more superficial vessel. Therefore, in all but one instance MAGNUM's data is in good quantitative and qualitative agreement with the literature.

C. Computational Effectiveness
MAGNUM was developed as a step toward a complex type of laser tissue interaction model-one which requires a great deal of computational power. In order to represent a sizable volume of tissue in three dimensions with high resolution and record the distribution of absorbed photons, it is necessary to use large material and energy deposition arrays. In the most extensive simulations performed (3)(4)(5)(6)(7), this involved matrices with dimensions of 140 140 100, for a total of nearly 2 million elements per array. The propagation of 10 000 000 photons required approximately 1000-1500 minutes of processor time on an IBM RS/6000 workstation.
When considering MAGNUM's computational effectiveness, the relationship between data resolution and noise level, inputs such as grid and voxel size, number of photons, and computational considerations such as memory and processor time, must be considered. Increasing the resolution of a model imposes greater demands on memory size. In the change from simulation 2 to 3, the number of elements in the material and absorption matrices went from 245 000 to 1 960 000, each. With increased resolution came increased statistical variation since there were fewer photons deposited per bin. In order to increase resolution and maintain accuracy, a proportional increase in the number of photons is needed, and, consequently, more processor time. Therefore, when setting input parameters relating to grid size and number of photons, it is worthwhile to consider the goal of the simulation and the limitations of the computer system. This study has presented energy deposition data for various tissue geometries. A thorough analysis of PWS laser treatment will require the addition of a finite difference or finite element method component for computation of heat conduction. The energy deposition (J/cm ) from the Monte Carlo model is used as the source term in thermal models that compute the transient temperature response. Such a model should have the capability of accounting for nonhomogeneous tissue, various types of boundary conditions, the cooling effect of perfusion, and changes in perfusion due to thermal damage [23], as well as dynamic changes in optical and thermal properties. Thermal damage can be predicted by including a temperature, time-dependent rate process [24].
It should be noted that accurate communication of four (or more) dimensions of data in a 2-D black and white (or color) format is difficult, if not impossible. The ability to view and graphically manipulate complex, multidimensional data sets would greatly enhance one's comprehension of the information. Until such a format is available to facilitate mass communication of data among researchers, the scientific community as a whole will remain at a significant disadvantage.

V. CONCLUSION
This augmented Monte Carlo technique will allow better representation of real biological tissues and lead to improved simulation of laser surgery procedures. Although the concept of modularization is not new, this is the first time, to the authors' knowledge, that it has been applied to optical modeling of complex, 3-D, anatomical structures. Good approximation of skin tissue with PWS can be accomplished, although care must be taken to use the proper resolution when specifying small, curved vessels. Benchmarking against data from the literature has helped to establish that MAGNUM can provide qualitatively and quantitatively reliable energy deposition information.
The cost of this improvement of the Monte Carlo technique is an increase in computational requirements. Access to a workstation of moderate speed and memory size is highly recommended when using MAGNUM with a moderate-to high-resolution grid. Future improvements in the model will most likely necessitate even greater computational power.
In order to raise MAGNUM to a level at which meaningful simulation is possible, two advancements are planned. First, real tissue geometries will be reconstructed by converting digitized pictures of histology or confocal and coherencedomain system images into material grids on a pixel-bypixel basis. Second, a thermal component that incorporates the capabilities mentioned in Section IV will be added to MAGNUM.