Pulsed photothermal profiling of hypervascular lesions: some recent advances

Pulsed photothermal radiometry (PPTR) can be used for non- invasive depth profiling of port wine stain (PWS) birthmarks, aimed towards optimizing laser therapy on an individual patient basis. Reconstruction of laser-induced temperature profile from the experimentally obtained radiometric signal involves the skin absorption coefficient in the infrared detection band. In the commonly used 3 - 5 micrometer detection band (InSb), the absorption coefficient varies by two orders of magnitude, while assumed to be constant in the reconstruction algorithms used thus far. We discuss the problem of choosing the effective absorption coefficient value to be used under such conditions. Next, we show how to account explicitly for the strong spectral variation of the infrared absorption coefficient in the image reconstruction algorithm. Performance of such improved algorithm is compared to that of the unaugmented version in a numerical simulation of photothermal profiling. Finally, we analyze implementation of a bandpass filter which limits the detection band to 4.5 - 5 micrometer. This reduces the absorption coefficient variation to a level that permits the use of unaugmented algorithm. An experimental test of the latter approach for in vivo characterization of the depth of PWS lesion and epidermal thickness will be presented, including a novel technique that uses two laser excitation wavelengths in order to separate the epidermal and vascular components of the radiometric signal.


INTRODUCTION
Port-wine stain birthmarks (PWS) are hypervascular lesions in human skin, which consist of an excess of ectatic blood vessels. These are usually fully contained within the most superficial millimeter of the skin. The exact depth varies from patient to patient, but on average, the highest fractional blood content is found 200-400 tm below the epidermal-dermal un' PWS are currently treated by selective photocoagulation ofthe ectatic vasculature using a pulsed green (KTP/YAG, 532 nm) or yellow/orange laser (flashlamp-pumped dye, 577nm, 585-600 nm). In order to optimize therapy on an individual patient basis, determination ofboth PWS depth and epidermal thickness is required, especially when cryogen spray cooling is applied.2'3'4 Pulsed photothennal radiometry (PPTR), which is based on time-resolved acquisition of infrared (IR) radiant emission following pulsed laser exposure, was recently introduced for assessment of laser-induced temperature profiles in PWS.5 '6 The ability ofPPTR to determine the depth of sub-surface chromophores has been demonstrated by profiling layered tissue phantoms, and by comparison with histological assessment ofPWS depth.7 Reconstruction of laser-induced temperature profiles from experimental radiometric signals involves the value of the tissue absorption coefficient at the infrared detection wavelength. While assumed to be constant in the reconstruction algorithms developed thus far, the absorption coefficient of water, which is the primary absorber in the commonly used 3-5.tm detection *Coespondence: 1002 Health Sciences Rd. East, Irvine; e-mail: majaron@bli.uci.edu, boris.majaron@ijs.si; fax: + 1-949-824-6969 In Lasers in Surgery: Advanced Characterization, Therapeutics, and Systems X, R. Rox Anderson, band (InSb), varies by two orders of magnitude (Fig. 1). 8 We discuss below, how to obtain the effective value of the JR absorption coefficient under such conditions. Next, we show how to account explicitly for such strong spectral variation of the absorption coefficient in the image reconstruction algorithm. Performance of such an improved algorithm is compared to that ofthe unaugmented version in a numerical simulation ofphotothermal profiling of PWS.
Alternatively, a bandpass filter can be implemented, which narrows the JR detection band to a level that permits the use of the unaugmented algorithm. An experimental test of the latte:r approach using a 4.5-5 tm acquisition band for in vivo characterization of PWS will be presented. Note that in lesions with PWS in close proximity to the epidermal-dermal junction, epidermal heating due to broad melanin absorption may prevent determination of PWS depth, owing to limited spatial resolution of PPTR.6'9 In order to overcome this problem, we use an approximation technique, which utilizes two laser excitation wavelengths to separate the epidermal and vascular components of the radiometric signal.1° In this way, the epidermal thickness and PWS depth can be assessed with adequate accuracy and reliability.

