Analysis of thermal relaxation during laser irradiation of tissue

BACKGROUND AND OBJECTIVE
Thermal relaxation time (tau(r)) is a commonly-used parameter for estimating the time required for heat to conduct away from a directly-heated tissue region. Previous studies have demonstrated that temperature superposition can occur during multiple-pulse irradiation, even if the interpulse time is considerably longer than tau(r). The objectives of this study were (1) to analyze tissue thermal relaxation following laser-induced heating, and (2) to calculate the time required for a laser-induced temperature rise to decrease to near-baseline values.


STUDY DESIGN/MATERIALS AND METHODS
One-dimensional (1-D) analytical and numerical and 2-D numerical models were designed and used for calculations of the time tau(eff) required for the peak temperature (T(peak)) to decrease to values slightly over baseline (DeltaT(base)). Temperature values included T(peak)=65 and 100 degrees C, and DeltaT(base) = 5, 10, and 20 degrees C. To generalize the calculations, a wide range of optical and thermal properties was incorporated into the models. Flattop and gaussian spatial beam profiles were also considered.


RESULTS
2-D model calculations of tau(eff) demonstrated that tau(eff) (2-D) was as much as 40 times longer than tau(r). For a given combination of T(peak) and DeltaT(base), a linear relationship was calculated between tau(eff) (1-D) and tau(r) and was independent of optical and thermal properties. A comparison of 1-D and 2-D models demonstrated that 1-D models generally predicted longer values of tau(eff) than those predicted with a 2-D geometry when the laser spot diameter was equal to or less than the optical penetration depth.


CONCLUSION
Relatively simple calculations can be performed to estimate tau(eff) for known values of tau(r), T(peak) and DeltaT(base). The parameter tau(eff) may be a better estimate than tau(r) of tissue thermal relaxation during multiple-pulse laser irradiation.


INTRODUCTION
Pulsed laser radiation is used to heat tissue and precisely ablate or coagulate tissue. In laser surgery, the general goal is to selectively heat a desired region of tissue while minimizing collateral thermal damage. The theory of selective photothermolysis proposed by Anderson and Parrish [1] is generally applied to a given situation to determine the appropriate laser pulse duration (t p ).
From the theoretical analysis of Anderson and Parrish [1], the concept of thermal relaxation time has become a popular term for selection of laser pulse duration. The thermal relaxation time (t r ) of a heated region of tissue is the time required for the peak temperature rise (Á T peak ) in a heated region of tissue to decrease to 37% of the total rise: where d is usually the penetration depth (cm). The appropriate value for d is discussed in more detail in the Discussion section. A common belief is that if t p is less than t r of the irradiated region, thermal damage is minimized. However, for procedures in which high peak temperatures (T peak ! 100 C) are expected, such as laser skin resurfacing, the temperature value may still be suf®ciently large to induce thermal damage after t r has elapsed. For example, assuming an initial temperature of 30 C and a laserinduced temperature rise to 100 C, T peak will be approximately (100 C À 30 C)* 0.37 30 C 56 C after t r has elapsed. If another laser pulse is applied to the tissue at this point, thermal effects may be more severe due to temperature superposition effects. Complete thermal relaxation of tissue may take several thermal relaxation times [2]. Also, t r is applicable only for speci®c tissue conditions. Tissue optical and thermal properties may change during the laser pulse and for consecutive pulses due to evaporation and tissue denaturation [3,4]; as these properties change, the value of t r changes as well. Several theoretical and experimental studies have shown that temperature superposition can occur even if the time between pulses is several orders of magnitude greater than the thermal relaxation time [5±8]. In previous measurements of skin radiometric surface temperatures during pulsed CO 2 laser irradiation (t r $700 ms) [5], it was demonstrated that elevated temperatures remained for 20±40 ms after the end of a laser pulse. These results indicate that complete tissue thermal relaxation to baseline temperatures is a relatively slow process.
The purpose of this study was to investigate theoretically the thermal response of soft tissue to laser-induced heating. The time required for T peak to decrease to values just above baseline is calculated using a 2-D numerical thermal model. Results of a 1-D analytical solution to the heat conduction equation are compared to those of the numerical model to demonstrate its applicability as a ruleof-thumb estimation for the effective relaxation time (t eff ) of tissue.

