1.

Introduction

Success of laser based therapeutic procedures depends critically on delivery of sufficient light to the target volume to induce the desired effect. A diagnostic measurement of the distribution of laser energy deposition and ensuing temperature rise can provide the clinician with necessary knowledge to customize treatment parameters and thus optimize the efficacy of laser treatment.

In our laboratory, we are interested in improving laser therapy of port wine stain (PWS) lesions, which are congenital vascular malformations present typically on the face and neck in 0.7 of the population.¹ Histopathological studies of PWS show a normal epidermis overlying an abnormal plexus of dilated blood vessels located in the dermis. Epidermal thickness (50–150 μm) and PWS blood vessel diameter (30–300 μm) and depth distribution (150–1000 μm) vary on an individual patient basis and even between different areas on the same patient.² Current therapy involves use of pulsed dye laser light (λ=585–595 nm), which is absorbed by hemoglobin constituents in blood and epidermal melanin, in conjunction with cryogen spray cooling.³ The goal of therapy is to deliver a sufficient quantity of light to destroy permanently PWS blood vessels and use cryogen spray cooling to minimize laser induced thermal injury to the epidermis. Clinical studies show promising results;⁴ however, multiple treatments are required, and complete blanching of the PWS lesion is achieved in only ∼10 to 20 of patients.^{5 6 7 8} A major limitation of current laser therapy is that treatment plans are based primarily on subjective criteria such as PWS color and physician experience. Numerous studies^{9 10 11 12 13 14} indicate that knowledge of specific PWS skin parameters, such as epidermal thickness (z _epi ), melanin absorption, and PWS vessel size and depth (z _PWS ), could result in optimization of laser therapy on an individual patient basis. For example, cryogen spray cooling parameters could be optimized with knowledge of z _epi and z _PWS . ¹⁴ Also, the maximum radiant exposure that can be used while safely avoiding epidermal thermal damage is dependent on ΔT _epi,max . Finally, selection of optimal laser wavelength may depend on both z _PWS and z _PWS,max . ^{15 16}

Pulsed photothermal radiometry (PPTR) involves measurement of time resolved changes in blackbody emission after irradiation of an object with a laser pulse.¹⁷ Superficial light absorption results in an immediate change in emission. Heat generated by light absorption in deeper structures will affect emission only after it diffuses towards the object surface. Previous studies detail capabilities and limitations of PPTR for depth profiling of chromophores in optically scattering media.^{18 19 20 21} PPTR has been employed for depth profiling of PWS skin.^{18 22} Since PWS blood vessels may be located just below the epidermal basal layer (∼60–100 μm below the skin surface),² high spatial resolution is required.

In this paper, we report on the development of a one-dimensional (1-D) numerical model for simulation of optical-thermal response of skin to pulsed dye laser irradiation. Information extracted from histological sections of PWS skin^{23 24} is incorporated in the model. An iterative, non-negatively constrained conjugate gradient algorithm¹⁸ is applied to computed PPTR signals to obtain depth profiles of the simulated PWS skin. In our analysis, reconstructed profiles are presented using a single PPTR signal (single wavelength excitation, SWE) or PPTR signals obtained using a dual wavelength excitation (DWE) approach.^{25 26} Comparison of the actual and reconstructed profiles from PPTR data revealed a good match. Results of this study indicate that PPTR is a viable approach for depth profiling of PWS skin.

2.

Optical-Thermal Model

To evaluate PPTR, a 1-D optical-thermal model incorporating PWS skin structure from a digitized biopsy²³ was developed. An overview of the overall model scheme is presented in Fig. 1.

Figure 1

Flow chart of the optical-thermal modeling technique using a biopsy-defined PWS skin geometry. Rectangles represent input or output, and ellipses represent data processing steps. Histology sections were first converted from 2-D cross-sectional slices to 1-D absorption profiles. A Monte Carlo model was applied to simulate light transport and energy deposition in each section. This procedure was repeated for each of the 42 histology sections used in this study. The 42 light distributions were converted to a corresponding number of initial temperature profiles immediately after pulsed laser irradiation. An average heat source profile was calculated and used as input in an explicit finite difference model to simulate ensuing heat transfer dynamics and compute a PPTR signal.