PPTR PROFILING
The basic theory of PPTR profiling will be reviewed first.5'6 The temperature field evolution following a pulsed laser exposure of skin can be treated by one-dimensional heat-diffusion equation (using thermal diffusivity value D = 1 .l• 10 7 m2/s). Such an approach is reasonable, since the laser irradiated spot is typically much larger (>5 mm) than the dermal depths considered in this study (<1 mm). Using the Green's function approach, the time evolution of temperature increase AT(z,t) can be calculated from the initial laser-induced temperature profile AT(z ',O) as: where KT(z',z,t) represents the thermal point spread function of the problem.
(1) Note that, owing to the finite penetration depth of the emitted infrared radiation, JR detectors do not measure the surface temperature. The radiometric signal obtained from a non-uniformly heated body is in general a sum of contributions from different depths z (measured from the sample surface inward): Wavelength (tim) 8 Figure 1 : Infrared absorption coefficient of water, displaying a variation of two orders of magnitude within the detection band of the lnSb detector (? = 3-5 im).8 S(t) = C IRJB(T(Zt)) e_1RZ dz . (2) Here, constant C accounts for skin emissivity, geometry and transmission of the collection optics, sensitivity of the detector, and similar factors, .tIR is the tissue absorption coefficient at the detection wavelength, and B(T) is given by Planck's law of radiation. For relatively small temperature increases encountered in our experiments, Planck's law can be linearized around the initial skin temperature I'0, resulting in S(t) = CB(T0) + C Bx(To) lRJAT(Zt) e_1RZ dz . (3) For the discussed case of pulsed laser heating, the dynamic component of the radiometric signal AS(t) can now be written by inserting the temperature increase z\T(z,t) (1) into (4), and reordering: After the point spread function of the problem K1(z',z,t) is established, integration over z can be performed, and it becomes clear that the PPTR signal (minus the background) AS(t) is linear in the initial temperature profile AT(z ',O): Defmition of the kernel function K(z',t) is evident from (4).

Experimental setup and methods
A PWS lesion on a volunteer patient is irradiated with a 1.5 ms long pulse from a flashlamp-pumped dye laser (ScleroPLUS, Candela, Wayland, MA), which allows a choice of four wavelengths (585 to 600 nm in increments of 5 nm). Energy density at the center of a 6 mm diameter laser spot is 5-6 J/cm2. The transient increase in infrared radiant emission from the central 1 .9x 1 .9 rmn2 area is recorded with an JR camera sensitive to a wavelength band of 3-5 tm (Gallileo, Raytheon, Dallas, TX), at a rate of 500 frames per second. PPTR signals AS(t) are obtained by calibrating the response of the camera's 64x64 array detector elements with a computer-controlled blackbody (BB7O1, Omega Engineering, Stamford, CT), averaging over the JR camera array, and subtracting the radiant emission level prior to the pulsed laser exposure (background).
Reconstruction of the initial temperature profile T(z ',O) from the experimental PPTR signals AS(t) is a severely ill-posed inverse problem. We solve it using a dedicated iterative algorithm based on a nonnegative-constrained conjugate-gradient method. 6 Since the PPTR signals, as well as the numerical solution in terms of the estimated initial temperature profile, are inevitably discretized, the linear mapping (5) of the initial temperature profile vector T (T = iT(z,O)) onto the signal vector S (S, = AS(t)) can be represented by a kernel matrix j: S = K T; K K(z'1,tj). (7) The computational load is somewhat reduced by the fact that the kernel matrix elements K have to be calculated only once for each image reconstruction example. Nevertheless, a characteristic inversion run involving a 400-element signal vector and a 64-element temperature profile vector requires a few minutes CPU time on a Sun Sparc II workstation. The computation is regularized by early termination.6 Rather than dwelling on the issue of determining the optimal number of iteration steps, we show in the discussed example several iterative solutions, which cover the whole range from smooth, blurred temperature profiles in under-iterated solutions, near optimal and often also over-iterated solutions. As discussed previously by Mimer et al.,6 the latter display characteristic oscillations, as the algorithm attempts to fit experimental noise and other mismatches between the experimental and theoretically predicted signals.

