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.


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][4][5][6].
Experimental cell survival results have accumulated for many cell lines in recent decades [1,[7][8][9][10][11][12][13][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 induceddamage, 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 easyto-use graphical tool to predict the cell survival upon both hyperthermia and anti-chaperone treatment.

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 D. 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

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: In the generic dynamical model of cell survival Equations (1) where t denotes the time in hours from the stress onset, T denotes the temperature in Celsius degrees, T ¼ T C for 0 t D 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 CÞ ¼ 1: The function f has been previously estimated to f ðTÞ ¼ Equation (1a) accounts for three biochemical processes: (1) the denaturation process, which creates misfolded proteins D from heat-labile proteins with a heat-dependent rate r f ðTÞ; (2) the degradation process, which destroys misfolded proteins without repair trials with a rate d D 2 HþD ; (3) the renaturation process, which restores back misfolded proteins into heat-labile proteins, with a rate c H D HþD : Equation (1a) includes four cell line-dependent parameters: r, the amount of damages created per hours at 37 C, H the abundance of molecular chaperones; d, the life time of damage in the high damage limit (where D ) H); and c the renaturation rate in the low damage limit (where D ( H). In the high damage limit, the renaturation rate is limited by the MCEL H and the maximal renaturation rate is c H: Alternatively, the maximal renaturation rate c H may be expressed as the renaturation rate for a so-called breakup temperature T B as c H ¼ r 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) describes the cell population S decaying with a rate a D S proportional to damages. The population equation (Equation 1b) does not include cell division process because the initial value of S is set to Sð0Þ ¼ exp ða D 37 t exp Þ where D 37 is the D value in a steady state at 37 C and t exp defines the experiment duration. Therefore, corresponds to the cell survival measured by cell colony assay, where Sðt exp Þ ¼ 1 for a constant exposure at 37 C. The parameter a 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) an extra equation which accounts for the unknown incubator or water bath rise time. So, the temperature is assumed to follow the equation dT dt ¼ ðT c ÀTÞ=s inc where T C denotes the setpoint temperature and s inc the incubator rise time. The parameter s inc is estimated together with the core model parameters.

Parameter estimation procedure
In the first setting, the dynamical model defined by Equation (1) displays five parameters (r, d, c, H, and a) which may vary among cell lines. This extensive parameterization (Equation 1) 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 ¼ r=d; all cells have therefore the same stationary levels of effective temperature-dependent damages D Ã =u ¼ f ðTÞ in the high damage regime (D ) H). We also assign the parameter c to its estimated value of c ¼ 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, j ¼ c=d the rescaled renaturation rate; and q ¼ a u: It is worth noting that the maximal renaturation rate defines a breakup temperature T B ¼ f À1 ðj hÞ by an inverse problem. Using breakup temperature T B conveniently determines the thermal resistance of a given cell line from parameter values.
Only three parameters retain to be estimated (j, h, and q, in addition to the incubator temperature rise time s 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 (j, h, q, s 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) between the experimental data and the numerical results (open circles and line of Figure 2).
where S exp and S 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 ðj opt , h opt , q opt , s opt inc Þ such that: 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 s inc value is always well determined with a narrow distribution around the typical time scale of a standard incubator ( 10 min); the estimated values are moreover similar for the different experimental results arising from the same lab. Regarding the free parameters (j, h, and q), 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 overfitting.
The systematic parameter estimation procedure reveals that, for most cell lines, j and q 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 q and j defined as global parameters i.e., with common value for all cell lines, whereas only h and s inc are optimized with cell-line specific values.
In practice, parameter non-identifiability is reduced by defining two global parameters ðj glob , q glob Þ: For each of ðj glob , q glob Þ parameter values, optimal parameter set ðh opt , s opt inc Þ is re-estimated using LMA: while the optimal values of the global parameters ðj glob, opt , q glob, opt Þ are also estimated using LMA: ðj glob, opt , q glob, opt Þ ¼ argmin where N exp ¼ 17 is the number of experiments (see Table 1). Parameter distributions for h and s inc are generated by MCMC-method through Metropolis-Hastings algorithm with a fixed value of the global parameters ðj glob, opt , q glob, opt Þ: This procedure ends up with q ¼ 0:831 h À1 and j ¼ 209 as optimal global parameters, while h and s 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 lM [22,23]), we get an estimation of u ¼ 10 l M and thus are able to define all parameter values: Table 2 summarizes the parameter values or expressions used in the original model (given by Equation 1) as well as their used estimation methods. The choice of uvalue does not affect the quality of the description but permit to fix a parameter set.
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.

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 D, thereby neglecting the normalization by the control population at 37 C during the exposure duration D.
At a constant temperature T, damage can reach its fixed point value D Ã , which is a stationary solution of Equation (1a), _ D ¼ 0: Analytical expression of the fixed point can be expressed in a simple setting by the use of the rescaled tem-  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.
How fast damage evolves to the fixed point depends on parameter 1=d (here 1=d ¼ 10 min). In the limit 1=d ( D, 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 D and high D Ã T , such that Sðt exp Þ ¼ exp ðÀa D Ã T DÞ can be rewritten as: where is a rescaled exposure duration at T B and defines the survival for one hour exposure at T B temperature. The rescaled exposure duration corresponds to the exposure duration (D r ¼ D) when the temperature corresponds to the breakup temperature ( For each cell line, we compute the CEM43 for each experimental data point by using the standard expression

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 D and setpoint temperature T C . The model given by Equation (1) 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.
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.
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).  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 0 in the CEM43 description and the MCEL H in the new description.

Rationalization of thermal injury quantification methods
In the here developed description D r is a new rescaled exposure duration at T B ; and D 0 defines the survival for one hour exposure at T B temperature. The rescaled exposure duration corresponds to the exposure duration (D r ¼ D) when the temperature corresponds to the breakup temperature (T C ¼ T B ). D r is thus similar to CEM43 whereas D r D0 is equivalent to the damage index. Detailed mathematical expressions of D r and D 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 log ðS exp Þ vs. CEM43 in the CEM43 description (Figure 4) or by plotting D 0 log ðS exp Þ vs. D 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. In more details, the CEM43 and D 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 D r appears when the value T 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 D r is highlighted in the Figure 2 of the Supplemental Material by plotting the quantity ÀX log ðS exp Þ=D where X refers to D 0 in the CEM43 description and D 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 : H ( H, H ' H, and H ) H (Table 4). Unfortunately, none of these experiments displays a sufficiently wide range of exposure temperature to discriminate between the two predictions.   Table 4. Thermal injury quantification tools.  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.  determination is maximal where the H isolines are well separated on the graph. 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 exp ¼ 37%, for T ¼ 42 C and D ¼ 1 h), Figure 5 predicts that a two-fold inhibition of molecular chaperones H ! 2 leads to a 4.5-fold increase of À log ðS exp Þ=D 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.

CEM43 Description
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 nonlethal 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.

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, stepdown 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) 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][37][38][39] and molecular chaperones inhibition has also been suggested and experimentally tested as a promising drug target for anticancer treatment [40][41][42][43][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), a death probability directly accounts for fractional killing and the surviving fraction thus is one of the model variables. The Hparameter 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 0 , 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 0 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 wellestimated H-parameter value used as a reference (named H ref ). Then, to design the hyperthermia treatment, the value of the H-parameter of the tumor considered can be estimated relatively to H 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 antichaperones 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).