2.2.

Conversion of 2D Histology Section to 1D Profile

With custom software written in Matlab™ (Version 6.1, The MathWorks, Natick, MA), each row in each two-dimensional (2-D) section was classified as one of the following (Fig. 2):

Figure 2

Example of 2-D histology section divided into different layers. In this example, the gray region is the epidermis, the black regions are cross sections of blood vessels, and the lightly shaded region is bloodless dermis. Each of the 229 rows in the image was classified into one of eight types (see figure legend). Adjacent rows of the same type were grouped into layers. Layers without blood vessels were analyzed row by row. Layers with blood vessels were analyzed as a single layer. In this example, two blood vessel layers (“h”) are present. Note that only six of the eight layer types listed in the legend are present in the presented example. Scale bar=50 μm=25 rows.

a. Air

b. Air/epidermis

c. Air/epidermis/dermis

d. Epidermis

e. Epidermis/dermis

f. Epidermis/dermis/blood

g. Dermis

h. Dermis/blood

Next, an absorption coefficient μ_a was assigned to each bloodless row as a weighted sum of epidermal and dermal μ_a values. The epidermal absorption coefficient (μ_a,epi ) was chosen to represent fair human skin. Within the reported range of values,^{8 27} μ_a,epi =13 cm ⁻¹ was selected at 585 nm, representative of fair-skinned human subjects.^{25 26} An epidermal melanin fraction f _mel associated with the selected μ_a,epi value at 585 nm was determined to compute μ_a,epi at 600 nm using the following equations:²⁷

Eq. (1a)

$${\mu }_{a,\text{epi}}={\mathit{f}}_{\text{mel}}{\mu }_{a,\text{mel}}+(1-{\mathit{f}}_{\text{mel}}){\mu }_{a,\text{baseline}},$$

Eq. (1b)

$${\mu }_{a,\text{mel}}=6.6\times {10}^{11}{\text{cm}}^{-1}{\left(\frac{\lambda }{\text{nm}}\right)}^{-3.33},$$

Eq. (1c)

$${\mu }_{a,\text{baseline}}=0.244{\text{cm}}^{-1}+(85.3{\text{cm}}^{-1})\times \text{exp}[-\frac{(\text{λ}-154\text{nm})}{66.2\text{nm}}],$$

where μ_a,mel is the melanosome absorption coefficient, μ_a,baseline is the average absorption coefficient of other skin structures, and λ is the laser wavelength (nm). In the present study, f _mel was ∼3, resulting in μ_a,epi =12 cm ⁻¹ at 600 nm (see Table 1).

Table 1

Tissue absorption coefficients μa used in the Monte Carlo optical model, for 585 and 600 nm light. Blood μa values represent 75 oxygen saturation (see Ref. 29).
Tissue	585 nm (cm−1)	600 nm (cm−1)
Epidermis	13	12
Dermis	0.37	0.37
Blood	177	32

Adjacent rows containing blood were grouped together to create layers containing entire blood vessels (Fig. 2); μ_a of each of these layers was calculated as

Eq. (2)

$${\mu }_{a,\text{layer}}={\mathit{f}}_{\text{epi}}{\mu }_{a,\text{epi}}+{\mathit{f}}_{\text{der}}{\mu }_{a,\text{der}}+\sum _i{\mathit{f}}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}{\mu }_{a,\text{bld}},$$

where f_x is the fractional area occupied by component “x”; and subscripts “epi,” “der,” and “bld” denote epidermis, dermis, and blood, respectively. For example, for a given row, f _epi was determined by counting the number of pixels representing epidermal tissue and then dividing this number by the total number of pixels in the same row. The term f_i is the area fractional content of vessel “i” in the layer, and C_i is a correction factor for blood vessel i. C_i was introduced by Verkruysse et al.;²⁸ to account for optical screening in PWS blood vessels. Due to high absorption of 585 and 600 nm light in blood, red blood cells towards the vessel center interact less with incident light than those at the periphery. This effect becomes more prominent as vessel diameter increases. C_i effectively reduces the blood fraction of each vessel, allowing the assumption of a homogeneous blood distribution in each layer. This simplifies computations yet yields a correct fluence distribution and energy deposition profile. The following equation was found to fit data from Verkruysse et al.;:²⁸