Response function R(z)
In order to analyze the influence of the strong spectral variation of the JR absorption coefficient IR(X), radiometric signals must be calculated by adding radiative contributions (2) within the detection band (?i to ?.h): xl z=O where we have also taken into account the spectral response ofthe detector D(?.). Equation (8) can be rewritten as: The response function R(z), which is obtained by integrating over the detection wavelengths, represents the radiometric signal resulting from a unity temperature increase in a delta-layer at a depth z. In general, its functional dependence deviates from the previously assumed exponential decrease with depth (3).
l9mm' R(z),calculated for tu9) variation of skin (70%water content), detector (InSb) spectral response proportional to wavelength, 3-5 tm detection band, and initial temperature of T0 =36°C (solid line). Three response functions obtained with various constant absorption coefficient values are also presented for comparison, displaying the exponential decay with depth (dashed and dotted lines; see the legend for values).

Effective absorption coefficient
In the next step, we attempt to determine which effective value of the JR absorption coefficient in should be used in PPTR image reconstruction algorithms that ignore its spectral dependence within the detection band. By introducing a spectrally constant effective value integration over 2. in (8) yields a constant, which can be taken out ofthe integral over depth z: dz .
By equating the response functions R(z) in (9) and (10), the effective absorption coefficient value tfl is given by the implicit equation: S D(2) Bx(To) The numerical solution of (1 1) for the 3-5 .tm detection band and spectral response of the InSb detector is presented in Figure 3 (solid curves). Note that the effective absorption coefficient value varies strongly with depth (z), and is doublevalued. While the latter effect is intuitively quite surprising, it does not present a computational problem, as any of the two branches can evidently be used to calculate the signal from (10) (or (4)). The strong depth dependence, however, indicates that there is no constant value of that would yield a correct time dependence of the PPTR signal under the discussed experimental conditions. Consequently, the reconstruction algorithm has to be augmented by taking into account the spectral dependence according to (9) or, alternatively, the detection bandwidth should be decreased to a level that permits the use of the unaugmented algorithm. Both approaches are discussed in sections 4 and 5, respectively.

