Mammalian cell sensitivity to hyperthermia in various cell lines: a new universal and predictive description

Abstract

Introduction

The Cumulative Equivalent Minute at 43 °C (CEM43) thermal dose model has been empirically derived more than 30 years ago and still serves as a benchmark for hyperthermia protocols despite the advent of regulatory network models. However, CEM43 suffers from several limitations regarding its inability to predict the effect of complex time varying profiles (thermotolerance, step-down heating), to predict synergistic effects with drug treatments or to explain the specificity of a cell line in thermal resistance.

Objective

Define a new generic predictive tool for thermal injury based on regulatory network models. Identify the biological parameters that account for the thermal resistance.

Materials

Comparative study of cell survival upon hyperthermia collected from literature (17 sets in 11 publications that cover 14 different cell lines from 8 different tissues).

Results

A dynamical model describes accurately cell survival according to the amplitude and duration of exposure but also molecular chaperone expression level. In the case of square shape hyperthermia, approximated analytical expression of the cell survival is derived from the dynamical model and compared to CEM43 description. The molecular chaperone expression level defines the thermal resistance of a given cell line and can be estimated from a single experimental result through an easy-to-use graphical tool.

Conclusion

The tools offered here can be useful for designing treatments combining hyperthermia and chemotherapy targeting molecular chaperones, but also for designing personalized hyperthermic treatment by prior biochemical screening of molecular chaperones. These tools could advantageously replace the description of CEM43.

Keywords:

• Cell survival; modelling; thermal dose; hyperthermia; chaperones

1. Introduction

Quantification of thermal injury of tissues is important for a wide variety of intentional and unintentional thermal exposures. In the oncological environment, thermal-based treatments are generally subdivided into two groups with respect to the targeted tissue temperature. Ablation therapy relies on the direct cytotoxic effect of temperature elevation (in general above 50 °C) [1]. In contrast, hyperthermia-based therapies rely on elevating the tissue temperature between 40 and 45 °C to either kill cells or enhance the efficiency of other treatments [1,2].

Two methods have proved valuable for quantifying the relationship between thermal exposure and damages. The damage index (or the 50 °C isothermal contour) is used by the ablation community, whereas the cumulative equivalent minutes at 43 °C (CEM43) is typically used by the hyperthermia community [1,2]. Both the damage index and CEM43 methods are based on the underlying premise that tissue damages follow an irreversible first-order chemical reaction, with the rate constant following an Arrhenius relationship. These methods have been used for the last decades, with some rationalizations between them reported in the published literature [3–6].

Experimental cell survival results have accumulated for many cell lines in recent decades [1,7–14]. They reveal that the survival response strongly varies among tissues but also among cell types of a given tissue. Breast cancer cells for instance are much more resistant than colon cancer cells: a 10 min exposure at 50 °C leads to 60% of MCF10A cell survival [14] but fewer than 5% of RWPE cell survival [13]. In the framework of phenomenological descriptions (CEM43 or Damage index), the variation of survival response with cell line is described by a simple parameter shift [15], without providing any interpretation regarding the biological origins of these differences. However, all cells share the same architecture of the heat shock response network and cell death. Such variability raises several questions: Can we identify the key biological parameters that account for the cell type-specific thermal sensitivity? How can cell sensitivity to thermal exposure be enhanced or inhibited?

To answer those questions, we perform a comparative study of cell survival upon hyperthermia based on the literature survey and mathematical modeling. Comparative study based on CEM43 parameters cannot help: CEM43 parameters are not related to any intracellular properties (heat induced-damage, repair or damage-dependent effect on division, or survival pathways). In contrast, a dynamical model has recently been proposed to accurately describe cell survival data in HeLa cell line while overcoming some CEM43 limits [16]. This model is derived from a data-driven biochemical network model of the heat-shock response in HeLa cells [17], and connects the cell survival to elementary biochemical processes like synthesis, repair, degradation and the overall protein abundance of molecular chaperones. This is the reason why we use this dynamic model to perform mathematical modeling in the comparative study; the fact that the model components can be specifically targeted and modulated by drug therapies are also a clear motivation for its use.

First, we collect cell survival datasets from the literature (17 sets in 11 publications that cover 14 different cell lines from 8 different tissues listed in Table 1 and replotted in Figure 2). These datasets are finely adjusted by a low-dimensional dynamical model where the Molecular Chaperone Expression Level (MCEL) is the prevailing biological parameter that varies among cell lines and explains their relative thermal resistance. Finally, we perform a mathematical analysis of the model to formulate a rescaled exposure time as a basis of comparison with CEM43, and also present an easy-to-use graphical tool to predict the cell survival upon both hyperthermia and anti-chaperone treatment.

Table 1. Collected experimental data sets.

Cell line	Organism	Tissue	Disease	Ref	Label
L1210	mouse	skin	lymphocytic leukemia	[9]	L1210
CHO-10B	chinese hamster	ovary		[10]	CHO^a
CHO-10B	chinese hamster	ovary		[1]	CHO^b
HTB-66	human	skin	malignant melanoma	[1]	HTB66
PC3	human	prostate	adenocarcinoma	[13]	PC3
RWPE-1	human	prostate	normal	[13]	RWPE
SW 620	human	colon	colorectal adenocarcinoma	[12]	SW620
WiDr	human	colon	colorectal adenocarcinoma	[12]	WiDr^a
WiDr	human	colon	colorectal adenocarcinoma	[11]	WiDr^b
HeLa	human	cervix	adenocarcinoma	[7]	HeLa^a
HeLa	human	cervix	adenocarcinoma	[8]	HeLa^b
MCF-7	human	breast	adenocarcinoma	[14]	MCF7
MDA-MB-231	human	breast	adenocarcinoma	[14]	MDAMB
MCF-10A	human	breast	fibrocystic disease	[14]	MCF10
U-87MG	human	brain	glioblastoma-astrocytoma	[11]	U87MG
CCD-18Lu	human	lung	normal fibroblast	[11]	18LU
A-549	human	lung	carcinoma	[11]	A549

