Modeling of thermal injury produced by laser ablation of biological tissue

A thermal model to predict the effects of laser parameters on the zone of thermal injury produced by laser ablation of biological tissue is presented. A dimensionless parameter based on the ablation velocity and the optical and thermal properties of the target is key in determining the resulting zone of thermal injury. The zone of thermal injury is minimized when this parameter, known as the Peclet number (Pe), is much larger than one. This occurs because the rapid movement of the ablation front prevents the diffusion of energy beyond the laser that absorbs the laser radiation. For Pe less than one, the slow movement of the ablation front allows for diffusion of energy away from the region of energy deposition and leads to larger zones of thermal injury. The model predictions are compared with data available in the literature. Deviations between the model predictions and published data are discussed and potential effects of pyrolysis, temporally varying pulse shapes and pulse repetition rates are explored.


Introduction
As the use of lasers for surgical P11P05e5 is becoming more prevalent. there is increasing interest in predicting the extent of therma.1 injury surrounding laser incisions. The photothermal effects produced by laser ablation are important since they result in hmotasis when a zone of tissue adjacent to the incision is coagulated and serves as a barrier to blood perfusion. This is valuable in many surgical applications dealing with soft and highly vascula.r organs ( the spleen, for example) which do not lend themselves easily to clamping, suturing or electrocoagulation {5}. However, this zone of coagulated tissue is no longer viable and often delays the healing process and leads to significant scarring [2,3,12]. In other applications, such as photorefractive keratectomy (PR.K). it is essentia.l to minimize the thermal effects which accompany the ablation process. As a result, it would be useful if the extent of the zone of thermally injured tissue could be predicted. In order to do this, two things must be determined: the temporal and spatial profile of temperatures developed within the tissue and the mechanism h vliich the laser-induced temperature changes j)rO(luCe therma.l injury.
This paper presents a model which estimates the zone of thermal injury produced by the laser cutting process. The model is not an attempt to describe the dynamics of the ablation process. Instead, the model focuses on aspects of the ablation process which are relevant to the production of thermally damaged tissue. We do this by determining the temperature profile established within tissue during and after laser ablation. The zone of thermal injury is found by calculating the maximum depth to which particular isotherms propagate within the tissue. The results of the model are then compared to data in the literature.
Many investigators have attempted to model the ablation process in biological tissue [1, 14,16,26]. These models assume that tissue has thermophysical properties similar to that of water and that the ablation process is essentially one of vaporization. We will take a similar approach in this model even though the process of tissue ablation is fundamentally different from and much more complex than a. simple phase change. This is because the solution of a. phase change 1)roblem if appropriately constructed. ca.n identify the parameters which control the ablation process and its effects. The model presented here hopes to identify the parameters which have the largest impact on determining the zone of thermal injury resulting from the ablation process. We make the following assumptions.
( a) The ablation process can be treated as a vaporization process.
(b) Explosive ejection of tissue without phase change is not significant.
( c) The optical and thermal properties of the tissue are similar to that of water and are constant throughout the ablation process.
( d) Energy losses from the tissue via convective and radiative transport are negligible.
( e) Beer's law accurately describes the distribution of laser radiation within the tissue.
(f) Absorption and scattering of the incident laser radiation by the ablation plume are negligible.
( g) There are no temporal variations in the irradiance of the laser pulse.