Eq. (3)

$$C_i=0.03895+0.48588exp\left(-\frac{{\mu }_{a,bld}{r}_i}{0.1928}\right)+0.46803exp\left(-\frac{{\mu }_{a,bld}{r}_i}{0.914 43}\right).$$

The effective radius r_i of the ith blood vessel was computed from the associated cross-sectional area (A_i) as (A_i/π)^1/2. Blood μ_a values represent 75 oxygen saturation²⁹ (Table 1).

Since the axial extent of each 2-D histology section was only 458 μm or less, additional PWS skin structure was added to the model to increase the total simulated skin depth to 2 mm. This step was necessary to minimize deep-edge artifacts and ensure that realistic light and temperature distributions were computed. The μ_a values at depths below the bottom of each slice were computed using Eqs. (1 2 3 4) and mean blood vessel sizes and blood fractions provided in Barsky et al.;² (Table 2).

Table 2

Mean blood vessel diameters d avg , blood fractions f bld , and absorption coefficient (μa) values at 585 and 600 nm for additional skin layers used to increase model geometry depth. Initial depth of the most superficial additional skin layer (e.g., 0.450 mm) was adjusted accordingly to correspond to the bottom of each 2-D biopsy-defined slice.
Depth (mm)	d avg  (mm)	f bld	μa,585 (cm −1)	μa,600 (cm −1)
0.450–0.699	0.037	0.080	4.4	2.2
0.700–0.899	0.034	0.038	2.4	1.3
0.900–1.099	0.027	0.025	2	1
1.100–2.000	0.024	0.020	1.7	0.88

3.

Results

Figure 3 presents a comparison between laser induced (λ=585 nm) temperature profiles, computed from the model using two sets of epidermal optical properties: (1) μ_a,epi =13 cm ⁻¹ and μ_s,epi =200 cm ⁻¹ (solid line); and (2) μ_a,epi =9 cm ⁻¹ and μ_s,epi =500 cm ⁻¹ (dashed line). In general, the two profiles were very similar in shape and demonstrate the complex structure of PWS skin. With higher μ_s,epi , the epidermal temperature rise was slightly higher, but PWS temperature profiles were nearly identical in shape and amplitude. Total diffuse reflectance values were 19 and 24 for μ_s,epi =200 and 500 cm⁻¹, respectively. A diffuse reflectance value of 19 was comparable to values (14–17) measured from fair-skinned PWS patients.³⁷ Furthermore, use of a higher μ_s,epi value of 500 cm⁻¹ (at a correspondingly reduced μ_a,epi =9 cm ⁻¹) had a minimal effect on the computed PWS temperature profile, suggesting that the ambiguity in μ_s,epi values^{8 16 32 33} does not detract from the interpretation of our modeling results. Due to these observations, we used μ_s,epi =200 cm ⁻¹ for the remainder of the study.

Figure 3

Comparison of computed initial temperature profiles at 585 nm using two different sets of epidermal optical properties: μ_a=13 cm ⁻¹, μ_s=200 cm ⁻¹ (solid line); and μ_a=9 cm ⁻¹, μ_s=500 cm ⁻¹ (dashed line). Epidermal profiles differed slightly, but PWS profiles were nearly identical in shape and amplitude.

Comparison of fluence profiles at 585 and 600 nm reveals that the fluence at 600 nm was consistently higher than that at 585 nm (Fig. 4). This trend was due to the substantially lower blood absorption of 600 nm light, resulting in deeper light penetration but also more backscattered light. The increased quantity of backscattered light resulted in a diffuse skin reflectance (30), considerably higher than at 585 nm (19).

Figure 4