^a,b

Different studies of the same cell line.

2. Material and methods

2.1. Literature survey

Experimental data were collected from literature (Table 1). Publications containing graphs depicting survival of cells as a function of hyperthermic stress duration for different temperatures were selected when the data points could be quantified precisely enough by data extraction software. Only single-dose rectangular-shaped temperature profiles were considered. Each experimental point consists of three coordinates: survival fraction or percentage S, hyperthermia temperature T_C, and duration of hyperthermia Δ. In addition, where available, time after exposure to hyperthermia before measuring survival has been taken into account, and reproduced in the numerical experiments. Table 1 contains the list of cell lines together with references and label identification code used in the remaining of this paper. Experimental cell survival results are re-plotted as open circles in Figure 2.

2.2. Dynamical model of the cell survival upon hyperthermia

The generic dynamical model of cell survival consists of two differential equations describing the time evolution of the damage D (the amount of misfolded protein), and the cell population S: (1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) (1b) $$\frac{\mathrm{d} S}{\mathrm{d}t}=-\alpha D S.$$(1b)

In the generic dynamical model of cell survival Equations (1)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) where t denotes the time in hours from the stress onset, T denotes the temperature in Celsius degrees, T = T_C for $0\le t\le \Delta$ and 37 °C otherwise. The temperature-dependent function f denotes the temperature-dependent production rate of damage D where f is normalized such that $f(T={37}^{°}\mathrm{C})=1.$ The function f has been previously estimated to $f(T)=\frac{10}{6}(1-0.4 e^{37-T})\times {1.4}^{T-37}$ [18].

Equation (1a)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) accounts for three biochemical processes: (1) the denaturation process, which creates misfolded proteins D from heat-labile proteins with a heat-dependent rate $\sigma f(T);$ (2) the degradation process, which destroys misfolded proteins without repair trials with a rate $\delta \frac{D^2}{H+D};$ (3) the renaturation process, which restores back misfolded proteins into heat-labile proteins, with a rate $\gamma \frac{H D}{H+D}.$

Equation (1a)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) includes four cell line-dependent parameters: σ, the amount of damages created per hours at 37 °C, H the abundance of molecular chaperones; δ, the life time of damage in the high damage limit (where $D\gg H$); and γ the renaturation rate in the low damage limit (where $D\ll H$). In the high damage limit, the renaturation rate is limited by the MCEL H and the maximal renaturation rate is $\gamma H.$ Alternatively, the maximal renaturation rate $\gamma H$ may be expressed as the renaturation rate for a so-called breakup temperature T_B as $\gamma H=\sigma f(T_B),$ defining T_B by an inverse problem. The breakup temperature roughly indicates the temperature above which the renaturation saturates and then the degradation dominates.

The second equation (Equation 1b(1b) $$\frac{\mathrm{d} S}{\mathrm{d}t}=-\alpha D S.$$(1b) ) describes the cell population S decaying with a rate $\alpha D S$ proportional to damages. The population equation (Equation 1b(1b) $$\frac{\mathrm{d} S}{\mathrm{d}t}=-\alpha D S.$$(1b) ) does not include cell division process because the initial value of S is set to $S(0)=\hspace{0.17em}\text{exp}\hspace{0.17em}(\alpha D_{37} t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})$ where D₃₇ is the D value in a steady state at 37 °C and $t_{\hspace{0.17em}\text{exp}\hspace{0.17em}}$ defines the experiment duration. Therefore, (2) $$S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=S(0)\times \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\alpha {\int }_0^{t_{\hspace{0.17em}\text{exp}\hspace{0.17em}}}D(t)dt\right)$$(2) corresponds to the cell survival measured by cell colony assay, where $S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=1$ for a constant exposure at 37 °C. The parameter α quantifies the cellular lethality in response to a given accumulated damage, which may differ from one cell line to another.

To quantitatively describe the cell survival response obtained over decades in different labs, we also add to Equation (1)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) an extra equation which accounts for the unknown incubator or water bath rise time. So, the temperature is assumed to follow the equation $\frac{\mathrm{d}T}{\mathrm{d}t}=(T_c-T)/{\tau }_{\text{inc}}$ where T_C denotes the setpoint temperature and ${\tau }_{\text{inc}}$ the incubator rise time. The parameter ${\tau }_{\text{inc}}$ is estimated together with the core model parameters.

2.3. Parameter estimation procedure

In the first setting, the dynamical model defined by Equation (1)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) displays five parameters (σ, δ, γ, H, and α) which may vary among cell lines. This extensive parameterization (Equation 1(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) ) is likely to suffer from a lack of parameter identifiability in the absence of quantitative data on the misfolded protein concentration induced by the hyperthermia protocols. To circumvent this issue, we perform a normalization procedure for damage and chaperone to express them in an adapted concentration unit $u=\sigma /\delta ;$ all cells have therefore the same stationary levels of effective temperature-dependent damages $D^{*}/u=f(T)$ in the high damage regime ($D\gg H$). We also assign the parameter γ to its estimated value of γ = 1200 h^–1 [19] in U2OS, as the renaturation rate is very unlikely to differ between cell lines.

