Theoretical evaluation of the impact of diverse treatment conditions by calculation of the tumor control probability (TCP) of simulated cervical cancer Hyperthermia-Radiotherapy (HT-RT) treatments in-silico

Abstract

Introduction

Hyperthermia (HT) induces various cellular biological processes, such as repair impairment and direct HT cell killing. In this context, in-silico biophysical models that translate deviations in the treatment conditions into clinical outcome variations may be used to study the extent of such processes and their influence on combined hyperthermia plus radiotherapy (HT + RT) treatments under varying conditions.

Methods

An extended linear-quadratic model calibrated for SiHa and HeLa cell lines (cervical cancer) was used to theoretically study the impact of varying HT treatment conditions on radiosensitization and direct HT cell killing effect. Simulated patients were generated to compute the Tumor Control Probability (TCP) under different HT conditions (number of HT sessions, temperature and time interval), which were randomly selected within margins based on reported patient data.

Results

Under the studied conditions, model-based simulations suggested a treatment improvement with a total CEM43 thermal dose of approximately 10 min. Additionally, for a given thermal dose, TCP increased with the number of HT sessions. Furthermore, in the simulations, we showed that the TCP dependence on the temperature/time interval is more correlated with the mean value than with the minimum/maximum value and that comparing the treatment outcome with the mean temperature can be an excellent strategy for studying the time interval effect.

Conclusion

The use of thermoradiobiological models allows us to theoretically study the impact of varying thermal conditions on HT + RT treatment outcomes. This approach can be used to optimize HT treatments, design clinical trials, and interpret patient data.

Keywords:

• Hyperthermia
• radiotherapy
• biological modeling
• treatment conditions
• tumor control probability

1. Introduction

Hyperthermia (HT) is a well-known radiosensitizer that aims to heat the tumors to 39–43 °C for a certain period (usually an hour). The benefits of combining HT with radiotherapy (RT) have been clinically proven for different cancer locations [1–7]. Furthermore, several in-vitro studies have demonstrated that HT induces various biological mechanisms that may explain the obtained clinical improvement, such as repair impairment [8–10], direct HT cell killing [10–13] or reoxygenation of hypoxic cells [14–16].

Nevertheless, in-vitro studies correspond to simplified scenarios (isolated 2D cell cultures) with well-controlled conditions (constant temperature and oxygen level, known number of initial cells, etc.). Patients constitute a more complex scenario with heterogeneous conditions and processes simultaneously occurring at the molecular [17–19], cellular [20, 21] and tissue level [22–24]. Moreover, patients treated with HT in combination with RT exhibit variations in the conditions (e.g., temperature, number of HT sessions or time interval between HT and RT) that may affect HT-related biological processes. Therefore, identifying how each biological process affects treatment outcomes in conventional clinical studies is challenging.

In addition to patient data analyses, biophysical models may investigate the translation of deviations in treatment conditions into outcome variations (e.g., Scheidegger et al. [25] presented a model to study the immune system response under different HT + RT schedules). These models can be used to study the theoretical impact of altered conditions on treatment outcome [25–29], design clinical trials [30–33] and to test the required model complexity [34] and model parameter resolution [35–37].

Diverse in-silico models have been proposed to study the cellular effects of HT + RT [38–43]. For example, van Leeuwen et al. [44] developed an extended linear-quadratic (LQ) model to optimize and evaluate the quality of different HT + RT treatment plans using Equivalent Dose (EQD) calculations. Additionally, several publications from this research group evaluate HT + RT plans for individual patients by performing 3D voxel-based EQD calculations to account for heterogeneous radiation dose and temperature distributions [44–49].

Our study primarily focused on inter-patient heterogeneities in treatment conditions instead of those among individual patients. Therefore, in a model-based simulation, we studied how variations in HT-related treatment conditions (number of HT sessions, temperature, and time interval) can influence the outcome of HT + RT cervical cancer treatments. To this end, the extended LQ model presented by van Leeuwen et al. [44] was used to perform Tumor Control Probability (TCP) calculations. First, four HT + RT schedules were simulated as examples to study the impact of low-/high-temperature HT sessions on the treatment outcome. Additionally, we examined the TCP dependence on the temperature, time interval, and total thermal dose (CEM43). Finally, the importance of including inter-patient heterogeneity was investigated.

2. Material and methods

2.1. In-silico radiosensitization and direct HT cell killing model

