2.1.

Optical Transport Model

Our simulation of light transport in PWS skin is based on the well-established weighted-photon MC technique.¹ In this approach, a large number of energy packets (“photons”) propagate through the tissue and deposit a fraction of their energy into specific volume elements, according to the local absorption properties. Stochastic relations are used to determine each photon's random trajectory, according to physical laws of light scattering, reflection, and refraction. Specifics of these relations depend on optical properties, i.e., absorption coefficient (μ _a), scattering coefficient (μ _s), anisotropy (g), and index of refraction (n) of the involved tissues. With an increasing number of simulated photons and ever finer spatial discretization of the treated volume, the resulting energy deposition map ideally converges toward the correct spatial distribution of the absorbed light energy.

Extending the multilayer MC model¹ to three dimensions allows treatment of skin inclusions of arbitrary complex shape and optical heterogeneity. In the presented examples, the VOI has dimensions of 1 × 1 mm² in the plane parallel to the skin surface (coordinates x and y) and 1.5 mm in depth (z). The spatial discretization step is 2 × 2 μm² in the (x, z) plane and 20 μm in the y direction, resulting in a 3D grid with 500×50×750 voxels.

A special care is devoted to realistic treatment of the VOI side boundaries. As indicated above, neither mirror, nor escape BC is suitable for high-resolution treatment of relatively small volumes in human skin, irradiated with spatially confined laser beams. In our approach, named extended layers, photons are launched according to the laser beam intensity profile and propagated in laterally infinite tissue layers until they deposit all their energy or escape into the air. Optical properties of these layers are inherited from the border VOI elements at the respective depths, and the majority rules whenever not all of them represent the same tissue. Perhaps most importantly, energy deposition is computed and recorded only within the VOI, which keeps memory requirements at a manageable level and greatly reduces the computation time.

The same principle is applied also at the bottom boundary of the VOI. Specifically, the entire semi-space below the VOI bottom boundary inherits the optical properties from the lowest VOI layer, which is occupied by adipose in all presented examples.

The MC code is implemented in Visual C++ (Visual Studio 2008, Microsoft, Redmond, Washington). Utilization of multiple threads helps speed up the computation.

2.2.

Skin Model Geometry and Optical Parameters

Our basic skin model is composed of three layers, which represent the epidermis (thickness: 60 μm), dermis (1.5 mm), and a semi-infinite layer of subcutis.

Melanin is assumed to be homogeneously distributed within the epidermis. The effective fractional content of melanin (f _mel) is calculated from absorption coefficient at 694 nm (μ _{a ,epi}) as determined from diffuse reflectance measurements by Svaasand ⁷ To this end, we use the well-known relations:⁸

The resulting melanin fractional contents for fair Caucasian, moderately pigmented Middle Eastern or Asian, and African skin are f _mel = 1.3%, 3.5%, and 9%, respectively. The corresponding epidermal absorption coefficients at 532, 585, and 595 nm, obtained by using the wavelength dependence relations (1), are collected in Table 1.

Table 1

Epidermal absorption coefficients at three common laser wavelengths, representative of three skin phototypes.

For blood, we assume a hematocrit of 40% and oxygen saturation of 70%. The absorption coefficients at three considered wavelengths are obtained by linear interpolation of the values for oxygenated and deoxygenated whole blood as reported by Meinke ⁹ (see Table 2). The absorption coefficients of subcutis (adipose) are from Bashkatov ¹⁰

Table 2

Absorption coefficient (μ a), scattering coefficient (μ s), anisotropy coefficient (g), and refractive index (n) values used in our MC model of optical transport in skin.

See Table 1.

A wide range of scattering parameter values for constituents of human skin can be found in the literature. We have tested several combinations of the reported values in MC simulations of both normal (assumed 1% blood in dermis) and PWS skin. For this study, we have selected those values of μ _s and g from recent reports, which resulted in good agreement between the predicted diffuse reflectance and reported experimental values. The value of μ _s for adipose was calculated from the reduced scattering coefficient ( [TeX:] $\mu_{s}^{\prime}$ ${\mu }_s^{\prime }$ ) as measured by Baskhatov,¹⁰ and the anisotropy coefficient (g) from Flock ¹¹