IMPROVED KERNEL MATRIX
In the following, we compare the performance of an image reconstruction algorithm that explicitly accounts for the spectral variation ofthe JR absorption coefficient jt1R(A) with the unaugmentedversion (using a constant value J-teff). This comparison is performed in a numerical simulation of photothermal profiling. First, an initial temperature distribution AT(z ',O) ("object") is chosen, and the corresponding PPTR signal AS(t) calculated using the kernel function K(z',t), which involves the "correct" response function R(z) (9). Then, this signal, with some noise superimposed onto it, is used as an input to the iterative image reconstruction algorithm. The reconstruction is performed first with the kernel matrix elements K (7) calculated using the simplified response function (10) with teff 22 1111Ti, Ifld then using the correct response function R(z) (9) (see Fig. 2).
Results from one such comparison are presented in Fig. 4. As the "object", we use a hyperGaussian temperature profile AT(z ',O) = M'0 exp(-2(x---xo)4/w4) with AT0 = 1°C, xo 300 pm, w = 100 tm (dashed curves in Fig. 4), which vaguely resembles a laser-heated PWS layer. The corresponding signal AS(t), which includes also a small amount of normal noise (SNR = 1000) is calculated at 250 equidistant time points (At =2 ms), and the solution (depth) space which represents the most superficial millimeter of the skin is divided into 64 intervals. As a matter of convenience, the simulation was performed on a PC (Pentium 11-350), using a slightly adapted reconstruction method (details to be presented elsewhere). The under-iterated solution obtained after 20 iteration steps (n = 20) is a blurred and deformed image of the actual temperature profile, and is very similar for both approaches. When using the simplified kernel function (Fig. 4a), increasing the number of iteration steps (n) seems to temporarily improve the result, but then very rapidly produces an oscillatory a n =20 U I,. Depth (mm) solution, which indicates the presence of an additional PWS layer starting at 120 imclearly a computational artifact. In contrast, the solution obtained using the corrected kernel matrix converges rapidly toward a significantly more representative image of the actual temperature profile (Fig. 4b). The position of the shallower edge of the profile, which is the most important parameter in the targeted application, is reproduced with remarkable accuracy. The peak temperature value is overestimated by only 8 %, compared to 43 % when using the simplified kernel function. The convergence of the reconstruction is also improved remarkably; when using the corrected matrix, the residual norm (which is minimized in the iteration) reaches a value of 9.6 i09 after 50 iteration steps, and changes only minimally (< 1%) by iteration step 200. In contrast, with the simplified kernel function, the residue norm after 200 iteration steps is still 2.7 106 (!) a 3 % decrease fromstep 100.
The difference between these two results is illustrated further in Fig. 5, where the radiometric signal calculated by each method from its last solution (n = 200) is compared to the signal used as the input to the reconstruction algorithm. Using the simplified kernel function (Fig. 5a), the mismatch between the input ("object") and predicted signal ("image") is particularly large in the first 40 ms, where the error exceeds 30 % of the signal value at that point, and 3 % of the maximalsignal value overall. In contrast, with the corrected kernel function, the two signals are undistinguishable (Fig. 5b), and even the 100-times expanded residue can be hardly observed (see the inset). Our general experience from other similar simulations is that higher levels of noise in the PPTR signal prevent successful convergence ofboth methods. The benefit ofthe corrected kernel function may thus not show up as clearly as in the presented example. This is particularly true for more complex initial temperature profiles. Nevertheless, we are confident that the performance of PPTR profiling should benefit from the suggested correction of the kernel matrix in many realistic experimental situations. In our PPTR profiling experiments, we have regularly observed a structure of the residue (similar to the one presented in the inset ofFig. 4a) well above the noise level.

BANDPASS FILTERING
As discussed above, an alternative approach to using the more accurate (and complex) formula in calculation of the kernel matrix elements could be a reduction of the detection bandwidth to a level that permits the use of the unaugmented algorithm. Using the equation (1 1), we have calculated t(z) for watery tissue for several detection bands. Based on the results of such an analysis, and wavelength dependence of Planck's factor B(T), we decided to reduce the detection bandwidth of our system optics of the InSb-array camera. As indicated by Fig. 3 (dashed curve), radiometric signals obtained with such a modified detector can be accurately described by the simplified formula (10) and effective absorption coefficient equal to teff 26.5 mm'. We have estimated that by using the additional filter, the amplitudes of the acquired signals AS(t) would decrease roughly by a factor of two. Provided that performance of the image reconstruction was limited primarily by error due to spectral variation of t(A), rather than the level of experimental noise, such modification should improve the overall performance ofPPTR profiling. . When using the reduced detection band, the camera integration time was increased from 0.5to 0.7 s, in order to compensate for the anticipated reduction of the radiometric signal. These two signals were used as input to the image reconstruction algorithm based on the simplified approach with a constant jtIR (4). In accordance with results presented in Figs. 2a, 3, the value of 1eff 22 mm was used to analyze the signal "a", and = 26 mm' for signal "b". The two families of iterative solutions (Figs. 7a, b) run from smooth, blurred temperature profiles in under-iterated solutions (see n = 2 in Fig. 7b), to over-iterated solutions. The latter display characteristic oscillatory artifacts, resulting from the illposedness of the inverse problem, combined with the experimental noise and other mismatches between the experimental and theoretically predicted signals. This is particularly pronounced in example "b", which corresponds well to the approximately two times lower SNR in the PPTR signal (Fig. 6).
In both examples, the most reliable solution seems to be found after approximately five iteration steps (solid curves).
It is evident, however, that with either approach, the temperature profile in this specific PWS lesion cannot be resolved from the epidermal temperature rise, resulting from melanin absorption. This may happen when the two layers lie in close proximity to each other, and determination of PWS depth and epidermal thickness, which are required to optimize laser treatment, is thus prevented due to the limited spatial resolution of PPTR profiling. To overcome this problem, we have developed an approximation technique, which achieves improved selectivity by utilizing two excitation

Two-wavelength excitation technique
This approximation technique allows us to differentiate between the PWS and epidermal contribution to the PPTR signal by exploring different spectral properties of the two chromophores (hemoglobin, melanin) at two excitation wavelengths. This enables detennination of PWS depth and epidermal thickness even when the two layers are in close proximity to each other. Moreover, it breaks down the relatively complex initial temperature profile into two simpler "objects", which are easier to reconstruct.
The method is based on the linearity between the PPTR signal AS(t) and the initial temperature distribution AT(z,O) (5). An experimental signal obtained, for example, following pulsed irradiation at 585 nm, can thus be expressed as a sum of two components, originating from the heated PWS layer [x(t)], and epidermis [v(t)]: At 600 nm, the absorption and scattering properties of the epidermis and dermis are essentially equivalent to those at 585 nm' As a result, the temperature profile in the epidermis and, consequently, the corresponding signal component, is practically the same with both excitation wavelengths. In contrast, light absorption by blood at 600 nm is significantly weaker than at 585 12 This results in a lower PWS contribution to the PPTR signal obtained after 600 nm irradiation, which can therefore be approximated to the first order by AS600(t) =ax(t) + 3y(t) .

122
Depth (mm) Depth (mm) Based on the above discussed spectral properties, we can expect a < 1 and 3 1.
From (12) and (13), the PWS component ofASs8s(t) can be derived as By using x(t) as input to the image reconstruction algorithm, the temperature profile resulting exclusively from blood absorption could be obtained, but values of a and f3 need to be determined first. We rely on the fact that the radiometnc signal immediately following pulsed irradiation results predominantly from temperature increase in the superficial layer of tissue -a few penetration depths of the detected JR radiation (40-50 trn). Since the blood vessels are located deeper in the skin, we can estimate f3 by requiring that the PWS signal component starts from x(O) =0. A signal 3 s585(t) -s600(t), proportional to x(t) (14), can then be calculated. The initial temperature profile calculated from such signal is proportional to the temperature profile in PWS, and shows its top boundary very clearly (Fig. 8).
Comparison of the two examples in Fig. 8 shows that with the reduced detection band (and thus more constant the iteration can be run longer before the non-physical oscillations appear in the solution (n = 20 in Fig. 8b vs. n = 10 in Fig. 8a).
Furthermore, the residual norm reaches a lower value (12.3 at steps 10 and 20, vs. 13 .8 in example "a"), despite the significantly higher noise level and absolute signal level in the corresponding input PPTR signal ("b" in Fig. 6). This indicates a closer match between the theoretically predicted signal shapes and the experimentally obtained PPTR signals when using the appropriately reduced detection band. Note that, whereas the PWS depth (most superficial boundary) is the same in both examples (within the estimated accuracy of the solution), the depths of the temperature peak are significantly different (220 im in Fig. 8b vs. 190 tm in Fig. 8a).
Similarly as above, the epidermal contribution to 1S585(t) can he derived as To determine the value of a, we use the knowledge that the pure epidermal (melanin) temperature profile must be zero at depths deeper than 150-200 tm into the skin. We therefore calculate signals y(t) (15) with increasing values of a, starting from 0, and observe the resulting temperature profiles. These usually differ little in the superficial layer, while temperatures deeper in the skin decrease monotonically with increasing a, and fmally disappear. At this value of a, the solution in general oscillates less, and the residue norm reaches a lower value, than at lower (and sometimes also higher) values. Thus obtained estimate of a can then be used to renormalize the previously calculated PWS profiles.
The final result, displaying both PWS and epidermal temperature profiles, is plotted in Figure 9. The depth of the PWS most superficial boundary, which is clearly distinguished from the epidermis, is determined with reasonable accuracy. The overlap of the two profiles at lOO tm most likely results from the uneven boundary between the epidermis and dermis, but could in part result also from broadening of one or both profiles by the reconstruction algorithm.
The epidermis profile deeper than 0.25 mm is an artifact of the method and should be disregarded. We attribute it to differences between the temperature profiles induced in PWS with 585vs. 600 nm irradiation. Since the latter penetrates deeper into the PWS lesion, the validity of proportionality between the two PWS temperature profiles, assumed in (13), deteriorates with depth. However, such a first order approximation definitely holds in the thin top layer of the PWS, where both temperature profiles (using 585 and 600 mu excitation) can be linearized. The solutions presented in Fig. 9 are thus representative up to a depth of '.O.25 mm. At greater depths, the PWS profile solution may differ from the actual temperature profile induced by the 585 nm excitation, as some contribution from deeper blood vessels is contained in signal y(t).
More experimental examples and a thorough discussion of the strengths and limitations of this approach will be presented elsewhere.

CONCLUSIONS
PPTR can be used for non-invasive depth profiling of PWS, with a potential for on-line guidance of laser therapy on an individual patient basis. In reconstruction of the laser-induced temperature profiles from the experimentally obtained