Numerical optimization of sequential cryogen spray cooling and laser irradiation for improved therapy of port wine stain

Despite application of cryogen spray (CS) precooling, customary treatment of port wine stain (PWS) birthmarks with a single laser pulse does not result in complete lesion blanching for a majority of patients. One obvious reason is nonselective absorption by epidermal melanin, which limits the maximal safe radiant exposure. Another possible reason for treatment failure is screening of laser light within large PWS vessels, which prevents uniform heating of the entire vessel lumen. Our aim is to identify the parameters of sequential CS cooling and laser irradiation that will allow optimal photocoagulation of various PWS blood vessels with minimal risk of epidermal thermal damage.


INTRODUCTION
Treatment of cutaneous vascular lesions usually employs pulsed lasers with wavelengths of 532 nm or 585-595 nm, according to the principle of selective photothermolysis [1]. Introduction of cryogen spray (CS) precooling [2] has allowed significant protection against unwanted thermal damage to the epidermis, thus enabling application of higher radiant exposures and improved treatment efficacy. Nevertheless, the response of individual port wine stain (PWS) lesion remains highly variable and unpredictable, with many patients' lesions fading only minimally [3,4]. In addition to nonselective absorption of laser light by epidermal melanin, which limits the maximal safe radiant exposure, one plausible reason for poor treatment outcome is optical screening in large PWS vessels, due to the limited penetration depth of the strongly absorbed laser light [5]. This effect prevents uniform coagulation of the entire vessel lumen with single short laser pulses, allowing the vessels to recover from such partial damage [6,7].
One possible solution is to divide radiant exposure into multiple laser pulses (MLP) while controlling the epidermal temperature with multiple cryogen spurts (MCS) applied before each laser pulse. Our recent numerical and animal model study showed that such MCS-MLP treatment can safely provide photocoagulation of much larger PWS blood vessels as compared to the customary single cryogen spurt single laser pulse (SCS-SLP) approach, without adverse side effects [8]. Moreover, since no characterization of PWS is available to guide therapy on an individual patient basis, development of a treatment protocol enabling complete and uniform coagulation of vessels with various diameters and depths would be of great clinical advantage.
Herein, we apply a custom numerical model, encompassing optical and thermal transport in PWS as well as tissuespecific coagulation kinetics, to analyze the effects of laser pulse number (n) and repetition rate (f) on coagulation of PWS blood vessels of different diameters and depths in two darker skin phototypes. For each combination of irradiation sequence and target parameters, we determine the threshold radiant exposures for complete vessel coagulation and for onset of epidermal damage, respectively, when using a Nd:YAG/KTP laser (wavelength: l ¼ 532 nm) with 1 ms long individual pulses. The Nd:YAG/KTP laser is considered (instead of, e.g., a pulsed dye laser at l ¼ 585 nm), because it is capable of generating the rather high repetition rates required for effective MCS-MLP treatment. Our specific aim is to identify the MCS-MLP parameters that will allow more effective and safer photocoagulation of PWS blood vessels than the customary SCS-SLP approach across a wide range of PWS vessel sizes, depths, and skin phototypes. Such treatment protocol might present a viable approach to improved treatment of PWS, especially for patients with larger PWS blood vessels and/or darker skin.

METHODS
Our numerical model has three main components. A custom Monte Carlo (MC) model simulates light transport in PWS in three dimensions (3D), to predict the laser energy deposition maps. The finite-element model (FEM) simulates heat transfer during and following the MCS-MLP treatment, to compute the temperature field evolution in PWS. The extent of blood coagulation and dermal and epidermal thermal damage are computed using the customary Arrhenius model of tissue denaturation kinetics. Each model component is described below in more detail.