Van Leeuwen et al. [44] defined a model to fit clonogenic data (i.e., measurements of the cell’s ability to duplicate after exposure to a specific stress [50]) at different temperatures and time intervals. This model redefines the parameters of the LQ model [51, 52] to describe RT cell killing including the HT radiosensitization effect (SF_RT): (1) $${\mathit{\text{SF}}}_{\mathit{\text{RT}}}\left(D,T,t_{\mathit{\text{int}}}\right)=e^{-\alpha \left(T,t_{\mathit{\text{int}}}\right)·D-\beta \left(T,t_{\mathit{\text{int}}}\right)·D^2}$$(1) with (2) $$\alpha \left(T,t_{\mathit{\text{int}}}\right)={\alpha }_{37}\cdot \mathit{\text{exp}}\left[\frac{T-37}{41-37}\cdot \mathit{\text{ln}}\left(\frac{{\alpha }_{41}}{{\alpha }_{37}}\right)\cdot \mathit{\text{exp}}\left(-\mu \cdot \left|t_{\mathit{\text{int}}}\right|\right)\right]$$(2) (3) $$\beta \left(T,t_{\mathit{\text{int}}}\right)={\beta }_{37}·\mathrm{\text{exp}}\left[\frac{T-37}{41-37}·\mathit{\text{ln}}\left(\frac{{\beta }_{41}}{{\beta }_{37}}\right)·\mathrm{\text{exp}}\left(-\mu ·\left|t_{\mathit{\text{int}}}\right|\right)\right]$$(3) where ${\alpha }_{37}=\alpha \left(37,0\right)$, ${\alpha }_{41}=\alpha \left(41,0\right)$, ${\beta }_{37}=\beta \left(37,0\right)$, ${\beta }_{41}=\beta \left(41,0\right)$, T (°C) is the temperature, $\mu$ (h⁻¹) is the rate at which the radiosensitization effect of HT disappears, $t_{\mathit{\text{int}}}$ (h) is the time interval between HT and RT and D (Gy) is the radiation dose per fraction.

Assuming an Arrhenius relationship [20] and a first-order reaction, direct HT cell killing (SF_HT) is given by (4) $$S{F}_{\mathit{\text{HT}}}\left(T,t_{\mathit{\text{HT}}}\right)=\mathit{\text{exp}}\left[-k\left(T\right)\cdot t_{\mathit{\text{HT}}}\right]$$(4) where t_HT (s) is the heating time, and k(T) is the rate of direct hyperthermic cytotoxicity, defined as (5) $$k\left(T\right)=2.05\cdot {10}^{10}\cdot \left(T+273.15\right)\cdot \mathit{\text{exp}}\left[\frac{\Delta S}{2}-\frac{\Delta H}{2\cdot \left(T+273.15\right)}\right]$$(5) where T (°C) is the temperature, ΔS (cal/°C/mol) is the entropy of inactivation and ΔH (cal/mol) is the inactivation energy of the critical rate-limiting molecules that cause cell lethality.

Therefore, the final survival after a HT + RT treatment (SF) is given by (6) $$\mathit{\text{SF}}={\mathit{\text{SF}}}_{\mathit{\text{RT}}}\left(D,T,t_{\mathit{\text{int}}}\right)·{\mathit{\text{SF}}}_{\mathit{\text{HT}}}\left(T,t_{\mathit{\text{HT}}}\right)$$(6)

2.2. Used model parameters

Van Leeuwen et al. [44] used the SiHa and HeLa cervical cancer cell lines to perform cell survival experiments at diverse temperatures (39–42 °C), time intervals (10–240 min) and treatment sequences (HT before or after RT) to calibrate the model parameters (Table 1). These cell lines were selected because they represent one of the most complete data-sets available in the literature, and considering that both cell lines correspond to cervical cancer, our study is representative of this tumor location.

Table 1. Model parameters used in the presented simulations corresponding to the values obtained by van Leeuwen et al. [44] for SiHa and HeLa cell lines. All the model parameters (except ΔH which was fixed) were obtained from a regression model analysis and, for simplicity, their confidence intervals were not considered in this study.

Parameter	SiHa cell line	HeLa cell line
${\alpha }_{37}\left(G{y}^{-1}\right)$	$0.386$	0.230
${\alpha }_{41}\left(G{y}^{-1}\right)$	$0.668$	0.398
${\beta }_{37}\mathrm{ }\left({\mathit{\text{Gy}}}^{-2}\right)$	$2.16\times {10}^{-2}$	$4.6\times {10}^{-2}$
${\beta }_{41}\left(G{y}^{-2}\right)$	$8.84\times {10}^{-3}$	$2.62\times {10}^{-2}$
$\mu \mathrm{ }\left(h^{-1}\right)$	$0.027$	0.053
$\Delta S\left(\mathit{\text{cal}}/°C/\mathit{\text{mol}}\right)$	$392.08$	432.36
$\Delta H\mathrm{ }\left(\mathit{\text{cal}}/\mathit{\text{mol}}\right)$	$147,907.8$	160 105.6

Additionally, van Leeuwen et al. [44] stated that the obtained decay constants ($\mu$) are relatively low in comparison with other in-vitro [53] and in-vivo [54] data, where time decay constants around 0.5 h⁻¹ were obtained. Therefore, a value of $\mu =0.5 h^{-1}$ was also considered to study the case where the time interval has a higher impact. Figure 1 and Supplementary Figure S1 present a graphical representation of the time interval effect at different temperatures.

Figure 1. Survival fraction of the SiHa cell line at diverse time intervals between HT and RT after a radiation fraction of 2 Gy and a HT treatment at 39 °C (blue line), 41 °C (orange line) or 43 °C (red line) and with a decay constant of $\mu =0.027 h^{-1}$ (left) or $\mu =0.5 h^{-1}$ (right). Additionally, the survival baseline without HT (black line) is included and the time interval is positive or negative if HT is given before or after RT, respectively.

