Mathematical model for the thermal enhancement of radiation response: thermodynamic approach

Abstract

Radiotherapy can effectively kill malignant cells, but the doses required to cure cancer patients may inflict severe collateral damage to adjacent healthy tissues. Recent technological advances in the clinical application has revitalized hyperthermia treatment (HT) as an option to improve radiotherapy (RT) outcomes. Understanding the synergistic effect of simultaneous thermoradiotherapy via mathematical modelling is essential for treatment planning. We here propose a theoretical model in which the thermal enhancement ratio (TER) relates to the cell fraction being radiosensitised by the infliction of sublethal damage through HT. Further damage finally kills the cell or abrogates its proliferative capacity in a non-reversible process. We suggest the TER to be proportional to the energy invested in the sensitisation, which is modelled as a simple rate process. Assuming protein denaturation as the main driver of HT-induced sublethal damage and considering the temperature dependence of the heat capacity of cellular proteins, the sensitisation rates were found to depend exponentially on temperature; in agreement with previous empirical observations. Our findings point towards an improved definition of thermal dose in concordance with the thermodynamics of protein denaturation. Our predictions well reproduce experimental in vitro and in vivo data, explaining the thermal modulation of cellular radioresponse for simultaneous thermoradiotherapy.

Introduction

Despite considerable efforts for decades towards the improvement of early diagnosis and therapy, cancer has remained a serious global health problem, with 18.1 million new cases and 9.6 million cancer deaths reported worldwide, just in 2018¹. Since the 1980s, mild hyperthermia (heating tumour tissue to 40.0–42.5 °C for ~ 1 h) is known to enhance the therapeutic outcomes in cancer patients, when combined with radio-, chemo- and/or immunotherapy^2,3. Technological improvements in precise medical heating, imaging and non-invasive thermometry over the past decade have revived hyperthermia treatment (HT) as a precision cancer therapy^3,4,5,6, particularly when used in simultaneous combination with ionizing radiation^7,8,9. The number of ongoing HT clinical trials, either alone or in combination with different treatment modalities, evidences the increasing use of therapeutic HT (467 still ongoing clinical trials out of 1198 since 2000)¹⁰. Radiotherapy (RT) is supposedly a curative treatment modality, but the radiation dose required to eradicate all cancer cell subpopulations in a tumour can often not be applied due to severe acute or long-term side effects, which include radiation-induced tissue fibrosis and second malignancies¹¹. Hyperthermia is known to be one of the most potent radiosensitisers^{12,13,14,15,16}, meaning that less radiation is required to achieve the same local tumour cell kill, thereby reducing the adverse effects of radiation in the adjacent normal tissues, e.g.^{17,18,19,20,21}.

The efficiency of combined HT+RT treatment clearly depends on the scheduled sequence of the two types of treatment and the time interval between them; best outcome for the patient is expected from a simultaneous application^13,22,23,24. However, several theoretical and practical problems still need to be overcome to implement simultaneous HT+RT approaches in routine clinical practice worldwide²⁵. Indeed, in practical terms simultaneous treatment has remained challenging because spatially precise hyperthermia delivery is required to avoid unspecific synergistic, cytotoxic effects of HT+RT on the surrounding normal tissue, which would critically limit the therapeutic benefit. Current standard clinical equipment is still not well-suited for such simultaneous and precise thermoradiotherapy. Accordingly there are only a few reports/completed trials, in which both forms of radiation were concomitantly applied in patients^{7,26,27,28,29}. From the theoretical perspective, mathematical models to predict the therapeutic outcome of various combinatorial treatment schemes are essential for a better understanding of the synergistic potential and therapeutic window of the two sources of energy (HT and RT), and are highly relevant to the design of adequate and individualized treatment planning in the clinical setting^30,31.

Nomenclature

The key biological terms used in this work have been specified as follows:

• Cell kill (“dead state”): From the radiobiological perspective, a cell is considered to be dead (killed) when it loses its proliferative capacity, i.e. is no longer able to divide (becomes replication-incompetent). This encompasses not only cells losing their membrane integrity and truly dying (by apoptosis, necrosis, or other), but also living cells undergoing terminal differentiation, permanent cell cycle arrest or senescence. This type of cell kill leads to control of the malignant disease, independent of the underlying process.

• Cell survival (“alive state”): A cell is considered to survive if it remains replication-competent, i.e. retains its proliferative capacity after the treatment.

• Cell damage: Any type of deterioration of the cellular processes, regardless of origin, that advances the cell towards the dead state.

• Radiological parameters $\alpha$ and $\beta$: They characterise the radiosensitivity of cells or tumours.

  • $\alpha$ - Initial slope of logarithmic survival curves. It is associated to the mean number of DNA double strand breaks produced with a single radiation event³².

  • $\beta$ - Shoulder of logarithmic survival curves. It is associated to the mean number of DNA double strand breaks produced with two radiation events, i.e. two independent single strand breaks in close proximity that lead to formation of a double strand break³².

  • $\alpha /\beta$ ratio - Quantifies radiation sensitivity of tissue. The higher the ratio, the lower the sensitivity.

• Thermal enhancement ratio (TER): Ratio between the radiation dose required to achieve a specific endpoint with ionizing radiation alone, and the radiation dose required to achieve the same endpoint in combination with hyperthermia.

