Objective: To personalize radiation therapy dose fractionation protocols, it will be necessary to first quantitatively describe tumor volume reduction dynamics and subsequently simulate the results of alternative fractionation schemes.
Methods and Materials: The proliferation saturation index (PSI) model of tumor volume dynamics was fit to weekly tumor volume data from computed tomography scans of n = 39 head and neck cancer patients who received 66 to 70 Gy in standard daily fractions or with accelerated fractionation. Using the outputs of this model, we additionally simulated how these patients would respond to hyperfractionation with 1.2 Gy fractions twice a day. We identified PSI values that would improve responses and outcomes to hyperfractionation compared to single daily fractions.
Results: The PSI model fit volumetric tumor response dynamics data with high accuracy (R2 = 0.92) using patient-specific PSI values. Simulations of an alternative fractionation protocol demonstrated that a subset of patients with intermediate PSI values could potentially have improved locoregional control by switching to hyperfractionation.
Conclusions: This is the first demonstration of the PSI model fitting data from head and neck cancer patients, and the results suggest a benefit from alternative fractionation schemes for a selected subset of patients.
Head and neck cancers (HNC) are commonly treated with definitive radiation therapy (RT) with or without systemic therapy, with RT being the single most common oncologic treatment.1,2 Few inroads have been made to synergize biological and quantitative approaches with radiation biology and radiation oncology methodologies to personalize and optimize RT. Standard treatment is 56 to 70 Gy in 1.8 to 2 Gy fractions once daily with or without concurrent chemotherapy,1 derived from maximum tolerable dose concepts and dose-escalation trials premised on a “one-size-fits-all” RT dose. Current radiation oncology clinical practice is that every patient receives the same treatment, planned strictly on pre-radiation American Joint Committee on Cancer (AJCC) TNM stage (Tumor size, lymph Node involvement, Metastasis presence), without regard to individual tumor dynamics that may influence RT outcome. Oropharyngeal cancer still represents significant heterogeneity within the same stage group. Thus, it remains unclear why some patients with similar clinical stage, treated with identical RT dose and dose fractionation, are cured while others are not. There is a strong need for personalized RT for individual patients.3
Adaptive RT delivery has been suggested primarily for anatomical and spatial adjustments during the course of treatment.4 Although alternative fractionation schemes have been tested in clinical trials with modest improvements from altered fractionation, it is possible that these benefits may be offset by increased toxicities across the cohort.5 This sets up the need to be able to identify specific patients who may benefit from alternative fractionation protocols. To rigorously evaluate every possible radiation dose and dose fractionation is experimentally and clinically unfeasible.6 However, the burgeoning field of integrated mathematical oncology may make such analyses possible. The integration of quantitative approaches could provide novel, reliable biomarkers based on mathematics and patient-specific disease dynamics to guide RT treatment personalization.7-9
Mathematical modeling in radiation oncology has a long history, with the linear-quadratic (LQ) model,10-14 biologically effective dose (BED),15-19 tumor control probability,20–22 and normal tissue complication probability (NTCP)23-25 being widespread in research and practice. Recent mathematical modeling and simulation studies of the canonical radiobiological principles have led to the concept of temporally feathered radiation therapy (TFRT). In TFRT, different treatment plans are developed for each weekday of fractionated radiation to spare organs at risk in fractionated radiation schedules.26,27 In close collaboration with experimental biologists, a mathematical model of glioma response to radiation has been calibrated and validated to develop and subsequently experimentally confirm radiation protocols that optimally counteract stem cell de-differentiation dynamics.28
One mainstay of mathematical modeling in radiation oncology is to simulate the volume regression profiles during RT,29-33 and to predict responses to a variety of dose and dose fractionation protocols.33-37 Previous analyses revealed that parameters of complex models may be impossible to determine with limited clinical data,38 suggesting that the number of patient-specific parameters must be kept to a minimum to avoid overfitting and model uncertainty. To this extent, the novel concept of a proliferation saturation index (PSI) has been previously introduced to simulate non-small cell lung cancer patient-specific response to RT with a single parameter.39,40 Mathematical analysis revealed that PSI can robustly describe clinical data for a wide range tumor growth models.41 Here we show that the PSI model can simulate head and neck cancer patient-specific responses to RT with data from two clinical cohorts from Moffitt Cancer Center and MD Anderson Cancer Center. We then use the model to run in silico comparisons of a hyperfractionation protocol to identify which patients may most benefit from hyperfractionated radiation.
A total of 39 head and neck cancer patients were treated with 66 to 70 Gy RT with and without concurrent chemotherapy. Seventeen patients were treated at Moffitt Cancer Center where they received a total of 66 to 70 Gy RT in 2 Gy weekday fractions with concurrent chemotherapy, and the remaining 22 patients were treated at MD Anderson Cancer Center where they received a total of 66 to 70 Gy RT (2 or 2.12 Gy weekday fractions or with accelerated fractionation) with half of the patients receiving concurrent chemotherapy. The patient cohort was comprised of different primary sites (oropharyngeal, laryngeal, nasopharyngeal), HPV status, and clinical stage (stage T1-T4). Patient characteristics for the two cohorts are detailed in Table 1.
Tumor volumes were delineated on weekly cone-beam computed tomography (CBCT) or CT-on-Rails system combining a GE Smart Gantry CT scanner (GE Healthcare) and a 2100EX linear accelerator (Varian) (256 total CT-scan-derived tumor volumes). It should be noted that contouring tumor volumes from CBCT scans is highly dependent on adequate image quality. Insufficient contrast, obstruction of the tumor by other patient anatomy, and other factors may prevent delineation of tumor volumes. We only included tumor volumes from patients who had contourable images. Nearly 50% of the patients considered for the study had inadequate image quality and were thus excluded. Locoregional control censored up to 5 years was abstracted as outcome measure and determined by biopsy confirmation or imaging sufficient to initiate additional treatments.
Tumor growth was modeled as logistic growth as governed by the following differential equation:
where λ is the intrinsic volumetric growth rate [day-1], V is the gross tumor volume [cc], and K is the tumor carrying capacity [cc], which is defined as the maximum tumor volume that the local tissue can support. We have previously proposed characterizing individual patient tumor growth rates with the PSI, as opposed to the patient-specific growth rates.41 PSI is defined as the ratio of the initial tumor volume at the start of RT, V0, to the tumor carrying capacity:
PSI is defined between 0 and 1 and PSI represents the fraction of nonproliferative cells in the tumor (Figure 1).
Response to radiation was modeled using the linear-quadratic (LQ) model:13
where SF(d) is the surviving fraction of cells after receiving a radiation dose, d [Gy]; and α [Gy-1] and β [Gy-2] are radiosensitivity coefficients. This is connected to the change in tumor volume by assuming that only proliferating cells are killed by radiation, in line with the Norton-Simon hypothesis that the rate of tumor regression under therapy is proportional to the tumor growth rate.42-44 The change in tumor volume with each radiation fraction is modeled as follows:
where γ is the volumetric death term per radiation fraction that is coupled to the LQ model as follows:
Custom scripts were written in MATLAB (Mathworks) to simulate volume trajectories given inputs of fractionation schedule, dose per fraction, λ, PSI, α, and α/β. To evaluate the ability of the model to fit the patient data, an optimization script was written using the particle swarm optimization function from the MATLAB Global Optimization Toolbox to find patient-specific PSI values, given a particular (λ,α) pair for the entire cohort, assuming α/β = 10 Gy. The optimal values for λ and α were found by performing a full grid search over λ ∈ (0.02, 0.18) day-1 with a step size of 0.025 and α ∈ (0.04, 0.25) Gy-1 with a step size of 0.01 to find the (λ,α) pair that minimized the mean square error to the patient data across the entire cohort for the first 4 weekly measurements.
The parameter optimization identified that the values of λ = 0.07 day-1 and α = 0.09 Gy-1 best fit the tumor volume dynamics across the entire patient cohort. Although the values of λ and α were the same for all patients, the patient-specific PSI values allowed the model to fit both shallow-volume regression dynamics (high PSI) and rapid-volume reduction responses (low PSI) with high accuracy (Figure 2A). Across the entire cohort, the model fit the weekly tumor volumes with high confidence; R2 = 0.92 (Figure 2B). Notably, while PSI was allowed to range from 0 to 1, only values > 0.5 were obtained from the fitting procedure using the optimized λ and α values, which indicates that the model was not overconstrained or overfitting the data (Figure 2C).
To assess the clinical significance of the model parameters, we compared patient-specific clinical measures with the fitted model parameters. Initial tumor size did not correlate with PSI (Figure 2D), but percent change in tumor volume after 4 weeks of RT showed a strong correlation with PSI (Figure 2E).
To translate tumor volume reduction to long-term patient outcomes, we tested midtreatment tumor volume reduction as a predictor for locoregional control (LRC). We found that the median tumor volume reduction at week 4 (-ΔV = 32.2%) perfectly separates out those patients who had early locoregional recurrence (Figure 3). This is comparable to the recent observation that midtreatment nodal decrease ≥ 43% in oropharyngeal cancer is prognostic for locoregional control.45
Given that the standard protocol with 2 Gy weekday fractions results in a subset of patients recurring locoregionally, we simulated an alternative hyperfractionation scheme of 1.2 Gy fractions 2 times a day (1.2 Gy BID). These simulations were performed using the previously optimized λ, α values across the cohort, and patient-specific PSI values. As expected, these simulations predicted that all the patients would have increased tumor volume reduction due to the higher total dose (Figure 4A). However, only a subset of patients with intermediate PSI ([0.835-0.91], Case II: 12/38 patients, 32%) is predicted to have sufficient marginal increase in tumor volume reduction to cross the previously determined volume-reduction threshold for LRC, indicating a potential to benefit from a switch to hyperfractionation (Figure 4B). On the other hand, patients with low PSI values have highly proliferative and radiation-sensitive disease (PSI < 0.835; Case I: 15/38, 39%), with both 2 Gy QD and 1.2 Gy BID yielding midtreatment tumor volume reductions > 32.2% indicating no additional benefit from hyperfractionation. On the other extreme, patients with high PSI values have less proliferative and more radioresistant tumors (PSI > 0.91; Case III: 12/38, 32%), and in this case neither fractionation protocol provides robust midtreatment tumor volume reduction.
Herein, we presented a simple mathematical model based on proliferation saturation index, PSI, that was able to characterize head and neck cancer patient-specific volume changes to fractionated RT with only 1 patient-specific parameter with high accuracy. Notably, the growth rate is assumed to be constant across the entire cohort, and all interpatient heterogeneity is captured in the PSI values. Of interest, initial tumor size, and thusly T stage, did not correlate with PSI and midtreatment response. The need to only determine one patient-specific parameter may facilitate future patient-specific modeling and predictions. It may be possible to estimate patient-specific PSI values using midtreatment volume reduction, which was shown to correlate with PSI. Furthermore, these modeling results demonstrate that treatment response dynamics may provide valuable insights to understanding the nature of the disease. These results complement previous analysis of pretreatment tumor volume dynamics associated with outcomes in patients with oropharyngeal cancer.46
Additionally, since the tumor volumetric death term was coupled to the LQ model, we were able to simulate the potential response to an alternative fractionation scheme of hyperfractionation with 1.2 Gy BID. Although hyperfractionation is not widely used in clinical practice to treat squamous cell head and neck cancers, large meta-analyses have shown an 8% difference in overall survival at 5 years post-treatment compared with standard fractionation.5,47 However, this benefit may be offset by long-term toxicities from the increased overall dose and logistical difficulties in delivering multiple fractions a day. Mathematical models of tumor volume dynamics may eventually serve as part of a framework to identify patients who would benefit from hyperfractionation over standard fractionation.
It is important to note that the α values identified in this study are for radiosensitivity for the entire tumor volume and are not directly comparable to radiosensitivity parameters derived from clonogenic assays in vitro. With the correlation of locoregional control with a volumetric reduction threshold at 4 weeks of RT, we used these simulations to identify which patients could benefit from hyperfractionation to achieve robust midtreatment tumor volume reductions that correlate with locoregional control. This study demonstrates mathematical-modeling-derived PSI as a radioresponse biomarker in head and neck cancer and supports previous findings in non-small cell lung cancer that patients with intermediate PSI values may benefit from hyperfractionation protocols.
Of interest mathematically, it is difficult to accurately identify a unique (λ,α) parameter pair that fits the entire cohort, as other parameter pairs with the same λ/α ratio would yield similar model fits. Thus, it is important to accurately identify either the volumetric growth rate before the start of therapy or tumor radiosensitivity. This can potentially be accomplished with just two pre-treatment CT scans spaced a few weeks apart, such as at diagnosis and at treatment planning or simulation, similar to a previous study of tumor volumetric growth velocity.46 Radiation sensitivity index, RSI, may be a candidate for radiosensitivity that will be explored in future studies.3 The alternative fractionation simulations are also limited by the canonical radiobiological assumptions of the LQ model. It will be important to verify whether this method of accounting for different dose sizes holds up in similar data with different doses per fraction.
Notably, our model cannot capture transient increases in tumor volume during the first few weeks of RT. These types of small increases in tumor volume, or pseudo-progression, in early weeks of treatment have been seen in the other studies. In a recent study of 44 node-positive oropharyngeal cancers,45 the authors saw changes in nodal volume ranging from a 74.3% volume increase to a 73.6% volume decrease at day 10 of chemoradiation; a 48% volume increase to a 94.9.% volume decrease at day 20; and a -18.1% volume increase to a 95.6% volume decrease at day 35. This can potentially be addressed by building a model with a delayed effect of RT.
We acknowledge that this analysis was performed on a small number of patients (n = 39), but we would like to note that our methodology of analyzing longitudinal tumor volumes (6 to 8 CT scans per patient) increases total data points to n = 256 CT scans, significantly strengthening the analysis. The usability of this model in making patient-specific predictions and clinical recommendations still requires independent parameter calibration and validation in independent datasets with larger and more homogenous cohorts.48 However, the results presented here are important steps in the mathematical modeling of RT response in this cancer type and demonstrate the potential for testing alternative treatment schemes in silico to inform the design of clinical trials for personalizing RT prescriptions.
Zahid MU, Mohamed ASR, Latifi K, Rishi A, Harrison LB, Fuller CD, Moros EG, Caudell JJ, Enderling H. Proliferation Saturation Index to Characterize Response to RT and Evaluate Altered Fractionation in Head and Neck Cancer. Appl Rad Oncol. 2021;18(1):18-25.
Dr. Zahid is a postdoctoral fellow, Department of Integrated Mathematical Oncology, H. Lee Moffitt Cancer Center & Research Institute, Tampa, FL. Dr. Mohamed is an instructor, and Dr. Fuller is an associate professor, both in the Department of Radiation Oncology, The University of Texas MD Anderson Cancer Center, Houston, TX. Dr. Latifi is an assistant member, Dr. Rishi is a resident physician, Dr. Harrison is the chair and a senior member, Dr. Moros is a senior member and chief of medical physics, and Dr. Caudell is an associate member, all in the Department of Radiation Oncology, H. Lee Moffitt Cancer Center & Research Institute, Tampa, FL. Dr. Enderling is an associate member, Department of Integrated Mathematical Oncology, and Department of Radiation Oncology, H. Lee Moffitt Cancer Center & Research Institute, Tampa, FL. Disclosure: Dr. Caudell reports grants, personal feesand other fees from Varian, outside of the submitted work. Dr. Enderling reports grants from NIH/NCI Physical Sciences in Oncology Network (1 U01 CA244100-01), while conducting the study. Dr. Zahid and Dr. Enderling have a patent “Personalized Radiation Therapy” (63/010,327) pending. Dr. Fuller reports grants from NIH, NIBIB, Elekta, NSF, University of Texas Health Science Center at San Antonio, and the American Association of Physicists in Medicine. Dr. Moros reports grants/payments from Varian outside of the submitted work. No other authors have conflicts of interest to disclose, and no part of this article has been previously published elsewhere.