Figure 1 and Supplementary Figure S1 show that direct HT cell killing prevents survival from returning to the RT baseline when the radiosensitization effect fades away. Additionally, because of the absolute value included in the time interval term of Equations (2)(2) $$\alpha \left(T,t_{\mathit{\text{int}}}\right)={\alpha }_{37}\cdot \mathit{\text{exp}}\left[\frac{T-37}{41-37}\cdot \mathit{\text{ln}}\left(\frac{{\alpha }_{41}}{{\alpha }_{37}}\right)\cdot \mathit{\text{exp}}\left(-\mu \cdot \left|t_{\mathit{\text{int}}}\right|\right)\right]$$(2) and (3)(3) $$\beta \left(T,t_{\mathit{\text{int}}}\right)={\beta }_{37}·\mathrm{\text{exp}}\left[\frac{T-37}{41-37}·\mathit{\text{ln}}\left(\frac{{\beta }_{41}}{{\beta }_{37}}\right)·\mathrm{\text{exp}}\left(-\mu ·\left|t_{\mathit{\text{int}}}\right|\right)\right]$$(3) , the survival is independent of the treatment sequence (i.e., HT before or after RT). This assumption was consistent with the experimental data presented by van Leeuwen et al. [44]. However, for other cell lines, applying HT before or after RT may have an impact on the radiosensitizing effect; therefore, this statement may not be generalizable to all cell lines.

2.3. Tumor control probability (TCP) calculation

Using the previously described thermoradiobiological model and assuming complete repair between consecutive irradiation fractions, the number of surviving cells after each fraction (N) was calculated as follows: (7) $$N=N_0·\mathit{\text{SF}}=N_0·{\mathit{\text{SF}}}_{\mathit{\text{HT}}}\left(T,t_{\mathit{\text{HT}}}\right)·e^{-\alpha \left(T,t_{\mathit{\text{int}}}\right)·D-\beta \left(T,t_{\mathit{\text{int}}}\right)·D^2}$$(7) where N₀ is the initial number of cells, D is the dose per fraction, and SF is the cell survival for a particular fraction calculated using the previously described model.

The final number of surviving cells (N_final) after treatment was calculated as follows: (8) $$N_{\mathit{\text{final}}}=N_0·\prod _{i=0}^{N_{\mathit{\text{fractions}}}}{\mathit{\text{SF}}}_i$$(8) where N_fractions is the total number of fractions, and SF_i is the survival of the i-th fraction.

Finally, the TCP was calculated following a Poisson distribution [34]: (9) $$\mathit{\text{TCP}}=e^{-N_{\mathit{\text{final}}}}$$(9)

2.4. Thermal dose CEM43 calculation

Sapareto et al. [55] derived an equation to estimate a thermal dose that translates different time–temperature treatments to an equivalent exposure time at 43 °C with the same HT-related biological effect. The thermal dose is based on the Arrhenius relationship [20] and is calculated as follows: (10) $$\mathit{\text{CEM}}43={\int }_0^t{R}^{\left(43-T\right)}\mathit{\text{dt}}$$(10) where T (°C) is the temperature and R is the factor to compensate for a 1 °C temperature change (R = 0.5 for T > 43 °C and R = 0.25 for T ≤ 43 °C).

2.5. Simulated treatment conditions

In the Cervical Cancer Dutch Deep Hyperthermia Trial [1], most patients had tumor diameters between 60 and 80 mm. Hence, considering a spherical tumor shape with a 70 mm diameter, a tumor volume of 180 cm³ was used, and assuming an average volume of 2$\times$10⁻⁹cm³ per cell, the initial number of cells (N₀) was equal to 9 × 10¹⁰.

Regarding the HT + RT treatments, two RT schedules consisting of 30 and 40 fractions of 2 Gy were assumed for the SiHa and HeLa cell lines, respectively. For each cell line, the same RT schedule was used for all simulated patients to obtain a baseline TCP (RT only) around 0.5. These RT schedules are clinically relevant when compared with the Cervical Cancer Dutch Deep Hyperthermia Trial [1] (cumulative pelvic side wall dose up to 60 Gy) and German guidelines for cervical cancer treatment [56] (at least an equivalent dose of 85 Gy should be achieved). For the HT treatment, the ranges of study for the HT treatment conditions (Table 2) were established based on the data reported in the reviews of Carrapiço-Seabra et al. [57] and Ademaj et al. [58]. These margins were defined by considering data (number of HT sessions, all temperature metrics, and time intervals) from diverse cancer locations treated with deep HT to cover most of the reported patients. Finally, a one-hour steady state was assumed for all HT sessions.

Table 2. Margins established for the thermal treatments based on patient data reported in the reviews of Carrapiço-Seabra et al. [57] and Ademaj et al. [58].

Number of HT sessions	HT temperatures	Time interval between HT and RT
3-6	39–43 °C	10–240 min

2.6. Impact of low/high temperature HT sessions

From the temperatures reported in Table 2, the minimum (39 °C), the mean (41 °C) and the maximum (43 °C) temperature values were used to define a low, medium and high temperature HT session, respectively. Four different schedules were used as examples to study the impact of low- and high-temperature HT sessions (Figure 2):

Figure 2. Graphical scheme of the four different schedules used to study the impact of low/high temperature HT sessions for a diverse number of HT sessions. Each circle corresponds to a HT session at 39 °C (blue), 41 °C (green) or 43 °C (red).