MATERIALS AND METHODS
Conduction heat¯ow is governed by the heat conduction equation. In three dimensions, conduction heat¯ux q (W/cm 2 ) is de®ned as: q ÀkrTY 2 where k is thermal conductivity (W/cm/K) and r is the 3-D del operator. In Cartesian coordinates, The negative sign in Equation (2) is due to the direction of heat conduction; heat always¯ows from a higher temperature region to a lower temperature region.
In cylindrical coordinates, Equation (2) becomes: where q gen is the laser source term (W/cm 3 ) and a is thermal diffusivity (cm 2 /s). If 1-D heat¯ow is assumed, Equation (2) takes on the following form:

2-D Finite Difference Model
A ®nite difference solution [9] to the heat conduction equation (Equation (2)) was derived. Cylindrical coordinates were used, resulting in a 2-D axisymmetric tissue geometry (Equation (3)). The tissue was assumed to be homogeneous and its dimensions were 1-cm thickness and 1-cm radial extent. The size of each grid element was 40 Â 20 mm.
Accurate thermal modeling of tissue subjected to relatively high temperatures is dif®cult due to potential dynamic changes in optical and thermal properties as a function of temperature and water content [3,4]. For the 2-D simulations, we assumed that tissue consisted of 70% water (representative of dermis [10]) and calculated thermal properties from the following equations [11]: where r is density (g/cm 3 ), c is speci®c heat (J/g/K), and W is water content (e.g., for 70% water, W 0.7). For the 1-D simulations (described below), other sets of thermal properties were considered that corresponded to different water contents. Absorption coef®cient (m a ) values of 10, 100, and 1000 cm À 1 were selected. These values were representative of those commonly encountered during laser-mediated tissue heating. Intermediate absorption coef®cients were used in the 1-D simulations (described below). Three different spot diameters (o o 500 mm, 1.5 mm, and 2.5 mm diameters) were used to take into account various spot size/penetration depth (d 1/m a ) combinations. Flattop and gaussian spatial beam pro®les were considered. For gaussian pro®les, o o represents the 1/e 2 diameter.
Pulse durations were chosen to correspond to 1/100th of the thermal relaxation time. This rule for pulse duration selection was chosen since it is reasonable to assume that such a value would generally be considered as suf®cient for limiting heat conduction during the pulse.
Radiant exposures were chosen so that the peak temperature T peak at the end of the laser pulse corresponded to either 65 or 100 C. The value of 65 C was chosen because it corresponds to a commonly-cited threshold value for collagen denaturation [7,12,13], although in reality thermal damage is a temperature-time phenomenon [14]. 100 C corresponds to water boiling at atmospheric pressure and is a good estimate of the residual peak temperature in tissue immediately after ablation. The effective relaxation time (t eff ) required for the peak temperature to fall to ÁT base 10 C above baseline (T baseline 22 C) was calculated for each model run.

1-D Analytical Expression
An advantage of using 2-D models is that they offer more realism than 1-D models. A major disadvantage is the considerably longer computation time required to arrive at a solution. A 1-D equation derived by Anvari et al. [15] from basic principles presented by Carslaw and Jaeger [16] describes 1-D heat¯ow in a semi-in®nite medium with an adiabatic surface boundary condition. This equation was slightly modi®ed to the following form: is the complementary error function, and erfcx(x) exp(x 2 ) erfc(x). Equation (5) is valid for an initial temperature distribution that is proportional to expÀm a z; in other words, heat conduction during the pulse is assumed to be negligible. Equation (5) was ®rst compared to results of a 1-D implicit ®nite difference model in which heat conduction during the pulse was assumed to determine the relative accuracy of this equation. Results of this equation were then compared to those obtained from the 2-D model to determine its suitability for predicting t eff .
Since the time required to implement Equation (5) was considerably shorter than that required with the 2-D ®nite difference equations, this equation was applied to a wider variety of cases. Tissue thermal properties were evaluated for four different tissue scenarios: (1) 100% water, (2) 70% water, (3) 0% water, and (4) pure carbon. These sets of thermal properties represented asymptotic limits of tissue states encountered during high-power laser irradiation of tissue and are compiled in Table 1. Equations (4a)±(4c) were used to calculate thermal properties for the ®rst three cases; pure carbon thermal properties were taken from various sources [9,17,18].
The following absorption coef®cients were used: 10, 20, 30, 40, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, and 1000 cm À 1 . Radiant exposures were chosen so that peak temperatures of 65 and 100 C were obtained by using the following equation: where ÁT peak T peak À T baseline is the peak temperature rise ( C). Equation (6) is simply a 1-D solution to the heat conduction equation at the surface and in the absence of conduction. The time required for the peak temperature to decrease to ÁT base values of either 5, 10, or 20 C above baseline were calculated. The calculations for 10 C above baseline were directly comparable to the 2-D model runs.