Upon these two assumptions, we obtain a generic dynamical model with only three parameters left: $h=H/u$ the MCEL, $\kappa =\gamma /\delta$ the rescaled renaturation rate; and $\rho =\alpha u.$ It is worth noting that the maximal renaturation rate defines a breakup temperature $T_{\mathrm{B}}=f^{-1}(\kappa h)$ by an inverse problem. Using breakup temperature $T_{\mathrm{B}}$ conveniently determines the thermal resistance of a given cell line from parameter values.

Only three parameters retain to be estimated (κ, h, and ρ, in addition to the incubator temperature rise time τ_inc), which is sufficient to adjust the cell survival experimental set for various cell lines (solid lines in Figure 2). The parameters of the model defined in Equation 1(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) (κ, h, ρ, ${\tau }_{\text{inc}}$) are optimized to adjust each experimental data-set through Levenberg-Marquardt optimization Algorithm (LMA) [20]. LMA minimizes the score function defined as a Root-Mean-Square Error (RMSE) (Equation 3(3) $$F_{sc}(\kappa ,h,\rho ,{\tau }_{\text{inc}})=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}-S_{\text{num}}\left(\kappa ,h,\rho ,{\tau }_{\text{inc}}\right)\right)}^2},$$(3) ) between the experimental data and the numerical results (open circles and line of Figure 2). (3) $$F_{sc}(\kappa ,h,\rho ,{\tau }_{\text{inc}})=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}-S_{\text{num}}\left(\kappa ,h,\rho ,{\tau }_{\text{inc}}\right)\right)}^2},$$(3) where $S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}$ and $S_{\text{num}}$ are experimental and numerical survival fractions, respectively. N is the total number of experimental points for a given experiment. Since LMA is a local optimization procedure and can strongly depend on the initial point, the method is applied to a hundred different initial conditions and the selected parameter set corresponds to the smallest score function. For each experimental data-set, we obtain an optimal parameter set $({\kappa }^{\text{opt}},h^{\text{opt}},{\rho }^{\text{opt}},{\tau }_{\text{inc}}^{\text{opt}})$ such that: (4) $$({\kappa }^{\text{opt}},h^{\text{opt}},{\rho }^{\text{opt}},{\tau }_{\text{inc}}^{\text{opt}})=\text{argmin}\left(F_{sc}\left(\kappa ,h,\rho ,{\tau }_{\text{inc}}\right)\right).$$(4)

In addition, the parameter distributions are generated by Markov Chain Monte Carlo (MCMC) method through Metropolis-Hastings algorithm [21] for each cell line; and represented as Box and whisker plots on Figure S1 of the Supplemental Material (displayed in gray) to visualize the confidence intervals.

The incubator rise time ${\tau }_{\text{inc}}$ value is always well determined with a narrow distribution around the typical time scale of a standard incubator ($\le 10$ min); the estimated values are moreover similar for the different experimental results arising from the same lab. Regarding the free parameters (κ, h, and ρ), we found a weak local identifiability, characterized by parameter distributions covering several orders of magnitude, for breast and prostate cancer and for CHO in particular. Weak local identifiability undoubtedly indicates over-fitting. We therefore use a reductionist approach to reduce the parameter space dimensions and avoid over-fitting.

The systematic parameter estimation procedure reveals that, for most cell lines, κ and ρ parameter distributions almost overlap from one cell line to another (Figures. S1b and S1e) whereas only h, and consequently T_B, vary with cell lines. Therefore, we restart the parameter estimation process with ρ and κ defined as global parameters i.e., with common value for all cell lines, whereas only h and ${\tau }_{\text{inc}}$ are optimized with cell-line specific values.

In practice, parameter non-identifiability is reduced by defining two global parameters $({\kappa }^{\text{glob}},{\rho }^{\text{glob}}).$ For each of $({\kappa }^{\text{glob}},{\rho }^{\text{glob}})$ parameter values, optimal parameter set $(h^{\text{opt}},{\tau }_{\text{inc}}^{\text{opt}})$ is re-estimated using LMA: (5) $$(h^{\text{opt}},{\tau }_{\text{inc}}^{\text{opt}})=\text{argmin }\left(F{\prime }_{sc}\left(h,{\tau }_{\text{inc}}\right)\right),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}F{\prime }_{sc}(h,{\tau }_{\text{inc}})=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}-S_{\text{num}}\left({\kappa }^{\text{glob}},h,\rho ,{\tau }_{\text{inc}}^{\text{glob}}\right)\right)}^2},$$(5) while the optimal values of the global parameters $({\kappa }^{\text{glob},\text{opt}},{\rho }^{\text{glob},\text{opt}})$ are also estimated using LMA: (6) $$({\kappa }^{\text{glob},\text{opt}},{\rho }^{\text{glob},\text{opt}})=\text{argmin}\left(\sqrt{\frac{1}{N_{\hspace{0.17em}\text{exp}\hspace{0.17em}}}\sum _{j=1}^{N_{\hspace{0.17em}\text{exp}\hspace{0.17em}}}F{\prime }_{sc}(h^{\text{opt}},{\tau }_{\text{inc}}^{\text{opt}})}\right),$$(6) where $N_{\hspace{0.17em}\text{exp}\hspace{0.17em}}=17$ is the number of experiments (see Table 1). Parameter distributions for h and ${\tau }_{\text{inc}}$ are generated by MCMC–method through Metropolis-Hastings algorithm with a fixed value of the global parameters $({\kappa }^{\text{glob},\text{opt}},{\rho }^{\text{glob},\text{opt}}).$

