1.1.

History of the Damage Integral (Ω)

In 1947, Henriques and Moritz published a series of systematic studies on the thermal damage kinetics of skin.^1–4 The authors applied heated water or oil on the surface of human and porcine skin using a special manifold, and looked for redness over a wide range of temperatures (44°C to 60°C) and exposure durations (3 s to 5 h). Their damage results showed an exponential relationship to temperature and a linear correlation to time, indicating a first order reaction for the “temperature-time history.” Adopting the first order rate law of Arrhenius,⁵ which describes the kinetics of isothermal (constant temperature) chemical reactions, Henriques and Moritz determined their energy of activation ($E_a$) was similar in value to the energy of denaturation for most purified proteins. To provide a measure of accumulated thermal damage ($\mathrm{\Omega }$), Henriques integrated the Arrhenius rate law over time ($d{t}^{\prime }$), using the isothermal peak temperature ($T_p$) in kelvins, to produce the isothermal damage integral shown as

The same issue came into play after the invention of the laser,⁶ when research on photothermal damage rate kinetics was accelerated. Laser exposure of absorbing tissue results in a rapid temperature rise that might seem instantaneous, but with proper temporal resolution, a ramp in temperature is evident. The nonisothermal nature of this initial photothermal temperature rise requires relatively long exposure durations to approximate an isothermal condition, just as was the case for Henriques and Moritz. Without an alternative, researchers in the field of laser–tissue interaction adopted the Arrhenius isothermal method described by Henriques decades earlier.

To obtain the Arrhenius rate parameters, $A$ and $E_a$ for calculating $\mathrm{\Omega }$ values, an accumulated damage level of unity ($\mathrm{\Omega }\text{value}=1.0$) is used to indicate the damaged state, and Eq. (1) can be linearized and rearranged to as

To illustrate the nonisothermal nature of short laser exposures we present, as reported in a previous paper,⁷ an averaged thermal profile (temperature versus time) for a 0.25-s exposure at $2\mu \mathrm{m}$ [Fig. 1(a)], and how that data were used to graph [Fig. 1(b)] Eq. (2). The dashed lines in Fig. 1(a) show the square function of temperature and exposure duration, and the arrow indicates the region of the isothermal assumption that causes an overestimate of the thermal dose delivered to the cells. The shape of the thermal profile dictates the degree of overestimation, and this area decreases as a percentage of the total area with extended exposure duration. From the straight line given by Eq. (2), the slope represents $\frac{E_a}{R}$ and the $y$-intercept gives $-\mathrm{Ln}{A}_{\mathrm{iso}}$ [Fig. 1(b)]. Figure 1(a) clearly shows why the Arrhenius integral, with constant temperature assumption, is suited for isothermal conditions, such as enzyme reactions where a small volume of enzyme is added to a reaction mixture that is already thermally equilibrated. Of course, the significance of this overestimation in photothermal damage can only be determined by devising an alternative approach that subtracts this overestimate, which is one purpose of this paper.

Fig. 1

Comparing isothermal and nonisothermal heating. (a) Average threshold temperature-rise profiles for a 0.25-s laser exposure ($2\text{-}\mu \mathrm{m}$ at 30°C ambient temperature) showing the square function (dashed lines) used to determine the Arrhenius rate parameters ($A$ and $E_a$). Arrow identifies region of overestimated temperature history used in typical Arrhenius method. Inset to (a) identifies the boundary of cell death, which is the source of the “threshold” thermal profile data. (b) Arrhenius plot to empirically determine $A_{\mathrm{iso}}$ (Exp[-$y$-intercept]) and $E_a$ (slope x $R$) using Eq. (2). Dashed arrow shows where the square function of panel (a) falls, as part of the average value on the plot of Eq. (2). Portions of (a) and all of (b) taken from Figs. 5(c) and 5(f) of Ref. 7, respectively. Other thermal profiles used in this study are provided in Fig. 2.