1. Baseline schedule with all HT sessions at a medium temperature (41 °C).

2. Schedule with one HT session at a low temperature (39 °C) and the rest at a medium temperature (41 °C).

3. Schedule with one HT session at a high temperature (43 °C) and the rest at a medium temperature (41 °C).

4. Schedule with a HT session at a high temperature (43 °C) and the rest at a low temperature (39 °C).

It is essential to highlight that the calculations presented in this manuscript do not depend on the frequency (e.g., once or twice a week) or order of HT sessions.

2.7. TCP dependencies

Diverse simulated patients were generated to study the TCP dependencies on the treatment conditions.

2.7.1. Temperature dependence

10⁴ simulated patients were included in the simulation. Each patient was treated with four HT sessions, and in each session, the temperature was randomly selected from a uniform distribution between 39 and 43 °C (Table 2). A fixed time interval between HT and RT of 30 min was considered for all simulated patients.

2.7.2. Time interval dependence

10⁴ simulated patients were included in the simulation. Each patient was treated with four HT sessions; in each session, the time interval was randomly selected from a uniform distribution between 10 and 240 min (Table 2). This procedure was repeated at three fixed temperatures (39, 41, and 43 °C).

2.7.3. Temperature and time interval dependence

Simulated patients were generated with 3–6 HT sessions, and in each session, the temperature and time intervals were randomly chosen from uniform distributions within the limits shown in Table 2. Ultimately, 10⁴ simulated patients were selected for each number of HT sessions.

2.7.4. Thermal dose dependence

Simulated patients were generated with 3–6 HT sessions, and in each session, the temperature was randomly selected from a uniform distribution between 37 and 43 °C. In this case, the minimum allowed temperature was set as 37 °C (instead of the 39 °C reported in Table 2) to obtain thermal doses closer to 0 and to better show the TCP rise at low thermal doses. Four cases were studied, considering fixed time intervals of 10, 60, 120, and 240 min. Finally, to ensure sufficient statistics for all the thermal doses,10⁶ simulated patients were generated for each number of HT sessions and time interval. The TCP and CEM43 values were calculated, and the simulated patients were divided into total thermal dose (CEM43) intervals to calculate the mean TCP value.

2.8. Inter-patient heterogeneity

Two sources of inter-patient heterogeneity were considered: tumor volume and radiosensitivity. For the former one, a Gaussian distribution (µ = 201.6 cm³ and σ = 119.7 cm³) was defined from the tumor diameters reported in the Cervical Cancer Dutch Deep Hyperthermia Trial [1] and considering a spherical shape to calculate the tumor volumes. Conversely, regarding radiosensitivity heterogeneities, higher stability is observed in model parameter ratios (e.g., α/β) than in the individual parameters (e.g., α or β) [59]. Hence, for simplicity, only a variation in the parameter ${\alpha }_{37}$ was considered and the other radiosensitivity-related parameters (i.e., ${\alpha }_{41}$, ${\beta }_{37}$ and ${\beta }_{41}$) were adjusted to maintain constant ratios concerning ${\alpha }_{37}$. Therefore, the values of ${\alpha }_{37}$ were randomly selected from a Gaussian distribution with a mean value as presented in Table 1 and a coefficient of variation ($c.v.=\frac{\mathit{\text{standard deviation}}}{\mathit{\text{mean}}}\times 100\%$) of 5% or 50% to simulate scenarios of low or high inter-patient heterogeneity, respectively. Additionally, from in-vitro experiments with cervical cancer cell lines [44, 60–67], ${\alpha }_{37}$ values between 0.2 and 0.74 Gy⁻¹ were obtained. Hence, only ${\alpha }_{37}$ values within this range were considered in this study.

Figure 3 demonstrates that the radiosensitivity heterogeneity had a greater impact on the dose-control relationship than tumor volume heterogeneity. Nevertheless, both inter-patient heterogeneities were considered to obtain more realistic results.

Figure 3. Obtained dose-control curves for simulated patients (without HT). A case without inter-patient heterogeneity (left) and with inter-patient heterogeneity (right) are presented. Two cases assume inter-heterogeneity only in the radiosensitivity (red line) and the tumor volume (green line). The SiHa cell line model parameters were considered for this analysis.

Finally, because of the lower radiosensitivity of the HeLa cell line, the dose-control curve is shifted to the right in comparison to the curve for SiHa cell line (Supplementary Figure S2). This is why a higher number of RT fractions should be delivered to HeLa cells to have a similar TCP.

2.9. Statistical analyses

The Spearman’s rank order correlation (ρ) and Kendall’s rank correlation coefficient (τ) were used to evaluate possible TCP dependencies on diverse treatment conditions (e.g., minimum, mean or maximum time interval/temperature). Finally, because of the large number of simulated patients, p-values smaller than 10⁻⁵ were obtained in all calculations; therefore, they are not reported in this manuscript.

3. Results

Various scenarios were evaluated to study the influence of different parameters on the TPC calculation for cervical cancer HT + RT treatments. First, the impact of low/high temperature HT sessions was simulated (Section 3.1). Second, the dependencies of TCP on the diverse treatment conditions and the total thermal dose (CEM43) were studied (Section 3.2). Finally, the inclusion of inter-patient heterogeneity (Section 3.3) was considered.

