Mathematical modeling of radiotherapy: impact of model selection on estimating minimum radiation dose for tumor control

Introduction Radiation therapy (RT) is one of the most common anticancer therapies. Yet, current radiation oncology practice does not adapt RT dose for individual patients, despite wide interpatient variability in radiosensitivity and accompanying treatment response. We have previously shown that mechanistic mathematical modeling of tumor volume dynamics can simulate volumetric response to RT for individual patients and estimation personalized RT dose for optimal tumor volume reduction. However, understanding the implications of the choice of the underlying RT response model is critical when calculating personalized RT dose. Methods In this study, we evaluate the mathematical implications and biological effects of 2 models of RT response on dose personalization: (1) cytotoxicity to cancer cells that lead to direct tumor volume reduction (DVR) and (2) radiation responses to the tumor microenvironment that lead to tumor carrying capacity reduction (CCR) and subsequent tumor shrinkage. Tumor growth was simulated as logistic growth with pre-treatment dynamics being described in the proliferation saturation index (PSI). The effect of RT was simulated according to each respective model for a standard schedule of fractionated RT with 2 Gy weekday fractions. Parameter sweeps were evaluated for the intrinsic tumor growth rate and the radiosensitivity parameter for both models to observe the qualitative impact of each model parameter. We then calculated the minimum RT dose required for locoregional tumor control (LRC) across all combinations of the full range of radiosensitvity and proliferation saturation values. Results Both models estimate that patients with higher radiosensitivity will require a lower RT dose to achieve LRC. However, the two models make opposite estimates on the impact of PSI on the minimum RT dose for LRC: the DVR model estimates that tumors with higher PSI values will require a higher RT dose to achieve LRC, while the CCR model estimates that higher PSI values will require a lower RT dose to achieve LRC. Discussion Ultimately, these results show the importance of understanding which model best describes tumor growth and treatment response in a particular setting, before using any such model to make estimates for personalized treatment recommendations.