Six decades after the invention of the laser, no alternatives have been described that can match the damage prediction capability of the isothermal damage integral. Better methods for measuring and recording temperature responses of biological tissues are now available, such as thermal cameras with reasonable spatial and temporal resolution. Although advancing computational physics models that accurately predict temperature rises based on first principles⁸ is important, by focusing less on empirical measurements to support or refine the quality of damage prediction by the isothermal damage integral, there remains no alternative. Hence, no new data are available that address not only the isothermal, but several other known drawbacks⁹ of the damage integral, including the assumptions that the temperature and exposure duration dependence of the frequency factor are negligible. For a methodical evaluation of the significance of these drawbacks, one needs an experimental system that allows the collection of thermal data that reproducibly generates the same accumulated damage ($\mathrm{\Omega }$ value), and a nonisothermal mathematical model that eliminates the overestimated portion of the thermal dose. The first requirement was described in a recent publication,⁷ and we will use that previous data to highlight our newly derived nonisothermal model.

To understand the significance of the data previously reported, it is important to understand the term minimum visible lesion, or MVL. The production of a photothermal MVL requires a minimal amount of laser energy to generate the threshold temperature rise leading to damage in a given cell or tissue type. However, due to differences in optical and thermal properties between similar samples, and even within the same sample, delivering the same “threshold” laser irradiance in repeated exposures will likely lead to three types of results; no damage, MVL damage, or damage that extended in size beyond the MVL stage. Thus, the term threshold laser dose, whether irradiance or radiant exposure, is less informative than a threshold temperature history. Only the MVL data are useful because of the desire to assign each exposure an $\mathrm{\Omega }$ value of 1 for establishing relationships between laser exposure duration and ambient temperature ($T_{\mathrm{amb}}$). This would lead to inefficient data sets. Here, we denote the accumulated thermal history leading to a minimum visible lesion as a threshold thermal dose. A threshold thermal profile defines a threshold thermal dose, as experienced by the tissue during an MVL exposure or equivalent. The threshold thermal dose can be quantified, such as by the Arrhenius integral as a function of the thermal profile ($\mathrm{\Omega }$) to arrive at an accumulated damage value. An $\mathrm{\Omega }$ value of unity ($\mathrm{\Omega }=1$) continues to be used routinely as the amount of accumulated damage barely (threshold) causing death.

1.2.

Damage Rate Processes

Laser-induced thermal (photothermal) damage is a complex, multiphasic cellular response to increasing temperature over a specified time. The progression of damage involves thermal transitions (unfolding) of macromolecules, especially proteins, to the point of losing biological function. As the accumulation of nonfunctional biomolecules continues, the ability of the cell to repair the thermally induced stress diminishes until there is irreversible damage. As shown for whey proteins in solution, the rate of heating can influence the structure and solubility of proteins.^10,11 Of particular interest is the notion that unfolded proteins have some chance of refolding, unless they begin to aggregate and become irreversibly “denatured.” The reversible protein unfolding process is considered a first-order reaction. deWit and Klarenbeek¹⁰ found that, although their whey protein began to unfold at 70°C for both heating rates of $1°\mathrm{C}{\min }^{-1}$ and $1°\mathrm{C}{\mathrm{s}}^{-1}$, the temperature at which the protein started aggregating (also 60% unfolded in both groups) was drastically different. The slower heating rate produced 60% unfolded protein at 73°C, and the faster heating rate produced the 60% unfolded fraction at 85°C. Thus, for temperature rises, there are two kinetic rates (unfolding and aggregation) existing concurrently, depending upon both the peak temperature achieved and the rate of heating.

This example provides kinetics data from slow processes (0.17 to $1°\mathrm{C}{\mathrm{s}}^{-1}$) relative to the rates of temperature rise produced by laser exposures. For instance, the exposure shown in Fig. 1(a) produced a heating rate of $120°\mathrm{C}{\mathrm{s}}^{-1}$ (30°C rise in 0.25 s). Considering all the $2\text{-}\mu \mathrm{m}$ laser heating data from Denton⁷ et al., the range of heating rate in their average threshold thermal profiles was from $7.5°\mathrm{C}{\mathrm{s}}^{-1}$ (40°C ambient to 55°C peak in the first 2-s of a 20-s exposure) to $800°\mathrm{C}{\mathrm{s}}^{-1}$ (20°C ambient to 60°C peak in 0.05 s). If we extrapolate the trend described by deWit¹¹ above to the photothermal data, the amount of denatured protein would be low by the end of these short exposures. Perhaps cells are sensitive to a small amount of inactivated protein, especially a protein with a critical function. However, the extremely high protein concentrations within living cells relative to the test tubes of deWit would significantly increase the rate of aggregation of proteins once unfolded.¹² This highlights that, although trends found in test tubes reveal important features for live cells, direct comparisons and interpretations require caution.