This procedure ends up with $\rho =0.831$ h⁻¹ and κ = 209 as optimal global parameters, while h and ${\tau }_{\text{inc}}$ values still vary significantly across cell lines with the same trend as previously.

If we assume that the h median value (around 0.1) corresponds to a typical concentration of HSP70 in mammalian cells (around 1 μM [22,23]), we get an estimation of $u=10 \mu$ M and thus are able to define all parameter values: $\alpha =\rho /u=0.0831$ μM^–1 h⁻¹; $\delta =\gamma /\kappa =5.74$ h⁻¹; and $\sigma =\delta u=57.4$ μM h⁻¹. Table 2 summarizes the parameter values or expressions used in the original model (given by Equation 1(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) ) as well as their used estimation methods. The choice of u-value does not affect the quality of the description but permit to fix a parameter set.

Table 2. Summary of model Equation (1)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) expression or parameter values and their estimation methods.

Param.	Value or expression	Description	Cell dependent	Estimation method
f(T)	$\frac{10}{6}(1-0.4e^{37-T}){1.4}^{T-37}$		No	Literature [18]
γ	1200 h^–1	Renaturation rate	No	Literature [19]
δ	5.74 h^–1	Degradation rate	No	Global fit
σ	57.4 μM.h^–1	Denaturation rate	No	Global fit
α	300 μM^–1.h^–1	Cell sensitivity to damage	No	Global fit
H	$[0.4-2.6]\mu \mathrm{M}$	MCEL	Yes	Local fit
$T_{\mathrm{B}}$	$[42.5-47.5]{}^{°}\mathrm{C}$	Breakup temperature	Yes	From H value

The concentration values are obtained by assuming that h = 0.1 corresponds to H = 1 μM e.g., by setting $u=10$ μM.

The relevance of this procedure of parameter space reduction is assessed by a good RMSE score adjustment for all cell lines except those that were already poorly described (Supplementary Figure S1A), and by a much stronger local identifiability (Supplementary Figure S1B-E). Even though two parameters are freely adjusted, the differences between cell lines can be essentially explained by the relative MCEL H in this cell survival upon hyperthermia model because the incubator time scale is only related to the experimental set-up.

2.4. Steady state analysis of survival curves

We aim to establish a detailed comparison between this new description and the standard CEM43 description. So, we first derive an analytical expression of the cell survival obtained in the generic dynamical model under two assumptions. First, we assume quasi-equilibrium for dynamical variables D and S. Second, we assume constant cell populations at 37 °C during the exposure duration Δ, thereby neglecting the normalization by the control population at 37 °C during the exposure duration Δ.

At a constant temperature T, damage can reach its fixed point value $D^{*},$ which is a stationary solution of Equation (1a)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) , $\dot{D}=0.$ Analytical expression of the fixed point can be expressed in a simple setting by the use of the rescaled temperature parameter $x=\frac{f(T)}{f({\mathrm{T}}_{\mathrm{B}})}$ where $f({\mathrm{T}}_{\mathrm{B}})=\frac{\gamma H}{\sigma }:$ (7) $$D^{*}(x)=\frac{\gamma H}{\delta } \frac{1}{2}\left((x-1)+\sqrt{{\left(x-1\right)}^2+4\frac{\delta x}{\gamma }}\right).$$(7)

Table 3 summarizes the characteristic values of damage concentration in steady states. With parameter values given in Table 2, one finds a 10% decrease of the control cell population (at 37 °C) over one cell cycle (24 h) due to the heat-induced damage at 37 °C.

Table 3. Steady state concentration of damages as function of temperature.

T	x	$D^{*}(x)$	$D^{*}(T)$
37 °C		$\frac{\delta }{\gamma }$	$\frac{\delta }{\gamma }$
$f(T)\ll f(T_B)$	$x\ll 1$	$x H$	$\frac{\delta }{\gamma } f(T)$
$T=T_{\mathrm{B}}$	x = 1	$\frac{\gamma }{\sqrt{\delta \gamma }} H$	$\frac{\sigma }{\sqrt{\delta \gamma }} f(T_{\mathrm{B}})$
$f(T)\gg f(T_B)$	$x\gg 1$	$\frac{\gamma }{\delta } x H$	$\frac{\sigma }{\delta } f(T)$