Introduction
More than 50% of all cancer patients will receive radiation therapy (RT) during the course of their cancer treatment, either given with curative intent as a single agent, concurrently with systemic therapies, or (neo-)adjuvant to other therapeutic approaches, or in the palliative setting (1).Even modest improvements in treatment outcomes and quality of life for cancer patients undergoing RT would yield benefits for a large patient cohort.However, current radiation oncology practice does not personalize or adapt RT dose for individual patients, despite variance in individual patient radiosensitivity.Thus, many patients are potentially receiving either too much or too little RT dose.Recent efforts include more strategic integration of basic science approaches into radiobiology and radiation oncology to help better understand the mechanisms of radiation response dynamics and to help predict how to best personalize radiation to individual patients.Genomic signatures (2)(3)(4)(5), imaging metrics (6)(7)(8), and burgeoning machine learning and artificial intelligence approaches are being retrospectively and prospectively evaluated as novel biomarkers for radiation response (9)(10)(11).
Simple mathematical approaches have a long history in radiobiology and radiation oncology.The linear-quadratic (LQ) model that describes the clonogenic survival of a cell population to increasing acute doses of radiation has been extensively used to identify cell-intrinsic radiosensitivities (12)(13)(14).Prominent developments of the LQ model include the concept of biologically effective dose (BED), tumor control probability (TCP), and normal tissue complication probability (NTCP) (15)(16)(17)(18)(19).Many conceptual studies have attempted to explain the biological underpinnings of linear-quadratic response dynamics.However, the non-linear tumor growth and treatment response dynamics require deployment of population dynamics models.
Mathematical oncology may hold the key to mechanistic understanding of the complex adaptive dynamic tumor system and its response to radiotherapy (20)(21)(22)(23)(24), with demonstrated feasibility of translation into prospective clinical trials.Using TCP and NTCP concepts combined with a logistic differential equation that describes the recovery of normal tissues from sublethal radiation-induced damage, Scott et al. pioneered the concept of temporally feathered radiation therapy (TFRT) that prioritizes and de-prioritizes organ-at-risk doses at different times during treatment.TFRT was subsequently shown to lead to increased doses to the radiation target, or reduced cumulative doses to organs at risk (25).Leder et al. combined experimental and differential equation models to identify novel radiation schedules to significantly improve radiation efficacy by taking advantage of the dynamic instability of radioresistance (26), which was recently demonstrated to be feasible and safe to administer to glioblastoma patients (27).Our group has introduced the concept of a patient specific 'carrying capacity' in a logistic growth model, called the proliferation saturation index (PSI) as a putative biomarker for radiosensitivity in head and neck cancer as well as non-small cell lung cancer (28)(29)(30)(31) that is currently being evaluated as trigger for personalized radiation dose fractionation (NCT03656133).
One of the advantages of using mechanistic mathematical models to simulate radiation responses is that if an appropriate model is calibrated, validated, and predictive power demonstrated, then it may be used to simulate potential alternative treatments (32,33).However, it is critical to examine the effects of the underlying models on these alternative treatment recommendations.While two models with different mechanstic mathematical formulations may be trained to fit longitudinal dynamics and predict individual responses equally well, they could have different implications for alternative radiation dose fractionations.
Poleszczuk et al. previously analyzed that clinical predictions are strongly dependent on the specific growth law assumed, and that the applicable growth law should be known to be utilized in clinical practice (29).The objective of this paper is to examine the impact of the underlying tumor volume dynamics models on the estimated optimal RT dose (Figure 1).In a previous study (34), our group used a pre-specified model of tumor volume response to RT to estimate the minimum RT dose required for locoregional control of head and neck tumors.Herein, we compare two different mathematical models of response to RT: (1) cytotoxicity to cancer cells that lead to direct tumor volume reduction (DVR) and (2) radiation damage to the tumor microenvironment that lead to tumor carrying capacity reduction (CCR) and subsequent tumor shrinkage.Both of these models have been shown capable of fitting longitudinal tumor volume data from head and neck cancer patients (31,35).The comparison is focused on evaluating the impacts of both of these models on RT dose personalization.

Tumor growth model
Tumor growth models are plentiful, ranging from incredibly simple exponential growth to highly different growth dynamics as the tumor volume changes, either relative to itself or its (static or dynamic) microenvironment (36).While classical investigations sought a universal tumor growth model (37)(38)(39), the search of tumor growth laws is still very much ongoing (40).The seminal study by Benzekry et al. demonstrated that the Gompertz growth model best captured pre-clinical in vivo breast and lung cancer growth dynamics but fell short of adequate forecasts beyond one subsequent measurement (41).More recently, Kather and his team provided the first such model comparison analysis in clinical data of 1,472 patients undergoing chemotherapy or cancer immunotherapy for solid tumors (42).Again, the Gompertz model provided the best balance between goodness-of-fit and number of parameters, but once more early treatment response was only moderately correlated with final treatment responses.
We have previously shown that logistic growth dynamics provide excellent fits to clinical data of head and neck and nonsmall cell lung cancer during fractionated radiotherapy, and demonstrated predictive power of final tumor volumes with sufficient patient-specific data (30,35).Logistic growth is described by the differential equation: where V is tumor volume (cc), l is the intrinsic tumor growth rate (day -1 ), and K is the carrying capacity of the tumor (cc), which is the maximum size tumor that the local tissue can support (28).In logistic growth, the tumor volume grows initially exponentially but growth monotonically decelerates as the volume approaches the defined carrying capacity, visually indicated by the horizontal asymptote (Figure 2A).In describing the distinct types of growth dynamics that this model can capture, our group has previously defined the proliferation saturation index (PSI) a measure of the effective tumor growth rate in the absence of RT (28).PSI is defined by the expression:  Study Overview and Hypothesis.The primary objective of this study is to assess how the selection of the underlying mathematical model of tumor volume dynamics affects estimates for optimal RT dose.This is premised on the idea that even when the same observations of tumor volume changes (both due to off-treatment growth and regression due to treatment effect) are used as inputs for differing models of response to RT there may be different estimates of the optimal RT dose.