2.1.

Derivation of Nonisothermal Damage Integral

Due to the complex physical changes in state and chemical aggregation leading to damage, approaching the kinetics as an overall composite system was an alternative starting point to the isothermal chemical reactions of Arrhenius.⁵ The currently accepted kinetic equation for a heterogeneous condensed phase system¹³ is

The rate of conversion of the measured variable, $\alpha$, was assumed to equal a separable mass-loss function, $f(\alpha )$, and a temperature dependent rate function, $k(T)$.^13,14 Think of $\alpha$ as the property of the material that was converted during the temperature dependent process. In laser bioeffects, this would be equivalent to the one or more proteins responsible for cellular damage when denatured. The degree of conversion of $\alpha$ (0.0 to 1.0) that causes cell damage can be arbitrary. Due to the complexity of the cellular components, including our lack of knowledge about the actual protein(s) involved in the damage rate process, $\alpha$ becomes impossible to measure. Therefore, assuming a separable $f(\alpha )$ and $k(T)$ functions, the heterogeneous condensed phase system provides a good starting point. Considering the familiar Arrhenius rate equation to describe the temperature dependent rate function yields

For variable heating rates, separation of variables and integration of Eq. (5) yields a measure of converted sample ($g(\alpha )$):

This formulation was an alternative derivation of the damage integral ($\mathrm{\Omega }$) [Eq. (1)] used by Henriques and Moritz and others.^4,7,15–21 Classically, threshold thermal damage occurs when the integration in Eq. (6), with respect to $\alpha$, generates a value of one at some time during a nonisothermal insult. The resulting equation is provided as

Operating Eq. (7) on thermal data shown in Fig. 1(a) entails the region below the line of the profile and represents the nonisothermal heating history. To quantify the region of overestimation shown in Fig. 1(a) (arrow), one can subtract the region under the profile from the region representing the square function (dashed box). To express the isothermal square function as Eq. (8), we replaced $T(t^{\prime })$ in the exponent of Eq. (7) with peak temperature ($T_p$), used the bulk time frame of τ, and set the equation equal to 1.0:

We correct for these assumptions using an exposure duration dependent correction factor $\delta c$ ($\tau$), that also takes into account the variable heating rate. In essence, we define ${\Omega }_{\mathrm{non}}(\tau )$ [Eq. (7)] as the false square function [Eq. (8)] minus a correction factor, $\delta c$ ($\tau$), as shown in Eq. (9):

However, researchers do not integrate the square function when calculating ${\mathrm{\Omega }}_{\mathrm{iso}}$. Instead, modelers predict damage by integrating nonisothermal profiles using the $A_{\mathrm{iso}}$ obtained from the isothermal function [Eq. (2), Fig. 1(b)], as depicted in Eq. (10):

Because of the fact that $E_a$ is the same value when solved via Eq. (2), rearrangement provides the simplified definition of ${\delta }_c(\tau )$

Equation (11) indicates that the correction factor is applicable to the common usage of $A_{\mathrm{iso}}$ to determine the damage integral used by the laser research community (${\mathrm{\Omega }}_{\mathrm{iso}}$) and does not correct the entire step function [Eq. (8)]. As such, ${\delta }_c(\tau )$ is considerably smaller than that if modelers used Eq. (9). Clearly, the correction in Eq. (10) is entirely a function of the frequency factor correction, where $A_{\mathrm{iso}}$ is smaller than $A_{\mathrm{non}}$, which consistently leads to an underestimate of damage (${\mathrm{\Omega }}_{\mathrm{iso}}$ values $<{\mathrm{\Omega }}_{\mathrm{non}}$) by the modelers. When using threshold thermal profiles, where the value for ${\mathrm{\Omega }}_{\mathrm{non}}$ is defined as 1, the correction factor simplifies further to ${\delta }_c(\tau )={\mathrm{\Omega }}_{\mathrm{iso}}-1$, reinforcing the fact that the correction factor as defined in Eq. (10) is a negative value.