Fluence profiles at 585 and 600 nm. The μ_s,epi =200 cm ⁻¹ was used in the calculation. The fluence was considerably higher at 600 nm due to the substantially lower blood absorption coefficient, resulting in a higher quantity of backscattered and deeply penetrating light.

Comparison of initial temperature profiles after irradiation at 585 and 600 nm illustrates the effect of a large difference in μ_a,bld between the two wavelengths [Fig. 5(a)]. Epidermal temperature rise at 600 nm (dashed line) is slightly higher than at 585 nm (solid line), despite the lower melanin absorption at the longer wavelength (Table 1). This is a direct consequence of the significantly higher epidermal fluence at 600 nm in comparison with 585 nm (Fig. 4). Corresponding PPTR signals computed with the FD model are presented in Fig. 5(b). Immediately after the end of the laser pulse, the radiometric temperature rise is nearly identical at both wavelengths due to the similar initial epidermal temperature rise [Fig. 5(a)]. At later times, as heat diffuses from the PWS towards the skin surface, the PPTR signal contains information pertinent to heating of increasingly deeper skin layers. Since the PWS temperature rise was considerably higher at 585 nm, the 585 nm PPTR signal is characterized by a higher delayed component than the 600 nm counterpart.

Figure 5

(a) Computed initial temperature profiles at the end of 1.5 ms laser pulses at 585 and 600 nm. (b) PPTR signals computed for incident 585 and 600 nm laser pulses. The FD model was used to calculate heat transfer dynamics with the initial temperature profiles in (a) used as input. Assumed signal-to-noise ratio was set to 1000.

These simulated PPTR signals were input into an iterative non-negatively constrained conjugate gradient algorithm to reconstruct initial temperature profiles. From the SWE-PPTR reconstructed profile [Fig. 6(a)], ΔT _epi,max and z _PWS,max can be determined (Table 4); however, z _epi and z _PWS cannot be determined. With separate epidermal and PWS profiles obtained from DWE-PPTR, z _epi and z _PWS can be determined in addition to ΔT _epi,max and ΔT _PWS,max [Fig. 6(b)], demonstrating that the DWE approach allows us to characterize features of PWS skin even when the PWS is in close proximity to the epidermal-dermal junction. DWE-PPTR parameters α and β [Eq. (8)] were determined as 0.47 and 1.07, respectively. The small error bars in Fig. 6 illustrate the stability of the iterative reconstruction process over a range of iterations in close proximity to the corner of the L curve.

Figure 6

Comparison of (a) SWE- and (b) DWE-PPTR reconstructed profiles with the actual 585 nm initial temperature profile. In (a), ΔS₅₈₅(t) was used as input, and “Actual” represents the actual initial profile and “Reconstructed” the profile reconstructed by applying the conjugate gradient algorithm to the SWE-PPTR signal. In (b), both signals from Fig. 5 were used as input. The actual 585 nm initial temperature profile (Actual) is shown in gray. X(z) and Y(z) represent profiles reconstructed from signal components x(t) and y(t), respectively. Error bars represent standard deviations of five iterative solutions. DWE parameters: α=0.47, β=1.07.

Table 4

A deep component was evident in the reconstructed epidermal profile [Fig. 6(b), dashed line], in contradiction with anatomical knowledge that no melanin is present in human skin at depths greater than 100–200 μm.³⁸ Majaron et al.;^{25 26} attributed this deep component to the approximate nature of Eq. (8). A 600 nm light penetrates more deeply in PWS skin than 585 nm light. However, Eq. 8(b) only includes a multiplicative factor α as the difference between the PWS components of the PPTR signals at the two wavelengths. As a result of this discrepancy, the solution Y(z) reconstructed from signal component y(t) includes a deep component which originates from PWS blood absorption.²⁵ This subcomponent must be added to the profile X(z) reconstructed from x(t) to constitute the actual PWS temperature profile [Fig. 6(b)]. After this modification was implemented, we used our model to compute separately initial temperature profiles due to epidermal and blood absorption [Figs. 7(a) and 7(b), respectively]. Comparison of actual and DWE-PPTR reconstructed epidermal profiles [Fig. 7(a)] revealed a good match in shape, as well as in z _epi and ΔT _epi,max values (Table 4). Similarly, the general shape of the actual PWS profile [Fig. 7(b)] was reproduced reasonably well.

