A new thermal dose model based on Vogel-Tammann-Fulcher behaviour in thermal damage processes

Abstract Thermal dose models are metrics that quantify the thermal effect on tissues based on the temperature and the time of exposure. These models are used to predict and control the outcome of hyperthermia (up to 45°C) treatments, and of thermal coagulation treatments at higher temperatures (>45°C). The validity and accuracy of the commonly used models (CEM43) are questionable when heating above the hyperthermia temperature range occurs, leading to an over-estimation of the accumulation of thermal damage. A new CEM43 dose model based on an Arrhenius-type, Vogel-Tammann-Fulcher, equation using published data, is introduced in this work. The new dose values for the same damage threshold that was produced at different in-vivo skin experiments were in the same order of magnitude, while the current dose values varied by two orders of magnitude. In addition, the dose values obtained using the new model for the same damage threshold in 6 lesions in ex-vivo liver experiments were more consistent than the current model dose values. The contribution of this work is to provide new modeling approaches to inform more robust thermal dosimetry for improved thermal therapy modeling, monitoring, and control.


Introduction
Quantifying the amount of thermal damage in biological systems based on thermal exposure (time and temperature of exposure) has been an active area of research for many decades [1][2][3][4]. Earlier, these studies focused mainly on skin burns [1,2,5]. The importance of heat-induced tissue damage studies increased with the use of heat in cancer treatments. The prediction of thermal damage incurred during or after treatment is essential for treatment planning, monitoring, and post-treatment assessment. The concept of a thermal dose is used for this purpose. The thermal dose is meant not merely to be a measure of the amount of energy delivered. Instead, it is meant to be a measure of the amount of tissue damage caused by heat. Heating has been investigated as an adjuvant to other forms of therapy, including radiation therapy and chemotherapy [6,7]. Controlled and mediated heating [8] can be used to activate drug release using thermosensitive liposomes [9]. Heating has also been used to speed up healing of injuries, and act as a surgical 'scalpel' to coagulate or ablate tissue [10]. Therefore an ideal thermal dose model would be applicable to various thermal processes and biological end-points.
Thermal damage, or heat-induced cell death, is believed to be caused by protein denaturation, for the range of temperatures used in thermal therapy [11][12][13]. The first thermal dose model known as the Arrhenius model was introduced by Henriques in 1947 [2] based on rate process analysis. In this model, thermal damage is quantified by a single parameter (denoted by X) which depends on the time-temperature profile and two tissue-dependent parameters that are experimentally determined for each type of tissue. Henriques found the parameters based on in vivo experiments for both porcine and human skin [1]. In those experiments, time-temperature combinations to achieve the same thermal damage end-point were obtained. The Arrhenius parameters reported for other tissue types from subsequent studies are significantly different [14], which is a challenge with using the Arrhenius dose model. Sapareto and Dewey [3] introduced another thermal dose model, which we refer to here as the CEM43 model, that quantified damage by the equivalent exposure time at some reference temperature. Although it was derived from empirical observations, they showed it to be a good approximation to the Arrhenius model in a limited temperature range. 43 C was chosen as the reference, where all thermal exposures are converted to an equivalent time, t 43, at this temperature. The formula that converts a time-temperature profile to a t 43 value was based on many experimental observations of different biological systems [1,15,16]. In most biological systems, a 1 model can be used for different tissue types without the need to obtain tissue-specific parameters; nevertheless, dosevalue thresholds to achieve a certain damage endpoint need to be known for different types of tissue. This simplicity, has led to its wide adoption for dose calculations in thermal therapy applications [17][18][19][20][21]. The current CEM43 model was originally developed for traditional hyperthermia applications (temperatures up to 45 C). The validity and accuracy of this model for high-temperature applications (above the hyperthermia temperature range) has remained uncertain. While the CEM43 model is simpler than the Arrhenius model because its unit of thermal dose is more clinically relevant and it does not require the determination of tissue-specific parameters to calculate the dose, its applicability to higher temperature exposures is uncertain. After introducing the CEM43 model, several publications presented alternative thermal dosimetry based on modified Arrhenius models or by linking Arrhenius and CEM43 models [22][23][24][25]. Nonetheless, thermal dosimetry is a challenging metric that requires further improvements. An examination of the accuracy of the CEM43 model at higher temperatures is presented here. Additionally, an alternative to the Arrhenius model, called the Arrhenius-VTF model is considered for thermal dose calculations. A mathematical relation between these Arrhenius-type models and the CEM43 model is presented, which leads to a new CEM43-type model, called the CEM43-VTF model. Finally, the performance of the new CEM43-VTF model is compared to the performance of the CEM43 and Arrhenius models for various biological applications of heat.