Upon additional rearrangements, we find the relationship between ${\delta }_c(\tau )$ and frequency factors, $A_{\mathrm{iso}}$ and $A_{\mathrm{non}}$ as

Substituting for $A_{\mathrm{iso}}$ in Eq. (2) yields

2.2.

Determination of $A_{\mathrm{non}}$ and ${\delta }_c(\tau )$ Using Individual Threshold Thermal Profiles

As a starting point for determining values for ${\mathrm{\Omega }}_{\mathrm{non}}$, we use the simple rearrangement of Eq. (7) to yield

2.3.

Source of Threshold Thermal Profiles to Calibrate Nonisothermal Damage Integral

As described in the previous section, once an $E_a$ value for a given $T_{\mathrm{amb}}$ is known, individual thermal profiles can provide values for $A_{\mathrm{non}}$ using Eq. (16). In contrast to the traditional plotting of Eq. (2) over a series of laser exposure durations, the use of individual thermal profiles enabled the determination of mean $A_{\mathrm{non}}$ and $A_{\mathrm{iso}}$ values with standard deviations within the same exposure duration. Ideally, “threshold” thermal profiles would be used to determine the Arrhenius rate parameters because then, by definition, the damage integral would solve to a value of unity. We have postulated that, as long as the damage area was smaller than that of the laser footprint, cells at the boundary of damaged and nondamaged cells in exposed cultures should represent the MVL condition [white line in Fig. 1(a), inset].⁷ If this were true, it would open the door for generating threshold thermal profiles using a wide range of laser irradiances and damage areas. Damaged monolayers would provide MVL data without being limited to a small minimum lesion in the center of an exposed area.

Our method for determining the temperature history at the boundaries of damage, termed microthermography,^7,15 uses high-speed and high-magnification thermal imaging, followed by an image overlay with postexposure fluorescence images indicating damaged areas. Spatially resolved temperature responses of nonpigmented retinal pigment epithelial (RPE) cells to $2\text{-}\mu \mathrm{m}$ laser irradiation were recorded at 800 fps with an InSb-based camera (ThermoVision SC6000, FLIR Systems, Inc., Boston, Massachusetts), with an effective pixel pitch of about $8\mu \mathrm{m}$ at the image plane. Custom LabVIEW image processing routines performed the overlays with the appropriate stretching, resizing, and orientation changes. Once overlaid correctly, our software identified the thermal pixels associated with the boundary by defining a single-pixel region of interest (ROI) around the perimeter of the damaged area [white line in Fig. 1(a), inset]. The final extraction of the thermal history of the boundary ROIs also used the suite of LabVIEW programs. Large damage areas, no $>90\%$ of the laser area, provided a larger number of useful thermal pixels per boundary ROI, which strengthened the statistical power of the method.

We chose 16 different sets of exposure parameters to investigate the relationship between threshold thermal dose of both laser exposure duration and the $T_{\mathrm{amb}}$ at which the cells were exposed. We investigated exposure durations of 0.05, 0.25, 1.0, and 20 s, which was much broader than our initial range of 0.1 to 1.0 s.¹⁵ Exposures at the four durations were delivered at four different ambient temperatures. The $T_{\mathrm{amb}}$ values chosen were 20°C, 25°C, 30°C, and 40°C. We ensured that the RPE cells were healthy throughout the duration ($<15\min$) of the experiments when held at these ambient temperatures. Multiple recordings (12 to 20) of thermal responses of cells in real-time with laser exposure provided the thermal data (Fig. 2 and Ref. 7) we used for each of the 16 different combinations of $T_{\mathrm{amb}}$ and $\tau$.

Fig. 2