How fast damage evolves to the fixed point depends on parameter $1/\delta$ (here $1/\delta =10$ min). In the limit $1/\delta \ll \Delta ,$ damage reaches the fixed point value fast enough compared to the duration of hyperthermic stress, so that the contribution to the integral of damage during the transient dynamics at the moment and after the stress is negligible compared to the contribution during the stress. The integral of damage can, then, be approximated by the rectangular function of duration Δ and high $D_T^{*},$ such that $S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=\hspace{0.17em}\text{exp}\hspace{0.17em}(-\alpha D_T^{*} \Delta )$ can be rewritten as: (8) $$S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\Delta }_r}{{\Delta }_0}\right),$$(8) where (9) $${\Delta }_r=\Delta \sqrt{\frac{\gamma }{\delta }} \frac{1}{2}\left((x-1)+\sqrt{{\left(x-1\right)}^2+4\frac{\delta x}{\gamma }}\right)$$(9) is a rescaled exposure duration at T_B and (10) $${\Delta }_0=\frac{\sqrt{\delta }}{\alpha \sqrt{\gamma }H}$$(10) defines the survival for one hour exposure at T_B temperature. The rescaled exposure duration corresponds to the exposure duration (${\Delta }_r=\Delta$) when the temperature corresponds to the breakup temperature (T_C = T_B).

For each cell line, we compute the CEM43 for each experimental data point by using the standard expression $\mathit{CEM}43=\Delta {(R_{\text{cem}})}^{T-43}$ with $R_{\text{cem}}=\{\begin{array}{cc}0.233\hfill & T\le 43\hfill \\ 0.428\hfill & T>43\hfill \end{array},$ and estimate the D₀ parameter by a linear regression analysis.

3. Results

3.1. Survival response of cell lines mostly depends on their chaperone concentrations

A generic model of cell survival, schematically drawn in Figure 1, accounts for the heat inhibition of the cell survival via the heat induced misfolded proteins. Three molecular biochemical processes are used: (1) the denaturation process, which creates misfolded proteins D from heat-labile proteins with a heat-dependent rate; (2) the degradation process, which destroys misfolded proteins without repair trials; (3) the renaturation process, which restores back misfolded proteins into heat-labile proteins. The rate of the last two processes vary with the MCEL H. In this model, it is assumed that misfolded proteins D inhibit the cell survival S. This study aims to describe at best the cell survival experimental data in various cell lines upon transient hyperthermia protocol of duration Δ and setpoint temperature T_C.

Figure 1. Cartoon of the generic model description of cell survival. Biochemical processes (orange arrows) control the misfolded protein concentration (D). Biochemical rates are modulated (red arrows) by heat (T) or by MCEL (H). The misfolded proteins inhibit the cell survival (S). The network globally describes a heat inhibited cell survival.

The model given by Equation (1)(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) describes suitably the overall collected cell survival data, but with varying degree of accuracy depending on cell lines and parameters. The best data adjustments are obtained for human breast, cervix and colon cancer, melanoma, and mouse leukemia cell lines, while the data adjustment of human prostate cancer and Chinese hamster ovary cell lines being less satisfactory. The available data are nevertheless not always enough to really constrain the optimal parameters, independently of the data description goodness.

By a reductionist approach based on a parameter estimation procedure, we identify the MCEL H as the key cell type dependent parameter, all other parameters of the model being fixed to consensus values. With this single free parameter, the generic dynamical model succeeds to describe the collected cell survival data (Figure 2). For each cell line, the estimated value of H with confidence interval (quartiles) is indicated directly on the corresponding panel in Figure 2 and reported in Table S1 of the Supplemental Material.

Figure 2. Generic model description of cell survival in various cell lines. Each panel corresponds to a cell survival experimental data-set (open circles) taken from the literature and displays the cell survival as a function of the exposure duration (Δ) at various maximum temperature T_C. The corresponding experiment is identified directly on the panel by its abbreviation code according to Table 1. A color code is used to indicate the maximum temperature (see the legend on the top). Best descriptions via the generic model are displayed as solid lines the estimated MCEL and the confidence interval (quartiles of the MCMC parameter distributions) are displayed on the top of each panel.

To ease the interpretation, the H value may alternatively be converted into a so–called breakup temperature T_B. The breakup temperature indicates the threshold temperature value above which the degradation process dominates the renaturation process in a steady state. In this framework, thermal resistant cell lines are thus characterized by high value of H and T_B (see Table S1 of the Supplemental Material).

Up to now, the H is a model parameter value, mathematically adjusted to fit the experimental data. To test our prediction and link the estimated H values with the HSP70 protein abundance for various cancer types, we compare for each tissue the estimated values with Reverse Phase Protein Arrays (RPPA) measurements obtained from Cancer Proteome Atlas (Figure 3). The good correlation between the two validates the here used modeling and reductionist approach.

Figure 3. Correlation between estimated H values and HSP70 protein abundance measured by Reverse Phase Protein Arrays (RPPA) (data are taken from The Cancer Proteome Atlas [24]). Each rectangle corresponds to a different human tissue. The width of a rectangle corresponds to the dispersion of the h parameter values (given in Figure 2) over the various cell lines coming from the same human tissue (or the confidence interval if the cell line is unique), the height of a rectangle corresponds to the interquartile range of the HSP70 protein abundance in cancer cells of the tissue obtained from The Cancer Proteome Atlas [24]. The color of a rectangle defines the tissue as indicated directly on the figure.

In the collected data sets, three cell lines (HeLa, WiDR, and CHO) include duplicate results obtained in different labs. The CHO data are consistent and are characterized in our framework by a similar value of H, while HeLa and WiDR data display significant discrepancies between the labs. In particular, the fact that R.J. Palzer et al. survival data (HeLa^a) display a shoulder, whereas E.G. Gerner et al. survival data (HeLa^b) do not, is here explained by a lower value of MCEL. Such a variability of the MCEL in the same cell line but different labs could arise from different cell culture conditions. It is more difficult to derive conclusion about WiDR because the data are poorly described by the reduced model. This is probably because the model does not take into account the slow (>3 h) transcriptional upregulation of molecular chaperones, while WiDr^b data consider much longer exposure (up to 180 h).

3.2. Rationalization of thermal injury quantification methods

As demonstrated in the previous result sections, the generic dynamical model is a suitable description of the cell survival upon transient hyperthermia in various cell lines. For those hyperthermia protocol, and under reasonable assumptions (see M&M section for details), one can derive an analytical expression of the cell survival obtained in the generic dynamical model so as to establish a detailed comparison between this new description and the standard CEM43 description. The two descriptions are summarized in Table 4 where the constant parameter values have been replaced by their values for clarity. In both descriptions, one free parameter remains to be determined: the dose parameter D₀ in the CEM43 description and the MCEL H in the new description.

Table 4. Thermal injury quantification tools.

CEM43 Description
$S=\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{\text{CEM}43}{D_0}\right)$	$\frac{\text{CEM}43}{D_0}=\Delta \frac{1}{D_0}{(R_{\text{cem}})}^{T-43}$
	$R_{\text{cem}}=\{\begin{array}{cc}0.233\hfill & T\le 43\hfill \\ 0.428\hfill & T>43\hfill \end{array}$
	D0 parameter to determine
Generic Model Description
$S=\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\Delta }_r}{{\Delta }_0}\right)$	$\frac{{\Delta }_r}{{\Delta }_0}=\Delta 10.3 \left(\mathcal{H}-H-\sqrt{{(\mathcal{H}-H)}^2+0.02 \mathcal{H} H}\right)$
	$\mathcal{H}=0.08(1-{0.4}^{37-T}) {1.4}^{T-37}$
	H parameter to determine