Several mathematical models for individual RT and HT have been proposed, but there is poor consensus when it comes to the efficacy of combined treatment regimes. For RT, the LQ-model is the most extensively used approach to predict the effect of irradiation on cell populations^33,34. This model describes the surviving fraction of cells as a function of the applied radiation dose $D_R$ by means of two main variables, called “radiological parameters” $\alpha$ and $\beta$³⁴. In the context of radiobiology, “survival” means the conservation of the cell’s proliferative capacity³⁵ (see definitions box). Regarding HT, there is considerable literature describing the impact of heat on different cellular components^36,37,38,39, and several models are aimed to predict the survival of cells under HT treatments, e.g. the Jung’s model⁴⁰, multiple-states models, Arrhenius models, biochemical models, stochastic models (reviewed in^41,42), among many others which include derivations of the LQ-model for RT^42,43 . For thermal-radiosensitisation using temperatures of 40–46 °C, there is a general agreement on a relevant role of DNA repair impairment by heat-induced protein denaturation in the processes of radiosensitisation^{12,33,36,37,39,44}. The majority of previous approaches to model the combined efficacy of hyperthermia and radiation on mammalian cells have implemented the thermal effects in the LQ-model by proposing empirical temperature dependencies for the radiological parameters^43,45,46, but the physical principles and the detailed mechanisms underlying this empirical dose-lowering concept are still elusive⁴⁴. The link between modelling concepts and plausible mechanistic explanations still needs to be established to serve as a more reliable framework for predictions.

Here, we describe a survival model for the simultaneous application of HT and RT that provides insights from a thermodynamic perspective. The intention of our work is not to predict the outcome of the HT or RT alone, but to propose a new model for the synergistic effect of the combined modality. This is based on the modulation of the radiological parameters in the LQ-model as a function of the HT temperature and treatment time. In our framework, this modulation arises directly from the definition of the thermal enhancement ratio (TER). It compares the radiation dose required to achieve a specific endpoint with ionizing radiation alone ($D_R$), e.g. surviving fraction of cells or tumour control probability, and the radiation dose required to achieve the same endpoint in combination with hyperthermia ($D_{R+H}$) $\text {TER}=\frac{D_R}{D_{R+H}}$⁴⁷. We propose the enhancement to be a rate limiting process, proportional to the energy invested in sensitising a cell to die. Our approach presents a theoretical basis to understand how hyperthermia results in radiosensitisation, a process that depends on treatment time and temperature. We show that our findings are consistent with previous experimental studies in the range of RT combined with hyperthermia between 40 and 46 °C, where irreversible protein coagulation and thermal ablation effects are not relevant⁴².

Methodology

Model considerations

In the context of our study, the cell killing process in thermoradiotherapy is composed of two main stages: sensitisation and damage fixation, as highlighted in Fig. 1. We, however, propose that in combination with RT, HT mainly affects the sensitisation stage, making the cell more vulnerable to die. According to our hypothesis, only a minor fraction of the thermal energy relates to damage fixation. The thermal enhancement of radiotherapy is usually explained by the reoxygenation of the perfused tissue and by protein denaturation (and coaggregation) occurring at different cellular compartments^39,42,48,49. The latter produces diverse cytotoxic effects; these are believed to mainly relate to the sensitisation in the HT regime. On the other hand, molecular oxygen favours the fixation of radiation-induced damage, and therefore, reoxygenation enhances the RT outcome via DNA damage fixation. Notably, reoxygenation only takes place in vivo where heat increases tumour perfusion, but is not reflected in current in vitro models. Furthermore, in culture medium molecular oxygen is less dissolved at higher temperatures⁵⁰. In order to develop a model that is valid for both, in vitro and in vivo conditions, our framework first focuses on HT-induced sensitisation by protein denaturation while the oxygen effect will be implemented at a later stage.

Figure 1

Left: Schematic survival probabilities for the three cases depicted on the right. (a) Cell killing as a single rate process with transition rate from alive (A) to dead (D) $\alpha$. (b) Two-step cell killing process in the LQ-model for radiation. A cell transits from the alive state (A) to the dead state (D) through two possible paths: $\alpha$ for direct killing (a single hit suffices to kill), and $\beta$ for indirect killing (when two hits are required to kill). (c) Combined HT+RT: HT-induced damage elevates cells from state (A) to an activated state ($A^\prime$), effectively reducing the $\alpha / \beta$ ratio. Since $\beta$ is more efficiently reduced, the direct path $\alpha$ dominates the killing process and consequently reduces the survival probability.

Thermal enhancement of radiotherapy

When a certain radiation dose $D_R$ is applied to a set of living cells, the reduction rate is proportional to the number of cells at the time of the treatment

$$\begin{aligned} \frac{dN}{dD_R}=-\alpha N. \end{aligned}$$

(1)

Therefore, the direct transition from the alive state of the cell to the dead state obeys an exponential behaviour $S=e^{-\alpha D_R}$, where $S=N/N_0$ is the survival fraction, and $\alpha$ defines the transition rate per dose, as depicted in Fig. 1a³⁴. If the killing effect is composed of a direct killing path $\alpha$, and a secondary path composed of two or more stages of sublethal damage (SLD) accumulation, the logarithmic survival curve acquires a shoulder, as depicted in Fig. 1b. In the particular case of the LQ-model, the exponent has a linear and a quadratic contribution, corresponding to direct killing and SLD accumulation respectively

$$\begin{aligned} -\ln (S)=\alpha D_R +\beta D_R^2 \text { .} \end{aligned}$$

(2)