Threshold thermal profiles (temperature rises) for nonpigmented RPE cells exposed at $2\mu \mathrm{m}$. Initial 2 s of thermal responses to $2\text{-}\mu \mathrm{m}$ laser exposure at ambient temperatures of (a) 20°C, (b) 25°C, (c) 30°C, and (d) 40°C. Panel (e) provides full thermal profiles for the duration of the 20-s exposures at each of the four ambient temperatures. Temperatures, taken from the boundaries of cell death across a wide range of laser irradiances and damage sizes, are represented as the mean value (with SEM bars) at each time point. Taken from Ref. 7.

In the prior paper, and using the typical isothermal method for determining $E_a$ and $A$, we verified that, within 10% experimental error, all 16 thermal profiles we extracted from the boundaries of cell death were equivalent to the same $\mathrm{\Omega }$ value.⁷ Due to the threshold nature of the thermal profiles, we assumed an $\mathrm{\Omega }$ value of 1.0 for each. Using empirically derived averaged values for $A$ and $E_a$, each of the 16 average threshold thermal profiles were corrected by scaling the temperature values at each time point by a small increment, until the $\mathrm{\Omega }$ value solved to 1.0. The scaling factors ranged from 0.6% to 10.3%. Therefore, the thermal profiles from our prior study are equivalent to those from an MVL exposure, allowing us to evaluate differences between isothermal and nonisothermal damage integrals in the current paper.

In preparing the 16 thermal profiles for use in this study, it became apparent that averaging the thermal profiles at each time point to generate one representative thermal profile, as shown in Fig. 2, was not optimum for data analysis. Here, we implemented the more rigorous strategy of a nearest-neighbor interpolation function, which determined numerical expressions of individual profiles. This generated only slight differences to values like $T_p$, critical temperature ($T_{\mathrm{crit}}$), and $E_a$ relative to the prior data, while providing superior statistical support to our new data. Here, we return to the fundamental aspects of the kinetic rate for a heterogeneous condensed phase system and derive a modified damage integral that removes the overestimated thermal dose. To highlight our results, we used the previously published thermal profiles to compare the differences in $E_a$ and $A$ values associated with nonisothermal treatment relative to the traditional isothermal method. Again, the thermal profiles provided threshold temperature history because they were derived from the boundary of cell death [Fig. 1(a), inset]. As such, we will set the nonisothermal damage integral (${\mathrm{\Omega }}_{\mathrm{non}}$) to unity. The variance from unity found previously for the isothermal damage integral (${\mathrm{\Omega }}_{\mathrm{iso}}$) values provided a measure of error. Another chief outcome of the nonisothermal derivation is that once temperature and exposure duration dependent $A$ and $E_a$ values are obtained for a given sample environment, one can predict threshold $T_p$ ($T_p^{\mathrm{thr}}$) and threshold $T_{\text{crit}}$ ($T_{\mathrm{crit}}^{\mathrm{thr}}$) without the need to integrate a thermal profile.

3.1.

Traditional (Isothermal) Convention for Determining Arrhenius Parameters

Traditionally, determining the Arrhenius rate parameters, $A_{\mathrm{iso}}$ and $E_a$, using the Arrhenius plot [Fig. 1(b)] of Eq. (2) requires measuring peak temperatures for multiple laser exposure durations. This gives one value for each parameter across the entire range of $\tau$ values within a given ambient temperature. Table 1 provides the values for $A_{\mathrm{iso}}$ and $E_a$ that were determined with an Arrhenius plot using the interpolated thermal profiles from Ref. 7. The laser exposure duration ($0.05\pm 0.002$, $0.250\pm 0.002$, $1.000\pm 0.002$, and $20\pm 0.1\mathrm{s}$) for each thermal profile was plotted versus the corresponding peak temperature [Eq. (2)]. The linearization of the data generated a single $A_{\mathrm{iso}}$ and $E_a$ parameter for each ambient temperature (Table 1). Both $A_{\mathrm{iso}}$ and $E_a$ values in Table 1 were greater than the values determined without interpolation.⁷ Differences due to interpolation had little effect (0.4% to 8.7%) on slope ($E_a$) values. However, due to slight changes in $T_p^{\mathrm{thr}}$ values, which propagated into slight changes in $y$-intercept values, there were large (three to five orders of magnitude) difference in the $A_{\mathrm{iso}}$ values between the prior and the current studies. These differences exemplify the large effect on $A_{\mathrm{iso}}$ with small changes in the $y$-intercept when taking the antilog function.