For our analysis of blood photocoagulation in laser therapy of PWS, we include one blood vessel with diameter d = 150 μm and a horizontal axis 400 μm below the skin surface (Fig. 1).

Fig. 1

Cross-sectional view of our PWS model geometry, including the epidermal layer (dark gray), six dermal layers with different absorption properties (gray), one horizontal blood vessel (white), and a semi-infinite adipose layer underneath. Dashed lines delineate the VOI.

In order to account for optical shading by other, not explicitly modeled PWS vessels,¹² the dermis is divided into five 0.2-mm thick layers and one layer with a thickness of 0.5 mm. The corresponding absorption coefficients μ _a,i (Table 3) are calculated by considering the reported fractional contents of blood (f _i) in each layer:¹³

Table 3

Mean blood vessel radii (r i), fractional blood contents (f i) (Ref. 11), and the resulting absorption coefficient values (μ a ,i) at λ = 532 nm for the six model PWS layers.

2.3.

Heat Transport and Thermal Damage

Temperature field evolution in PWS skin is computed using the two-dimensional heat diffusion equation

Two-dimensional treatment of heat transport in the central plane perpendicular to the vessel axis (y = 0.5 mm) will be sufficient because—as we will demonstrate below—deposition of laser energy in the discussed examples is nearly homogeneous along the y direction (i.e., parallel to the vessel axis). Consequently, heat flow in this direction is negligible.

Heat exchange due to tissue perfusion can be safely neglected during the short time periods (below 0.2 s) considered in this study.¹⁷ In addition, due to the low blood velocity in the ectatic PWS venules^{18, 19} and their tortuous shape, only a small fraction of heated blood can leave the VOI within the studied time interval. Moreover, since our VOI represents only a small part of the irradiated volume, this blood will be replaced by blood with a similarly elevated temperature. For these reasons, we neglect heat loss due to blood flow in the present demonstration example.

Thermal properties of most skin constituents are the same as used in our recent report.²⁰ This includes the temperature dependent specific heat of bloodless skin,²¹ and the way we account for the latent heat of water vaporization (in agreement with other similar studies).^{22, 23} In the present study, we additionally take into account the temperature dependence of tissue thermal conductivity, k(T), which was found to have an observable effect on predicted temperature dynamics, and consequently the extent of epidermal, blood, and perivascular thermal damage.²⁰ This temperature dependence is estimated by considering the data for water²⁴ and accounting for the protein content in blood and other skin layers.

For the duration of CS cooling, we apply a convective boundary condition at the skin surface, with a heat transfer coefficient H = 5000 W/m²K and outer temperature of −50 °C.^{25, 26} At all other times, and for the examples without active cooling, we apply a typical natural convection value (H = 10 W/m²K) and ambient temperature (25°C). Adiabatic boundary condition is applied at the lateral and bottom boundaries of the VOI and the initial skin temperature is set to 35°C.

Equation 3 is solved using the finite element method within the software package FEMLAB (COMSOL, California).

The laser-induced thermal damage is assessed within the customary Arrhenius model of protein denaturation kinetics, by computing the thermal damage coefficient,

All calculations cover an interval of 100 ms after the laser pulse, to ensure that the value of Ω has stopped increasing.

3.1.

Diffuse Reflectivity

As a first test of our MC optical model, we compute diffuse reflectance of healthy skin (with 1 vol.% of blood distributed evenly throughout the dermis) at 532, 585, and 595 nm. The epidermal layer thickness is varied from 60 to 80 and 100 μm, which are typical values for human skin,^{32, 33} in order to obtain a realistic range of diffuse reflectance (DR) values.