PSI ≡
V 0 K 0 where V 0 and K 0 are the initial tumor volume and tumor carrying capacity before treatment, respectively.K 0 can be prospectively calculated using the following expression: where V 0 is the current tumor volume before radiation (usually obtained at radiation simulation), V Dx is the tumor volume measure at diagnosis, and Dt is the time between the two volume measurements (usually a few weeks). 5PSI is defined between 0 and 1.As PSI approaches 0 tumor growth approaches exponential growth, which indicates a tumor microenvironment capable of sustaining a much larger tumor than what currently exists.In contrast, as PSI approaches 1, tumor growth approaches its carrying capacity, which corresponds to high tumor proliferation saturation in the constraints of the tumor microenvironment limiting further proliferation of the tumor.

Modeling response to radiotherapy
We simulate the effect of RT with 2 different models (1): Direct Tumor Volume Reduction (DVR) and ( 2) Tumor Carrying Capacity Reduction (CCR).

Direct tumor volume reduction model
In the DVR model (Figure 2B), we simulate the effect of an RT fraction as an instantaneous reduction in proliferating tumor volume due to cancer cell death: where V + is the tumor volume after the RT fraction; V -is the tumor volume before the RT fraction; g is the cancer cell death rate; and K is the tumor carrying capacity.The parameter g is derived from the linear-quadratic model (12,43): where d [Gy] is the RT dose per fraction, and a [Gy -1 ] and b [Gy -2 ]  are the LQ radiation sensitivity parameters, respectively.In this study, we set the ratio a b = 10 Gy, as seen in many cancer types that are treated with fractionated RT, including head and neck cancer (44).Of note, here we model the effect of radiation on tumor volume and not individual cells in a clonogenic assay.Therefore, the absolute value of a may not be directly comparable to the preclinical radiobiology literature.

Tumor carrying capacity reduction model
In the CCR model (Figure 2B), we model the effect of RT as an instantaneous reduction in the tumor carrying capacity: where K + is the tumor carrying capacity after each RT fraction; K -is the tumor carrying capacity before the RT fraction; and d is the proportion that the carrying capacity is reduced with each RT fraction, ranging from 0 to 1. Modeling the effect of RT as a reduction in the tumor carrying capacity is motivated by observations of how RT alters components of the tumor microenvironment, such tumor vasculature (45) or the release of tumor-specific antigens and damage-associated molecular patterns (DAMPs) that stimulate antitumor immunity (46), that may reduce the tumor carrying capacity.As of yet, the actual dose dependency of radiation-induced carrying capacity reduction is unknown.Thus, we limit this study to the effect of the total dose given in 2 Gy fractions, without consideration of alternative dose fractionations.

Simulating tumor volume dynamics during RT
All simulations were done using custom scripts developed in Java and subsequent analyses were done in Python.The code for both is available at the following Github repository: https://github.com/akutuva21/SPARK-Project.All simulations of RT were performed using a common schedule for fractionated RT, where d = 2 Gy fractions are delivered daily Monday-Friday with no RT delivered on Saturday and Sunday.All simulations of RT were performed using schedules and doses routinely used in treating head and neck cancer patients.Between-fraction tumor volume changes were simulated with the logistic growth model using a 1-hour time resolution.