In the here developed description Δ_r is a new rescaled exposure duration at ${\mathrm{T}}_{\mathbf{B}};$ and ${\Delta }_0$ defines the survival for one hour exposure at T_B temperature. The rescaled exposure duration corresponds to the exposure duration (${\Delta }_r=\Delta$) when the temperature corresponds to the breakup temperature (${\mathrm{T}}_{\mathbf{C}}={\mathrm{T}}_{\mathbf{B}}$). Δ_r is thus similar to CEM43 whereas $\frac{{\Delta }_r}{{\Delta }_0}$ is equivalent to the damage index. Detailed mathematical expressions of Δ_r and ${\Delta }_0$ are presented in the Material and Methods section only.

According to those descriptions all experiments may be recast on one simple curve by plotting $D_0\hspace{0.17em}\text{log}\hspace{0.17em}(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}})$ vs. CEM43 in the CEM43 description (Figure 4) or by plotting ${\Delta }_0\hspace{0.17em}\text{log}\hspace{0.17em}(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}})$ vs. Δ_r in the here developed description. As shown in Figure 4, if both descriptions give satisfactory results, the description developed here slightly increases the quality of the description.

Figure 4. Survival and time rescaling. The dose survival response for various cell-line experiments (see the legend for color code) plotted as a function of cumulative equivalent minutes at 43 °C according (upper panel) or as function of equivalent minute at $T_{\mathbf{B}}$ according to Equation (8)(8) $$S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\Delta }_r}{{\Delta }_0}\right),$$(8) (lower panel). D₀ and CEM43 (upper panel) or ${\Delta }_0$ and Δ_r (lower panel) are estimated from regression.

In more details, the CEM43 and Δ_r are always positively correlated even if their values for the same experiment may significantly differ (see Figure S2 of the Supplemental Material). A difference between CEM43 and Δ_r appears when the value ${\mathrm{T}}_{\mathrm{B}}$ found by the optimization process moves away from the 43 °C of CEM43; this is especially clear for the breast (MCF10A, MCF7, and MDA) and prostate (PC3 and RWPE) cell lines (see Table 1 of the Supplemental Material). The temperature dependency of CEM43 and Δ_r is highlighted in the Figure 2 of the Supplemental Material by plotting the quantity $-X \hspace{0.17em}\text{log}\hspace{0.17em}(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}})/\Delta$ where X refers to D₀ in the CEM43 description and ${\Delta }_0$ in the generic model description (Table 4). In the CEM43 description, two different regimes are predicted: below and above the breakpoint (43 °C here), whereas the Generic Model Description predicts three different regimes: $\mathcal{H}\ll H,$ $\mathcal{H}\simeq H,$ and $\mathcal{H}\gg H$ (Table 4). Unfortunately, none of these experiments displays a sufficiently wide range of exposure temperature to discriminate between the two predictions.

3.3. Graphical determination of molecular chaperones concentration from a single experimental result

Based on the analytical results derived in the previous section, one can build a universal graphical tool (Figure 5) to estimate the H parameter value from one single cell survival experiment. Let us assume that one measures the cell survival $S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}$ in another cell line for one transient thermal exposure at temperature T for an exposure duration Δ, for instance 37% of cell survival for 1 h at 42 °C. Then by computing $-\hspace{0.17em}\text{log}\hspace{0.17em}(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}})/\Delta$ (1 in the previous example), and using Equation 8(8) $$S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\Delta }_r}{{\Delta }_0}\right),$$(8) , one can estimate the H parameter value from the overlapping of the experimental result and the H-isolines ($H\simeq 0.4$ in the previous example), without having to implement fitting procedures. Obviously, the validity of this graphical approach depends on the conditions previously defined: the exposure duration must be much longer than the incubator rise time, and the precision of the H value determination is maximal where the H isolines are well separated on the graph.