Model Development
We begin by determining the temperature distribution established during the ablation process. To do this the steady-state ablation equation in one dimension is considered [18]: where T is the temperature, z the axial position within the tissue, Vabi the velocity of the ablation front, a the thermal diffusivity of the tissue, q"(z) the distribution of energy deposited by the laser, and k the thermal conductivity of the tissue. The following dimensionless variables are introduced into eqn.
After introducing a dimensionless heat flux, ,b = q"/kji(T -T) (Pe -1), the solution to eqns. (2) and (3) is: In order to apply eqn. (4), the Péclèt numl)er must be determined. From energy conservation, it is easily seen that Vabi = q"/p [c (Tm -I) + h19}. This leads to the following expression for the Péclèt number: bSf Pe=5f Sf 1' where Sf = Cp(Tm -T)/hj9, the Stefan number. If the laser pulse is sufficiently long, eqn. (4) gives the temperature distribution within the tissue just prior to the end of the laser pulse {22]. After the end of the laser pulse, the temperature profile decays. Since it is assumed that a negligible amount of energy is lost via convection or radiation through the tissue surface, the subsequent behavior of the temperature profile is given by the solution to the one dimensional heat diffusion equation, shown below in dimensionless form: d20 dO d(2 = dr' (6) where r is a dimensionless time variable, r = cit. The appropriate initial condition is given by eqn.
Since we assume that cooling at the tissue surface is negligible, the following boundary conditions apply: (( = 0, r) = 0 and O(C ' oo T) ' 0.
where erfc(x) represents the complementary error function. Eqn. (9) gives the behavior of the temperature distribution for all positions and times after the laser irradiation is completed.
The foregoing analysis assumes that the tissue heating occurs isochorically. That is, we assume the tissue undergoes no volume expansion prior to vaporization. Alternatively, one could assume that the tissue expands isobarically as it is heated. This would result in the formation of a layer of mixed phase at the tissue surface a.s has been postulated by some investigators [26,27,29]. Both formulations represent important limiting cases since some tissue expansion is expected during the ablation process and has been observed prior to the onset of material remova.l [6,13,23]. However, this expansion is never sufficient to achieve an isobaric process and the actual ablation process is somewhere in between these two formulations. The analysis of the isobaric case is complicated by the formation of the vapor-liquid layer and has been considered in detail [22].

Determination of the Zone of Thermal Injury
Thermal injury of tissues can occur by several means: damage to the genetic apparatus of cells, enzyme inactivation or protein denaturation [4]. Studies of thermal injury resulting from laser ablation have focused on collagen denaturation in tissues such as skin [24,25,28], cornea [17] and cartilage [15]. Although collagen denaturation is a thermally activated rate process, there exists a set of conditions where the process will proceed when a threshold temperature is reached and where the threshold temperature is nearly independent of the duration of the thermal transient. This approach has been taken by several investigators and will also be adopted in this work {22, 24,28]. The conditions under which this approximation is valid will now be investigated.
The process of thermal injury can be thought of as one which transforms a complex Ato another complex B by means of a chemical reaction that passes through an activated state A* which lies at a Gibbs free energy larger than both states A and B. Since the temperature affects the equilibrium concentrations of A and A*, the process A -+ B is temperature dependent. In order to determine the zone of thermal injury we wish to determine the amount of complex A which is converted to B at any point within the tissue when given the temperature history at that point. Using the theory of absolute reaction rates as proposed by Eyring one can derive the following expression (known as the Arrhenius damage integral) for the concentration of complex A, x, remaining at any time during the temperature transient [8]: t kT(t') Eqn. (11) indicates that both enthalpic aid entropic contributions are important when considering the kinetics of such rate processes. For protein denaturation both the heat and entropy of activation are quite large [8]. In this limiting case iH -TAS becomes the difference between two large numbers and the exponential term in eqn. ( 1 1 ) becomes a very strong function of temperature. This permits the definition of a critical temperature, T, such that the region of denaturation can be taken as all points where this temperature is reached [11]. In this limiting case, T is approximated by the implicit relation [11]: R1n(kTt/h)+iS (12) where te S the time over which the region is exposed to T. When considering thermal injury in skin, it has been found that when the Arrhenius damage integral is written in the form Q(t) = A f exp [-LE/RT(t')] dt', A = 3 x 10 and LE = 6.3 x i05 J rno11 [10]. From these values we can extract an activation enthalpy and entropy of LH = 6.3 x i05 J rnol and LS = 1.64 x i0 Jmo11 K1 , respectively. From these values we can calculate the critical temperature T versus the exposure time, te . The result of this calculation is shown in Figure 1 . The result of this calculation indicates that the critical temperature is indeed weakly dependent on the exposure time. Thus, studies which have investigated the thermal injury in skin by subjecting them to temperatures in a heat bath for relatively long times ('-' 102 s) should yield critical temperatures which are also applicable on the time scales relevant for these ablation processes. Such studies indicate that this critical temperature lies in the range of 65 and 73°C [24,28] which corresponds to ® = 0.55-0.64 if an initial temperature of 22 °C is assumed. Thus, for the temperature profiles that have been calculated, the depth of thermal injury will be taken as the maximum dimensionless depth ( to which the O isotherm has propagated.  . Note that an increase in Péclèt number increases the slope of the temperature profile near the surface which is necessary to supply the flux of energy for vaporization. As the Péclèt number decreases the slower movement of the ablation front allows for significant diffusion of energy beyond the optical penetration depth resulting in deeper penetration of the temperature profile. In addition, as the Péclèt number decreases, the amount of energy stored in the temperature profile increases. This is seen by integrating eqn. (4 ) to produce a dimensionless measure of the energy contained in the temperature profile. Doing this we get: 0 Pe This result shows that as Pe ' the energy contained in the temperature profile approaches a minimum, normalized value of 1. As the Péclèt number is reduced, the amount of energy in the tissue bulk increases due to diffusion. Pc = 1 represents the point where half of the energy in the tissue bulk arrived via direct energy deposition and the other half via diffusional transport. Therefore, under steady-state conditions, both spatial confinement of the temperature profile and minimization of the energy remaining in the tissue bulk is achieved at high Péclèt numbers. Thus, one would expect the zone of thermal injury would be minimized at high Péclèt numbers. Figure 3 plots the model prediction for the dimensionless zone of thermal injury as a function of Péclèt number. The figure also plots data from the literature [7. 9, 15. 17. 19. 21, 24] for zones of thermal injury produced by severa.1 pulsed and continuous wave laser sources in various tissues whose primary constituent is collagen. The laser parameters used in the studies from which the data were taken are listed in Table 1. Several important observations can be made when comparing the model predictions with the data. First, for Pe 1 the model overestimates the damage produced by both CO2 lasers and erbium lasers. Secondly, in this same range of Péclèt number the model considerably overpredicts the damage zones created by Ho3+ :YAG and Tm3+ :YAG lasers. Lastly, for Pe 1 the model underestimates the zone of thermal injury.
In general, the lack of agreement between the experimental data and the model predictions are probably due to effects not considered by the model. These effects can be classified into four categories: (a) assumptions within the model relating to the geometry and boundary conditions of the process, (b) scattering of the incident laser radiation which may significantly alter the distribution of laser energy deposition, (c) pyrolysis which would violate the premise of the vaporization analysis and (d) laser pulse shapes whose irradiance varies SPIE Vol. 1882 Laser-Tissue Interaction IV(1993) / 7  Table 1 according to the number assigned to each data point.

Effect of Model Geometry and Boundary Conditions
For Pe ;: 1 the model overpredicts the zone of thermal injury in almost all cases. In the model, we consider the ablation process and the transport of energy in the tissue bulk in one dimension rather than considering the effects of a finite laser beam. These effects include thermal diffusion in two dimensions (if cylindrical symmetry is assumed) and the geometric changes associated with crater formation. The simplification to one dimension has the effect of reducing the temperature gradients thereby increasing the distance over which the energy is dissipated and increasing the predicted zone of thermal injury. Further, when considering the evolution of the temperature profile following ablation we enforce an adiabatic boundary condition at the tissue surface. That is, we do not consider the possibility of energy transfer by convection (e.g., by air, saline or blood) or by mass transfer (i.e., post-pulse material removal). This ensures that all the energy contained within the temperature distribution is dissipated via conduction into the tissue bulk. This maximizes the extent to which this energy can affect thermal injury. Thus, both the geometry and the boundary conditions employed in the model tend to overestimate the zone of thermal injury for Pe 1.

Effect of Optical Scattering
For 11o3+:YAG and Tm3+:YAG lasers it may be inappropriate to consider ;' as the characteristic optical penetration depth. At these wavelengths scattering of incident laser radiation within tissue may be significant. Further, the effect of scattering may be enhanced because the spot. sizes used in the experiments are comparable to ç1 . It is known that for a given value of /La, an increase in the scattering coefficient reduces the average depth to which laser radiation is deposited [20]. If sufficient scattering is present , a coefficient representing the inverse of the optical penetration depth, say Peff , should be used in place of ,ua S the length against which the variables should be normalized against. If scattering is significant at the Ho3+:YAG and Tm3+:YAG wavelengths, the value of the I-eff would be larger than jua. This would increase the values of the dimensionless zone of thermal injury found through experiment and reduce the error between the experimental observations and the model predictions. Quantification of this effect on the model predictions would require a measurement of the optical properties of tissue at the appropriate wavelengths, (2.01 and 2.12 urn) followed by a reworking of all the model calculations using a modified mathematical expression for the distribution of the volumetric heat generation.

Effect of Pyrolysis
The observed zone of thermal injury is greater than the model prediction for ablation events in which char was formed. The presence of char indicates that pyrolysis occurred during the ablation process. Pyrolysis should result in damage zones that are larger than the model predictions for two reasons. First, with pyrolysis the tissue is subjected to significantly higher temperatures (up to ''1OOO°C) than those during vaporization [23] which greatly increases the magnitude of the temperature profile computed earlier. Secondly, pyrolysis is a less efficient mechanism for material removal than vaporization and reduces the ablation velocity, which increases the amount of energy which remains in the tissue bulk and fails to achieve the confinement of the temperature profile which occurs at large Péclèt numbers (shown in Figure 2). These effects are not considered in the model and may result in underestimates of the zone of thermal injury relative to these experimental observations.

Effect of Temporal Pulse Shape
The zones of thermal injury reported by Ren et al. [17] for the free running Er3+:YAG and Er3+:YSGG laser sources at Pe < 1 are also higher than predicted by the model. In these experiments the temporal pulse shape (200,us, full width at half maximum {FWHM]) i not constant. Because the irradiance varies during the laser pulse, the ablation velocity and the Péclèt number are not constant during the ablation event. The calculated Péclèt number is obtained using an 'average irradia.nce" where the pulse energy is divided by the FWHM pulse duration and the spot area. In general, the difference between the true, time varying Péclèt number and the average Péclèt number must be considered only when the following criteria are satisfied: S The pulse duration, t, is of the same order or larger than the characteristic thermal diffusion time i.e., a/atp (9(1).
. The average Péclèt number is small enough that changes in Péclèt number during the pulse would lead to substantia.l changes in the predicted zone of therma.1 injury as shown by Figure 3. That is, the computed Péclèt number should be (9(1).
In the case of Ren's experiments both criteria. are satisfied and the temperature distribution given by eqn. (4) is inaccurate because no steady state temperature distribution is established during the pulse. This is significant because the characteristic decay length of the temperature distribution at the end of the laser pulse is not known. If the tail of the laser pulse is on the same order of the characteristic thermal diffusion time and the Péclèt number during the tail of the laser pulse is lower than the average Péclèt number during the pulse, the actual penetration of the temperature profile (and resulting zone of thermal injury) would exceed the prediction made by the model. This results in a predicted zone of thermal injury which underestimates that which is seen experimentally [17].

Effect of Pulse Repetition Rate
The introduction of "superpulsed" CO2 lasers having the ability to deliver a succession of short high intensity pulses at varying repetition rates, has stimulated interest into understanding how pulse repetition rate may affect the zone of thermal injury. The repetition rate at which cumulative thermal effects can be expected should correspond roughly with the reciprocal of the time at which the 0 = 0.55 isotherm propagates to the damage depth, -rd. For Péclèt numbers much larger than 1, Td = 32. For a CO2 laser, this corresponds to a time of 40 ms assuming ,ua 77 000m1 and c = 1.3 x iO m2 s1. This in turn predicts a pulse repetition rate of 25 Hz at which cumulative thermal effects should be seen. Experimental work [21] confirms that increases in thermal damage are first seen at repetition rates near 20 Hz and continues to increase until a repetition rate of roughly 200 Hz is reached.

Conclusions
A thermal model has been presented to predict the effect of various laser parameters on the zone of thermal injury produced during the laser ablation of tissue. The model predicts that for Péclèt numbers larger than one, the temperature distribution is confined to a thin layer in the tissue and minimizes the zone of thermal injury. The model is in qualitative agreement with experimental data.
There is much room for improvement with respect to both the determination of the thermal transients which accompany the ablation process and modeling of the processes thought to produce the thermal injury. Concerning the issue of thermal transients, future models should consider ablation geometry and boundary conditions and also account for the effects of pyrolysis, optical scattering of the incident laser radiation, and complex temporal pulse shapes which no doubt play a significant role in affecting the temperature profile in the tissue bulk. This can only be done once the physical mechanisms controlling the ablation process are determined and understood. With respect to modeling the damage process efforts should be made to determine the enthalpy and entropy of activation for thermal damage in tissues such as cornea and cartilage and existing numerical codes should be employed to determine the zones of thermal injury once the temperature profiles generated from these ablation processes are determined4