The LQ-model was originally employed as an empirical approach⁵¹; later Chadwick and Leenhouts³² proposed a molecular interpretation based on a statistical approach. In their interpretation cell death occurs due to double-strand breaks (DSB) of DNA, such that $\alpha$ and $\beta$ account for the probability of producing irreparable DSB as a consequence of one or two photon/particle hits, respectively. As a consequence of the sensitisation effect of HT, the radiation dose $D_{R+H}$ required to produce the same surviving fraction is reduced. This reduction implies in Eq. (2), that $\alpha$ and $\beta$ are increased to $\alpha ^*$ and $\beta ^*$ (in order to obtain the same therapeutic outcome), assessing the increased sensitivity of the cells as a consequence of heat

$$\begin{aligned} -\ln (S)=\alpha ^* D_{R+H} +\beta ^* D_{R+H}^2 \text { .} \end{aligned}$$

(3)

This radiosensitising effect of hyperthermia is quantified by the thermal enhancement ratio TER. It is defined as the ratio between the radiation dose required to achieve a specific endpoint with ionizing radiation alone ($D_R$), and the radiation dose resulting in the same endpoint value when combined with hyperthermia ($D_{R+H}$):

$$\begin{aligned} \text {TER}=\frac{D_R}{D_{R+H}}, \end{aligned}$$

(4)

with $D_R>0$ and $D_{R+H}>0$. The new linear and quadratic coefficients of the LQ-model are obtained by replacing $D_R$ with $D_{R+H}\text {TER}$ in Eq. (2):

$$\begin{aligned} -\ln (S)=\alpha \text {TER} D_{R+H} +\beta (\text {TER})^2 D_{R+H}^2 \text { .} \end{aligned}$$

(5)

Comparing Eqs. (3) and (5) shows how the radiobiological parameters are effectively rescaled by hyperthermia to $\alpha ^*= \alpha \text {TER}$ and $\beta ^*=\beta \text {TER}^2$. Notably TER has a stronger effect on $\beta ^*$, bending the survival curves to lower survival values, in accordance with previous empirical data from experimental and clinical values studies^30,52, bending the survival curves to lower survival values. We propose a model for TER as a function of HT parameters, namely temperature and time, which is incorporated into the LQ-model to predict the survival probability of RT combined with HT. As detailed in the results section, TER is assumed to be proportional to the energy absorbed in the transition from the live state (A) to the more vulnerable state ($A^\prime$) $E_{A \rightarrow A^\prime }$, which in turn is defined as a rate-limited process

$$\begin{aligned} TER \propto E_{A \rightarrow A^\prime } = c_1 + c_2 k(T)t \text { ,} \end{aligned}$$

(6)

where $c_1$ is the baseline of TER, and $c_2$ accounts for the cell-line specific radio- and thermal sensitivity. In the absence of hyperthermia $TER=1$, resulting in $c_1=1$. The transition rate from (A) to ($A^\prime$) k(T) is modelled assuming protein denaturation as the mechanism responsible for heat-induced cell damage, as described in the next section.

Temperature dependence of the transition rate

The temperature dependency of the transition rate k(T) is modelled by means of the Eyring’s transition state theory⁵³:

$$\begin{aligned} k(T)=\left( \frac{K_B}{h_p}\right) T e^{-\frac{\Delta G(T)}{k_B T}}, \end{aligned}$$

(7)

where $k_B$ and $h_p$ are the Boltzmann and Planck constants, respectively, and T is the temperature in Kelvin. We next introduce a suitable model for the change in Gibbs energy $\Delta G(T)$ consistent with protein denaturation.

All conformation changes during protein denaturation arise from the competition between formation and breakage of chemical bonds. Protein denaturation becomes thermodynamically more favourable with increasing temperature. The dynamics of protein bonds is quantified by the standard heat of reaction $\Delta H_0$ and the thermal work function $\Delta W(T)$ respectively. We model the mixture of proteins sensitive to hyperthermia, as an average equivalent protein⁵⁴. Its overall heat capacity changes as a result of the state changes of individual proteins within the mixture. All the $\Delta$ symbols refer to changes in the thermodynamic properties of this “equivalent” protein, before and after the transformation.