Figure 5. Graphical determination of H value from a single experimental result. The H parameter value can be graphically determined from one cell survival experimental result by using this abacus. For instance, a 37% cell survival ($S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}=e^{-1}$) for a 1 h exposure at 42 °C corresponds to $H\simeq 0.4$ (brown continuous line). Continuous lines are obtained by using Equation (8)(8) $$S(t_{\hspace{0.17em}\text{exp}\hspace{0.17em}})=\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\Delta }_r}{{\Delta }_0}\right),$$(8) with κ = 209 and $\rho =0.831,$ the color code defines the H value (see the legend on the top).

The Figure 5 may also be used to predict graphically the effect of a drug treatment targeting molecular chaperones. In the framework of the previous example ($S_{\hspace{0.17em}\text{exp}\hspace{0.17em}}=37$%, for $T={42}^{°}$C and Δ = 1 h), Figure 5 predicts that a two-fold inhibition of molecular chaperones $H\to 2$ leads to a 4.5-fold increase of $-\hspace{0.17em}\text{log}\hspace{0.17em}(S_{\hspace{0.17em}\text{exp}\hspace{0.17em}})/\Delta$ and thus 1% of cell survival. The Figure 5 clearly identifies the parameter region where molecular chaperone inhibition strongly decreases the cell survival upon hyperthermia e.g., the parameter area where a significant vertical shift is obtained by moving from one H isoline to another without changing the temperature.

To test the graphical determination method of H, we search in the literature for cell survival experimental results on a given cell line with quantified varying MCELs and we found a complete experimental data set obtained in the field of thermotolerance study. It is well known that a first non-lethal heat treatment induce the HSP70 synthesis and thus increases the cell survival obtained in a second treatment. Long time ago, G. Li and J. Mak contributed to establish this knowledge by studying the induction of Heat Shock Protein synthesis during the development of thermotolerance with squamous cell Carcinoma implanted in the flank of C3H mice ([25]). These authors applied a first priming heat treatment in vivo (41 °C for 60 min, 42 °C for 40 min, 43 °C for 15 min, or 44.5 °C for 10 min) followed by graded second treatments at 45 °C in vitro where they performed cell survival assay. The various thermal protocols are schematically displayed in Figure 6. In parallel they measure the relative rate of HSP70 synthesis induced by the first priming heat treatment. As shown in Figure 6, the H parameter value that is estimated from survival assay in each thermal protocol by using Figure 5, is highly correlated with the HSP70 synthesis rate measured by G. Li and J. Mak. Therefore, the here presented model is a valuable tool to predict the thermotolerance once the molecular chaperone induction is known.

Figure 6. Thermotolerance induced by HSP70 induction. The upper panel schematically exhibits the thermal protocols used to study the induction of HSP70 during the development of thermotolerance. The first priming treatment corresponds to no priming treatment (black), 41 °C for 60 min (red), 42 °C for 40 min (light blue), 43 °C for 15 min (orange), 44.5 °C for 10 min (dark blue). For each first priming treatment, survival assay is performed through three different second treatments: (1) no second treatment, (2) 45 °C for 30 min, or (3) 45 °C for 45 min. Three cases are considered: (1) no priming treatment (filled square), (2) no recovery between the two treatments (filled circles), and (3) a 2 h recovery between the two treatments (filled stars). The lower panel highlights the correlation between the rate of HSP70 synthesis and the estimated H value from survival assay.

4. Discussion

Thermal dose tools quantify hyperthermia lethality, whether caused by intentional exposure or not, and assist treatment design and side damage prevention. The widely used CEM43 description suffers from numerous limitations acknowledged for long, thereby restricts elaboration of new treatment strategies. In particular, CEM43 leaves out thermotolerance, step-down heating effect, and more generally lacks of accuracy in case of complex time varying hyperthermia [16,26]. Furthermore, the CEM43 phenomenological origin disconnects parameters values from fundamental biological or biochemical processes, and thus leaves the thermal sensitivity variation with cell line unexplained [27]. Such a disconnectedness harms new therapies design as treatments personalization. So, beyond the simple rational satisfaction, unveiling such links would indeed assist chemotherapeutic hyperthermia design, either by designing the drug that specifically targets the pathway responsible for thermal resistance, or by designing hyperthermic protocols that reckon on the chemotherapy-induced downregulation of thermal sensitivity. Unveiling such links would also assist treatments personalization in which identified biomarkers would quantify thermal resistance of tumors based on biopsy [28].

To overcome those CEM43 limitations while keeping its accuracy, we propose here a dynamic and mechanistic modeling framework as a CEM43 alternative. The new description may, thus, be used by either performing numerical integration of the dynamical model (Equation 1(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) ) or by using the approximate analytical description (Table 4 and spreadsheet provided in the Supplemental Material) or by using the graphical determination methods (Figure 5). However, some limitations of CEM43 still remain in this new thermal injury quantification tool: both tools focus on cell survival and are not suitable to describe nor predict the improved perfusion and oxygenation in tumor after a mild hyperthermia [29,30]. The new description uses MCEL as a master regulator of thermal sensitivity, in line with parameter analysis and experimental cross checking. The new description therefore predicts cell survival in the event of hyperthermia in the first place, but also in the event of MCEL inhibition. Because the HSP70 chaperone function prevents heat-induced apoptosis [31] and lysosomal cell death [32,33], HSP70s indeed negatively influence antiproliferative and apoptosis-inducing therapies [34,35]. HSP70 has therefore been suggested as a prognosis biomarker [36–39] and molecular chaperones inhibition has also been suggested and experimentally tested as a promising drug target for anticancer treatment [40–44]. The new description could for instance be used to optimize hyperthermic intraperitoneal chemotherapy design [45].