Optical Transport Model
Simulation of light transport in PWS skin is based on the weighted-photon MC technique [9]. In this model, a large number of energy packets (''photons'') propagate through the tissue and deposit a fraction of their energy into specific volume elements, according to local tissue absorption properties. Stochastic relations are used to simulate each photon's random trajectory, according to physical laws of light scattering, reflection, and refraction.
Our MC code is an extension of the basic MC light propagation technique [9], to allow treatment of arbitrary 3D tissue structures with different optical properties. Special care was devoted also to implementation of realistic boundary conditions. In earlier documented 3D MC models, either the photon escape [10] or reflective side boundaries [11] were applied. Clearly, neither approach is appropriate for treatment of a relatively small volume inside human skin with arbitrary inclusions and irradiated with a spatially confined laser beam. To overcome this problem, we propagate the photons exiting the finely discretized volume of interest (VOI) in laterally infinite tissue layers, until they return to the VOI, escape into the air, or lose all their energy. The same principle is applied also at the bottom of the VOI.
Our skin model consists of laterally infinite layers, representing the epidermis (thickness: 60 mm), dermis (1.5 mm), and a semi-infinite adipose layer underneath. In the dermis, we model one horizontal blood vessel with diameter varied from 30 to 150 mm, with an axis located 200, 400, or 600 mm below the skin surface (Fig. 1). The spatial discretization step is 2 mm Â 2 mm in the plane perpendicular to the vessel axis, and 20 mm along the vessel axis. The photon launching pattern simulates an incident laser beam with a top-hat profile (diameter: 5 mm), centered above the vessel axis.
In the model, melanin is homogeneously distributed throughout the epidermis. The epidermal absorption coefficient, m a,epi , is calculated using the following equations [12]: To account for optical shading by other, not explicitly modeled PWS vessels [5], the dermis is divided into five 200-mm thick layers and one layer with a thickness of 500 mm. The absorption coefficients for each layer (m a,i ) are calculated by considering the respective fractional blood volumes (f i ) and mean blood vessel sizes according to PWS histology data [14] (Table 1), as: Here, m a,der is the absorption coefficient of bloodless dermis (0.053 mm À1 ) [15] and C i is the correction factor to account for optical screening in a PWS blood vessel of radius r i [16,17]: For blood, we assume an oxygen saturation of 70%. We calculate the absorption coefficient by linear interpolation of the values for oxy-and deoxygenated red blood cells at 532 nm from Meinke et al. [18]. Since their values were reported for a 33.2% hematocrit, we have multiplied the interpolated value by 0.40/0.332 to obtain m a,bld ¼ 19.55 mm À1 for a hematocrit of 40%.
The absorption coefficient in the adipose is m a,adi ¼ 0.075 mm À1 [19].
A rather wide range of scattering coefficients (m s ) and anisotropy factor values (g) can be found in the literature for each human skin constituent [9][10][11][12][13][14][15][16][17][18][19]. We have tested a few combinations of the reported scattering parameters with absorption values as described above in MC simulations of both normal and PWS skin. Herein, we have selected the m s and g values that resulted in best agreement between the simulated diffuse reflectance and reported experimental values in both normal skin of different phototypes, and PWS (Table 2). (More details will be presented elsewhere.) The scattering coefficient in the adipose was calculated from the reduced scattering coefficient value reported by Bashkatov et al. [19] and the anisotropy g from Flock et al. [20].