Figure 7

Comparison of (a) epidermal and (b) PWS actual initial temperature profiles and the DWE-PPTR reconstructed profiles.

4.

Discussion

In this study, a biopsy-defined PWS skin geometry²³ was employed. Using serial histological sections of PWS skin as input, the resulting laser induced temperature profile [Fig. 5(a)] was considerably more complex than those used in previous studies.^{18 19 35}

If the most superficial PWS blood vessels are distinctly separated from the epidermis, it is usually possible to extract clinically relevant PWS characteristics from SWE-PPTR reconstructed temperature profiles. However, in most PWS patients, the most superficial blood vessels are located within the upper 200 μm of the dermis,² and thus SWE-PPTR reconstructed profiles oftentimes resemble the profile in Fig. 6(a), from which only ΔT _epi,max and z _PWS,max can be extracted.

In DWE-PPTR reconstructed temperature profiles of experimental data,^{25 26} an overlap of the epidermal and PWS components was often observed. Proposed explanations were undulations in the epidermal basal layer and broadening of one or both reconstructed profiles due to the limited spatial resolution of PPTR.^{19 25 39} From the actual epidermal and PWS temperature profiles (Fig. 7), an overlap of the epidermal and PWS profiles was indeed present. Furthermore, in the presented simulated example of DWE-PPTR profiling, separation of the reconstructed temperature profile into epidermal and PWS contributions yielded a ∼50 μm overlap between the two solutions [Fig. 6(b)]. Comparison of the reconstructed PWS profile to the actual initial profile [Fig. 7(b)] revealed that the former was indeed broader. An increased sampling rate of ΔS(t) may reduce this effect but may increase the noise level in ΔS(t), potentially affecting the stability of the iterative reconstruction process.

DWE parameters α and β were determined as 0.47 and 1.07, respectively. These values are in good agreement with those obtained from experimental data of fair-skinned patients.^{25 26 35} Majaron et al.;²⁵ also determined that the sum of the integrals of the epidermal and PWS reconstructed temperature profiles matches the integral of the SWE-PPTR reconstructed profile within 1–2. From our results, the difference between the two integrals is only ∼2, supporting the selected values of α and β.

A qualitative comparison demonstrates that DWE-PPTR can provide accurate temperature profiles of the epidermis (Fig. 7(a), Table 4). However, in the reconstructed PWS profile, high spatial frequency components in the profile were not reproduced, in agreement with previous theoretical and experimental studies.^{18 19} The area of skin included in the histology sections was approximately 0.25 mm×1.6 mm, which was ten times smaller than the 3.7 mm² region typically interrogated in our experiments.^{25 26} If more histological sections were utilized in the simulations, the PWS profile probably would have contained fewer high spatial frequency components, thereby improving our ability to reconstruct the initial temperature profile. Nevertheless, location and shape of the superficial part of the PWS temperature profile was reconstructed reasonably well, allowing determination of z _PWS within 20 μm and z _PWS,max within 15 μm (Table 4). Using DWE-PPTR, all four clinically relevant PWS skin characteristics can be determined with an accuracy sufficient to guide selection of optimal laser treatment parameters on an individual patient basis (Table 4).

5.

Conclusions

A 1-D optical-thermal model was developed to evaluate PPTR depth profiling of PWS skin. A biopsy-defined PWS skin geometry was used to add realism to the model. From a preliminary sensitivity analysis, the computed initial temperature profile was minimally affected by large differences in optical properties. From SWE-PPTR reconstructed profiles, ΔT _epi,max and z _PWS,max were determined accurately, but z _epi and z _PWS could not be determined. Comparison of the actual and DWE-PPTR reconstructed profiles revealed a good match for all four PWS skin characteristics of clinical interest. Results of this study indicate that DWE-PPTR is a viable approach for depth profiling of PWS skin, for the purpose of laser therapy optimization.