Parameter sweep analysis
To understand the impact of the model parameters on tumor volume dynamics, we conducted parameter sweeps of both the radiosensitivity parameters (a for the DVR model and d for the CCR model) and the intrinsic tumor growth rate, l.For the radiosensitivity parameter sweeps, we set l = 0.1 day -1 and PSI = 0.9 for all simulations.These parameters are arbitrarily chosen to investigate qualitative response dynamics.Dynamics for different growth rate and PSI parameters are comparable and intuitively derivable from the below results.For the DVR model, we tested a ∈ (0, 0.20) Gy -1 , with a step size of 0.01 Gy -1 .For the CCR model, we tested d ∈ (0, 0.20), with a step of size of 0.01.
For the intrinsic growth rate (l) sweeps, we tested l ∈ (0, 0.10) day -1 for both models.For the other parameters, we set a PSI = 0.7 for both models for rich model dynamics, and a = 0.1 Gy -1 for the DVR model and d = 0.1 for the CCR model.These sweeps were done using simulations of standard 6 weeks of RT with standard weekday fractionation.Parameter ranges were informed by previous studies fitting these models to longitudinal tumor volume data from head and neck cancer patients that received fractionated RT (31,34,35).However, herein we focus on qualitatively demonstrating response dynamics without emphasis on actual values for a specific cancer type.

Estimating minimum RT dose for locoregional tumor control
In head and neck cancer, mid-treatment volumetric responses correlate with outcome (31,47).Patients with greater than 32.2% tumor volume reduction after 4 weeks of RT were 100% locoregionally controlled (LRC) at a mean follow-up time of 20 months.We have previously used this tumor volume reduction threshold to estimate patient-specific RT doses to achieve locoregional control (LRC) using the CCR model (34).While it is conceivable that a greater tumor volume reduction would not jeopardize tumor control, higher RT doses are c o r r e l a t e d w i t h h i g h e r n o r m a l t i s s u e c o m p l i c a t i o n probability (NTCP).
Here, we will use the 32.2% tumor volume reduction threshold to estimate the minimum cumulative RT dose required to achieve locoregional control (LRC) in both DVR and CCR models.We simulate RT up to 8 weeks (allowing consideration of modest dose escalation) using the same fractionation schedule and dose/ fraction described above and then finding the minimum cumulative dose (D min ) where the tumor volume shrinks below the volume reduction threshold.These simulations were done over the following parameter ranges: PSI ∈ (0.6,1.0), a ∈ (0.06, 0.14) Gy -1 for the DVR model, and d ∈ (0.01, 0.09) for the CCR model.

Parameter sweep analysis
Intuitively, when the radiosensitivity parameters (a for the DVR model, d for the CCR model) increase, the reduction in tumor volume increased for both models of RT response (Figure 3).However, the effect of the intrinsic growth rate, l, on tumor volume reduction were opposite in the two models.In the DVR model, as l decreases the net tumor volume reduction increases (Figure 4A).This is because with lower l values there is less repopulation between RT fractions.In the CCR model, however, higher l values result in higher net tumor volume reduction (Figure 4B).This counterintuitive result comes from the fact that in the CCR model tumor volume only decreases when V > K, which makes dV dt < 0 and results in l becoming the rate at which the tumor volume approaches the current carrying capacity from above.Of interest, this contrasts with the response dynamics during the first week of RT (Figure 3B, inset).During early radiation fractions, the carrying capacity remains greater than the tumor volume, which results in dV dt > 0 and continued tumor growth, albeit slower with each fraction as V approaches K. Thus, initially, higher l values yield higher transient tumor volumes, followed by steeper volume reduction.

Minimum cumulative dose estimation
In the DVR model, higher radiosensitivity (a) leads to lower estimated D min for LRC, while higher PSI values lead to higher estimated D min (Figures 5A, B).Similarly, in the CCR model higher radiosensitivity (d) leads to lower estimated D min for LRC.However, increasing PSI in the CCR model leads to lower estimated D min for LRC (Figures 5E, F).This is, again, due to tumor reduction being achieved only when the carrying capacity drops below the current tumor volume, and consequently dV dt < 0. The closer the tumor volume is to its carrying capacity (i.e., higher PSI), the faster radiation can reduce the carrying capacity below the current value.The implications of these analyses can be more clearly seen by looking at D min for LRC as a function of PSI and the radiosensitivity parameters (Figures 5C, D, G, H 1-4).In the DVR model, higher PSI always yields a higher estimate for D min regardless of the value of a (Figure 5C).The opposite is true in the Effect of radiation sensitivity parameter in DVR and CCR models.