The analysis in this study focused on two cervical cancer cell lines, SiHa and HeLa. In light of the similar results obtained for both cell lines, the analysis for SiHa is presented in the main text and, most of the analysis for HeLa cells is included in the Supplementary Materials.

3.1. Impact of low/high temperature HT sessions

The schedules described in Section 2.6 were simulated to study the effects of the low-/high-temperature HT sessions.

As expected, the best schedule was one that reached a high temperature (43 °C) in one HT session while maintaining reasonable temperatures (41 °C) in the rest of the sessions (red bars in Figure 4 and Supplementary Figure S3). Nevertheless, similar results were usually obtained by a schedule with a medium temperature (41 °C) in all HT sessions (green bars in Figure 4 and in Supplementary Figure S3). The main difference between these two schedules is the direct HT cell killing produced during the high temperature HT session. This is why the highest differences between both schedules correspond to cases with fewer HT sessions (i.e., individual HT sessions have a higher impact) or cases considering direct HT cell killing and a long time interval (i.e., when the radiosensitization effect has a smaller contribution). Thus, under these conditions, a schedule with a high-temperature HT session and low-temperatures in the remaining HT sessions (orange bars in Figure 4 and Supplementary Figure S3) can obtain a reasonably good outcome. However, when the number of low-temperature HT sessions was increased, the importance of the high temperature HT session is reduced, and this schedule resulted in a smaller TCP value than the other schedules studied. Finally, it was demonstrated that good treatments could still be performed with a low-temperature HT session as long as at least a medium temperature was obtained in the rest of the sessions (blue bars in Figure 4 and Supplementary Figure S3).

Figure 4. TCP obtained for the diverse number of HT sessions considering (right column) or not (left column) direct HT cell killing and with a time interval of 10 (upper row) or 120 (lower row) min. Four different schedules were considered: a schedule with one session at 39 °C and the rest of the sessions at 41 °C (blue bars), a schedule with all the sessions at 41 °C (green bars), a schedule with a session at 43 °C and the rest of the sessions at 39 °C (orange bars) and a schedule with a session at 43 °C and the rest of the sessions at 41 °C (red bars). Finally, in the legend the parameter T stands for the temperature. The SiHa cell line model parameters were considered for this analysis.

3.2. TCP dependencies

3.2.1. Temperature dependence

Simulated patients were generated under the conditions described in Section 2.7.1. to study the dependence of the TCP on the minimum, mean, and maximum temperatures when direct HT cell killing was considered or not (Table 3, Figure 5, Supplementary Table S1, and Figure S4).

Figure 5. Obtained TCP values for 10⁴ simulated patients (each point in the plots corresponds to a patient) after four HT sessions with 30 min time intervals, achieving a temperature randomly selected from a uniform distribution between 39 and 43 °C for each HT session. The minimum (left column), the mean (middle column) and the maximum (right column) temperatures are shown for each patient plotted when direct HT cell killing is considered (lower row) or not (upper row). The SiHa cell line model parameters were considered for this analysis.

Table 3. Correlation coefficients obtained when comparing the TCP with the minimum, mean and maximum temperature. The values presented in this table correspond to when only radiosensitization is considered (outside parentheses) and direct HT cell kill is also included (inside parentheses). The SiHa cell line model parameters were considered for this analysis.

Temperature	ρ/τ
Minimum	0.626/0.451 (0.533/0.375)
Mean	0.998/0.962 (0.972/0.858)
Maximum	0.707/0.520 (0.804/0.610)

A robust correlation was obtained when the TCP was compared with the mean temperature of all HT sessions, whereas a moderate correlation was obtained for the minimum/maximum temperature. This means that the treatment outcome was more closely to the overall treatment (i.e., the combined effect of all HT sessions) than to a single HT session. Furthermore, this correlation usually worsened when direct HT cell killing was included.

3.2.2. Time interval dependence

Simulated patients were generated under the conditions described in Section 2.7.2. to study the dependence of the TCP on the minimum, mean, and maximum time interval (Table 4 and Supplementary Table S2).

Table 4. Correlation coefficients obtained when comparing the TCP with the minimum, mean and maximum time interval and considering diverse temperatures and decay constants.

Variable	Temperature	ρ/τ for $\mu =0.027\hspace{0.17em}{h}^{-1}$	ρ/τ for $\mu =0.5\hspace{0.17em}{h}^{-1}$
Minimum time interval	39 °C	−0.673/−0.491	−0.814/−0.622
	41 °C	−0.675/−0.492	−0.829/−0.639
	43 °C	−0.677/−0.494	−0.843/−0.654
Mean time interval	39 °C	−0.9998/−0.989	−0.963/−0.836
	41 °C	−0.9998/−0.987	−0.954/−0.818
	43 °C	−0.9997/−0.985	−0.946/−0.803
Maximum time interval	39 °C	−0.662/−0.481	−0.519/−0.365
	41 °C	−0.660/−0.479	−0.503/−0.353
	43 °C	−0.658/−0.478	−0.490/−0.343

The correlation coefficients are independent of the direct HT cell killing because the temperatures are fixed.