Theory
The CEM43 model and the Arrhenius model are described below followed by the Arrhenius-VTF model and a mathematically unification of the three models. This leads to a thermal dose equation for the new CEM43-VTF model. Finally a reformulation of the three models for the special case of constant temperature heating is provided, which allows them to be bench-marked against published heating results.

The CEM43 model
The most commonly used thermal dose model, CEM43, is based on the Sapareto-Dewey relationship [3] that quantifies thermal dose in terms of the equivalent amount of time needed at a constant 43 C to get the same thermal effect. This is an approximation to the Arrhenius model and is given by where T is the time-varying temperature (in C) during a treatment of duration s. C ¼ 4 when T<43 C and C ¼ 2 for T ! 43 C:

The Arrhenius model
A commonly used dose mode is based on Arrhenius analysis of reactions [26], is the Arrhenius model [1,2]. In this model thermal damage is modeled as a first-order rate process in which tissue constituents transform from native state to damaged state with a reaction velocity (or rate constant), k [2,11]. The dependence of k on temperature can be described by the Arrhenius equation [26]: Using this assumption the thermal dose, XðsÞ, is where NðsÞ is the concentration of the native tissue at a heating duration of s (making N(0) the concentration of native tissue prior to heating), A is a frequency factor, E a is the activation energy, and R is the universal gas constant.

The Vogel-Tammann-Fulcher model
Herein, we propose and investigate the use of a modified (and more general) Arrhenius-like model, known as the Vogel-Tammann-Fulcher [27][28][29] (Arrhenius-VTF) model for thermal dose, where where A 0 is a frequency factor, a is a parameter related to a temperature-dependent activation energy [30][31][32], and T 0 is the temperature where the damage process can no longer be thermally activated. This means that the rate constant k ¼ 0 for temperatures below T 0 ; therefore, there is no dose accumulation below this temperature.

A Unified thermal dose equation
The equations in all three dose models can be generalized to a form where the amount of thermal damage, D, depends exponentially and cumulatively on temperature, The pre-integral factor for the CEM43 model is a ¼ 1, while it is a ¼ A in the Arrhenius model, and a ¼ A 0 for the Arrhenius-VTF mode. The exponent for each model is, where bð43 CÞ is the value of b at T ¼ 43 C: This yields a generalized CEM43-type equation Substituting (6c) in (8), yields the CEM43-type version of the Arrhenius-VTF model: t 0 43 is used instead of t 43 to differentiate the dose values calculated using the two CEM43-type models. Equation (9) is hereafter referred to as the CEM43-VTF model.

Parameter fitting plots
The constant parameters in the dose models can be obtained using well controlled experiments that heat at various constant temperatures for various durations. This is done using special temperature-time plots, that are slightly different for CEM43 models compared to Arrhenius models.

CEM43 plots
For the special case of a treatment performed for a duration, s, at some constant temperature, T, and if the amount of damage at this temperature is the same as the damage incurred at a constant 43 C after a duration of s 43 , then, ae bðTÞ s ¼ ae bð43 CÞ s 43 : This results in the relationship, ln A plot of ln ðs 43 =sÞ versus T can be used to fit each model to experimental data using Equations (6a), (6b), or (6c).