Discussion
In this study we have shown that different mathematical models of RT response yield different clinical implications.Thus, the choice of RT response model is critical when estimating RT doses that best shrink tumor volumes based on intrinsic model parameters.Although the two models studied herein both estimate that patients with higher radiosensitivity will require a lower RT dose to achieve LRC, they make opposite estimates on the impact of pretreatment tumor growth dynamics biomarker, PSI.Fast-growing tumors have lower PSI values and slow-growing tumors have higher PSI values.The DVR model, which assumes that the effect of RT comes from the direct radiation-induced death of tumor cells, estimates that tumors with higher PSI values, i.e. lower pretreatment proliferation, will require a higher RT dose to achieve LRC.However, the CCR model, which assumes that the effect of RT comes from a reduction in the tumor carrying capacity of the local tissue, estimates that higher PSI values will require a lower RT dose to achieve LRC.It is therefore of utmost importance to know which mathematical model best describes the radiobiology that underlies the observable radiation response dynamics.It is encouraging, however, that both radiosensitivity and PSI could be measured or estimated in the clinic: PSI can be calculated by using two temporally separated tumor volume measurements before the start of treatment (28) and radiosensitivity via genomic measures such as the radiosensitivity index (48).Estimates of radiosensitivity may increase in accuracy by using serial measurements of the tumor volume during RT to dynamically update estimates of tumor radiation response (35).
The opposing estimates of the effect of PSI on the required minimum RT dose for LRC stem from the fact that tumor cell death is modulated by PSI in the DVR model as the model assumes that only proliferative cells are killed by RT.This means that as PSI increases, higher and higher doses will be required to achieve the same reduction in tumor volume.On the other hand, in the CCR model tumor volume reduction only occurs once the carrying capacity is less than the tumor volume.This means that tumors with higher PSI values require less RT dose for the carrying capacity to drop below the tumor volume, as the initial values for V and K are already relatively close to each other.This result suggests that by determining which model more accurately describes on-treatment tumor volume dynamics in a particular scenario it may be possible to determine which effect of RT is more dominant.
Herein, we focused on two particular models of radiation response mechanisms in the logistic growth model and PSI framework that have been previously presentedand studied which of the two mechanisms has the predominant effect tumor volume dynamics (28,35).While it is conceivable that both mechanismsdirect cancer cell kill and modulation of the tumor microenvironment via carrying capacity reductioncontribute to the clinically observed tumor responses, the current framework unable to combine both models without significant adjustments to the underlying mathematics.Combining both models (SI Methods) leads to scenarios where V/K > 1 resulting in numerical artifacts of unrealistic spikes and large oscillations in tumor volume (SI Figure 1) for the majority of tested model parameters (SI Figure 2).It may be possible to prevent this issue by simulating DVR and CCR on different timescales, but the required mathematics are beyond the scope of this study and left for future analysis.
Furthermore, we have constrained the models to only simulate changes in tumor volume immediately before and during an RT treatment course.There are documented delayed and cumulative effects of RT that manifest in the months following radiation (49)(50)(51).The discussed models, however, are specifically trained to simulate on treatment tumor response dynamics.For tumor decay following radiation, for example due to activated immune responses or clearance of necrotic debris, more complex models will be required that do not contribute to the implications of the herein discussed results.Additionally, the chosen models require measurable tumor volume.Any dynamics below the limit of detection where stochastic effects may dominate need to be simulated with different mathematical approaches (52,53).
Ultimately, it is critical to understand if a model appropriately describes tumor growth and treatment response for specific cancer subtypes, or individual patients, before using any such model to make estimates for personalized treatment recommendations.It may eventually be possible to select appropriate patient-specific or tumor site-specific models, but this will require further study.One route for studying which models of response to RT are best fit for different contexts will be increased acquisition of both direct and indirect measurements of tumor burden during the course of RT.This may be enabled by emerging techniques such as RT with MRIguided linear accelerators (MR-LinAc) (54, 55) and liquid biopsies to measure circulating tumor DNA (ctDNA) and cell-free DNA (cfDNA) (56)(57)(58).In the absence of sufficient evidence for cancerspecific or patient-specific model selection, however, the more prudent approach may be to consider model ensembles for making prediction or treatment recommendations.Ensemble modeling is commonly utilized in weather forecasting, transport modeling, ecology, or financial forecasting to account for model biases, measurement uncertainties, and forecast uncertainty (59)(60)(61)(62)(63).In the context of modeling tumor response to RT, if using an ensemble of models, one might only recommend a change from standard treatment or dosing when a sufficient number of models in the ensemble agree on the direction or magnitude of the estimated treatment personalization.
Tumor dynamics models.(A) Simulated example of tumor growth modeled as logistic growth.The blue curve shows tumor volume, V, over time and the orange line the tumor carrying capacity, K (B) Simulated examples of response to RT, which is modeled by either direct tumor volume reduction (DVR, left) or tumor carrying capacity reduction (CCR, right).Timing of RT fractions, simulated as a standard fractionated RT course with weekday fractions, is shown by the arrows above each plot.