Table 1

Arrhenius rate coefficients from Arrhenius plot.

3.2.

Nonisothermal Convention for Damage Integral

3.2.2.

Threshold heating rate is reciprocally related to laser exposure duration

As indicated in previous sections, the frequency factor in a nonisothermal process ($A_{\mathrm{non}}$) changes as a function of temperature. Thus, it is not a constant parameter throughout the exposure duration, and it increases as a function of time and is expected to remain constant when the temperature reaches steady state. The values of $A_{\mathrm{non}}$ computed here were mean values for each exposure duration using threshold thermal profiles so they, along with the appropriate $E_a$ values, are suited for calculating damage rates (Arrhenius) and accumulated damage (nonisothermal damage integral) for the laser-induced damage processes in our cell culture model.

Rather than studying the complex relationships between $A_{\mathrm{non}}$ and both ambient temperature and laser exposure duration shown in Table 2, a natural extension is to examine the average heating rates (${\beta }_{\mathrm{avg}}=\frac{1}{\tau }{\int }_{T_{\mathrm{amb}}}^{T_p^{\mathrm{thr}}}\mathrm{d}T$). As expected, the heating rate was inversely proportional to both laser exposure duration and ambient temperature. Figure 3(a) shows that the heating rate needed to cause damage (threshold thermal profiles) depended greatly on $T_{\mathrm{amb}}$ and laser exposure duration. Average heating rates across $T_{\mathrm{amb}}$ (same $\tau$ value) were linear (lines not shown). Rather than providing equations for the ${\beta }_{\mathrm{ave}}$ versus $T_{\mathrm{amb}}$ in Fig. 3(a), we summarize the $y$ intercepts and slopes, with correlation coefficients, in Table 3. Both Fig. 3(a) and Table 3 indicate that ${\beta }_{\mathrm{avg}}$ must necessarily be larger at lower ambient temperatures to obtain the $\tau$-dependent peak temperature for threshold damage. Note that the dependence on $T_{\mathrm{amb}}$ (slopes, Table 3) was greatest for the shortest exposure durations.

Fig. 3

Threshold average heating rates in RPE cells exposed to $2\text{-}\mu \mathrm{m}$ laser light. (a) Bar graph showing the relationship between both ambient temperature and laser exposure duration with average heating rate. Bars of yellow, gray, orange, and blue represent exposure durations of 0.05, 0.25, 1.0, and 20 s, respectively. (b) Power function trends between average heating rate and exposure duration. Slopes (exponents of power functions) indicate a reciprocal relationship (${\tau }^{-1}$) between ${\beta }_{\mathrm{ave}}$ and $\tau$. Note that $T_{\mathrm{amb}}$ + prefactor = $T_{\mathrm{crit}}^{\mathrm{thr}}$. All correlation coefficients ($R^2$) values were unity.

The average heating rates within the $T_{\mathrm{amb}}$ groups [Fig. 3(b)] had a similar power function relationship to laser exposure duration, where each slope (exponent) indicated that ${\beta }_{\mathrm{avg}}$ is not just inversely proportional but reciprocally related to laser exposure duration (${\beta }_{\mathrm{avg}}\alpha$ ${\tau }^{-1}$). It is interesting to note that this reciprocal relationship extends to all four durations, regardless of whether or not the threshold thermal profile achieved thermal steady state [Fig. 3(b)].