Thermal Transport Model
Temperature field evolution in PWS skin is computed from two-dimensional (2D) heat-diffusion equation rc p @Tðx; z; tÞ @t ¼ kr 2 Tðx; z; tÞ þ Sðx; z; tÞ where r is density, c p is specific heat, and k is thermal conductivity. The heat-source term S(x,z,t) is considered only during the laser exposure, and assumes the spatial distribution of energy deposition in the central plane perpendicular to the vessel axis as obtained from the MC simulation. Evaluating the heat transport in 2D is legitimate because heat deposition within the VOI is quite homogeneous in the direction along the vessel axis and, therefore heat flow in that direction is negligible [27]. This greatly reduces the memory requirements and computation time for the thermal transport model. Heat exchange due to blood perfusion can be neglected during the relatively short time periods considered in this study (below 0.6 second) [28]. Moreover, since blood velocity in PWS venules is so low [29,30] and our VOI represents a small part of the irradiated volume, the small amount of blood leaving the VOI within the studied time is replaced by blood with similar temperature.
Since temperature in blood can transiently reach boiling conditions, the latent heat of water vaporization is accounted for by augmenting the specific heat term, that is, c(T) ¼ c p þ qD(T), where q is the latent heat of evaporation and D(T) is a normalized Gaussian function [31]: Assuming that structural non-homogeneity, compounded by rather high shear velocity and strong convection (upon nucleation boiling) will prevent significant overheating of blood, we set the function parameters to T B ¼ 1018C and DT ¼ 18C.
All thermal properties of skin components are listed in Table 3. We take into account the temperature dependence of heatcapacity forbloodless skin (ranging from3.1 to 3.9 kJ/ kg K) due to protein denaturation, as measured by differential scanning calorimetry [24].  [19] At the skin surface (z ¼ 0), we apply the convective boundary condition with a heat transfer coefficient h ¼ 5000 W/m 2 K and outer temperature of À508C during the CS period [32,33]. During the laser irradiation and after the cooling/irradiation sequence has ended, we apply a typical natural convection value (h ¼ 10 W/m 2 K) and ambient temperature (258C). The initial skin temperature is assumed to be uniform, at 358C, and adiabatic boundary conditions are applied at the lateral and bottom boundaries of the VOI.
Finally, Eq. 4 is solved with the finite element method (FEM), using a commercial software package (FEMLAB TM by COMSOL, Burlington, MA).

Thermal Damage Kinetics
Thermal damage induced in the modeled tissues is estimated by computing the thermal damage coefficient, V, according to the Arrhenius model: where R is the universal gas constant. Calculations of V are always performed for 100 ms after the last laser pulse to ensurethat the accumulation of thermal damage has stopped. The frequency factor (A) and activation energy (E a ) values for the included tissues are compiled in Table 3. For the epidermis and dermis, we use the values from Gaylor [23]. In a recent experimental study [35], these were reported to better match histological data than, for example, the bulk skin coefficients proposed earlier by Henriques [34]. Indeed, when using the latter in our initial simulation runs, the results featured extensive perivascular damage even at radiant exposures, which did not induce vessel coagulation. No such unrealistic dermal damage is obtained when applying the coefficients by Gaylor (Fig. 4).

Simulation Protocol
All simulated MCS-MLP sequences consist of an initial 50 ms CS spurt followed by 1 ms laser pulse, followed in turn by several pairs of CS (duration: 33-150 ms) and laser irradiation. For each combination of treatment sequence parameters (number of laser pulses: n ¼ 1-5; repetition rate: f ¼ 7-30 Hz), target vessel depth, and skin phototype, we determine the minimal radiant exposure per pulse at which all simulated PWS vessels are completely coagulated (H PWS ). Analogously, we determine the maximal radiant exposure per pulse at which no epidermal damage occurs (H epi ).
To determine these thresholds, we simulate the MCS-MLP sequence for each vessel diameter at a selected radiant exposure, H. Subsequently, we calculate V at the epidermal-dermal junction (z ¼ 60 mm; above the vessel center) and along the vertical line through the vessel axis. The radiant exposure is then increased in increments of 0.1 J/cm 2 and the ensuing V values evaluated, until the epidermal (criterion: V epi 1) and vascular damage thresholds (V PWS ! 1) are determined.
In addition, we estimate the extent of perivascular damage by evaluating V in dermis above and below the vessel axis, at a distance of 50 mm from the vessel wall. The vessel wall thickness is set to 10% of the vessel radius.
In addition to the skin heat capacity, several other optical and thermal properties of the involved tissues are known to change dynamically upon pulsed laser irradiation. Because  including all of them in our model would make the present parametric study prohibitively extensive, we present in the discussion an analysis of their impact on the study results.