The energy source for bonds to break is the thermal content $\int _{T_0}^T C_p(T')dT'$, which refers to the heat absorbed during the process of protein unfolding while temperature increases. Only a part of the absorbed heat can be converted into bond-breaking work, as restricted by the second law of thermodynamics. The unused proportion of thermal content goes into entropy –the thermal work– and is proportional to the absorbed heat and the relative temperature increment. The expressions for enthalpy and work content read as⁵⁵:

$$\begin{aligned} \Delta H(T)&= \Delta H_0+\int _{T_0}^T C_p(T')dT' \text {,}\nonumber \\ \Delta W(T)&= \int _{T_0}^T C_p\frac{(T-T')}{T'} dT', \end{aligned}$$

(8)

where $\Delta H$ is the enthalpy of the reaction, containing the bond forming energy $\Delta H_0$ and the sum of isothermal transfers of heat $\int _{T_0}^T C_p(T')dT'$. Here $C_p(T)$ is the heat capacity, which might vary with temperature, according to the third law of thermodynamics. The net driving energy is then given by the Gibbs free energy

$$\begin{aligned} \Delta G(T)= \Delta H_0-\Delta W(T)=\Delta H-T \Delta S(T) \text {,} \end{aligned}$$

(9)

where $\Delta S(T)=\int _0^T \frac{C_p(T')}{T'}dT'=\Delta S_0+\int _{T_0}^T \frac{C_p(T')}{T'}dT'$ is the entropy change, with $\Delta S_0$ as reference value. Accordingly, the Gibbs energy is expressed as

$$\begin{aligned} \Delta G(T)= \Delta G_0 + \int ^T_{T_0} dT' C_p(T')\left[ 1-\frac{T}{T'}\right] \text { ,} \end{aligned}$$

(10)

where $\Delta G_0= \Delta H_0-T \Delta S_0$. The reference temperature can be chosen so that $\Delta H(T_0=T_h)=0$, $\Delta S(T_0=T_S)=0$, or $\Delta G(T_0=T_g)=0$. $T_g$ and $T_s$ are of particular interest since they define the melting and maximal stability temperatures of the protein, respectively. When bond formation and breakage reach a balanced state ($\Delta G(T_g)=0$) the reaction does not progress anymore. The melting temperature is defined as the temperature at which the half of the proteins are denatured⁵⁶. Due to the importance of protein denaturation, the melting point is used as the reference temperature from now on.

The next challenge is to model the heat capacity in aqueous solutions above physiological temperatures. The heat capacity is expected to increase with temperature before approaching the vicinity of the melting point, as a result of ongoing protein reconfigurations. Beyond the transition, exothermic co-aggregations of proteins occur and $C_p$ is expected to decrease due to the reduced degrees of freedom of more rigid proteins. With these arguments, we propose to consider the next order by introducing the heat capacity change as a linear function of $(T-T_g)$⁵⁷, $C_p(T)=A -B\vert T-T_g \vert$, which is the same as $C_p(T)=A +B(T-T_g)$ for $T \le T_g$, leading to

$$\begin{aligned} \Delta G(T)= \Delta G_c -\frac{B}{2}(T^2-T_g^2)+BT T_g\ln \left( \frac{T}{Tg}\right) . \end{aligned}$$

(11)

Here $\Delta G_c$, is the usual Gibbs energy resulting from the assumption of constant heat capacity change. By introducing Eq. (11) in Eq. (7), the transition rate for denaturation becomes

$$\begin{aligned} k(T)=\left( \frac{K_BT}{h_p}\right) e^{-\frac{\Delta G_c}{K_B T}} e^{\frac{B}{2K_B}(T-T_g) \left( 1+\frac{Tg}{T} \right) }. \end{aligned}$$

(12)

where the last term in Eq. (11) should vanish, because $T_g/T$ is about one in the Kelvin scale for the hyperthermia temperature range (40–50 °C). The first two factors of Eq. (12) slightly change ($\sim \pm 2.5\%$) in these regimes, and then the transition rate is dominated by the exponential behaviour. Based on these considerations, k(T) can be described as

$$\begin{aligned} k(T) \approx c\, e^{b(T-T_g)}, \end{aligned}$$

(13)

with $c=\left( \frac{K_BT}{h_p}\right) e^{-\frac{\Delta G_c}{K_B T}}$ and $b=\frac{B}{K_B}$ as –cell dependent– adjustable parameters of the model.

Results and discussion: Mathematical model for the outcome of simultaneous HT+RT

In the following we describe our theoretical model and its correspondence with different types of experimental data in the range of HT between 40 and 46 °C, derived from mammalian cell models. Numerous in vitro and in vivo studies reveal the successful and promising combinatorial application of HT and RT for anticancer treatment (see e.g. Refs^26,27,52). However, the majority of the documented data is quite limited or incomplete and thus insufficient to test our model. To this end, we chose three rather dated seminal studies because, to our knowledge, they are the only ones which compile complete sets of thermal enhancement ratios, systematically obtained for several temperatures or treatment times in the HT regime. Two of these data sets were collected in course of in vitro 2D culture experiments using Chinese hamster ovary (CHO)⁵⁸ and murine mammary carcinoma (M8013) cells⁵⁹, respectively. Another set of data comes from an animal study with C3H murine mammary carcinoma experimental tumours⁶⁰.

Hyperthermia affects the radiation dose-response curve

The LQ-model for radiotherapy predicts the surviving fraction of cells as an exponential function of the radiation dose, $S(D_R)=\exp \left\{ -(\alpha D_R + \beta D_R^2)\right\}$³⁴. When HT is applied in combination with RT the parameters $\alpha$ and $\beta$ are modulated by both the temperature T, and the application time t of heat^52,59,61,62. As a result, the sensitivity of cells to RT is increased and the radiation dose $D_{R+H}$ required to produce the same surviving fraction is lower. HT affects the survival probability curves in three ways: 1. the curves are shifted down as a consequence of cell killing from HT itself (offset at $D_R=0$), 2. there is a steeper initial slope ($\alpha$), and 3. the shoulder of the curve ($\beta$) is changed as illustrated in Fig. 1. In this work, the term “cell kill” is defined as the complete loss of proliferative capacity of a cell, regardless its membrane integrity.

The most accepted hypothesis for the radiosensitising effect of HT assumes the heat-induced denaturation of repair proteins impairs the DNA repair process upon irradiation^12,33,44. In the LQ-model hyperthermia mainly affects $\beta$, which is supposedly related to repairable DNA single-strand breaks (SSB), and the HT-induced sensitisation is generally associated with inhibition of DNA repair³⁰. Nevertheless, this description is incomplete because the change in $\alpha$ is not negligible. Given that $\beta$ is not exclusively related to pairs of SSB but also to clusters of DNA lesions³⁴, we propose to differentiate between repairable and sublethal DNA damage, which are not necessarily the same. We suggest to extend the hypothesis of repair inhibition to a more general explanation based on sublethal damage accumulation (whether reversible or not), to better understand the synergy between radiation and thermal energy when applied to biological tissue.

Modulation of $\alpha$ and $\beta$ by HT as a function of TER

We propose that the radiosensitising portion of the energy is invested in the accumulation of sublethal damage, facilitating radiation-induced cell death. In our model hyperthermia causes the cells to advance from an original undamaged state (A) to a more damaged state ($A^\prime$) in the sequence of sublethal damage (SLD) accumulation, as is illustrated in Fig. 1c. Starting from ($A^\prime$) instead of (A), the radiation energy required to produce lethal and sublethal transitions is reduced, and hence, $\alpha$ and $\beta$ are effectively rescaled to $\alpha ^*$ and $\beta ^*$. Further, we assume that this modulation comes directly from the definition of the TER, in such a way that the new parameter ($\alpha ^*$ and $\beta ^*$) become treatment-time and temperature dependent (see Methods section for details)

$$\begin{aligned} \alpha ^*(T,t)&=\alpha * \text {TER} \end{aligned}$$

(14)

$$\begin{aligned} \beta ^*(T,t)&=\beta * \text {TER}^2. \end{aligned}$$

(15)

Thermodynamic basis of TER

TER is expected to be proportional to the thermal energy absorbed by the cell, which is invested in the transition from (A) to ($A^\prime$) (transition towards “dead state”). We propose this energy to increase linearly with the time of heat exposure t, and with the rate of energy absorption $k_E(T)$. In a simplified version of the SLD accumulation induced by hyperthermia, the step from (A) to ($A^\prime$) is represented by a single rate process, with a net rate of transition k(T) (proportional to the rate of energy absorption) as depicted in Fig. 1c and expressed by

$$\begin{aligned} \text {TER}= o +a\,tk(T). \end{aligned}$$

(16)

Here, o is the onset of the thermal enhancement ratio, which should converge to one for no HT treatment ($t=0$); a is a parameter that accounts for the tumour size (or the amount of malignant cells) and the intrinsic sensitivity of the cells to RT and HT, and k(T) is the temperature-dependent rate of the sensitisation process. Based on the thermodynamics of protein denaturation, we found the transition rate of this process to increase exponentially with increasing temperature $k(T)=c\,e^{b(T-T_g)}$ (see methods section for details of the model). Such exponential behaviour with $(T-T_g)$ has been observed in previous works, but could not be explained^43,59,60. In this equation, c and b are cell-type dependent parameters and $T_g$ is the dominant transition temperature, i.e. the average melting point of cellular proteins undergoing denaturation. We achieve this theoretical prediction by considering the change of the heat capacity of the proteins as a linear function of the temperature, and not as a constant value as usually assumed in Arrhenius kinetics. The heat capacity of cellular proteins displays a Lorentzian-type function of the temperature^36,37, which can be approximated at first order as linear functions in the vicinity of the melting point^54,57. Remarkably, the melting point $T_g$ in both cases has good correspondence to the calorimetry studies performed by Lepock and collaborators^36,37 where they found the melting point in the hyperthermia treatment to be in the range of 45–48 °C for different mammalian cells. Plugging the obtained transition rate into Eq. (16), the TER reads

$$\begin{aligned} \text {TER}=o+a^\prime t\,e^{b(T-T_g)}, \end{aligned}$$

(17)

with $a^\prime =ac$ for simplicity. This model predicts exponential increase of $\alpha *$ and $\beta *$ with temperature, which is much more pronounced for $\beta *$. These predictions are consistent with experimental results in cell cultures^58,59, and data from human clinical trials^30,52.

Radiosensitising effects are also reflected and quantified by reductions in the $\alpha /\beta$ ratio, which is basically higher for intrinsically more radioresistant cells^30,34. For the combined RT+HT scheme the $\alpha /\beta$ ratio is reduced as a consequence of the enhancement of the sublethal damage over the direct damage. The ratio for the combined treatment then reads $\alpha ^*/\beta ^*=\frac{\alpha /\beta }{\text {TER}}$.

Predictions of experimental data from literature

We tested the performance of our model (Eq. 17) on three experimental data sets that document thermal enhancement values in simultaneous HT+RT treatments for different temperatures. The data comes from three murine biological models, which are helpful in the study of non-cancer epithelial cells (CHO) and mammary carcinoma (C3H and M8013) respectively. In these referenced experimental studies, the heat source for HT was a precision-controlled water-bath. In two of them TER was measured for different treatment times and temperatures, and the third data set presents $\alpha$ and $\beta$ values obtained for various temperatures but just one treatment time.

Thermal enhancement ratio: The first data set was recorded in in vitro 2D cell culture experiments (CHO cell line)⁵⁸, and the second one derived from an in vivo animal study (C3H mammary carcinoma tumour mouse model)⁶⁰. For both datasets we found that our model well predicts the outcome of these studies. Our model (Eq. 17) predicts a linear dependence of TER as HT time t increases for a fixed temperature. As shown in Fig. 2a,c, both datasets display this linear dependency for all tested temperatures, indicating a rate-dependent nature of the TER function. For each temperature, our model anticipates a temperature-dependend slope, i.e. the rate, which is exponential. As can be seen in Fig. 2b,d, CHO cells in vitro and C3H mammary carcinoma tumours in vivo exhibit this exponential behaviour. The parameters and the respective coefficients of determination $R^2$ are summarized in Table 1 for both examples.

Figure 2

(a) and (c) show the linear dependency of the thermal enhancement ratio (TER) on time of exposure for CHO cells in vitro and C3H mammary carcinoma cells in mice tumours in vivo, respectively. The slope of the linear fitting clearly depends on the temperature of the hyperthermia treatment, and the natural logarithm of the slope was plotted as a function of temperature for both datasets in (b) and (d). The linear trend lines show the exponential behaviour of the temperature dependent rate k(T) according to Eq. (6). The data for CHO cells (a, b) and C3H mammary carcinoma (c, d) was extracted from⁵⁸ and⁶⁰, respectively.

Table 1 Parameters of the TER model Eq. (17), obtained from CHO and C3H cell models^58,60.

Thermal modulation of $\alpha$ and $\beta$: The third data set was documented in cultured M8013 murine mammary carcinoma cells⁵⁹. In this study, the normal (non-thermotolerant) cell-line was compared with a thermotolerant modification. The authors determined the radiobiological parameters $\alpha$ and $\beta$ for cells irradiated halfway through a 30 min hyperthermia treatment (temperatures from 42 to 46 °C). In this case, we calculated the thermal enhancement of $\alpha$ and $\beta$ from Eqs. (14) and (15) to test our model:

$$\begin{aligned} TER_{\alpha }=\frac{\alpha (T)}{\alpha } \quad \text {and} \quad TER_{\beta }=\sqrt{\frac{\beta (T)}{\beta }}. \end{aligned}$$

(18)

Here, $\alpha$ and $\beta$ are the radiobiological parameters without HT. To assess the behaviour of the temperature-dependent rate k(T), we calculated $TER_{\alpha (\text { or }\beta )}-1$ for every data point to compare with the result of Eq. (17), which was rearranged for this purpose as follows:

$$\begin{aligned} \text {TER}-1=atk(T)=a^\prime t\,e^{b(T-T_g)}. \end{aligned}$$

(19)

The results shown in Fig. 3 display an exponential dependency of k(T) with temperature in both cases - as predicted by our model. The model parameters ($a^\prime$, b, $T_g$) and the coefficients of determination are presented in Table 2. Notably, the melting temperatures are quite similar for the two sublines, but the main difference comes from the slope of the calorimetry function $b=B/2k_B$, reflecting a possible slower denaturation of cellular proteins in the thermotolerant subline in response to heat. The parameters were adjusted for all TER values obtained from Eqs. (18). However, it must be noted that the authors of this study reported problematic deviations in the measurements of $\alpha$⁵⁹, which may explain the low coefficients of determination shown in Table 2. When the adjustment is made using only the $TER_{\beta }$ experimental points, it improves to $R^2$ = 0.986 and 0.951 for thermotolerant and non-thermotolerant M8013-cells, respectively.

Figure 3

Thermal enhancement $(TER-1)$ as function of the relative temperature $(T-T_g)$ for M8013 mouse mammary carcinoma cells in vitro⁵⁹. (a) Thermotolerant modification of the cell line and (b) Non-thermotolerant cells. Vertical axes displayed in logarithmic scale. The lines are exponential fittings of the $TER_{\alpha }-1$ and $TER_{\beta }-1$ points together.

Table 2 Parameters of the TER model Eq. (19), obtained from M8013 mouse mammary carcinoma in vitro⁵⁹.

We must stress that this linear model is valid in the regime of non-ablative HT (40–46 °C), which is used for radiosensitisation purposes at which heat-induced damage is primarily sublethal^37,40. This sublethal damage infliction and the corresponding cellular responses are complex biological processes that may involve other variables than protein denaturation, especially at low temperatures. Indeed, many recent studies have shown DNA repair mechanisms, cell cycle redistribution and genomic expression alteration to be involved in hyperthermia-induced radiosensitisation^26,42,44,63. Nonetheless, these stress responses are also triggered by temperature-induced chemical reactions. Hence, independent of whether or not these alteration relate to protein destabilisation/denaturation, they are covered by our mathematical function because the thermodynamic approach (based on Eyring’s transition theory with the temperature-dependent heat capacity of the substrate) is also valid and useful to describe the effect of temperature on non-protein denaturation chemical reactions. In this case, the parameters of the model represent the thermodynamic properties of an effective relevant chemical substrate instead of an “average protein” undergoing thermal transformation⁵⁴. This interpretation is supported by the good fitting of our model to three different datasets, even at low temperatures. We therefore claim that the proposed thermodynamic approach is valid and general enough to cover the different types of chemical reactions relevant in HT-induced radiosensitisation.

As can be seen in Tables 1 and 2, the parameter related to the number of cells ($a^\prime = ac$) is set to one for all the tumour cell models. Doing so, one could speculate that the slope in Eq. (17) is completely modelled by the exponential factor which is solely a function of the thermodynamic quantities describing the heat capacity, namely the melting point $T_g$ and the slope of the calorimetry peak. Remarkably, this calorimetry peak is also very similar for the three non-thermotolerant tumour cell models, but lower for the thermotolerant one. This result may indicate the calorimetry peak as a possible marker for cellular thermotolerance. Moreover, we obtained different melting points for CHO (in vitro), M8013 (in vitro), and C3H (in vivo), which are in the range of 46–49 °C. This result is consistent with the findings from Lepock and collaborators that show different melting points for distinct cell types³⁷ in that temperature range. This should be verified by meticulous future experimental work. Calorimetry assays together with systematic TER measurements in various tumour cell models will be particularly relevant in this context, and can lead to a considerable reduction in the number of adjustable parameters. Nonetheless, our model already quite well predicts and explains the modulation of radioresponse caused by HT treatment from thermodynamic principles with at most three adjustable parameters.

It is recognized that tumour tissue is quite complex and heterogeneous with respect to histomorporphology as well as (local) response to treatment. Temperature inhomogeneities within the tissue can constitute a significant difficulty when applying predictive mathematical models in real tumours²⁶. To overcome this problem, the temperature distribution is nowadays simulated in homogeneous tissue subzones pre-identified from CT scans by state-of-the-art treatment-planning software^64,65. During the treatment application, precision thermometry then monitors, verifies and controls the temperature distribution in real time²⁶. The strength of the presented thermodynamic-based mathematical approach is that it well models the TER not only in homogeneous 2D cultures but also in an experimental in vivo tumour. Furthermore, although the underlying data all came from experiments using precision water baths for heating, our mathematical model does not depend on the source of heat since it is based on HT-induced, temperature and time-dependent protein denaturation. For this reason, we believe that our model may be equally valid for other heating techniques using, for example, photo-induced/magnetic nanoparticles, focused ultrasound, microwaves, radio frequencies, or lasers.

Simultaneous thermoradiotherapy supposedly leads to higher TERs than sequential HT + RT modalities^13,22,23,24. This is true for both tumour and normal tissue. Consequently, normal tissue has to be spared to achieve therapeutic gain, which requires precise application of heat to the tumour. This has remained challenging, due to blood flow, re-oxygenation and heat dispersion, despite the fact that precise real-time temperature control and monitoring techniques are already in place³¹. Today, various technologies, which are based on ultrasound, micro or electromagnetic waves, are available for simultaneous HT+RT treatment in the clinical setting and have for example been applied to treat breast cancer (reviewed in⁷) and different types of superficial malignancies²⁹. Others are under further development and/or envisioned for the treatment of different types of surface and deep-tissue tumours (reviewed in³¹). In this vein, new methodologies such as magnetic or photo-induced nanoparticles for HT^{66,67,68,69,70} and the design of precise simultaneous applicators^7,8,9 have opened attractive prospects for implementing precise simultaneous thermoradiotherapy in standard clinical practice. Comprehensive mathematical modelling to better predict treatment outcome - as documented in the present article - will critically contribute to this process towards clinical routine.

Thermal dose

The current standard concept thermal dose unit refers to “equivalent minutes at 43 °C” (CEM43°C). It was proposed by Sapareto and Dewey^71,72 more than 40 years ago based on the empirical Arrhenius activation theory. Given an exposure time t at temperature T, the proposed function estimates the equivalent exposure time necessary to obtain the same biological response at the reference temperature 43 °C. In this approach, the biological response is modelled as chemical reaction, whose rate is empirically proposed as $Ae^{\frac{-E_{a}}{K_B T}}$. Here $E_{a}$ is the activation energy for the transition. Equating the chemical products at 43 °C and the temperature of interest T, and solving the function for the treatment time at 43 °C, the equivalent minutes at 43 °C are:

$$\begin{aligned} CEM43^\circ :=t_{43^{\circ } \text {C}} =t R^{T-43^{\circ } \text {C}}=t e^{(T-43^{\circ } \text {C})\ln (R)}, \end{aligned}$$

(20)

where the factor A is assumed to be independent of time and temperature, and $R=\exp \lbrace \frac{-E_{a}}{K_B T(43+273.15)} \rbrace$ is approximated to $R=4$ for T < 43 °C and to $R=2$ for T > 43 °C. If the temperature profile is not constant, the doses are added over n subintervals at constant temperature $t_{43^\circ \text {C}} =\sum ^n_{i=1} t_i R^{T_i-43^{\circ }\text {C}}$.

The CEM43°C concept is routinely used in the clinical context, although it presents several theoretical and practical problems, particularly when different heating rates are used^73,74,75,76. Firstly, the Arrhenius approach presumes a constant Gibbs energy for transitions to occur at 43 °C; this is not necessarily true in different cell types and tissues. Indeed, in our study, we found strong temperature dependences for the heat capacity, and therefore for the Gibbs energy (see “Temperature dependence of the transition rate” section). In addition, the heating rates affect the calorimetry proles, and thus the thermodynamic properties of cellular components⁷⁷. As described in “Methodology” section, the reaction responsible for radiosensitisation in our model is protein denaturation. This assumption, together with thermodynamic calculations, lead to the transition rate described in Eq. (13). Applying the same approach of CEM43°C, we propose the thermal dose to be defined as the isoeffective time at the melting temperature, i.e. CEM$T_g$. Using Eq. (13) to equate the chemical products ($t_{T_g}k(T_g)=t k(T)$), and solving for the time of the treatment at $T_g$ (reference temperature) we get:

$$\begin{aligned} D_T&:=t_{T_g}=t \frac{T}{T_g} \exp \left\{ \frac{\Delta G_c(T-T_g)}{K_B T T_g} \right\} \exp \left\{ b(T-T_g)\right\} \nonumber \\&\simeq t \exp \left\{ b(T-T_g)\right\} . \end{aligned}$$

(21)

Since all the temperatures are in Kelvin scale, the fraction $T/T_g\simeq 1$ and the factor $\frac{\Delta G_c(T-T_g)}{K_B T T_g}\simeq 1$ in the range between 40 and 50 °C. Comparing Eqs. (20) and (21) we find that our approach and CEM43°C are equivalent when the responsible reaction is protein denaturation, the transition temperature is the melting point $T_g=43^\circ$ C, and R is associated with the slope of the calorimetry curve $b=\ln (R)$. Defining the quantities in this way, our temperature-dependent transition rate multiplied by the treatment time serves as an alternative thermal dose. In combination with properly defined chemical potential for the denaturated proteins, this thermal dose has the prospective to be translated into Gray units (Gy=J/Kg), namely absorbed energy per mass of tissue. Within this proposition, the thermal enhancement depends linearly with the new defined thermal dose, that depends not only on the time-temperature combination, but on the thermodynamic properties on the specific cell line/type.

We tested our thermal dose concept ($D_T$) on the datasets highlighted in the previous section (see Fig. 2)^58,60 and on additional HT monotreatment data (without radiation), also extracted from the study of Dikomey and Jung⁵⁸. To determine the performance of our dosimetry approach, we rescaled the time axis for the different temperatures according to Eq. (20) ($\hbox {CEM43}^\circ \hbox {C}$) and 21 (Dt) and compared the fitting results. Overall, the two concepts reflect the different types of in vitro and in vivo data similarly well, with an advantage for one or the other concept in distinct datasets (see coefficient of determination $R^2$ in Figs. 4 and 5). Both dosimetry concepts presented outliers in some cases. In our $D_T$ approach the outliers correspond to the in vitro data points at 43 °C. In the original work of Dikomey and Jung⁵⁸, they observed a biphasic behaviour in treatment response with a transition at this temperature, which was accounted for in CEM43°C by a different value for the parameter R in Eq. (20) for temperatures below and above 43 °C. This problem could be similarly solved in the $D_T$ approach for the respective datasets by defining two values of b in Eq. (21). However, the necessity and stability of such a general cutoff temperature for modelling therapeutic outcome remains to be further elucidated and proven because our basic concept performs as well as CEM43°C without a transition temperature, and even more proficiently describes the TER data in vivo.

Figure 4

Survival curves of cells after HT treatment. Data extracted from reference⁵⁸ (a) $-\ln (S)$ as function of HT treatment time for different temperatures. The lines correspond to quadratic fittings. (b) $-\ln (S)$ as function of the thermal dose as defined in Eq. (20). (c) $-\ln (S)$ as function of the thermal dose as defined in Eq. (21).

Figure 5

Thermal enhancement TER as function of the thermal dose as defined by Eq. (20) ((a) and (c)) and Eq. (21) ((b) and (d)). (a) and (b) show the performance of both dosimetry concepts for CHO cells in vitro⁵⁸. (c) and (d) show the performance for C3H mammary carcinoma cells in mice tumours in vivo⁶⁰. The filled (orange) symbols in (b) correspond to the TERs obtained at 43 °C.

Perspectives

The model introduced here for simultaneous treatment is based on the modulation of the radiobiological parameters of the LQ-model. It is suitable to reproduce clonogenic survival curves, and very well reproduces TERs of in vitro and in vivo experiments. However, translation into more relevant clinical outcomes is still required, i.e. tumour control probabilities or control doses, where “disease control” means the long-term extinction of replication-competent tumour cells in vivo after completion of the treatment³⁵. This translation is usually done by means of simple logistic functions. Since the extraction of systematic radiobiological parameters from humans is difficult, even unfeasible for different treatment times and temperatures, the underlying parameters often come from classical 2D cell cultures. However, these have been found to insufficiently estimate radiation response^78,79. Indeed, more accurate translations require more elaborated approaches to reflect the treatment response in a more realistic and complex in vivo-like environment. Examples of factors that might be affected by HT include cell-cell and cell-matrix interactions, oxygen distributions, proliferative activity and cell cycle progression in a 3-D cellular context. We therefore intend to combine the present theoretical findings with cellular automaton simulations in a next step, to model the treatment outcome in 3D tumour cell models such as in silico multicellular tumour spheroids.

Our encouraging findings implicating a putative predictive power of our model is the basis for the targeted implementation of more complex scenarios, i.e., with respect to (1), the heterogeneous tumour cellular environment and micromilieu in a 3D geometry, and (2) the often transient nature of oxygen distribution in tumour tissues. In this context, heterogeneous blood flow, thermal washout, as well as cyclic hypoxia (H) and reoxygenation (R) phenomena need to be incorporated into the experimental design as modulators of local thermotolerance and radioresponse. Tools to systematically alter concentrations of oxygen along with variable durations and frequencies of H–R periods, thereby generating dynamic pathophysiological conditions that better mimic the in vivo situation, are increasingly employed for in vitro research, but are still in their infancy^80,81. Such approaches, especially when combined with sophisticated 3D culturing, will be key to better address the challenge of clinical translation. Moreover, different treatment schedules including fractionated regimes but also sequential treatments with different HT-RT/RT-HT intervals are to be considered. The latter are known to be easier in clinical handling and thus of high practical relevance. The present work paves the ground for a more elaborate unied mathematical model, which is in the focus of our ongoing work, with the aim of describing the individual treatments and their sequential application from common general principles.

Conclusion

Taken together, our model interprets the enhancement of radiotherapy by hyperthermia to result from an increased vulnerability of a cell in the temperature regime between $\sim 40-46^\circ$C. Radiosensitisation is achieved by the accumulation of sublethal damage either repairable or not due to protein denaturation. The increased vulnerability then in turn reduces the survival probability of cells undergoing radiation, and therefore, a larger proportion of the tumour is controlled in the combined scheme. In this model, the synergistic effect quantified by TER is proportional to the energy invested to induce the damage. It is proposed to be a rate-dependent effect that increases linearly with the time of HT and exponentially with temperature in the aforementioned temperature range. Our model offers a thermodynamics-based approach to explain previous experimental observations. Despite the tumour heterogeneity in vivo and the complex cellular response to thermoradiotherapy, the present work shows that thermodynamic principles of chemical reactions, including but not limeted to protein denaturation, can explain to a good extent the TER both in vitro and in vivo. It constitutes a crucial step for implementing more complex scenarios, in which thermodynamic reversible/repairable effects need to be considered for subsequent treatment planning.