2-D Finite-Difference Model
A total of 36 simulations were run.  chosen so ÁT peak at the end of a 17.6-ms pulse was 65 C. The time required for T peak to decrease to 32 C was 9.57 seconds. For the experimental conditions simulated in this case, t r was 1.76 seconds. Thus, in this example, t eff was 5.44 times longer than t r . A summary of model results is provided in Table 2. The dimensionless variables À and t represent the ratios o o /d and t eff /t r , respectively. Note that for all but one model run, t is greater than unity, indicating that the times required for T peak to decrease to 32 C were longer than t r . Figure 2 shows a plot of t eff calculated for gaussian beam pro®les as a function of t eff for¯attop beam pro®les. Each coordinate pair represents identical simulated experimental conditions except for the spatial beam pro®le ( Table 2). The solid line represents the boundary at which t eff,gaussian t eff,¯attop . All data points fall below this line, indicating that for a given set of experimental conditions, the value of t eff required for gaussian temperature distributions to thermally relax is shorter than that for uniform distributions. The discrepancy is larger for longer values of t eff .

1-D Analytical Equation
Equation (5) was applied to a total of 456 different experimental scenarios. Calculations of t eff for relatively short pulses (t p 0.01 t r ) were compared to those obtained from 1-D implicit ®nite difference calculations and agreed to within 1%. Results of all model runs are shown in Figure 3. Note that all calculations of t eff fall above the solid line, indicating that the times required for T peak

Model run
Beam type m a (cm À 1 ) decreasing to values 5±20 C above baseline are longer than t r . In Figure 3, it is evident that for given combinations of peak temperature rise (T peak 43 or 78 C) and ®nal peak temperature above baseline (ÁT base 5, 10 or 20 C), a plot of t eff vs. t r yields data that fall along straight lines. Average slopes of linear ®ts are listed in Table 3. The slope of each line dictates the offset of each line on the plot; a larger slope indicates a higher offset. Note that for a given combination of ÁT peak and ÁT base , the linear slopes are identical regardless of the combination of tissue optical and thermal properties used in the model.  (2-D model). In general, the 1-D model overpredicts t eff (e.g., y-axis values greater than unity). This overprediction is associated with the ratio (À) of spot size (o o ) to penetration depth (d). When o o ) d (e.g., for large À), the 1-D solution is accurate at the center of the laser spot, and thus the values of t eff calculated with the 1-D and 2-D models are comparable (Fig. 5).