RESULTS
Moderately Pigmented Skin Figure 2a showsamap of deposited energy obtained by MC simulation in a PWS model for moderately pigmented skin and a blood vessel (d ¼ 150 mm) located at depth z 0 ¼ 400 mm. The radiant exposure, H ¼ 5.4 J/cm 2 , matches the epidermal damage threshold for the case of single-pulse irradiation. The energy deposition within the vessel is nonuniform with the largest values near the top of the vessel, despite rather strong scattering of light in the dermis. The resulting temperature distribution in the same PWS vessel at the end of the 1 ms laser pulse preceded by 50 ms CS precooling is presented in Figure 2b. Temperature at the bottom of the vessel is approximately 608C, which is significantly lower than at the top ($1008C). Figure 3a shows temperature evolution at the epidermaldermal junction for a single laser pulse (n ¼ 1; left panel) and a MCS-MLP sequence (n ¼ 5, f ¼ 20 Hz; right). In each example, the radiant exposure per pulse is equal to the epidermal damage threshold, H epi (i.e., 5.4 and 3.7 J/cm 2 for the first and second case, respectively). The initial CS reduces the epidermal temperature to 108C. Thereafter, the single laser pulse heats the epidermis to the maximum temperature of 618C. In contrast, the maximum  temperature of 578C is reached with the MCS-MLP sequence, 48C less than with the single laser pulse.
In Figure 3b, we present the temperature evolution at the top, center, and bottom of the vessel with diameter d ¼ 30 mm located at depth z 0 ¼ 400 mm. For the single pulse, the maximum temperature of 1008C is obtained in the center of the vessel (thicker, black line). For the MCS-MLP, the maximum temperature at the vessel center gradually increases from 908C obtained after the first CS laser pulse pair to 988C at the end of the sequence. Significantly lower temperatures are obtained at the top and bottom of the vessel lumen, due to thermal relaxation during laser irradiation.
Temperature evolution for the vessel with diameter d ¼ 150 mm is presented in Figure 3c. Here, with the single laser pulse, the highest temperature is obtained at the top of the vessel (858C; thinner, red line), while reaching only 568C near the bottom (blue line). For the MCS-MLP sequence, the increase of temperature on the vessel axis (thicker, black line) by each subsequent pulse is more prominent as compared to the smaller vessel (Fig. 3b). The maximal temperatures are 93 and 728C at the top and bottom of the vessel, respectively. Hence, due to significant heat accumulation within the larger vessel, the MCS-MLP sequence results in higher and much more uniform temperatures as compared to SCS-SLP.
The extent of resulting coagulation is illustrated in Figure 4, showing the cross-sectional distributions of the damage parameter, V. Note that V ! 1 indicates blood coagulation or tissue necrosis, while the white and light gray areas (0.5 < V < 1) represent healthy and partly compromised tissues, respectively. We assume that coagulation of the entire vessel lumen is required to induce irreversible damage [7]. The 30 mm vessel is completely coagulated by both the single laser pulse (Fig. 4a) and MCS-MLP sequence (n ¼ 5, f ¼ 20 Hz ; Fig. 4b). In contrast, the 150 mm vessel is only partly coagulated by SCS-SLP (Fig. 4c), but the MCS-MLP sequence results in complete coagulation of the lumen (Fig. 4d). Since both radiant exposures are just below the respective epidermal damage thresholds, no epidermal thermal damage is evident in either case. The longer sequence causes some thermal damage to perivascular dermis, evident as a band of coagulated tissue extending $25 mm above the top of the vessel lumen. Figure 5 shows the effect of increasing the number of laser pulses (n) on tissue coagulation for vessels with diameters of 30 and 150 mm. At each n, the radiant exposure is set to the respective epidermal damage threshold, H epi . The smallest vessel (Fig. 5a) is completely coagulated for any number of laser pulses, and V at the characteristic points within the vessel does not vary significantly with n. The perivascular damage (diamonds) does not exceed 1.
In contrast, coagulation inside the largest vessel (Fig. 5b) depends strongly on n. Complete coagulation of the vessel lumen is obtained only with 3 or more laser pulses, while one and two pulses can not provide coagulation of the bottom part of the vessel. The perivascular damage increases with n and exceeds V ¼ 1 at n ! 3. Yet, the difference between V at the bottom of the vessel (down-triangles) and perivascular area (diamonds) increases with n. Therefore, perivascular damage could be avoided by applying MCS-MLP with suitably reduced radiant exposure.
Values of tissue damage for PWS vessels with different diameters are presented in Figure 6. For n ¼ 1 (Fig. 6a), damage at the bottom of the vessel decreases with increasing diameter, resulting in incomplete coagulation of the lumen for d ! 120 mm. For n ¼ 5 (Fig. 6b), the differences between the V values achieved in the smaller versus the larger vessels are much smaller as compared to the previous example, especially at the bottom of the vessel lumen (down-triangles). In addition, damage within larger vessels (d ! 100 mm) is more homogeneous as compared to n ¼ 1. By using the MCS-MLP with n ¼ 5, safe coagulation of all modeled vessels is achieved, albeit at the expense of some dermal damage (diamonds) around the larger vessels.
In Figure 7, minimal radiant exposures required for coagulation of all model PWS vessels at z 0 ¼ 400 mm (H PWS ; circles) and maximal permissible radiant exposures (H epi ;squares),arecomparedforrepetitionratesf ¼ 7-30 Hz. For single LP (n ¼ 1), complete coagulation of all PWS vessels can not be achieved, since H epi < H PWS . Using the MCS-MLP with three laser pulses enables successful coagulation of all simulated PWS (H epi > H PWS ). At n ¼ 4-5, the difference between H epi and H PWS becomes even larger. Therefore, the laser treatment is potentially safer, because the radiant exposure can be selected from a wider range of effective and safe radiant exposures.
In general, lower radiant exposures per pulse (H PWS ) are required to completely coagulate the PWS vessels with increasing n. The lowest H PWS is obtained at a 30 Hz repetition rate (Fig. 7, bottom panel). However, the epidermal damage threshold (H epi ) is reduced to almost the same value, because less heat can be extracted from the epidermal layer during the shorter CS periods. Therefore, such high repetition rates are not advantageous for therapy. In contrast, H epi does not vary significantly with increasing n at f ¼ 7 Hz (Fig. 7, top panel), indicating that the epidermis thermally relaxes between successive laser pulses.
The optimal repetition rates are evident from Figure 8. For MCS-MLP sequence with n ¼ 5, the difference between H epi and H PWS reaches 1.6 and 1.4 J/cm 2 for repetition rates of 10 and 14 Hz, respectively. At n ¼ 3 (blue symbols, dashed lines) the safety margin is somewhat smaller, but the optimal repetition rate remains around 10 Hz.   Figure 9 presents the threshold radiant exposures for shallower PWS vessels (z 0 ¼ 200 mm). These vessels can be coagulated while preserving the epidermal layer by all simulated MCS-MLP sequences. Nevertheless, the MCS-MLP sequences result in significantly larger differences between H epi and H PWS , thereby providing safer therapy than the SCS-SLP approach (Fig. 9a). The difference between H epi and H PWS is largest at lower repetition rates (Fig. 9b).
In contrast, none of the simulated MCS-MLP sequences can provide safe coagulation of all PWS vessels located at depth z 0 ¼ 600 mm (Fig. 10a). However, since using more laser pulses reduces the difference between H PWS and H epi , adding a few additional CS-LP pairs might enable successful treatment.
Considering the repetition rate, f ¼ 10 Hz provides the smallest difference between the two thresholds at n ¼ 5 (Fig. 10b). As Figure 11 demonstrates, such MCS-MLP sequence can completely coagulate all vessels at depth z 0 ¼ 600 mm, except the bottom part of the largest vessel (d ¼ 150 mm). Some perivascular damage is predicted for vessels with diameters above $100 mm (V > 1, diamonds).