Arrhenius plots
In order to use this model for a certain tissue type, one needs to experimentally find the coefficients A and E a for that tissue. This can be done using a series of heating experiments that produce the same amount of thermal damage (which is defined to correspond to X ¼ 1 [2,[33][34][35]). Since the heating is at a constant temperature, equation (2) can simplified to A plot of ln ð1=sÞ versus 1=ðT þ 273 CÞ can be used to obtain the activation energy, E a , from the slope of the plot, and A its y-intercept for the basic Arrhenius model. Similarly, for the Arrhenius-VTF one can reach an Arrhenius-VTF analogue: From a plot of ln ð1=sÞ versus 1=ðTÀT 0 Þ, one can find a and A 0 from the slope and the y-intercept, respectively. It should be noted that the process-based temperature, T 0 , that replaces the absolute zero in the Arrhenius model, has to be known in order to do this fitting.

Fitting the models to published results
Four published data sets from four different experiments were considered. Each of them provides the exposure time, s, at different constant temperatures, T, needed to achieve a certain isothermal effect. Three of the experiments were carried out in vitro by different groups [36][37][38] and reported by Henle and Dethlefsen [16]. The fourth was done in vivo [1]. The equations to be fit to the results all require the heating time at 43 C, s 43 : Since the fourth experiment starts at 44 C, s 43 was extrapolated from the other data points. The other three experiments had measurements at 43 C heating.

Epidermal injury
In their work titled 'Studies of thermal injury' in 1947, Moritz and Henriques tried to find information about the rate at which skin burning occurred at any given surface temperature. They gained time-temperature thresholds resulting in different levels of skin injuries. They did their experiments in vivo for both porcine and human skin. Experiments were done at constant temperatures. The skin surface was immediately brought to and maintained at a predetermined constant temperature using an applicator by which a running stream of temperature-controlled water was brought in direct contact with the skin [1]. In introducing the Arrhenius based dose calculations, Henriques used two of the epidermal injury thresholds: 'threshold A' and 'threshold B' [2]. Threshold A was used to produce one of the data sets used in Figure 1. This threshold represents the minimum exposure   time resulting in complete transepidermal necrosis. Threshold B represents the maximum exposure time resulting in reversible epidermal injury. Threshold B was used to test the CEM43 model as shown in Figure 2 in the results.

Optical properties and dose
Increase in tissue optical scattering due to damage accumulation can be significant [39]. These changes are clear in many examples; like the change of egg white from transparent (low scattering coefficient) to opaque (high scattering coefficient) after heat exposure. Thermal damage, with its related protein denaturation, manifests as an increase in optical scattering [40]. Since the scattering coefficient is correlated to thermal damage, the optical scattering is expected to be an injective function of thermal dose. We used previously published experimental results on the changes in the optical scattering due to heating of tissue at different temperatures [39] to study the ability of the models to predict changes in optical scattering. The normalized reduced scattering coefficient, l 0 s ðtÞ=l 0 s ð0Þ, is used and compared to the thermal does values obtained using the CEM43 and CEM43-VTF models. Here, l 0 s ðtÞ is the reduced scattering coefficient at a certain time during the heating and l 0 s ð0Þ is the the native reduced scattering coefficient. Henceforth, we will use 'normalized l 0 s ' as a short form of the normalized reduced scattering coefficient.

Laser-induced heating
Thermal dose models are needed for clinical applications like cancer treatment where temperatures vary spatially and temporally during treatment. Ex-vivo liver heating experiments were performed to test the performance of the thermal dose models under such realistic conditions. In this experiment, six thermal lesions were created on the surface of an adult bovine liver sample using laser energy. The sample was originally frozen, and was brought quickly to room temperature before the start of the experiment. The tissue was allowed to cool to room temperature after the creation of each lesion. A multimode optical fiber with a 1 mm core diameter and 0.37 numerical aperture (ThorLabs, Newton, NJ, USA), coupled to a Diomed 60, 810 nm diode laser (Diomed, Cambridge, UK) was placed 13 mm away from the tissue surface. A constant laser power of 1.3 W À 1.6 W for 1 min. to 6 min. exposures was used to produce lesions of varying size and temperature histories. The surface temperature map for each lesion was collected during heating using a thermal camera (FLIR ThermaCAM SC2000, FLIR Systems, Burlington, On, Canada) with a Close-up lens (34 mm horizontal field-of-view and 80 mm focal length) at a rate of one frame per second. The thermal camera continued to collect data two minutes after the laser was turned off since damage keeps accumulating after heating has stopped. Temperatures at the end of this period were too low to accumulate more damage significantly.
Since thermal damage in liver tissue is correlated with change in tissues optical properties, observed as tissue whitening [39], an arbitrary brightness threshold was chosen in an optical colored image acquired after the experiment. This was done by converting the colored image to a grayscale image and then choosing a certain image intensity (which corresponds to choosing a certain damage threshold) to create a binary image of six white areas, wherein the image intensity is higher than the selected threshold, in the black background where the image intensity in less than that threshold. Then, the white areas inside these threshold boundaries were calculated for all six lesions (A1 to A6).
Thermal dose maps were calculated using the thermal camera temperature data using both the CEM43 and CEM43-VTF models. The two-minute cool-down period was included in the calculations. The calculated dose values were used for each lesion to find the dose threshold that encloses the same area obtained earlier from the binary image for that lesion to determine the thermal dose threshold's consistency between lesions. Figure 1 shows the fit of the three existing models to four published sets of experimental results. The results support the universality of equivalent-minute calculations since data from different experiments in different biological systems can be fit using a single function. It also shows that the relationship used in the current CEM43 dose model (the black dashed line) fits the data reasonably well for hyperthermia temperatures range, but it diverges increasingly for higher temperatures. The results demonstrate that the CEM43 model over-estimates dose accumulation at higher temperatures. The results also show that the Arrhenius-VTF based relationship (the red solid line) fits the data much better than the basic Arrhenius based relationship (the blue dotted line). The fit to Equation (11)

Epidermal injury
An ideal thermal dose model would predict the same dose value for each combination of the time-temperature combinations for the threshold B isothermal effect in the epidermal injury studies [1,2]. Figure 2 shows that the CEM43-VTF dose values (t 0 43 ) are of the same order of magnitude and are more consistent than the CEM43 dose values (t 43) , which have 2 orders of magnitude difference at higher temperatures. More specifically, the coefficient of variation of the t 0 43 values is 0.312, while for the t 43 dose values it is, 2.72. This shows that the new CEM43-VTF model is more consistent in predicting thermal damage than the current CEM43 model.

Optical properties and dose
Using previously published results of the changes of optical scattering due to heating [39] the normalized l 0 s plotted against dose is shown in Figure 3. There is a more consistent relationship between the normalized l 0 s and t 0 43 than there is between l 0 s and t 43 as temperature is varied, which shows that the CEM43-VTF model is capable of quantifying biological end-points other than cell death. The use of the CEM43-VTF model in simulations of laser thermal therapy to include the dynamic changes of the optical scattering coefficient of tissue due to heating will result in more accurate simulation results.
As mentioned in the methods (subsection 3.3), the optical scattering is expected to be an injective function of the thermal dose; i.e., the normalized scattering must be, ideally, a one-to-one function of the thermal dose value. The closer the relation to this ideal case, the more truthful the dose model, is in quantifying thermal-induced changes in tissues. While the results in Figure 3(b) are not perfect, they are significantly closer to a one-to-one relation than the results in (a).

Laser-induced heating
The optical image that was acquired after the experiment is shown in Figure 4, along with the corresponding binary image that was produced as descried is subsection 3.4. Figure 5 shows the CEM43 and CEM43-VTF threshold dose values such that the liver surface areas that received thermal dose higher than or equal to these thresholds match the thermal damage areas (A1-A6) for the six thermal lesions  shown in Figure 4(b). Ideally the dose threshold for all the lesions would be the same. The dose threshold values for the six lesions using the CEM43 model is shown in Figure  5(a). They have an average of 3:4 Â 10 6 eq. min. with a standard deviation of 3:3 Â 10 6 eq. min. The coefficient of variation was 0.99. The CEM43-VTF thresholds ( Figure 5(b)) have an average, of 3:0 Â 10 4 eq. min. with a standard deviation of 1:1 Â 10 4 eq. min. and a coefficient of variation of 0.37. Similar results were obtained when the same procedure was repeated for a second liver sample. For the second sample, the coefficient of variation for the CEM43 values (0.86) was again higher than the one for the CEM43-VTF values (0.40). These results demonstrate the consistency of the new CEM43-VTF dose model in predicting thermal damage in comparison with the current CEM43 dose model.

Discussion
Thermal dose models are empirically derived formulae to be used in clinical applications. The models quantify thermal damage by finding relationships between exposure time, temperature, and a biological end-point. These relationships have been studied for many biological systems both in vivo and in vitro. The studies show that thermal damage rises exponentially with temperature for a given exposure time and rises linearly with exposure time for a given temperature. The current CEM43 thermal dose model assumes that the slope of the curve in a CEM43 plot (section 2.5.1), like the one in Figure 1, has a break at 43 C in all biological systems. The values of the slope is ln ðCÞ, which is obtained by substituting Equation (6a) into Equation (11). While the value C should vary continuously with temperature [3], using only two values is a good approximation for traditional hyperthermia applications (up to 45 C). On the other hand, this approximation is not suitable for high-temperature applications, above the hyperthermia range as evident in Figures 1  and 2. The CEM43-VTF model has a continuously varying slope for all temperatures (Figure 1), making it a more flexible dose model.
The Arrhenius-VTF model is a generalized form of the basic Arrhenius equation. In Arrhenius-VTF there is a temperature, T 0 , below which the reaction is not thermally    activated [30,31,41], whereas the reaction is always thermally activated above absolute zero (-273 C) in the basic Arrhenius case. The benefit of using the Arrhenius-VTF formulation for modeling thermal damage can be seen in the Arrhenius plots (section 2.5.2) shown in Figure 6 for various thermal damage experiments. The Arrhenius-VTF plot was generated using the activation temperature T 0 ¼ 23:5 C: The activation energy in the Arrhenius model, E a and the activation energyrelated parameter, a, in the Arrhenius-VTF model, are obtained from the slope in a linear regression fit to the data in these plots. The Arrhenius-VTF plot ( Figure 6(b)) shows that the data from various experiments [2,36,38,42] are more linear than the points in the basic Arrhenius plot ( Figure  6(a)) for the same measurements. This implies that E a in the basic Arrhenius model is temperature-dependent as has been observed in some studies [4,13,14,43]. Using the Arrhenius-VTF model essentially eliminates this dependence. Other publications [24,25] have used a modified-Arrhenius model that has a similar mathematical form as the Arrhenius-VTF equation. However, we used the established Arrhenius-VTF model based on a physical interpretation of the T 0 parameter being an activation temperature (below which the process is not thermally activated). While the T 0 ¼47.5 C value in these works was arbitrarily chosen, the T 0 ¼23.5 C value was obtained by fitting to experimental results. In addition, we derived the VTF-CEM43 equation from the Arrhenius-VTF equation. A recent publication [22] has linked and combined both the Arrhenius and the CEM43 models, but they used the original Arrhenius model, not the improved VTF-Arrhenius model.
The CEM43-VTF model assumes that all tissue types have the same slope in a Arrhenius-VTF plot, but this is not the case. Equation (14) was derived using a ¼ 382 C: This is an acceptable value for various tissue types since the tissue-specific parameter ranged between a ¼ 371 C and a ¼ 415 C in Figure 6(b). The goodness of fit to these tissues in Figure 1 was R 2 ¼ 0:99 for a ¼ 382 C: Specific tissue can be used to fit a and A 0 from the Arrhenius-VTF plots (assuming T 0 ¼ (4)) if a more precise, tissue-specific, is desired. Thermal dose models use relative metrics that are positive-valued increasing functions of the thermal damage. One is the dimensionless X in the Arrhenius models, and others are the t 43 or t 0 43 in the CEM43 models. These dose metrics are different in the way they quantify thermal damage. The metrics t 43 or t 0 43 are independent of the tissue type, but the threshold for damage is tissue specific. The Arrhenius metric, X, is tissue-specific, and the threshold X ¼ 1 quantifies the same damage end-point, independent of tissue type. Although this study introduces two new models for thermal dose, the Arrhenius-VTF and CEM43-VTF models, we propose the adoption of the CEM43-VTF model since it is simpler and more universal. If an Arrhenius-type dose model is preferred, then the Arrhenius-VTF model should be used instead of the basic Arrhenius model. Since Arrhenius-VTF is new to thermal dosimetry studies, new parameters a and A are needed to be derived for different tissue types using the T 0 value that was established in this work. Furthermore, since the relation between the CEM43 models and the Arrhenius-type models is established here, if the parameter a for a certain type of tissue is determined, it can be used for CEM43-VTF calculations as well.

23:5 C in equation
One of the motivations for this work was the development of a thermal dose model that could be used for a broad range of thermal applications, including those that target cell death as the end-point, changes in the optical properties of tissues, or the amount of drug released from thermosensitive liposomes. The results from the two laser heating experiments (sections 4.3 and 4.4) show that the CE43-VTF model was able to more accurately predict the amount thermal damage, more specifically, the reduced optical scattering coefficient (directly in section 4.3 and indirectly through tissue whitening in 4.4) in compared to the CEM43 model. Another motivation was to develop a thermal dose model that was able to predict cell death when heat was used as an adjuvant to other therapies (such as chemo-  and radio-sensitization). This was not investigated in this work and this goal may never be realized for the case of radio-sensitization with heat. That is because one of the mechanisms of radio-sensitization with heat is the reduction of hypoxia due to increased blood perfusion caused by heating. Unfortunately there is no functional relationship between the response of blood perfusion and temperature. For example an elevated temperature can cause an initial increase in perfusion but later produce a reduction in perfusion due to vascular shut-down [44].

Summary and conclusion
Accurate quantification of thermal damage is critical for clinical applications of thermal therapy but the damage process in biological systems is complicated, and understanding it is not trivial. To improve thermal dosimetry, experimental results for a wide range of temperatures in different types of tissue along with thermal damage models that use these experimental results more effectively are needed. A mathematical formulation for dosimetry is provided that makes it easier to develop new models using a variety of experimental results. Using this formulation, a new model for CEM43type dose calculations is introduced that improves on the existing CEM43 model in terms of accuracy and consistency over a wider temperature range. This consistency is evident in the epidermal injury in-vivo results (Figure 2), where the coefficient of variation for the dose values calculated using the old model was 2.72, is almost an order of magnitude higher than that calculated using the new dose model (0.312). Even for the range of temperature below 50 C, where the old model performs well, the coefficient of variation of 0.316 is higher than the coefficient of variation of 0.255 from the new model. While a new dose model has been introduced and fit to various experiments, this work introduces the foundation for a new dose model. While the use of this new dose model will lead to more robust thermal dosimetry and improvements in thermal therapy modeling, monitoring, and control, further experiments in more clinically relevant scenarios are needed to better fit and tailor the model parameters. In addition, new studies are needed to find thresholds that represent the minimum dose (as equivalent time) needed to cause irreversible damage to tumor tissues that are commonly targeted by heat treatments, such as prostate, liver, and brain.The CEM43-VTF dose thresholds that correspond to the maximum dose that tissues of some vital organs can receive without getting irreversibly damaged need to be determined. The suitability of the CEM43-VTF model for predicting other end-points, such as drug release from thermosensitive liposomes, is also needed.Such work will be very important for the clinical use of thermal therapy.