DISCUSSION
The main advantage of t r is that it is a simple quantity to calculate. It can also be used to adequately predict thermal effects from single-pulse irradiation procedures, such as Port Wine Stain treatment [19]. For cases in which multiple pulses are delivered, such as laser skin resurfacing or laser-assisted blepharoplasties, our results demonstrate that t r is not a good indicator for estimating thermal effects ( Fig. 3 and Table 2). Model calculations predict that t eff (2-D) can be as much as 40 times longer than t r for peak temperature relaxation from 100 to 32 C for a 22 C baseline temperature; other researchers have estimated comparable values of t eff [7,8].
Since t eff can be on the order of 30 seconds (e.g., model run #6 in Table 2), it is necessary to consider the effects of blood perfusion on t eff . Welch et al. [20] estimated that perfusion has a signi®cant effect on temperatures in normal skin after approximately 100 seconds of laser heating. Thus, the values of t eff calculated in this study are not expected to differ signi®cantly if blood¯ow were considered.
van Gemert and Welch [21] demonstrated that the value of d in Equation (1) depends on À o o /d. For large values of À, d is the penetration depth of light. For small values of À, d is better represented by o o . Many laser-based therapeutic procedures involve conditions in which À is large, and so d is the penetration depth. To minimize confusion, in this study we assumed that d is the penetration depth for all values of À. Fig. 2. Log-log plot of effective thermal relaxation time for a gaussian beam pro®le vs. effective relaxation time t eff for a¯attop beam pro®le. All points fall below the solid line, indicating that t eff for a gaussian beam is less than that for a¯attop beam. For this plot, t eff is time required for peak temperature to decrease to 10 C above baseline. Figure 2 suggest that for a given parameter set, thermal relaxation is a faster process for incident gaussian beams than for¯attop beams. This can be explained by a discussion on radial temperature gradients. Normalized, representative temperature pro®les at the end of laser irradiation are shown in Figure 6 for incident gaussian and¯attop beams that are 2.5 mm in radius. Initially, at r 0, no net radial heat conduction will occur for thē attop case because the temperatures at adjacent radial positions are equal in magnitude. On the other hand, for the gaussian example, net heat¯ow away from the central axis occurs almost immediately because adjacent regions (e.g., r b 0) are at lower temperatures. The discrepancy between t eff for gaussian and¯attop beams increases as À gets larger because, for a given m a , as spot size increases, the region of T T peak during¯attop-beam irradiation increases as well. For a gaussian beam, T peak occurs only at the central axis and its location is insensitive to changes in spot size.

Data in
Schomacker et al. (7) observed that heat required for ablation during CO 2 laser ablation of in vitro guinea pig skin increased from 4.8 to 6.6 kJ/g as o o decreased from 680 to 250 mm. Our model results ( Table 2) show that for a given m a , t eff (2-D) decreases as o o gets smaller, suggesting that heat loss from the irradiated region is faster for smaller beam diameters. This hypothesis offers a reasonable explanation for why a decrease in o o leads to an increase in heat required for ablation.
A plot of t eff (1-D) calculated using Equation (5) as a function of t r (Fig. 3) shows that a 1-D thermal model  predicts t eff (1-D) b t r for the range of t r considered in this study. The linear relationship between t eff (1-D) and t r suggests that a simple multiplicative factor can be used to calculate t eff (1-D). The value of this factor depends solely on ÁT peak and ÁT base and was calculated for six combinations of ÁT peak and ÁT base (Table 3). Further derivation of Equation (5) demonstrates the reason for this linear relationship. For all times t b 0, the peak temperature induced by laser light absorption occurs at the tissue surface because the following conditions were assumed during model runs: (1) light scattering was negligible and 2) convective cooling was not considered. Thus, at z 0 and at t t eff , Equation (5) becomes In terms of the models, ÁT(0,t eff ) is equivalent to ÁT base and the quantity H o m a /rc equivalent to ÁT peak . The term erfcx(b) takes into consideration the dynamic change in peak temperature over time. Thus, Equation (7) takes on the following form ÁT base ÁT peak erfc xbX 8 Using Equation (1) and t t eff , b can be rewritten as b t t eff m a a t eff 1a2 1 2 Thus, b represents the ratio t t eff /t r . Since erfcx (b) is a monotonically-decreasing function as b increases, for a given set of values for ÁT base and ÁT peak , there exists a single value for b that satis®es Equation (8). In this situation, b is a constant, so the ratio t eff /t r is also constant. This relationship is true for any combination of optical and thermal properties. Table 3 provides examples of the linear relationship existing between t eff and t r for six combinations of ÁT base and ÁT peak . A comparison between t eff values calculated from the 1-D and 2-D models (Fig. 4) indicates that the 1-D model is comparable to or overestimates t eff (2-D). The amount of overestimation depends on À (Fig. 5). As described above, for larger spot sizes, radial heat conduction is not as signi®cant as for smaller spot sizes. For À b 10, t eff is equivalent for 1-D and 2-D model calculations; this condition is generally satis®ed for incisional or ablation procedures. These results show that the 1-D model can at least provide a conservative estimate of the time required between pulses to minimize potentially deleterious effects of temperature superposition. A worst-case value for t eff can be calculated using a wide range of thermal properties (Equations 4 (a±c), Table 1). Equation (5) is more complex than the equation for t r (Equation (1)), but it can still be solved relatively easily using conventional spreadsheet programs.
Since t eff is linearly related to t r , the accuracy of both values suffers from the same limitation in that accurate knowledge of local optical and thermal properties is required. In particular, to calculate b, values of m a and a must be known. Since t eff is longer than t r , the relative error associated with the uncertainty in a given set of optical and thermal properties is smaller for calculation of t eff than for t r . It is necessary to keep in mind that both t eff and t r are guidelines. The primary purpose of t eff is to provide an estimate of the time required for heated tissue to cool from T peak to a value ÁT base above baseline, whereas t r represents the time required to cool from T peak to a 1/e value.

CONCLUSIONS
Theoretical calculations showed that the concept of thermal relaxation time does not adequately provide a measure of time that can be used during multiple-pulse irradiation of tissue. It is necessary to consider the slow temperature decay that occurs for times longer than t r . The following equation can be used to estimate t eff .
where m is the slope of the line ®tting t eff (calculated with Equation (5)) as a function of t r . Calculated values of m for six combinations of ÁT peak and ÁT base are provided in Table 3; it is straightforward to calculate m with Equation (5) for other combinations.