Darkly Pigmented Skin
Threshold radiant exposures for PWS vessels in darkly pigmented skin at depth of z 0 ¼ 200 mm are presented in Figure 12. At f ¼ 14 Hz (Fig. 12a), safe coagulation of all PWS blood vessels is predicted for n ¼ 4 and 5. However, the difference between H epi and H PWS is much smaller than in moderately pigmented skin, only $0.2 J/cm 2 . The optimal repetition rates are around f ¼ 14 Hz (Fig. 12b).
At greater depth (z 0 ¼ 400 mm; Fig. 13), only PWS vessels with d ¼ 30-70 mm are completely coagulated at n ¼ 5 and f ¼ 10 Hz, while the bottom part of larger vessels (d ! 100 mm) remains uncoagulated. Nevertheless, the difference between H epi and H PWS is reduced by increasing the number of laser pulses from 1 to 5, and is smallest for repetition rates of 10-14 Hz (at n ¼ 5; not presented). Both trends are very similar to the results for z 0 ¼ 600 mm in moderately pigmented skin (Fig. 10a,b, respectively).

DISCUSSION
A recent numerical study [8] indicated that a MCS-MLP sequence with two laser pulses at repetition rate of 20 Hz is more effective for coagulation of larger PWS vessels (d > 90 mm) as compared to the customary SCS-SLP approach.
In this study, we analyze the effects of laser pulse number (n) and repetition rate (f) on coagulation of various PWS blood vessels, in search of the optimal MCS-MLP parameters. In order to isolate the effects of sequential CS cooling and laser pulse irradiation, a single irradiation wavelength is considered for all cases. A more elaborate optical model of a PWS lesion, with several dermal sublayers according to the blood vessel size and depth distribution from PWS histology data [14,17], is used to account for optical shading [5]. By updating the tissue optical parameters based on most recent reports [15,18,19], the resulting optical model correctly predicts diffuse reflectance values for skin of different phototypes [27].
For moderately pigmented skin, our simulation results indicate that the customary SCS-SLP approach is ineffective for coagulation of larger PWS vessels located at z 0 ¼ 400 mm or deeper. Yet, MCS-MLP sequences exist, where complete coagulation of PWS vessels of all diameters can be obtained with no adverse side effects (Figs. 5-8). Only the largest vessels at depth z 0 ¼ 600 mm can not be coagulated by 5 laser pulses (Fig. 11). Nevertheless, the trend indicated in Figure 10a suggests that complete coagulation of 150 mm vessels at such depth might be possible by using somewhat longer MCS-MLP sequences. Note, however, that blood vessels of such large diameter are rarely present in typical PWS [14,36]. Moreover, considering the pronounced scattering of visible light in the epidermis and dermis, and abundance of PWS vessels between 200-400 mm [14,36], the effect of such deep vessels on visual appearance of PWS might be rather small.
Very similar trends were found also for darkly pigmented skin. Complete coagulation of all PWS vessels at z 0 ¼ 200 mm is obtained with MCS-MLP sequences at n ¼ 4 and 5, and the optimal repetition rate is around f ¼ 14 Hz.  At z 0 ¼ 400 mm, only vessels of smaller caliber (d < 100 mm) can be coagulated. However, since such vessels prevail in a typical PWS lesion [14,36], the MCS-MLP approach appears promising for treatment of patients with darker skin.
Parameters of the applicable and near-optimal MCS-MLP sequences according to this study are collected in Table 4. It is evident that the clinically most demanding cases included in our study (i.e., involving PWS vessels with d ¼ 150 mm at z 0 ¼ 600 mm in moderately pigmented, and z 0 ¼ 400 mm in darkly pigmented skin) require the longest MCS-MLP sequences (n > 5) and allow the narrowest range of applicable repetition rates. Note, however, that PWS with vessels of such diameter are indicated as untreatable by SCS-SLP for all examples, except for the very shallowest depth in moderately pigmented skin (Figs. 7, 10a, and 12a). This is in good agreement with clinical experience, if the differences in pulse duration, wavelength, and cooling technology between specific vascular laser systems are considered.
Perhaps the most important finding, however, is that the optimal MCS-MLP parameters (f ¼ 10 Hz, n ¼ 5) indicated for the two most demanding examples (stated above) lie within the therapeutically effective ranges for the remaining, less demanding cases. Moreover, ultimate optimization of the sequence for the latter examples appears unnecessary, since the safety margin (indicated by the difference H epi À H PWS ) is comfortably large and does not vary strongly with small changes of f or n. Simply, the MCS-MLP sequence with f around 10 Hz and n ¼ 5 not only works best for the most demanding PWS cases included in our analysis, but works very well also for all less demanding examples (involving lighter skin phototypes, shallower PWS vessels, and/or smaller vessel diameters).
In fact, such MCS-MLP treatment appears superior to the customary SCS-SLP approach (n ¼ 1) for every PWS geometry under study, as long as vessels with d > 50 mm are involved, and efficiently equalizes the temperature rise and consequent coagulative effect achieved in the smallest versus the larger PWS vessels (Figs. 3-6). Consequently, MCS-MLP may turn out to be a treatment modality that alleviates the need for optimization of therapy on an individual patient basis-an elusive goal in absence of suitable lesion characterization tools.
One word of caution is in order, however, with respect to the predicted perivascular damage at larger PWS vessel diameters (Fig. 6). This effect increases with pulse number n (Figs. 6b, 11, and 13), since increasingly more heat accumulates within the target vessels and subsequently diffuses into the dermis. The clinical implications of such perivascular damage are unclear at this point and require further study.
In present study, we have taken a conservative approach, applying mostly temperature-independent optical and thermal properties of the involved tissues to avoid uncontrolled biasing of the results by potentially inaccurate modeling of the involved dynamical processes. Therefore, we analyze below in one example, central to the study, the sensitivity of our results to inclusion of most relevant dynamic variations of thermal (in addition to those already included: c p (T) for skin; plasma boiling) and optical properties in the involved tissues.
First, we estimate the temperature dependence of thermal conductivity k(T) in blood and bloodless skin by considering the data for water [37] and accounting for the respective protein contents. Using these dynamic values (and all other tissue parameters kept the same as in the core study) we assess the epidermal damage threshold for a MCS-MLP sequence with f ¼ 10 Hz and n ¼ 5, that is, H epi ¼ 4.9 J/cm 2 . The corresponding axial profile after the end of the MCS-MLP sequence (Fig. 14, dashed blue line) shows a $3% lower temperature rise within the vessel as compared to that obtained using the constant k (thin red line), especially in the lower part of the lumen.
In the subsequent step, we consider also conversion of oxyhemoglobin to met-hemoglobin, which occurs nearly simultaneously with hemoglobin denaturation [39]. To model this effect, we analyze the distribution of thermal damage inside the vessel just before each laser pulse in the MCS-MLP sequence, and modify the absorption coefficient at all positions where blood coagulation threshold has been exceeded (V ! 1). The respective absorption coefficient is computed by replacing all oxyhemoglobin with met-hemoglobin [40] while keeping (unconverted) deoxyhemoglobin. This yields m a,coag ¼ 14.1 mm À1 , which is 28% below the initial value for blood. In this way, an updated energy deposition map is obtained from a unique MC simulation run for each laser pulse in the MCS-MLP sequence, accounting for the extent of coagulation caused by preceding laser pulses.
The axial temperature profile obtained by including both k(T) and met-hemoglobin formation is presented in Figure 14 (heavy black line). In comparison with the profile that considers only k(T) (dashed blue line), the additional effect of met-hemoglobin formation is very small. Only a minor increase of temperature at the vessel center and in epidermal layer can be observed.
Using the latter approach, we have computed also threshold radiant exposures H PWS and H epi as a function of repetition rate for n ¼ 5 (Fig. 15). The main difference from the original result (Fig. 8) is a general increase of both threshold values, due to more effective thermal relaxation of both the target vessel and epidermis between consecutive laser pulses. The two effects are quantitatively different, however, resulting in an increased therapeutic safety margin at all tested frequencies. The corrected threshold values for the SCS-SLP approach are H epi ¼ 5.2 J/cm 2 (red open square at f ¼ 0) and H PWS ¼ 7.4 J/cm 2 (not plotted), almost identical to those obtained using the original model (Fig. 7).
Besides the increased thresholds, however, the indicated trends are the same as obtained with the former approach (Fig. 8). Frequencies between f ¼ 7 Hz and 20 Hz appear very safe, with the optimum indicated around 10 Hz-same value as before. This suggests that the optimal MCS-MLP frequencies are dictated primarily by thermal relaxation of the target vessels and epidermis, rather than specifics of optimal transport. Note that our approximation of a constant heat transfer coefficient h is not utterly unrealistic (e.g., experimental data in Fig. 8 of ref. [32,38]). Moreover, due to inherent properties of heat diffusion, the cooling effect at the epidermal basal layer and target vessel depth will depend primarily on the average rate of heat extraction at the skin surface, rather than specifics of its dynamics. Nevertheless, we present also the results for one CSC setup with a highly variable h(t) (Fig. 4 in ref. 38; Fig. 15, smaller red symbols). Evidently, the lower average heat extraction rate as compared to the former example has resulted in diminished H epi values, with a minor effect on H PWS . The optimal range of MCS-MLP frequencies is, however, only marginally affected, with the safety margin being the same from 10 to 20 Hz. The slight shift of the optimal range toward higher frequencies can be attributed to the fact that the applied h(t) peaks after t $ 20 ms, so the average h is largest for CS durations around 50 ms.
Because the 3D MC runs are computationally very intensive, the discussed augmentation of our optical/thermal model has slowed down the simulations by a large factor. At the available processing power and modeling software, accounting for dynamic changes of tissue properties throughout the parametric study would be prohibitive in terms of the required computation time. Even so, the described augmentation of the model ignores dynamic changes occurring during the individual 1 ms laser pulses. A more rigorous approach would require dividing the laser pulse into multiple shorter intervals and performing the above described steps for each of these intervals. Only then could one treat, for example, the bathochromic shift of the hemoglobin absorption line [39], which is significant only during the laser pulse. Nevertheless, we expect that its impact would be even smaller than that of met-hemoglobin formation, because the associated decrease of m a,bl at 532 nm at a 708C temperature increase is below 6%. We therefore believe that the present model encompasses most relevant physical phenomena to a degree where the indicated optimal MCS-MLP sequence parameters are relevant. Preliminary clinical tests involving MCS-MLP treatment at n and f values similar to the values presented in this study significantly improved clearing of PWS lesions with no adverse side effects.

CONCLUSIONS
The results of our numerical simulation indicate that MCS-MLP sequences exist, where complete and safe coagulation of PWS vessels is achieved in several cases, where the customary single laser pulse approach is not effective. Optimal sequences involve MLP (n ¼ 3-5) at repetition rates of 10-15 Hz. The presented approach could improve therapeutic outcome for patients with resistant PWS lesions, especially those involving very large blood vessels and darker skin phototypes.