Figure 2 presents the obtained DR values for the three epidermal thicknesses (see the legend) in different skin phototypes. The DR spectra as measured on forearms of healthy volunteers by Svaasand ⁷ are presented for comparison (dashed lines).

Fig. 2

Diffuse reflectance values at 532, 585, and 595 nm for healthy skin of: (a) fair Caucasian (open symbols), (b) moderately pigmented Asian (gray), and African phototype (black). The epidermal thickness is set to 60, 80, or 100 μm (see the legend). Dashed lines indicate the corresponding diffuse reflectance spectra in vivo (Ref. 5).

For fair Caucasian skin [Fig. 2a], our diffuse reflectance values (open symbols) fall completely within the range of the DR spectra measured on the volar and dorsal side of the forearm, respectively (dashed lines). In Asian skin [Fig. 2b, gray symbols], the predicted range of DR values matches the reported spectrum rather well. The predicted values at λ = 532 nm are a little higher, but are matched with the experimental report by Randeberg (R = 0.20).³⁴ An excellent match between the predicted DR values and experimental spectrum in vivo is observed for African skin, specifically for epidermal thickness of 80 μm (black squares).

Our model results illustrate how an increase in epidermal thickness leads to lower diffuse reflectivity of skin, due to increased melanin absorption. Moreover, this effect is enhanced with increasing melanin content (f _m). Consequently, we can observe a substantial change of DR with the epidermal thickness in Asian and African skin phototypes. In contrast, the effect is rather small, and probably often not discernable from experimental data, in fair Caucasian skin.

At all three phototypes, however, the largest influence of epidermal thickness on DR is predicted at the longest wavelength, 595 nm.

3.2.

Influence of Boundary Condition

Before embarking on simulations of PWS treatment, we analyze the influence of BC implemented at the sides of the VOI. For this analysis, we augment our model of healthy skin (fair Caucasian; f _m = 1.3%), as presented above, with one cylindrical blood vessel in the dermal layer. The vessel (diameter: 120 μm) lies horizontally under the laser beam axis, at a subsurface depth of 0.25 or 0.50 mm.

In addition, we consider the case where both vessels are present simultaneously. This enables quantitative assessment of the influence of one vessel on energy uptake in another vessel in close proximity (a.k.a. optical shading). The same geometry was used in one of the first studies of this effect, as well as of optical screening within individual blood vessels, using 3D MC modeling.² For even closer correspondence, we also set the epidermal thickness to 60 μm, laser beam diameter to 1 mm, and radiant exposure to 1 J/cm². For each model geometry, we calculate energy deposition maps upon irradiation at 585 nm when applying the mirror or escape BC, and compare them with our extended layer (EL) approach.

The resulting energy deposition profiles along the vertical beam axis (passing through both vessels’ axes) are presented in Fig. 3. Clearly, the choice of BC is reflected primarily in the result for the deeper vessel (z ₀ = 0.50 mm). For this vessel, using the mirror BC [fig. 3a] results in a 73% larger energy deposition (integrated over the entire vessel lumen) as compared to the EL result [Fig. 3c; see Table 4]. The escape BC [Fig. 3b], on the other hand, results in 16% lower energy deposition in the same vessel than the EL approach.

Fig. 3

Energy deposition profiles through the axes of two horizontal vessels located at subsurface depth of 0.25 and 0.5 mm (solid black line), compared with the result for a single vessel located at 0.25 mm (thinner, blue), or at 0.5 mm (dashed, red). The boundary condition applied at the sides of the VOI is: (a) mirror, (b) escape, and (c) extended layers.

Table 4

Energy deposition per unit length within a single vessel (z 0 = 0.25 or 0.50 mm) and with both vessels present simultaneously (see Fig. 3). Results for three boundary conditions are presented (λ = 585 nm).

For the shallower vessel (z ₀ = 0.25 mm), the influence of BC is somewhat less pronounced. Application of the mirror BC results in an overestimation of the energy uptake by 32%, while the escape BC leads to a 7% lower energy deposition estimate as compared to the EL result.