The trends shown in Fig. 3(b) provide more information than just a reciprocal relationship. The prefix values for the power functions, by definition, are the values of ${\beta }_{\mathrm{avg}}$ at 1 s. When the ambient temperature (Table 1 mean value) was added to the ${\beta }_{\mathrm{avg}}$ at 1 s, we get our expected critical temperature, $T_{\mathrm{crit}}^{\mathrm{thr}}$ (threshold peak temperature at 1 s). In this way, the $T_{\mathrm{crit}}^{\mathrm{thr}}$ values derived from ${\beta }_{\mathrm{avg}}$ are 54.8°C, 55.6°C, 54.8°C, and 58.4°C for 20.2°C, 24.8°C, 29.6°C, and 40.0°C ambient temperature (Table 1), respectively. The $T_{\mathrm{crit}}^{\mathrm{thr}}$ values from our ${\beta }_{\mathrm{avg}}$ analysis from 20°C to 30°C $T_{\mathrm{amb}}$ were in line with the values (within the standard deviations) in Table 2 (peak temperature for 1-s exposures). The $T_{\mathrm{crit}}^{\mathrm{thr}}$ value of 58.4°C was only 0.3°C above the standard deviation range for the empirical value in Table 2. This result suggests that if one can empirically determine the $T_{\mathrm{crit}}^{\mathrm{thr}}$ value in any experimental system (cell culture, skin, retina) using threshold thermal profiles determined for 1-s exposures at a given $T_{\mathrm{amb}}$ (e.g., body temperature), one can determine the heating rate required for that 1-s exposure, as well as any other exposure durations by applying the reciprocal relationship to $\tau$. More experimental data are needed to confirm this relationship.

3.2.3.

Nonisothermal predictive model

Threshold peak temperatures

Methods to compare the predictive consequences of using $A_{\mathrm{iso}}$ rather than $A_{\mathrm{non}}$ include comparing ${\mathrm{\Omega }}_{\mathrm{non}}$ (e.g., 1.0) to ${\mathrm{\Omega }}_{\mathrm{iso}}$ values calculated from either threshold thermal profiles or the values for ${\delta }_c(\tau )+1$. As shown in a previous section, the $A_{\mathrm{iso}}$ is a valid approximation under steady state conditions only. However, the $A_{\mathrm{non}}$ is valid for all exposure durations and is thus applicable to isothermal and nonisothermal heating equally well.

Another metric to examine the validity of correcting for non-isothermal heating rates is the calculation of threshold peak temperatures ($T_p^{\mathrm{thr}}$), which is analogous to the “$T_{\mathrm{TH}}$” value derived by Pearce and Thomsen.¹⁹ Table 2 compares differences between empirical threshold $T_p$ values relative to the $T_p^{\mathrm{thr}}$ values calculated using $A_{\mathrm{non}}$ in the following equation, derived from Eq. (14):

Eq. (17)

$$T_p^{\mathrm{thr}}=\frac{E_a}{R\mathrm{Ln}(\tau ·A_{\mathrm{non}}({\delta }_c(\tau )+1))}.$$

Notice Eq. (17) depends upon the laser exposure duration. The $E_a$ value (energy requirement) drives the damage process to account for the $T_{\mathrm{amb}}$ by requiring the proper temperature rise. In addition, the $A_{\mathrm{non}}$ values influence $T_p^{\mathrm{thr}}$ values due to both $T_{\mathrm{amb}}$ and $\tau$ trends (Table 2). Similar dependencies are true for calculating the threshold critical temperature ($T_{\mathrm{crit}}^{\mathrm{thr}}$), which is Eq. (17) solved with an exposure duration of 1 s, similar to that derived by Pearce and Thomsen.¹⁹ Equation (18) provides the simplified expression for the critical temperature:

Eq. (18)

$$T_{\mathrm{crit}}^{\mathrm{thr}}=\frac{E_a}{R\mathrm{Ln}(A_{\mathrm{non}}({\delta }_c(\tau )+1))}.$$

Notice how Eqs. (17) and (18) allowed the prediction of these two important threshold temperature values without a thermal profile. Once the proper threshold $E_a$ and $A_{\mathrm{non}}$ values are obtained, using a statistically significant number of threshold thermal profiles in a given system, $T_p^{\mathrm{thr}}$ and $T_{\mathrm{crit}}^{\mathrm{thr}}$ values can be predicted for any given laser exposure duration and ambient temperature for that system. At least this seems true in our in vitro retinal model.