For time-interval dependence, it can also be concluded that the treatment outcome is more closely related to the mean time interval of all HT sessions (overall treatment) than to the minimum/maximum value (single HT sessions). Therefore, only TCP dependence on the mean time interval is presented for low and high decay constants when direct HT cell killing is considered or not (Figure 6, Table 5, Supplementary Figure S5 and Table S3).

Figure 6. Obtained TCP values for 10⁴ simulated patients (each point in the plots corresponds to a patient) after four HT sessions at 39 (blue points), 41 °C (orange points) and 43 °C (red points) with time intervals randomly selected from a uniform distribution between 10 and 240 min for each HT session. The dependence on the mean time interval for each patient is plotted using a decay constant of $\mu =0.027 h^{-1}$ (left column) or decay constant of $\mu =0.5 h^{-1}$ (right column) and considering (lower row) or not (upper row) direct HT cell killing. The SiHa cell line model parameters were considered for this analysis.

Table 5. TCP variations (ΔTCP = TCP_max – TCP_min) obtained for diverse temperatures and decay constants. The values presented in this table correspond to when only radiosensitization is considered (outside parentheses) and direct HT cell kill is also included (inside parentheses). The SiHa cell line model parameters were considered for this analysis.

Temperature	ΔTCP for $\mu =0.027 h^{-1}$	ΔTCP for $\mu =0.5 h^{-1}$
39 °C	2.0 × 10^−2 (1.8 × 10^−2)	1.8 × 10^−1 (1.6 × 10^−1)
41 °C	2.1 × 10^−2 (1.1 × 10^−2)	2.9 × 10^−1 (1.6 × 10^−1)
43 °C	1.0 × 10^−2 (4.0 × 10^−4)	3.3 × 10^−1 (1.5 × 10^−2)

The obtained results present a high dependence on the decay constant (Table 5 and Supplementary Table S3): TCP variations of the order of 10⁻⁶ − 10⁻² are observed for low decay constants, while a higher decay constant ($\mu =0.5 h^{-1}$) leads to variations one order of magnitude higher. Furthermore, the time interval effect is reduced when direct HT cell killing is considered because this biological effect is independent of the time interval (Equation (4)(4) $$S{F}_{\mathit{\text{HT}}}\left(T,t_{\mathit{\text{HT}}}\right)=\mathit{\text{exp}}\left[-k\left(T\right)\cdot t_{\mathit{\text{HT}}}\right]$$(4) ). Additionally, the time interval effect is usually more pronounced at higher temperatures, independent of the decay constant, which agrees with the results reported by Kok et al. [49]. However, this statement is not true when direct HT cell killing is considered at 43 °C because, at this high temperature, most of the cells are directly killed by HT (TCP > 0.97 when radiosensitization is not considered) and, therefore, under these conditions, the time interval effect can only produce small TCP variations (<3 × 10⁻²).

3.2.3. Random patient results

Simulated patients were generated under the conditions described in Section 2.7.3. to study the dependence of the TCP on the temperature and time interval (Table 6, Figures 7 and 8, Supplementary Table S4 and Figures S6 and S7).

Figure 7. Obtained TCP values for 4 × 10⁴ simulated patients (each point in the plots corresponds to a patient) with different number of HT sessions (3–6). In each HT session, time interval and temperature were randomly selected from uniform distributions within the limits shown in Table 2. The dependence on the mean temperature (left column) and the mean time interval (right column) for each patient are plotted when direct HT cell killing is considered (low row) or not (up row). For these simulations, a decay constant of $\mu =0.027 h^{-1}$ is considered. The SiHa cell line model parameters were considered for this analysis.

Figure 8. Obtained TCP values for 4 × 10⁴ simulated patients (each point in the plots corresponds to a patient) with different number of HT sessions (3–6). In each HT session, time interval and temperature randomly selected from uniform distributions within the limits shown in Table 2. The dependence on the mean temperature (left column) and the mean time interval (right column) for each patient are plotted when direct HT cell killing is considered (low row) or not (up row). For these simulations, a decay constant of $\mu =0.5 h^{-1}$ is considered. The SiHa cell line model parameters were considered for this analysis.

Table 6. Correlation coefficient ranges obtained when comparing the TCP for a given number of HT sessions with the mean temperature and time interval. The values presented in this table correspond to when only radiosensitization is considered (outside parentheses) and when direct HT cell killing is also included (inside parentheses). The SiHa cell line model parameters were considered for this analysis.

Variable	Decay constant	ρ	τ
Mean Temperature	$\mu =0.027 h^{-1}$	[0.9897, 0.990] ([0.972, 0.974])	[0.913, 0.916] ([0.857, 0.866])
	$\mu =0.5 h^{-1}$	[0.435, 0.446] ([0.871, 0.882])	[0.297, 0.306] ([0.686, 0.700])
Mean Time interval	$\mu =0.027 h^{-1}$	[−0.114, −0.098] ([−0.060, −0.044])	[−0.076, −0.065] ([−0.040, −0.029])
	$\mu =0.5 h^{-1}$	[−0.814, −0.788] ([−0.355, −0.341])	[−0.620, −0.593] ([−0.240, −0.231])

As expected, no correlation was observed with the mean time interval when the decay constant was low (Table 6 and Supplementary Table S4). Additionally, this correlation became moderate/strong with a higher decay constant. However, a weak correlation was obtained when direct HT cell killing was included because this biological process does not depend on the time interval.