The energy deposition results in Table 4 also allow quantitative assessment of optical shading, i.e., the relative change of energy uptake by the presence of a nearby vessel. We can see that the influence of the deeper vessel on the laser energy deposition in the shallower one—albeit not very large—varies strongly with the choice of BC. Specifically, by applying the mirror BC at the side boundaries of the VOI, a reduction of deposited energy by 6.3% is obtained upon introduction of the deeper vessel. This prediction, however, is more than twice as large as the result obtained with the more realistic EL approach (i.e., −2.7%). Applying the escape BC, in contrast, yields to underestimation of the shading effect by ∼30%.

The optical shading effect in the deeper vessel, while obviously significantly larger in magnitude (−28%), is affected by selection of either mirror or escape BC to a significantly lesser degree.

Figure 4 presents energy deposition profiles through a single vessel (z ₀ = 0.25 mm) with diameter of 60 or 120 μm, upon laser irradiation at 532, 585, and 595 nm, when the EL approach is applied. At the wavelength of 595 nm, with the lowest μ _a value in blood, the near-central part of the larger vessel (lighter; red solid line) absorbs slightly more laser energy than at 532 and 585 nm (thicker; black and thinner; blue, respectively). This is not the case for the 60 μm vessel (dashed lines), which is almost optically thin (μ_a d < 1) at the considered wavelengths.

Fig. 4

Energy deposition profiles through the axis of a single vessel with diameter of 120 μm (solid lines) or 60 μm (dashed), irradiated at 532 nm (thicker black line), 585 nm (thinner; blue), and 595 nm (lighter; red). All results were obtained using the extended layer approach.

Energy depositions within the vessel lumen from the result in Fig. 4 are compared with predictions obtained when using the mirror and escape BC in Table 5. The largest discrepancy between the mirror BC and EL approach is observed at λ = 595 nm. The mirror BC yields 56% and 48% larger energy depositions in 60 and 120 μm vessels, respectively, as compared to the EL results. The escape BC, in contrast, leads to 11% lower predictions than the EL approach at both vessel diameters. At 532 and 585 nm, similar trends but somewhat smaller relative differences can be observed.

Table 5

Total energy uptake by vessels of 60 and 120 μm diameter located at 0.25 mm depth, upon irradiation at 532, 585, and 595 nm. Results for all three boundary conditions are presented.

3.3.

Deposition of Laser Energy in a Model PWS Blood Vessel

Figure 5 illustrates spatial distribution of deposited energy in a model PWS vessel, as obtained by our MC simulation upon irradiation with a 5-mm diameter hat-top laser beam at λ = 532 nm in moderately pigmented Middle Eastern or Asian skin (f _m = 3.5%). The radiant exposure is H = 7.6 J/cm², which is just below the epidermal damage threshold (see Secs. 3.5 and 3.6). In the central plane perpendicular to the vessel axis [y = 0.5 mm; Fig. 5a], the energy deposition within the vessel is nonuniform, with the largest values near the top of the vessel lumen. Figure 5b presents the energy deposition in the vertical plane through the vessel axis.

Fig. 5

Cross-sectional maps of laser energy deposited in a model PWS blood vessel (d = 150 μm, z0 = 400 μm) upon irradiation with a 5-mm diameter hat-top laser beam (l = 532 nm) and radiant exposure H = 7.6 J/cm² in the central planes: (a) perpendicular to the vessel axis (y = 0.5 mm) and (b) through the vessel axis (x = 0.5 mm).

Figure 6a presents the corresponding lateral energy deposition profile through the vessel axis (y = z = 0.5 mm). The result indicates a dip in perivascular dermis, caused by reduction of light fluence due to strong absorption in the blood vessel. Figure 6b presents the corresponding energy deposition profiles in the direction parallel to the vessel axis (x = 0.5 mm). Energy deposition along the vessel axis (solid black line) is practically constant throughout the VOI, due to a comparatively larger laser beam diameter (d _l = 5 mm). Energy deposition along the top of the PWS vessel (lighter, blue line) is significantly larger and shows a slight decrease toward the VOI borders. At the epidermal–dermal (ED) junction, energy deposition is again nearly constant, and larger than in the vessel center (dashed line).