Over a population or a tumor, MCEL cell-to-cell variability deeply influences thermal cytotoxicity: cells having a higher MCEL being more thermal resistant [46,47]. This known result differs from the one derived here: average MCEL of a cell line explains thermal cytotoxicity. Highlighting those two facts consistency requests to consider two kinds of MCEL variation. The first one to consider is the MCEL variation in a cell population (cell-to-cell variability). The MCEL lognormal distribution induces the cell-to-cell variability (phenotypic heterogeneity) of the heat shock network dynamic response [48] and of the thermal cytotoxicity [46,47]. In our opinion, the MCEL cell-to-cell variability explains the fractional killing: only cells having higher values of the MCEL survive. The second one to consider is the variation of the average MCEL from one cell type to another (cell line variability) e.g., a variation of the center of the MCEL distribution. This variation defines thermal sensitivity of a given cell line. In the cell population model given by Equations (1a,b(1a) $$\frac{\mathrm{d} D}{\mathrm{d}t}=\sigma f(T)-\delta \frac{D^2}{H+D}-\gamma \frac{H D}{H+D},$$(1a) ), a death probability directly accounts for fractional killing and the surviving fraction thus is one of the model variables. The H-parameter therefore corresponds to the average value of MCEL. In contrast, a single cell model would have given as an output a binary death (dead or alive) and the cell-to-cell variability of MCEL would have given fractional killing.

The CEM43 description is an isoeffect tool: CEM43 defines a relation between amplitude and exposure time to obtain the same cell survival fraction. Performing numerous colony formation tests, aiming to determine the parameter d₀, is indeed an essential prerequisite for cell survival quantification from CEM43. Consequently, in the laboratory, with each new cell line, with each new condition of laboratory culture, but also in the clinic, with each new tumor or each new patient, the parameter d₀ must be recalibrated beforehand to target a cell survival setpoint by tuning an hyperthermic protocol. This is obviously very tedious and almost impossible to implement. We believe that these new tools could constitute a valuable step toward personalized hyperthermia protocols. Let’s imagine having, in the lab, a control cell culture with well-established cell survival assay and thus well-estimated H-parameter value used as a reference (named $H_{\text{ref}}$). Then, to design the hyperthermia treatment, the value of the H-parameter of the tumor considered can be estimated relatively to $H_{\text{ref}},$ by a biochemical measurement (immunohistochemistry, western blot or PCR) targeting HSP70. To clarify, let us take the example of a measurement by ELISA indicating an over-expression of 20% of the HSPs in the tumor biopsy whereas the reference line is characterized by an H value of 0.4; in this case the tumor will be characterized by a value of the parameter H of 0.48. This protocol is much faster than colony assays and could therefore be implemented routinely. Once the H value is calibrated, the proposed abacus (Figure 5) for instance can be used to design hyperthermic protocols to achieve the targeted cell kill.

The reliability and accuracy of survival model predictions are essential expectations for designing optimal personalized hyperthermic treatments. As shown here, the predictions of the mathematical model are in good quantitative agreement with the experimental data, three possible refinements may nevertheless be pointed out. First, the reductionist approach used, which keeps only one free parameter to capture cell line-specific properties, can be considered as an extreme. Undoubtedly, some other differences between different cell lines contribute in some ways to thermal sensitivity, as extracellular pH or nutrient levels. However, the MCEL may still be the master regulator of thermal sensitivity, for instance the MCEL increases with pH [49] while the cell survival upon hyperthermia decreases whith pH [50], in consistency with the model predictions. Clear quantification of these secondary factors would require, for each cell line, joint measurements of chaperone concentration and cell survival to track discrepancies. Second, the denaturation rate considered has, however, been established by indirect fitting procedure rather than direct measurement [18] and, up to our knowledge, no information is available concerning the precision of its formulation. Detailed investigation on the denaturation function could, thus, be an important source of model refinement. Finally, third, the regulated synthesis of molecular chaperones, which is responsible for the phenomenon of thermotolerance [51], is missing in the present model. This simplifying assumption was convenient to reduce the number of model parameters, and thus avoid overfitting given the experimental data available. As we have shown, the model nevertheless depicts this phenomenon provided the concentration of molecular chaperone is correctly and manually adjusted for each pulse of hyperthermia. Yet, the model could be extended to implement active transcription of molecular chaperones, which would require a well-supplied set of experimental data to reach a survival prediction model for fractional treatments, which would be a major advance compared at CEM43.

To conclude, we present a new tool to define the cell survival in various cell lines upon hyperthermia and/or anti-chaperones treatment that challenges the CEM43 description accuracy while overcoming some CEM43 limitations and which offers refinement perspectives through extended experimental studies.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work has been supported by the LABEX CEMPI [Grant/Award Number: ANR-11-LABX-0007] and by the Ministry of Higher Education and Research, Hauts de France council, French Ministry of Higher Education and Research and European Regional Development Fund (ERDF) through the Contrat de Projets Etat-Region (CPER Photonics for Society P4S).