FIGURE 1
FIGURE 1 ). D min (PSI), D min (a), and D min (d) were fit to exponential functions with the form D min = a • e bx + c , where x is PSI, a, or d depending on the respective context (fitted coefficient values in SI Tables (A) Tumor volume trajectories simulated using the DVR model with values of a ∈ (0,0.2) Gy-1 , where larger a values lead to greater decrease in tumor volume.(B) Tumor volume trajectories simulated using the CCR model with values of d ∈ (0,0.2),where larger d values lead to greater decrease in tumor volume.For both models, l = 0.1 day -1 and PSI = 0.9, and the value of the respective radiation sensitivity parameters are indicated by the color bar.All simulations have an arbitrary initial tumor volume with a fractionated RT regimen, where treatment is applied every weekday for a total of five weeks of treatment.CCR model, where higher PSI always yields a lower estimate for D min , regardless of the value of d (Figure5G).Additionally, for both the DVR and CCR model, higher values of the radiosensitivity parameter yield lower estimates for D min (Figures5D, H).Overall, in the DVR model the highest estimated D min values are found at high PSI values and low a values, while in the CCR model the highest estimated D min values are found at low d values and low PSI values.

5 FIGURE 4
FIGURE 5 Minimum cumulative dose (D min ) for LRC estimates for DVR and CCR models.(A, E) Sample volume trajectories for representative parameter pairs across the parameter range, where the bold symbols indicate the location on the heatmap in (B, F). Green curves are the tumor volume plotted as function of cumulative dose, which increases linearly with treatment time; horizontal dashed line indicates the 32.2% volume reduction cutoff used to calculate D min ; vertical red dashed line indicates the calculated D min with the specific value of D min indicated on the x-axis.Patient-specific parameters for each simulation are indicated on the corresponding plots.For all simulations, l = 0.07 day -1 .(B, F) Heatmaps of D min over the clinically relevant range for the radiosensitivity parameter (a or d) and PSI.All simulations have an arbitrary initial tumor volume with a fractionated RT regimen, where treatment is applied every weekday for a total of five weeks of treatment.Black curves indicate "iso-dose" levels with the number of RT fractions required for the indicated dose.White areas indicate parameter regions where sufficient volume reduction was not achieved in the 8 weeks of simulated RT. (C, G) Plots of the radiosensitivity parameters (a or d) against D min for PSI = 0.7, 0.8, 0.9.Colored dots are data points sampled from the heatmaps in (B, F); the corresponding solid lines are exponential fits to the data (fitted coefficients in SI).(D, H) Plots of PSI against D min for 3 different values of the radiosensitivity parameters.Colored dots are data points sampled from the heatmaps in (B, F); the corresponding solid lines are exponential fits to the data (fitted coefficients in SI Tables1-4).