However, when the decay constant was small, comparing the TCP with the mean temperature resulted in a robust correlation that was slightly worsened if direct HT cell killing was included in the calculations. Furthermore, when the decay constant increased, the correlation with the mean temperature weakened; however, this correlation increased when direct HT cell was considered.

3.2.4. CEM43 dependence

This procedure, as explained in Section 2.7.4, was used to study the dependence of the TCP on the total thermal dose.

Figure 9 and Supplementary Figure S8 show that direct HT cell killing did not drastically change the TCP value when the time interval was short. Also, the TCP rapidly rises at small CEM43 values, showing a clear improvement with a total CEM43 of 10 min (this value can be obtained with 6 HT sessions at 40.5 °C). However, when the time interval increases (i.e., the radiosensitization effect is reduced), the direct HT cell killing becomes more important. A clear example is the case of the 240 min time interval: radiosensitization alone does not improve the outcome but, when direct HT cell killing is included, a clear improvement is observed for high total CEM43 values (i.e., high temperatures). Finally, schedules with more HT sessions had higher mean TCP values for similar thermal doses.

Figure 9. Mean TCP values obtained for diverse total CEM43 values when direct HT cell killing is considered (solid line) or not (dashed line). 10⁶ simulated patients with different number of HT sessions were generated with temperatures randomly selected from a uniform distribution between 37 and 43 °C. In each plot the time interval was fixed for all the simulated patients and a decay constant of $\mu =0.5 h^{-1}$ was considered. The SiHa cell line model parameters were considered for this analysis.

3.3. Inter-patient heterogeneity

First, a coefficient of variation of c.v. = 5% was assumed to study a case with small inter-patient heterogeneity, which might correspond to a case with patients at a similar cancer stage (Figure 10B and Supplementary Figure S9B). In contrast, a c.v. = 50% was assumed to represent the case in which most patients were included (Figure 10C and Supplementary Figure S9C).

Figure 10. Obtained TCP values for 10⁴ simulated patients (each point in the plots corresponds to a patient) after four HT sessions with 30 min of time interval achieving a temperature randomly selected for each HT session (considering direct HT cell killing). A case without inter-patient heterogeneity (left column) and with inter-patient heterogeneity with a c.v. of 5% (middle column) and 50% (right column) are shown. The SiHa cell line model parameters were considered for this analysis.

The correlation obtained without heterogeneity (Figure 10A and Supplementary Figure S9A) worsened when a small inter-patient heterogeneity was considered (Figure 10B and Supplementary Figure S9B) and disappeared when the heterogeneity was further increased (Figure 10C and Supplementary Figure S9C). In the last case (Figure 10C and Supplementary Figure S9C), the TCP values are too spread and no tendencies can be observed because the TCP calculation is quite sensitive to variations of ${\alpha }_{37}$. Hence, if the same RT schedule is performed for all the simulated patients, many will have TCP values of approximately 0 or 1. Thus, the TCP of a single simulated patient cannot be used to study the dependencies/correlations. However, the mean TCP value can be used as an alternative.

Figure 11 shows that the main difference between the two cell lines is that the HeLa cells have a significantly higher baseline mean TCP (RT only) when c.v.=50%. This is because the HeLa radiosensitivity (0.230 Gy⁻¹) is close to the lower limit of the radiosensitivity range imposed in this manuscript (0.2–0.74 Gy⁻¹) and, therefore, most of the selected radiosensitivities have higher values, increasing the mean TCP baseline.

Figure 11. Obtained TCP mean values for 10⁶ simulated patients after four HT sessions with 30 min time intervals, achieving a temperature randomly selected from a uniform distribution between 39 and 43 °C for each HT session (considering direct HT cell killing). A case without inter-patient heterogeneity (left column) and with inter-patient heterogeneity with a c.v. of 5% (middle column) and 50% (right column) are shown. The SiHa (upper row) and HeLa (lower row) cell line model parameters were considered for this analysis.

Nonetheless, a clear temperature dependence was observed for both cell lines. However, the improvement produced by HT decreases with inter-patient heterogeneity, which may result in an underestimation of the capacity of HT to improve RT treatments.

4. Discussion

This study demonstrated the utility of in-silico thermoradiobiological models in elucidating clinical findings, investigating the impact of variations in HT-related treatment conditions on treatment outcomes, and evaluating the efficacy of different HT + RT schedules. Furthermore, although the results obtained in this study pertain to the treatment of cervical cancer, the proposed methodology can be generalized to various tumor locations.

In terms of interpretation of clinical results, our findings demonstrate that the TCP dependence on temperature/time interval exhibit a stronger correlation with the mean value of all HT sessions than with the minimum/maximum value (Sections 3.2.1 and 3.2.2). Additionally, employing the mean temperature as a benchmark proves to be a valuable strategy for investigating the time interval effect (Section 3.2.3). This approach aligns with studies such as the one published by Kroesen et al. [68], in which the authors examined the impact of time interval by comparing the mean time interval with the clinical outcomes. Nonetheless, it would be interesting to explore the relationship between treatment outcome dependence and mean temperature, and compare it with the correlations observed in this study. Therefore, in-silico models offer a novel perspective for comprehending clinical data, enabling not only the investigation of the importance of a specific biological process (e.g., the time interval effect) is important or not, but also to get an idea on how significant it might be.