The closer a damage rate process model can predict $T_p^{\mathrm{thr}}$ and $T_{\mathrm{crit}}^{\mathrm{thr}}$ values relative to empirical values, the more optimal the model. All the $T_p^{\mathrm{thr}}$ and $T_{\mathrm{crit}}^{\mathrm{thr}}$ ($\tau =1$) reported in Table 2 fall within the standard deviations of the empirical threshold $T_p$ values, also shown in Table 2. This validates the predictive power of Eqs. (17) and (18) that uses $A_{\mathrm{non}}$ values. Most importantly, $A_{\mathrm{non}}$ not only accurately predicts the $T_p^{\mathrm{thr}}$ values, it maintains the required value of unity for ${\mathrm{\Omega }}_{\mathrm{non}}$ to confirm damage. This is the great advantage of taking the additional work needed to get the $A_{\mathrm{non}}$ value.

Average heating rate

For processes carried out at a constant heating rate ($\beta$), $=\frac{dT}{dt}$, Eq. (5) can be reorganized to Eq. (19):

Eq. (19)

$$\frac{d\alpha }{dT}=f(\alpha )\frac{A_{(t^{\prime })}}{\beta }{e}^{\frac{-E_a}{RT(t)}}.$$

Notice that assuming constant temperature (isothermal approach), $\beta =0$, generates an invalid assumption in the kinetic equation, eliminating its usage for analyzing nonisothermal laser-induced damage to tissue. The temperature history at the boundary location of cell death provides an estimate of the average heating rate required for threshold damage, (${\beta }_{\mathrm{ave}}^{\mathrm{thr}}(\tau )$), as a function of exposure duration. As previously indicated, threshold thermal damage is reached by means of a peak temperature and a heating rate together. In a thermal process with a heating rate below threshold, reaching the threshold peak temperature is not likely to result in the anticipated damage. Calculations of the threshold average heating rate for each exposure duration and ambient temperature used all 198 thermal profiles collected at the boundary of the cell death were tabulated (data not shown). Using the empirical data, we established a two-dimensional (2D) fit of the heating rate $\beta (\tau ,T_B)$ as a function of $T_{\mathrm{amb}}$ and $\tau$ at the boundary of the cell death was determined to be

Eq. (20)

$$\beta (\tau ,T_{\mathrm{amb}})=e^{-\left(\frac{\ln (\tau )+0.0284*T_{\mathrm{amb}}-3.8828}{0.93}\right)}.$$

Equation (20) is an empirical formula, but it was not a simple “best-fit” analysis. One can convey the analysis as that being a method by which the complex relationships between ${\beta }_{\mathrm{ave}}^{\mathrm{thr}}$ and $T_{\mathrm{amb}}$ and $\tau$, shown in Fig. 3, were iteratively related until the formula shown was found to consistently relate the two variables to the correct value for $\beta (\tau ,T_{\mathrm{amb}})$.

Figure 4 shows the plots of Eq. (20) for varying laser exposure durations at the same four ambient temperatures shown in Fig. 3(b). Notice the similarities between Figs. 3(b) and 4, which indicates a robustness represented using the empirical data to formulate Eq. (20). When the data points for the 16 exposure types ($T_{\mathrm{amb}}$ and exposure durations) were used to obtain power functions in Microsoft Excel [as was done in Fig. 3(b)], we obtained similar equations to those in Fig. 3(b). The equations provided in Fig. 4 indicate the same reciprocal relationship and very similar $T_{\mathrm{crit}}^{\mathrm{thr}}$ values [Eq. (18)]. This agreement between the empirical heating rates and the best-fit 2D function validates the predictability of Eq. (20) in our cell culture system using only the $A_{\mathrm{non}}$ and $E_a$ values.

Fig. 4

Plots of Eq. (20) at 20°C, 25°C, 30°C, and 40°C ambient temperature (solid lines). Calculated data points at 0.05, 0.25, 1.0, and 20 s are related to power functions (insets). The equations for the calculated values are essentially identical to those determined empirically [Fig. 3(b)]. Again, the power functions indicate reciprocity (${\tau }^{-1}$) and the $T_{\mathrm{amb}}$ + prefactor = $T_{\mathrm{crit}}^{\mathrm{thr}}$. All correlation coefficients ($R^2$) values were unity.