Fig. 6

Horizontal energy deposition profiles from the result in Fig. 5: (a) through the vessel axis (y = z = 0.5 mm) and (b) along the vessel axis (solid black line), along the top of the vessel lumen (thinner, blue line), and at the epidermal-dermal junction (z = 60 μm; dashed).

The lower noise level in Figs. 5a, 6a, as compared to Figs. 5b, 6b results primarily from the fact that the former correspond to a 10 times thicker tissue slice than the latter (i.e., Δy = 20 μm versus Δx = 2 μm). In addition, we have averaged 5 centrally located slices to further suppress photon noise in Fig. 5a and obtain a smooth heat source term for Eq. 4. Obviously, the same could not be done in the case of Fig. 5b, leading to rather noisy appearance of the profiles in Fig. 6b. This is not an issue, however, since the latter serves only as an illustration of the model behavior and is not used in the subsequent modeling steps.

The central depth profile of energy deposition is very similar to that shown in Fig. 4 (for λ = 532 nm, d = 120 μm), and is therefore not presented.

3.4.

Temperature and Coagulation Maps in Laser-Irradiated PWS Blood Vessel

Figure 7a presents temperature evolution in the model PWS blood vessel from Fig. 5, irradiated with a 1 ms laser pulse at 532 nm after 50 ms long CS cooling, as obtained using our thermal transport model 3. At radiant exposure of H _epi = 7.6 J/cm², just below the epidermal damage threshold, the temperature at the top of the vessel lumen (thinner, blue line) reaches 100°C, but remains below 70°C at the vessel bottom (dashed line).

Fig. 7

Temperature evolution in laser irradiated PWS: (a) at the top (thinner, blue line), center (black), and bottom (short dashed) of the vessel lumen; and (b) at the epidermal-dermal junction. Radiant exposure is H = 7.6 J/cm²; see text for other details.

Temperature evolution at the ED junction, meanwhile, shows an initial drop to 10°C upon CS precooling, followed by a jump to the maximal temperature of 68°C [Fig. 7b].

Figure 8a presents temperature distribution in the central cross-sectional plane of the model PWS geometry at the end of 1 ms laser pulse irradiation, preceded by 50 ms CS cooling. The temperature at the bottom of the vessel lumen is considerably lower than at the top.

Fig. 8

(a) Temperature field in the central cross-section of the model PWS geometry (y = 0.5 mm) immediately after pulsed laser irradiation preceded by CS cooling. (b) Temperature depth profiles through the vessel center obtained with (solid line) and without CS precooling (dashed). The radiant exposure values for CS cooling (7.6 J/cm²) and no cooling (4.6 J/cm²) match the respective epidermal damage thresholds, H _epi.

The corresponding temperature depth profile through the vessel center is presented in Fig. 8b (solid line) and compared with the result obtained for the case without CS precooling (dashed). Note, however, that these profiles correspond to different radiant exposure values, selected to match the respective epidermal damage thresholds (as determined further below): H _epi = 7.6 J/cm² with CS cooling versus 4.6 J/cm² without cooling. We can see that the two radiant exposure values indeed result in comparable maximal temperatures within the epidermal layer (T ≈ 80 °C). At the same time, CS cooling allows significantly larger temperatures to be achieved in the lower part of the vessel lumen, as compared to irradiation without precooling.