Furthermore, several authors have underscored the significance of inter-patient heterogeneities in clinical outcomes [69–71], a finding consistent with our results (Section 3.3). Consequently, the contribution of inter-patient heterogeneity should be meticulously considered when analyzing clinical data and comparing theoretical in-silico results with real patient data. Moreover, while intra-patient heterogeneities (i.e., differences within a patient) appear to have a lesser impact [72–75], they should be considered alongside heterogeneities in radiation dose and temperature distributions (e.g., through voxel-based TCP calculation) to achieve greater precision in results. First, RT-related treatment conditions, such as dose-rate, dose per fraction, number of fractions, or impact of the brachytherapy boost, may influence final cell survival and should thus be carefully considered into calculations. However, precise monitoring of temperature distribution within the tumor is challenging, complicating the feasibility of conducting voxel-based TCP calculations accurately. Therefore, caution should be exercised when applying this approach to clinical data.

Regarding the impact of variations in HT treatment conditions, our findings indicate that single HT sessions usually do not exert a significant impact on treatment outcomes. A HT session with poor heating can be compensated by subsequent sessions, while a HT session achieving high temperatures only enhances the treatment outcomes if the remaining sessions are reasonably good. Additionally, we observed a discernible improvement in treatment outcomes with total CEM43 thermal doses of approximately 10 min. Furthermore, for a given thermal dose, the TCP increased with the total number of HT sessions (Section 3.2.4). These findings align with those by Ohguri et al. [76], who demonstrated that CEM43 values higher than 1 min significantly increased the complete response rates in cervical cancer patients. Similarly, Franckena et al. [77] showed that a high median CEM43 level was associated with tumor control and improved survival in a large retrospective data analysis of cervical cancer patients.

Concerning the evaluation of HT + RT schedules, the most effective treatment corresponds to the shortest time interval and highest temperature. However, patient characteristics strongly influence the achieved temperature and may be constrained by the presence of hotspots (i.e., high-temperature regions), which may lead to normal tissue toxicities [78–81] as well as patient discomfort or pain [82, 83]. To address these concerns, additional HT sessions may be incorporated into the schedule to compensate for the lower temperatures, thereby enhancing patient comfort. However, variations in RT doses/schedules for different cancer types directly impact the HT schedule (e.g., a three/six-week RT schedule implies treating twice/once per week with HT to deliver a total of six sessions). Furthermore, thermotolerance, characterized by the temporary heat resistance of cells following a prior HT session, imposes limitations on the maximum number of HT sessions that can be administered to a patient within a given week [84–86]. Consequently, the feasibility of adding new HT sessions to a predefined HT + RT schedule warrants careful consideration.

In addition to the omission of the thermotolerance effect, one of the limitations of our calculations was the incapability to discern chronological differences in the treatment schedule. Furthermore, our calculations solely accounted for the effect of the HT therapeutic window (steady state) and assumed a constant heating duration of one hour. However, actual patient scenarios exhibit diverse temperature profiles, with some patients failing to reach a steady state, and the pre-heating part of the HT session may significantly contribute to HT-induced biological effects. In this regard, the integration of a dynamical model that considers the evolution of biological processes would be advantageous in addressing these aspects. Additionally, HT elicits a myriad of biological effects, including perfusion or immune reactions, which were not considered in our study. These additional factors may mask the results presented here. Thus, conducting an analysis of clinical patient data becomes essential to study the importance of such processes on the outcomes of HT + RT treatments under varying conditions and to validate the in-silico thermoradiobiological models. This validation represents a necessary step for the future development of biological treatment planning systems capable of identifying the ideal HT + RT treatment for each patient.

Finally, the variations in HT treatment conditions assumed in the presented results were modeled using uniform probability distributions. However, actual probability distributions may differ, potentially affecting the final results. Therefore, two Gaussian distributions were defined for temperature (µ = 40.46 °C and σ = 0.68 °C) and time interval (µ = 61 min and σ = 45 min) based on data reported in the reviews of Carrapiço-Seabra et al. [57] and Ademaj et al. [58] from diverse cancer types treated with deep HT. Repeating the calculations with these probability distributions resulted in a slight reduction (<10%) in the TCP ranges (i.e., the difference between the maximum and minimum TCP) due to the less extreme treatment conditions. Thus, while probability distributions appear to have an impact on the results, uniform distributions may serve as a reasonable approximation when actual probability distributions are unknown, as they produce similar results and cover a broader range of scenarios.

In conclusion, this study demonstrates how in-silico thermoradiobiological models can be used to theoretically analyze the impact of varying thermal conditions on the HT + RT treatment outcomes. Furthermore, biophysical models can help comprehend the importance of diverse biological processes, design new clinical trials, and eventually optimize HT + RT treatments.

Disclosure statement

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

Data availability statement

The methods described in this manuscript were implemented in Python (version 3.9.7), and plots were created using R (version 4.1.3). All the codes used for simulating and plotting the presented results are available in the GitHub repository at https://github.com/mingzhaw/TCP_paper.git

Additional information

Funding

This project (Hyperboost; www.Hyperboost-h2020.eu (accessed on 10 November 2021)) has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 955625.