The extent of tissue coagulation, resulting from either treatment approach, is compared in Fig. 9 by presenting the cross-sectional distributions of the damage parameter Ω 4. Coagulation of blood and dermis is indicated by darker gray shades (Ω ≥ 1), while white and magenta areas indicate healthy and partly compromised tissue, respectively. In moderately pigmented Asian skin (f _m = 3.5%), used in this model example, the model PWS vessel can be only partly coagulated at the maximal safe radiant exposure allowed for laser treatment without CS cooling [Fig. 9a]. The CS cooling approach [Fig. 9b], in contrast, allows safe application of significantly larger radiant exposure. This enables almost complete coagulation of the model PWS blood vessel, except for the very bottom part of the vessel lumen (white and magenta regions).

Fig. 9

Cross-sectional maps of the coagulation parameter Ω in PWS model geometry (f _m = 3.5%): (a) immediately after irradiation with 1 ms laser pulse after natural convection cooling; (b) after irradiation with 1 ms laser pulse preceded by CS cooling. The radiant exposure in either example is just below the epidermal damage threshold, i.e., 4.6 J/cm² and 7.6 J/cm² for (a) and (b), respectively.

3.5.

Radiant Exposure Thresholds

For three laser wavelengths (λ = 532, 585, and 595 nm), two skin phototypes (fair Caucasian and moderately pigmented Asian), and two cooling approaches (no cooling and CS cooling), we determine the minimal radiant exposure, at which the entire lumen of the model PWS vessel is coagulated (H _PWS), and the epidermal damage threshold, H _epi, as described before. Provided that H _epi is larger than H _PWS, the interval between the two thresholds indicates the range of safe and therapeutically effective radiant exposure values for the selected skin and treatment conditions.

In Fig. 10a, we compare the obtained radiant exposure thresholds for fair Caucasian skin. Clearly, a complete coagulation of the model PWS vessel can be safely achieved at all tested conditions (see the labels).

Fig. 10

Numerically predicted threshold radiant exposures HPWS and H _epi in (a) fair Caucasian skin (f _m = 1.3%); and (b) moderately pigmented Asian skin (f _m = 3.5%) for laser irradiation at 532, 585, and 595 nm, without (left) and with application of CS cooling (right; note the labels).

Specifically, the smallest safety margin is predicted for PWS treatment at 532 nm without precooling of skin (H _epi – H _PWS = 2.1 J/cm²). By addition of CS cooling, the epidermal damage threshold (H _epi) is markedly increased with almost no effect on H _PWS, thus increasing the safety margin to a much more comfortable value of 7.5 J/cm². These results agree very well with clinical experience. A good fraction of PWS lesions in fair skin could be treated with no adverse side effects using early pulsed 532 nm lasers without CS cooling. However, introduction of the latter has dramatically improved the safety record of PWS therapy. Moreover, by allowing safe use of significantly higher radiant exposures, it has enabled treatment of many “resistant” PWS lesions, featuring deeper and/or larger blood vessels.

At the two longer laser wavelengths (585 and 595 nm), with lower absorption in epidermal melanin as compared to 532 nm, the difference H _epi – H _PWS is considerably larger. For our model geometry, the predicted range of safe yet therapeutically effective radiant exposure values thus amounts to ∼6 J/cm² without cooling, and raises to ∼13 J/cm² with application of 50 ms CS cooling. This is also in good agreement with clinical experience, which made such pulsed dye lasers (PDLs) the preferred treatment modality for the majority of PWS patients.

In contrast with the above, our results predict that complete coagulation of the model PWS vessel in moderately pigmented Asian skin cannot be obtained safely without application of CS cooling, as H _epi < H _PWS at all tested wavelengths [Fig. 10b; left part]. With CS cooling (right), however, safe treatment of PWS appears possible at 585 and 595 nm, but with only a marginal difference between H _epi and H _PWS (1.8 to 2.0 J/cm²). From clinical practice, it is well known that treatment of PWS lesions in patients with moderately and darkly pigmented skin is very challenging, and often impossible, even with up-to-date PDL lasers equipped with CS cooling devices.

A more detailed comparison of our predicted therapeutic thresholds with radiant exposures from several clinical reports can be found in Sec. 4.