Optimal integration of molecularly targeted therapies, such as tyrosine kinase inhibitors (TKI), with concurrent chemotherapy and radiation (CRT) to improve outcomes in genotype-defined cancers remains a current challenge in clinical settings. Important questions regarding optimal scheduling and length of induction period for neoadjuvant use of targeted agents remain unsolved and vary among clinical trial protocols. Here, we develop and validate a biomathematical framework encompassing drug resistance and radiobiology to simulate patterns of local versus distant recurrences in a non–small cell lung cancer (NSCLC) population with mutated EGFR receiving TKIs and CRT. Our model predicted that targeted induction before CRT, an approach currently being tested in clinical trials, may render adjuvant targeted therapy less effective due to proliferation of drug-resistant cancer cells when using very long induction periods. Furthermore, simulations not only demonstrated the competing effects of drug-resistant cell expansion versus overall tumor regression as a function of induction length, but also directly estimated the probability of observing an improvement in progression-free survival at a given cohort size. We thus demonstrate that such stochastic biological simulations have the potential to quantitatively inform the design of multimodality clinical trials in genotype-defined cancers.

Significance:

A biomathematical framework based on fundamental principles of evolution and radiobiology for in silico clinical trial design allows clinicians to optimize administration of TKIs before chemoradiotherapy in oncogene-driven NSCLC.

Therapies targeting specific oncogenic driver mutations have dramatically improved survival in patients with metastatic cancer across several disease sites. Most significantly, tyrosine kinase inhibitors (TKI) have become a standard-of-care treatment for metastatic non–small cell lung cancer (NSCLC) harboring mutations in the EGFR (1, 2). Unfortunately, disease inevitably recurs as both preexisting drug-resistant cells and acquired de novo drug-resistant cells, mutating from drug-tolerant “persistent” cells during the course of treatment, outcompete the drug-sensitive cells, thus manifesting the evolution of drug-resistant tumors (3–6). In contrast to metastatic EGFR-mutant NSCLC, the role of TKIs in locally advanced, nonmetastatic NSCLC remains unknown. But, the addition of TKIs to standard-of-care concurrent chemotherapy and radiation (CRT) in patients with EGFR-mutant cancers has the potential to improve long-term survival in this cohort. The recent PACIFIC trial has demonstrated the dramatic benefit that integration of new systemic agents can yield in locally advanced NSCLC (7), using adjuvant durvalumab, an anti-PD-L1 checkpoint immunotherapeutic agent, and standard-of-care CRT improved progression-free survival (PFS) and overall survival significantly. However, there was no benefit to immunotherapy over chemotherapy in the subset of EGFR-mutant patients (8). Therefore, the PACIFIC trial exemplifies not only the possible benefits of including successful systemic agents into earlier stages of NSCLC treatments, but also the considerable, but underutilized, potential TKIs may play in the increasing number of successful treatment options for EGFR-mutant locally advanced NSCLC.

However, it is entirely unknown how to optimally administer TKIs with CRT to minimize the risk of acquired drug resistance and improve the efficacy of CRT. Consequently, the design of these multimodality therapy trials is largely empirical with great variability in treatment protocols and a “one-size fits all” approach. In addition, a large part of the treatment design space is not explored because of the infeasibility of running randomized controlled trials for every possible combination of treatments.

One useful method to sample treatment design space is bio-mathematical modeling, which enables the quantitation of abstract, interconnected phenomena making it a powerful tool in the field of translational oncology for both hypothesis testing and generation (9, 10). Specifically, mathematical modeling of tumor evolution has significantly impacted our understanding and interpretation of acquired resistance to targeted cancer therapies (11, 12) and has helped formulate the idea of collateral sensitivities to sequential drug regimens (13). In addition, combing evolutionary modeling of tumors with mechanistic biological models of radiation cell killing and plasma level drug concentrations has enabled novel dosing schedules of radiotherapy in glioblastoma (14), intercalated administration of multiple targeted and chemotherapy agents for melanoma (11), and novel pulsed injections of TKI therapy in EGFR-mutant NSCLC (15). This TKI delivery protocol designed using mathematical modeling was successfully translated into a clinical trial (15, 16).

In this study, we present a generalized biomathematical model to optimize TKI plus CRT multimodal regimens in EGFR-mutant locally advanced NSCLC, with a focus on the relative effect of therapy on local versus occult distant disease sites. In metastatic EGFR-mutant NSCLC populations, TKIs result in significant improvements in overall response and survival compared with chemotherapy, but acquired drug resistance inevitably leads to disease recurrence (1–6, 17). Optimal combination regimens with EGFR TKIs and CRT for patients with locally advanced NSCLC have not yet been studied. A major complicating factor is that failure after CRT can be local or distant, which may be differentially impacted by TKIs. The goal of this study was to inform the design of TKI-CRT multimodal clinical trials through mechanistic, biomathematical modeling, which encompasses fundamental principles of evolutionary targeted drug resistance, radiation biology, and comorbidities. The model aims at quantitative, in silico analyses of varying treatment design parameters to predict effect sizes in heterogeneous patient populations, and also the probability of observing such effects at a given sample size. We achieved this by defining distributions for model parameters to simulate interpatient heterogeneity, calibrating our model using several institutional datasets, and subsequently demonstrate accurate model predictions of recurrence rates in independent validation datasets of recent multicenter clinical trials and meta-analyses. Next, we exhaustively explored the multimodal treatment design space and showed that a TKI induction period of 2–3 months, as used in the recent clinical trials NCT01553942 and NCT01822496, may reduce the effectiveness of adjuvant TKI therapy due to the proliferation of TKI-resistant clones during the induction period. Instead, we propose and provide a quantitative rationale for an individualized induction period tailored to resistance evolution. This mathematical framework can not only inform the design of future multimodal TKI-CRT trials, but can also be extended to other oncogene-driven cancers to design more precise and personalized cancer treatment regimens.

Local and distant NSCLC progression model

Our tumor progression model is an advancement of our previously published biomathematical models of CRT in NSCLC (18) and evolutionary TKI resistance in advanced EGFR-mutant NSCLC (12). An exponential vector system (Eqs. A–N) tracked the number of TKI-resistant, drug-tolerant persistent, and sensitive clonogenic cells (clonogens). While TKI-sensitive clonogens are killed by the drug and TKI-resistant clonogens multiply irrespective of drug exposure, drug-tolerant persistent clonogens neither proliferate or succumb to the drug, but do mutate into drug-resistant clonogens (4, 5). The tumor burden was separated into local versus distant compartments to deconvolve in-field locoregional failures and out-of-field distant failures. Thus, time to local and distant failures (see Supplementary Eqs. S5 and S6) was directly modeled from the local and distant cell compartments, respectively, rather than modeling overall survival as a function of total tumor burden. This was an important distinction, as postprogression treatment in the era of targeted therapy can be highly patient specific with varying degrees of response (19). Comorbidity-related deaths were implemented into our model with a Monte Carlo Russian Roulette formalism (see Supplementary Eqs. S1–S4) and estimated from an analysis of the Surveillance Epidemiology and End Results (SEER) program for a regional lung cancer population (20, 21). While a recent SEER analysis demonstrated that patients with lung cancer receiving TKIs were less likely to be smokers, TKI recipients did not have a statistically significant difference in comorbidity as measured by the Charlson comorbidity index (0 vs. 1+; ref. 22). Progression was defined as the first occurring local failure, distant failure, or comorbidity-related death. A pictorial summary conceptualizing the model and its main endpoints is shown in Fig. 1.

Figure 1.

Graphical illustration of the distant and local tumor progression model. TKI-sensitive, -persistent, and -resistant cells are shaded in blue, purple, and red, respectively. Top, the distant tumor compartment is only susceptible to systemic agents, in this case, chemotherapy (chemo) and TKI therapy. Bottom, the local tumor compartment is additionally affected by radiotherapy. Note that TKI therapy is not only specific to sensitive cells, but also leads to slowed growth of persistent cells, as depicted by the downward purple arrow in the illustration. Chemoradiation, however, is assumed to have equal effects on each cell subpopulation. The initial number of distant cells is assumed to be a constant fraction of the cells in the initial primary tumor.

Figure 1.

Graphical illustration of the distant and local tumor progression model. TKI-sensitive, -persistent, and -resistant cells are shaded in blue, purple, and red, respectively. Top, the distant tumor compartment is only susceptible to systemic agents, in this case, chemotherapy (chemo) and TKI therapy. Bottom, the local tumor compartment is additionally affected by radiotherapy. Note that TKI therapy is not only specific to sensitive cells, but also leads to slowed growth of persistent cells, as depicted by the downward purple arrow in the illustration. Chemoradiation, however, is assumed to have equal effects on each cell subpopulation. The initial number of distant cells is assumed to be a constant fraction of the cells in the initial primary tumor.

Close modal

Model calibration to recurrence dynamics of NSCLC

While the tumor progression model tracked a single tumor volume trajectory, a general patient population with heterogeneity in their presentation, response to treatment, and recurrence patterns was simulated by creating truncated normal and truncated log-normal distributions for the model parameters as was done previously (18). These distributions were then randomly sampled for each run (i.e., simulated patient) of the tumor progression model, yielding a histogram of times to events, which in turn was used to create simulated Kaplan–Meier curves for the distributed population (see Supplementary Data). The model parameter distributions were then optimized such that the model-predicted Kaplan–Meier curves and matched clinically reported Kaplan–Meier curves, as was done by Geng and colleagues (18). In this work, the growth and radiosensitivity distributions were fitted to literature Kaplan–Meier curves of freedom from local and distant failure (FFLF and FFDF) in wild-type (WT) and EGFR-mutant locally advanced NSCLC populations receiving definitive concurrent CRT (19, 23, 24). The TKI model parameters were derived from a model-based analysis of patients with advanced EGFR-mutant NSCLC (see Supplementary Data). A full summary of the optimized model parameters along with the source of the calibration data is shown in Table 1.

Table 1.

Parameter overview: overview of key model parameters along with data sources used for model calibration and testing.

Key model parameters—estimated from literature and previous studies
ParameterDescriptionValueReference
|{{\rm{N}}_{{\rm{cells}}}}( 0 )$| Initial number of local clonogens [9.1 × 106–6.7 × 1011] [cells] Geng and colleagues (18) 
|K$| Gompertzian carrying capacity 8.2 × 1012 [cells] Geng and colleagues (18) 
|\eqalign{ {{\rm{V}}_{\rm{p}}}( 0 ) = {{\rm{N}}_{{\rm{persistent}}}}( 0 )/{{\rm{N}}_{{\rm{cells}}}}( 0 ){\rm{\ }}\cr {{\rm{V}}_r}( 0 ) = {{\rm{N}}_{{\rm{persistent}}}}( 0 )/{{\rm{N}}_{{\rm{cells}}}}( 0 ){\rm{\ }}$| Initial fraction of TKI-persistent and TKI-resistant clonogens [0.05–0.5] [-]
[10−4–10−1] [-] 
Grassberger and colleagues (12) 
|{\rm{\mu }}$| Resistant mutation probability 10−7 [-] Grassberger and colleagues (12) 
|{{\rm{\beta }}_{{\rm{TKI}}}}$| TKI cell kill parameter trunc-norm {μ = 2, σ = 7} [ml/μg]; (>1) Grassberger and colleagues (12) 
|{\rm{\beta }}$|c Chemotherapy cell kill parameter trunc-norm {μ = 0.028, σ = 6.8 × 10−4} [m2/mg]; (>0) Geng and colleagues (18) 
|{t_{1/2}}$| Chemotherapy plasma half life 24[hr] Geng and colleagues (18) 
|\kappa $| TKI plasma decay factor 0.0465 [hr−1Grassberger and colleagues (12), Foo and colleagues (30) 
|{{\rm{r}}_{{\rm{death}}}}$| Comorbidity death rate for regional lung cancer 0.55 [%/mo] SEER (20, 21) 
Key model parameters—estimated from literature and previous studies
ParameterDescriptionValueReference
|{{\rm{N}}_{{\rm{cells}}}}( 0 )$| Initial number of local clonogens [9.1 × 106–6.7 × 1011] [cells] Geng and colleagues (18) 
|K$| Gompertzian carrying capacity 8.2 × 1012 [cells] Geng and colleagues (18) 
|\eqalign{ {{\rm{V}}_{\rm{p}}}( 0 ) = {{\rm{N}}_{{\rm{persistent}}}}( 0 )/{{\rm{N}}_{{\rm{cells}}}}( 0 ){\rm{\ }}\cr {{\rm{V}}_r}( 0 ) = {{\rm{N}}_{{\rm{persistent}}}}( 0 )/{{\rm{N}}_{{\rm{cells}}}}( 0 ){\rm{\ }}$| Initial fraction of TKI-persistent and TKI-resistant clonogens [0.05–0.5] [-]
[10−4–10−1] [-] 
Grassberger and colleagues (12) 
|{\rm{\mu }}$| Resistant mutation probability 10−7 [-] Grassberger and colleagues (12) 
|{{\rm{\beta }}_{{\rm{TKI}}}}$| TKI cell kill parameter trunc-norm {μ = 2, σ = 7} [ml/μg]; (>1) Grassberger and colleagues (12) 
|{\rm{\beta }}$|c Chemotherapy cell kill parameter trunc-norm {μ = 0.028, σ = 6.8 × 10−4} [m2/mg]; (>0) Geng and colleagues (18) 
|{t_{1/2}}$| Chemotherapy plasma half life 24[hr] Geng and colleagues (18) 
|\kappa $| TKI plasma decay factor 0.0465 [hr−1Grassberger and colleagues (12), Foo and colleagues (30) 
|{{\rm{r}}_{{\rm{death}}}}$| Comorbidity death rate for regional lung cancer 0.55 [%/mo] SEER (20, 21) 
Key model parameters—derived in this study
ParameterDescriptionValueReference
|{{\rm{f}}_{{\rm{met}}}}$| Initial number of metastatic clonogens as a fraction of |{{\rm{N}}_{{\rm{cells}}}}( 0 )$| 2 × 10−6 (95% CI: 7.5 × 10−7 - 4 × 10−6) [-] This work 
|{\rm{\rho }}$| NSCLC growth rate trunc-norm {μ = 7 × 10−5, σ = 0.0055
(95% CI: 0.005-0.007)} [day−1]; (>0) 
This work 
|{\rm{\alpha }}$| WT NSCLC radiosensitivity trunc-norm {μ = 0.10 (95% CI: 0.04-0.18), σ = 0.17 (95% CI: 0.11-0.2)} [Gy−1]; (>0) This work 
|{{\rm{\alpha }}_{EGFR}}$| EGFR-mutant NSCLC radiosensitivity trunc-norm {μ = 0.16 (95% CI: 0.10-0.26), σ = 0.32 (95% CI: 0.20-0.44)} [Gy−1]; (>0) This work 
Key model parameters—derived in this study
ParameterDescriptionValueReference
|{{\rm{f}}_{{\rm{met}}}}$| Initial number of metastatic clonogens as a fraction of |{{\rm{N}}_{{\rm{cells}}}}( 0 )$| 2 × 10−6 (95% CI: 7.5 × 10−7 - 4 × 10−6) [-] This work 
|{\rm{\rho }}$| NSCLC growth rate trunc-norm {μ = 7 × 10−5, σ = 0.0055
(95% CI: 0.005-0.007)} [day−1]; (>0) 
This work 
|{\rm{\alpha }}$| WT NSCLC radiosensitivity trunc-norm {μ = 0.10 (95% CI: 0.04-0.18), σ = 0.17 (95% CI: 0.11-0.2)} [Gy−1]; (>0) This work 
|{{\rm{\alpha }}_{EGFR}}$| EGFR-mutant NSCLC radiosensitivity trunc-norm {μ = 0.16 (95% CI: 0.10-0.26), σ = 0.32 (95% CI: 0.20-0.44)} [Gy−1]; (>0) This work 
Clinical data used for calibrating and testing model this study
DatasetDescriptionEndpointReference
Calibration (N = 118, United States) Institutional LF and DF rates in EGFR-mutant/WT LA NSCLC FFLF & FFDF [mo] Mak and colleagues (24) 
Calibration (N = 95, Japan) Institutional LF and DF rates in EGFR-mutant/WT LA NSCLC FFLF & FFDF [mo] Lim and colleagues (23) 
Calibration (N = 185, South Korea) Institutional LF and DF rates in EGFR-mutant/WT LA NSCLC FFLF & FFDF [mo] Yagishita and colleagues (19) 
Testing (N = 86, France, Italy, Spain) Phase III clinical trial of TKIs in metastatic NSCLC (NCT00446225) PFS [mo] EURTAC Trial, Rosell and colleagues (1) 
Testing (N = 598, Worldwide) Phase III clinical trial of concurrent chemoradiation in LA NSCLC (NCT00686959) PFS [mo] PROCLAIM Trial, Senan and colleagues (35) 
Testing (N = 6 clinical trials) Meta-analysis of six clinical trials of sequential versus concurrent chemoradiotherapy in LA NSCLC FFLF & FFDF [mo] Auperin and colleagues (36) 
Clinical data used for calibrating and testing model this study
DatasetDescriptionEndpointReference
Calibration (N = 118, United States) Institutional LF and DF rates in EGFR-mutant/WT LA NSCLC FFLF & FFDF [mo] Mak and colleagues (24) 
Calibration (N = 95, Japan) Institutional LF and DF rates in EGFR-mutant/WT LA NSCLC FFLF & FFDF [mo] Lim and colleagues (23) 
Calibration (N = 185, South Korea) Institutional LF and DF rates in EGFR-mutant/WT LA NSCLC FFLF & FFDF [mo] Yagishita and colleagues (19) 
Testing (N = 86, France, Italy, Spain) Phase III clinical trial of TKIs in metastatic NSCLC (NCT00446225) PFS [mo] EURTAC Trial, Rosell and colleagues (1) 
Testing (N = 598, Worldwide) Phase III clinical trial of concurrent chemoradiation in LA NSCLC (NCT00686959) PFS [mo] PROCLAIM Trial, Senan and colleagues (35) 
Testing (N = 6 clinical trials) Meta-analysis of six clinical trials of sequential versus concurrent chemoradiotherapy in LA NSCLC FFLF & FFDF [mo] Auperin and colleagues (36) 

Note: Units in brackets, distributions in curly brackets, and 95% CIs in parentheses for parameters fitted in this work.

Abbreviations: DF, distant failure; hr, hour; LA, locally advanced; LF, local failure; mo, month.

While the 2-year FFLF for WT NSCLC has been shown to be relatively poor, in the range of 54%–63%, EGFR-mutant tumors have a noticeably better response with a 2-year FFLF in the range of 70%–87% (19, 23, 24). This differential response in FFLF was observed after CRT alone and the results were not confounded by TKI administration. Furthermore, in vitro studies have demonstrated an enhanced radiosensitivity of EGFR-mutant NSCLC compared with WT (25), and so in our model, unique radiosensitivity distributions were defined for WT and EGFR-mutant populations separately and optimized against FFLF. Distant metastases have been observed to be the most common form of recurrence with a 2-year FFDF in the range of 31%–43% with no statistically significant relationship to EGFR status (19, 23, 24). Therefore, a common metastatic fraction and growth rate distribution were defined for both populations and optimized against FFDF, as shown in Fig. 2. For each calibration step, two parameters were optimized simultaneously using brute force, and while each of the two parameters exhibited a correlative relationship, a global solution was determined (Supplementary Fig. S1A–S1F). The bootstrapped confidence intervals (CI) encompassed these correlative regions, and model-predicted failure rates matched the calibration data over these intervals (Supplementary Fig. S1G–S1L).

Figure 2.

Model calibration. A–C, Simulated freedom from distant (A) and local (WT, B; EGFR mutant, C) failure Kaplan–Meier curves using the calibrated model parameters. D–F, Corresponding model-predicted versus literature-reported failure rates at 1-, 2-, 3-, and 4-year timepoints. The solid line represents the model-predicted versus the weighted mean of literature-reported failures, with corresponding linear regression and summary statistics. The black dashed line represents unity. DF, distant failure; LF, local failure.

Figure 2.

Model calibration. A–C, Simulated freedom from distant (A) and local (WT, B; EGFR mutant, C) failure Kaplan–Meier curves using the calibrated model parameters. D–F, Corresponding model-predicted versus literature-reported failure rates at 1-, 2-, 3-, and 4-year timepoints. The solid line represents the model-predicted versus the weighted mean of literature-reported failures, with corresponding linear regression and summary statistics. The black dashed line represents unity. DF, distant failure; LF, local failure.

Close modal

The model-predicted FFDF curve for the optimal growth and metastatic fraction parameters is shown in Fig. 2A and exhibited a strong correlation to the weighted mean of the literature-reported values (r2 > 0.99; P = 0.002; Fig. 2D). Similarly, the model-predicted FFLF curves with the optimized WT and EGFR-mutant radiation sensitivity distributions matched the literature-reported values, as shown in Fig. 2B and C, respectively. For both populations, the model-predicted FFLFs were similarly strongly correlated with reported values (r2 > 0.97; P < 0.02; Fig. 2E and F). In addition, the model-predicted PFS was calculated with the optimized parameters for WT and EGFR-mutant populations, which yielded only a minor increase in PFS for the more radiosensitive EGFR-mutant population (HR, 1.06; 95% CI, 1.02–1.09; Supplementary Fig. S2). This was in accordance with literature-reported observations of a lack of statistically significant difference in time to first recurrence between the WT and EGFR-mutant populations (19, 24).

Mathematical implementation of local versus distant tumor progression model

We described the evolution of the number of TKI-resistant (⁠|{N_R}$|⁠), -persistent (⁠|{N_P}$|⁠), or -sensitive (⁠|{N_S}$|⁠) clonogens in the distant (⁠|{N_D}$|⁠) or local (⁠|{N_L}$|⁠) compartments (Eqs. E–J) with exponential factors accounting for the differential terms of Gompertzian cell growth, Norton–Simon cell kill by TKIs, log cell kill by chemotherapy, and linear-quadratic radiation cell kill (with the assumption of |\alpha /\beta $| = 10). While older reports suggested |\alpha /\beta $| could be higher than 10 for lung cancer (26), our results were insensitive to the exact value of |\alpha /\beta $| as different fractionation schema were not considered. In addition, the efficacy of chemoradiation was assumed to be independent of the proximity to TKI administration, that is, the model parameters associated with chemoradiation cell kill were considered to be constant throughout the treatment. The initial local cell number was based on tumor volume distributions for each stage of NSCLC and an assumed cell density of 5.8 × 108 cells/cm3 estimated from previous model-based work of NSCLC tumor growth (18, 27). The initial distant cell number was assumed to be a scalar fraction (⁠|{f_{met}}$|⁠) of the initial local cell number (Eq. B), which was fitted during the modeled calibration procedure. Time of distant failure was defined as the distant compartment reaching a PET-detectable volume of 1 cm3 (28), that is, when |{{\rm{N}}_{{{\rm{D}}_{{\rm{total}}}}}}( {\rm{t}} ) > 1$|cc (Supplementary Eq. S6). Time of local failure was defined as growth of the local compartment past its size at the start of treatment, that is, when |{{\rm{N}}_{{{\rm{L}}_{{\rm{total}}}}}}( {\rm{t}} ) > {{\rm{N}}_{{{\rm{L}}_{{\rm{total}}}}}}[ 0 ]{\rm{\ }}$|(Supplementary Eq. S5). By heuristically simulating a patient population, yielding histograms of time to failures, time-varying rates of FFLFs and FFDFs were modeled. Throughout treatment and regrowth, there was also a transitional term to account for acquired mutations from TKI-persistent to -resistant cells with the mutation rate |\mu $| assumed to be 10−7 based upon in vitro studies of TKI resistance (29). The model also assumed that the growth of the persistent compartment was slowed during TKI administration as seen clinically (Eq. M; refs. 5, 12). The first-order pharmacokinetic model of chemotherapy (⁠|{C_C}$|⁠, Eq. L) and TKI plasma concentrations (⁠|{C_{TKI}}$|⁠, Eq. N) used in our model was fully described in previous publications (12, 18, 30). The TKI dosing regimen of a single daily dose was implemented for both induction and maintenance, as is done clinically (NCT01553942). The state equations were implemented over 1 day time periods with 0.1 day resolution (⁠|\Delta t$|⁠), and with the differential effect of each term calculated independently. Discrete treatment events, such as radiation delivery or a bolus injection of chemotherapy, were assumed to occur instantaneously at a single timepoint (Eqs. K and L). At a discrete time step, if a cell compartment went below one cell, it was assumed to be controlled and set to 0. A total of five iterations of 1,024 patients (n = 5,120) were simulated for the calibration of model parameters, while a total of 1,024 patients having 12 unique combinations of initial resistant and persistent fractions (n = 12,288) were simulated for the model validation and multimodal treatment outcome predictions. The effect of simulated population size on variability of outcomes was examined in the simulated in silico induction trial section of the Results.

Software implementation

The simulations were implemented in Python v3.7 using both the SciPy v1.3.1 (31) and NumPy v1.17.2 (32) libraries. The scipy.odeint routine was used to implement TKI therapy and growth terms, while CRT was implemented exponentially. The seaborn v0.9.0 python library (33) was used to present data, and the lifelines v0.24 python library (34) was used for the survival analysis (see Supplementary Data). The source code for the simulations and resulting data are available at github.com/bomcclatchy/Modeling-Multimodal-TKI-ChemoRadiation.

formula

Retrospective model validation against clinical trial outcomes

With all the model parameters calibrated and fixed, the tumor progression model was then validated against independent datasets from multi-institutional phase III trials of either TKI alone or CRT alone (see Supplementary Data). First, the model-predicted PFS in an advanced stage IV population was compared against the TKI arm of the EURTAC trial, which compared chemotherapy versus TKI in advanced EGFR-mutant NSCLC (Fig. 3A; ref. 1). A recent meta-analysis of first-generation TKIs in EGFR-mutant advanced NSCLC found the EURTAC trial to have the most similar effect to the median of the six analyzed trials (HR EURTAC, 0.42; 95% CI, 0.27–0.64 and HR median, 0.37; 95% CI, 0.27–0.52), making it an appropriate benchmarking dataset (17). The model-predicted PFS was similar to the trial-reported PFS (r2 = 0.92; P < 1e-5; Fig. 3B), validating the rate of progression during TKI administration as predicted by the model. A waterfall plot of the maximum change in tumor volume from baseline is shown for a simulated patient population (n = 256) in Fig. 3C, with each bar color coded by degree of initial TKI resistance. From an analogous waterfall plot reported in the EURTAC trial, we see that the trial and the model-simulated populations have similar tumor response dynamics with a median tumor volume decrease of 50%–75% and <10% of patients exhibiting progressive disease, demonstrating the validity of the derived TKI model parameters. Furthermore, the capacity of the model to stochastically embody a heterogeneous population is shown in Fig. 3C, as tumor shrinkage was not simply a function of preexisting resistance, but was modulated by other factors such as growth rate.

Figure 3.

Model validation. A, Model-predicted and trial reported PFS Kaplan–Meier curves for stage IV, EGFR-mutant NSCLC populations receiving TKIs until progression. The 12 thin Kaplan–Meier curves represent each initial starting persistent and resistant fraction combination (see Supplementary Data), while the full line is the aggregate response. B, Plot of model-predicted versus trial reported PFS over 3 month intervals corresponding to A with linear regression summary statistics and the black dashed line representing unity. C, Waterfall plot of best model-predicted volume change from baseline for stage IV, EGFR-mutant NSCLC populations receiving TKIs until progression for a random subset of 256 modeled patients, color coded by the 12 initial TKI persistence/resistance conditions. D, Model-predicted and trial reported PFS Kaplan–Meier curves for stage III NSCLC populations receiving definitive CRT, with a corresponding plot of model-predicted versus trial reported PFS over 3 month intervals and summary statistics in E.

Figure 3.

Model validation. A, Model-predicted and trial reported PFS Kaplan–Meier curves for stage IV, EGFR-mutant NSCLC populations receiving TKIs until progression. The 12 thin Kaplan–Meier curves represent each initial starting persistent and resistant fraction combination (see Supplementary Data), while the full line is the aggregate response. B, Plot of model-predicted versus trial reported PFS over 3 month intervals corresponding to A with linear regression summary statistics and the black dashed line representing unity. C, Waterfall plot of best model-predicted volume change from baseline for stage IV, EGFR-mutant NSCLC populations receiving TKIs until progression for a random subset of 256 modeled patients, color coded by the 12 initial TKI persistence/resistance conditions. D, Model-predicted and trial reported PFS Kaplan–Meier curves for stage III NSCLC populations receiving definitive CRT, with a corresponding plot of model-predicted versus trial reported PFS over 3 month intervals and summary statistics in E.

Close modal

Next, model-predicted PFS for a WT locally advanced NSCLC population receiving concurrent CRT was compared against the results of the PROCLAIM trial (35), representing the pre-PACIFIC standard of care for CRT alone (Fig. 3D). Model-predicted PFS for concurrent CRT correlated very strongly to the aggregated reported PFS of the trial (r2 > 0.99; P < 1e-21; Fig. 3E), validating the absolute rate of recurrences during concurrent CRT as predicted by the model. In addition, the model-predicted local and distant failure dynamics for sequential versus concurrent CRT were in accordance with the results of a meta-analysis of six randomized trials comparing chemotherapy scheduling in locally advanced NSCLC (36). The model accurately predicted a significant benefit in local failure rates for concurrent CRT with a modeled HR of 0.79 (95% CI, 0.77–0.82), while the meta-analysis reported a HR of 0.77 (95% CI, 0.62–0.95; P = 0.01; Supplementary Fig. S3A). Furthermore, the model predicted no difference in distant failure rates for sequential versus concurrent CRT (HR, 1.00; 95% CI, 0.97–1.03), which was consistent with the findings of the meta-analysis that there was no statistically significant effect of chemotherapy scheduling on distant failure rates (HR, 1.04; 95% CI, 0.86–1.25; P = 0.69; Supplementary Fig. S3B). Together, these results provide strong evidence for the validity of the model-predicted local and distant recurrence patterns during various treatment schedules as compared with current NSCLC multicenter clinical trials.

Estimation of improved outcomes for TKI induction and maintenance

With the same simulated patient population and treatment parameters used during model validation (histograms shown in Supplementary Fig. S4), the expected recurrence dynamics of combining CRT and TKI therapy were explored. Two main treatment designs were simulated for locally advanced EGFR-mutant NSCLC: TKI induction with daily administration up to 16 weeks, followed by definitive concurrent CRT, followed with or without adjuvant TKI maintenance. These two treatment schemes were chosen to approximate the format of ongoing combined TKI + CRT trial protocols (see Supplementary Data). TKI therapy concurrent with CRT was not modeled as initial clinical experience suggests the potential for increased toxicity and also nonsynergistic efficacy (37, 38). The hypothesis explaining these results was that TKIs cause G1-phase cell-cycle arrest stunting cell replication antagonizing both chemotherapy and radiotherapy (39, 40), despite both agents having cytotoxic effects regardless of cell division.

The predicted 2-, 3-, and 5-year FFLF, FFDF, and PFS as a function of TKI induction length with and without adjuvant TKI maintenance are plotted in Supplementary Fig. S5A–S5C, respectively, while the endpoints over the entire design space are tabulated in Supplementary Fig. S6. For TKI induction without maintenance, the greatest predicted benefit for FFLF, FFDF, and PFS over CRT alone occurred with 2 weeks of induction (Δ5-year FFLF, 1.5%; Δ5-year FFDF, 11.1%; and Δ5-year PFS, 6.5%), with a decreasing benefit as the induction length increased, resulting in a similar predicted outcome to CRT alone with 16 weeks of induction. For each endpoint, there was a predicted additive benefit for including adjuvant TKI maintenance to induction regardless of the timepoint or duration of the induction period (Supplementary Fig. S5). But, the greatest benefit of including TKI maintenance was seen at shortest induction periods (Δ5-year FFLF, 5.8%; Δ5-year FFDF, 23.7%; and Δ5-year PFS, 15.8% with 0 weeks induction and Δ5-year FFLF, 3.5%; Δ5-year FFDF, 7.7%; and Δ5-year PFS, 6.7% with 2 weeks induction compared with CRT alone), with outcomes monotonically worsening as a function of induction length for each endpoint. The local failure rate was the least sensitive to the length of induction, given the high radiosensitivity and local control rates of EGFR-mutant NSCLC. The most dramatic effects were seen in the distant failure rates, which have an enhanced sensitivity to the evolutionary dynamics during induction as the subsequent chemotherapy was the only modeled therapeutic able to target the occult TKI-resistant subpopulation in the distant compartment.

Longer induction is predicted to induce TKI resistance

The maximum benefit observed with 2 weeks of TKI induction when maintenance was not administered was due to the fact that this timepoint presents a balance in benefit for both responders and nonresponders, if a patient does not respond at all there is not much additional growth that early, while the responders will have at least derived some benefit through volume shrinkage. Furthermore, the benefit of delaying progression by means of extending the induction length (≥4 weeks) was gradually outweighed by the proliferation of the resistant population and resulting tumor volume increase prior to CRT. But surprisingly, when adjuvant maintenance was administered, the shortest induction periods were predicted to have the lowest local and distant failure rates. This suggests that the relative benefits of shrinking the tumor volume before CRT with TKIs were outweighed by the cost of acquired TKI resistance caused by targeted evolution during induction, which dramatically reduced the efficacy of TKI maintenance. This tradeoff in competing advantages of tumor size and TKI sensitivity is exemplified in the Kaplan–Meier analysis and simulated tumor volume trajectories displayed in Fig. 4.

Figure 4.

Model predictions for variable induction periods. The local and distant tumor volume trajectory of the median simulated patient with an initial persistent fraction of 0.1 and initial resistant fraction of 0.01, stratified by TKI response cell subtypes, are shown in A and B, respectively. D–F, Model-predicted FFLF (D), FFDF (E), and PFS (F) Kaplan–Meier curves for various induction lengths. The simulated treatment regimen is TKI induction, chemoradiotherapy, and adjuvant TKI maintenance. HRs corresponding to the Kaplan–Meier curves in D–F are shown in C. DF, distant failure; LF, local failure; wks, weeks.

Figure 4.

Model predictions for variable induction periods. The local and distant tumor volume trajectory of the median simulated patient with an initial persistent fraction of 0.1 and initial resistant fraction of 0.01, stratified by TKI response cell subtypes, are shown in A and B, respectively. D–F, Model-predicted FFLF (D), FFDF (E), and PFS (F) Kaplan–Meier curves for various induction lengths. The simulated treatment regimen is TKI induction, chemoradiotherapy, and adjuvant TKI maintenance. HRs corresponding to the Kaplan–Meier curves in D–F are shown in C. DF, distant failure; LF, local failure; wks, weeks.

Close modal

The median tumor volume trajectories stratified by TKI-sensitive, -persistent, and -resistant cell subtypes receiving 2 and 12 weeks of TKI induction with adjuvant maintenance are plotted for the local compartment in Fig. 4A and the distant compartment in Fig. 4B. While both treatment schedules resulted in local control (Fig. 4A), a 2-week induction resulted in control of the initially microscopic resistant subpopulation of the distant compartment with chemotherapy during CRT and a slow regrowth of the persistent subpopulation during maintenance (Fig. 4B). But, over the course of the longer 12 weeks of induction, the resistant subpopulation of the distant compartment out competed and outgrew the slowly proliferating persistent subpopulation, resulting in a rapid, resistant distant recurrence after approximately 1 year of adjuvant TKI therapy (Fig. 4B), controlling for the same efficacy of chemotherapy during CRT. In addition, the resistant growth during induction was accelerated in the distant compartment compared with the local compartment because of the resource advantage at lower cell numbers inherently modeled with Gompertzian growth (12).

Kaplan–Meier curves for the simulated FFLF, FFDF, and PFS for various TKI induction and maintenance schedules show the long-term predicted benefit over CRT alone in Fig. 4DF, respectively. While TKI induction and maintenance were predicted to have modest FFLF benefit over CRT alone, there was little stratification between different induction lengths (Fig. 4D). But for FFDF (Fig. 4E) and PFS (Fig. 4F), there was a much stronger effect with greater stratification between induction lengths. However, stratification between induction lengths was not observed until approximately 1 year after the start of treatment, because of the time needed for distant cells to proliferate to an observable threshold after CRT. In this analysis, a schedule of no induction with CRT and adjuvant TKI maintenance was not considered despite having the best predicted outcomes, as some amount of induction with a measurable response in tumor volume was needed to warrant daily TKI administration after CRT until progression.

Calculated HRs for 2, 8, and 12 weeks of TKI induction with adjuvant maintenance quantified the relative effect size of each multimodal treatment schedule compared with CRT alone, independent of the simulated sample size. As shown in Fig. 4C, FFDF had the most dramatic effect size, having an HR of 1.89 (95% CI, 1.83–1.96) for 2-week induction dropping to 1.38 (95% CI, 1.34–1.43) for 12-week induction. The expected statistical significance of these proposed treatment schemes at clinically relevant sample sizes was stochastically investigated, as described in the next section.

Statistical significance and power in a model-based trial design

Evaluating the statistical significance between the predicted failure rates of different treatment arms was a nontrivial task as an arbitrarily large number of patients can be simulated, which could result in statistical significance even when HRs were very close to one. Therefore, we determined the probability of reaching a certain level of statistical significance as a function of the number of trial patients. We utilized the stochastic nature of the model by randomly sampling multiple iterations of simulated patients to yield an estimate of the false negative rate at lower sample sizes and the corresponding statistical power (see Supplementary Data). By doing so, our model-based analysis not only yielded an expected magnitude of effect between two multimodal treatments, but also quantified the probability of observing the effect in a population of a given size.

A simulated two-armed clinical trial of 2 versus 12 weeks of TKI induction with CRT and adjuvant TKI maintenance was modeled with FFDF as the endpoint, which has both the strongest effect (Fig. 4F) and is particularly relevant for future targeted therapy trials given the relatively high rates of local control in NSCLC with CRT for EGFR-mutant patients. An induction length of 2 weeks was chosen for arm 1 as it was predicted to improve outcomes and also would allow to screen for initial TKI response. For arm 2, 12 weeks of induction was chosen to mimic an induction length currently being investigated (NCT 01822496). A depiction of the simulated evolution of TKI resistance illustrating the hypothesized mechanism of benefit to a shorter induction period is shown in Fig. 5A. Modeled FFDF Kaplan–Meier curves displayed a constant magnitude of effect but wider CIs with decreased sample size (Fig. 5B). A heatmap of the log-rank P values testing the significance between the two arms over 1,000 iterations revealed the variability in detecting the effect with fewer patients in Fig. 5C. This is further quantified in Fig. 5D, where the distributions of median FFDF in each arm have a consistent mean, but increasing variance with decreasing sample size. Note that the simulation ran for 5 years and so the peak at 60 months can be attributed to iterations where the median FFDF was not yet reached. The distribution of log-rank P values between the arms for the 1,000 iterations is shown in Fig. 5E along with the median P values. The fraction of iterations reaching statistical significance of 0.05 was estimated to be the statistical power, and is plotted as function of sample size in Fig. 5F. Thus, this analysis has projected that a clinical trial would need 256 patients per arm for a power of 79%.

Figure 5.

Simulated in silico induction trial. A, Illustration of the differential evolution of the TKI-resistant and -sensitive populations for a given tumor burden between a short 2-week and long 12-week (wk) induction (ind) length. When CRT is done after a significant TKI induction period, the tumor shrinks with the targeted drug killing the TKI-sensitive cells (blue), but with more TKI-resistant (red) cells at the time of CRT, increasing the chance of a late TKI-resistant recurrence if CRT is not curative. B, Simulated FFDF Kaplan–Meier curves for 2 versus 12 weeks of induction lengths with increasing number of simulated patients. Each curve corresponds to the iteration with the median log-rank P value. C, A heatmap of log rank P values testing statistical difference between the 2- versus 12-week FFDF Kaplan–Meier curves for 1,000 iterations of the simulation at each sample size. D, Histograms of the median FFDF for the 1,000 iterations of the 2- and 12-week induction simulations at each sample size. E, Histogram of the log rank P value between the 2- versus 12-week induction simulations at each sample size. F, Estimated statistical power as a function of sample size. Here, statistical power was estimated as the fraction of iterations resulting in a P < 0.05. main, maintenance; mo, months.

Figure 5.

Simulated in silico induction trial. A, Illustration of the differential evolution of the TKI-resistant and -sensitive populations for a given tumor burden between a short 2-week and long 12-week (wk) induction (ind) length. When CRT is done after a significant TKI induction period, the tumor shrinks with the targeted drug killing the TKI-sensitive cells (blue), but with more TKI-resistant (red) cells at the time of CRT, increasing the chance of a late TKI-resistant recurrence if CRT is not curative. B, Simulated FFDF Kaplan–Meier curves for 2 versus 12 weeks of induction lengths with increasing number of simulated patients. Each curve corresponds to the iteration with the median log-rank P value. C, A heatmap of log rank P values testing statistical difference between the 2- versus 12-week FFDF Kaplan–Meier curves for 1,000 iterations of the simulation at each sample size. D, Histograms of the median FFDF for the 1,000 iterations of the 2- and 12-week induction simulations at each sample size. E, Histogram of the log rank P value between the 2- versus 12-week induction simulations at each sample size. F, Estimated statistical power as a function of sample size. Here, statistical power was estimated as the fraction of iterations resulting in a P < 0.05. main, maintenance; mo, months.

Close modal

In addition, this benefit to a shorter induction period was seen regardless of the initial TKI sensitivity. When simulated patients were stratified by initial resistant or persistent fractions (see Supplementary Data), and the trial was rerun for each combination of |{\rm{Vr}}( 0 )$| and |{\rm{Vp}}( 0 )$|⁠, the 2-week arm always benefitted as quantified by the HR (Supplementary Fig. S7). While this benefit appeared to be short-term for patients with higher levels of preexisting persistence and lower levels of preexisting resistance, a durable benefit was observed for higher levels of preexisting resistance and lower levels of preexisting persistence.

A corresponding analysis using PFS as the endpoint is displayed in Supplementary Fig. 8A–S8E, from which it was estimated that 512 patients per arm are needed to reach a power more than 80%. While these analyses assumed a general population, the modeled population can be readily stratified into clinically distinguishable groups. For instance, modeled populations above and below the median initial tumor size (7.2 cm) will exhibit diverging distributions of predicted failure rates, which will consequently decrease and increase the required number of patients in each subpopulation needed for an adequately powered clinical trial.

We have developed a mathematical framework to forecast the relative effectiveness of user-defined regimens of combined targeted chemoradiation therapy and have applied this methodology to optimize the multimodal administration of TKIs in locally advanced NSCLC. This framework integrates fundamental principles of tumor growth, radiation biology, and acquired drug resistance, with model parameters calibrated to institutional outcomes and predicted recurrence rates independently validated against clinical trial outcomes. While our previous modeling work has exclusively investigated targeted (12) or cytotoxic (18) therapy, the current framework not only integrates both treatment modalities, but also deconvolves local versus distant failures, creates population-based model parameters specific for EGFR-mutant patients, and includes stochastic noncancerous failures. We predicted local and distant recurrence rates for various lengths of TKI induction before CRT and estimated the expected number of patients needed to observe a clinical benefit. In contrast to the current design of multimodal clinical trials, which prescribe 2–3 month long induction periods, we discovered an inverse relationship between progression and length of induction. In this sense no induction period at all would be optimal. However, some TKI induction is indicated because tumor response to the TKI needs to be verified before proposing an adjuvant TKI regimen. Thus, we propose 2 weeks because a volumetric reduction at this point can be observed when measured by consistent volumetric segmentation.

After 2 weeks of TKI exposure, there may not be maximal macroscopic gross tumor volume change; however, at these early timepoints, the occult TKI-resistant and -persistent populations will still be microscopic and with greater potential for control by cytotoxic therapy. Even though we made the assumption that preexisting resistant cells exist at therapy initiation, our result does not entirely depend on this assumption. Even assuming no preexisting cells would favor shorter induction periods, as this reduces the probability of cells acquiring mutations that confer resistance. The sooner CRT is administered, the sooner the pool of possible persister cells that could acquire resistance is diminished, increasing the efficacy of the TKI maintenance regimen. Thus, the benefit to shorter induction periods holds true regardless of the absolute effect of CRT, as long as the CRT has the ability to control nonsensitive TKI cells. This should also hold true for second- and third-generation TKIs as these drugs still result in acquired drug resistance through a combination of preexisting and acquired mutations (3), despite overcoming the EGFR T790M mutation, which has indeed resulted in superior PFS compared with first-generation TKIs (2). Finally, the predicted benefit of a shorter induction period was robust to any assumed level of preexisting TKI sensitivity in our model, but a greater benefit was seen in tumors with higher levels of preexisting TKI resistance.

Our result that shorter induction times lead to better outcomes contradicts expectations in current multimodal trial protocols for locally advanced NSCLC, but reflects the increasing risk of TKI resistance as a function of TKI exposure. Induction before CRT has several benefits, most notably decreasing tumor size, which has been shown to improve overall survival when treated with CRT (41) and may also lead to surgical candidacy as outlined in the ASCENT trial protocol (NCT 01553942). But, upfront TKI exposure in advanced staged patients has resulted in a 50%–75% average volume change (1), which is similar to the expected surviving fraction after a single 2 Gy radiation fraction in NSCLC (42). In addition, upfront TKI exposure will extend the duration of treatment compared with CRT alone, which will ostensibly delay disease progression. But, as patients systemically become resistant over the course of months, adjuvant TKI maintenance therapy may be rendered ineffective. Conversely, with upfront chemoradiotherapy and the ability to control TKI-resistant cells, the adjuvant TKI maintenance has the potential to dramatically slow tumor regrowth, and in fact, a retrospective analysis of the treatment of brain metastases originating from EGFR-mutant NSCLC demonstrated that upfront radiotherapy with adjuvant TKIs nearly doubled overall survival compared with upfront TKIs until progression and subsequent radiotherapy (34.1 vs. 19.4 months; P = 0.01; ref. 43). But in the locally advanced setting, some initial TKI exposure is needed to demonstrate tumor response before maintenance therapy can be considered. Furthermore, there is evidence to suggest that TKI exposure may prime the cells for radiotherapy by inducing senescence (44).

Genetic heterogeneity present throughout a patients' tumor predicates the clinically observable acquired treatment resistance (45). There are two main models explaining the mechanism by which intratumoral heterogeneity and consequently treatment resistance arises: (i) branched clonogenic evolution where subpopulations arise from an accumulation of mutations and (ii) cancer stem cells (CSC) that are predominantly dormant and treatment resistant, but uniquely give rise to the diversity of differentiated cancer cells (9, 46). In addition, models have been proposed unifying CSCs and clonogenic evolution, suggesting that CSCs themselves can acquire favorable mutations leading to competing subpopulations (47). In this work, a form of clonogenic evolution was assumed with CSCs not explicitly modeled, as the focus of this analysis was acquired TKI resistance, which has been shown to occur through specific genetic mutations (5, 6). As such, our model assumes that TKI-resistant cells have the same chemo- and radiosensitivity as TKI-sensitive cells, which is in-line with experimental evidence (25).

While our model does not explicitly account for the constricting effect of TKIs on vasculature, resource deprivation is implicitly modeled through the slowed growth of the TKI-persistent subpopulations resulting in apparent stunted tumor growth. However, modeling-based in vitro and in vivo studies have shown that alterations to microenvironment can modulate the growth of resistant and persistent subclonogens independently, which in turn can lead to frequency-dependent tumor evolution (48, 49). Future work could incorporate the non-cell autonomous effects of vasculature modulation into the model, by both EGFR-TKIs and also antiangiogenic VEGF inhibitors, which have shown great synergy when administered in combination (50). While in vitro studies have suggested that TKI exposure can modulate DNA repair pathways affecting chemoradiation sensitivity (51, 52), clinical studies have shown chemotherapy intercalated with TKIs to be effective (53), and recent reports have demonstrated effectiveness of radiotherapy after TKIs in oligometastatic disease (54). Consequently, the effects of sequential administration of TKIs and CRT were conservatively assumed to be independent. However, TKIs and CRT were not simulated concurrently as uncertainty regarding the interplay and potential toxicity of concurrent TKIs + CRT remains (37, 38), but could be investigated in future studies as robust clinical data become available. Furthermore, while we parameterized our model using the chemotherapy regimen in RTOG 9410, clinical experience has demonstrated similar outcomes with a variety of chemotherapy regimens (55).

On the basis of our framework, future modeling studies might aim at stratifying simulated populations (e.g., high or low radiosensitivity |\alpha $| and tumor growth |\rho $|⁠) and determine which patient populations experience the most benefit from shorter versus longer induction periods, thus further personalizing therapy. Potential biomarkers for NSCLC radiosensitivity demonstrated to be associated with survival outcomes are rad51 expression (56) and the polygenic radiosensitivity index (57, 58). In addition, ki67 expression and fluorothymidine uptake imaged by PET have both been shown to be related to tumor cell proliferation and could be potential biomarkers of tumor growth rate (59). In addition, future modeling work could focus on predicting organ-specific distant failure rates in locally advanced NSCLC, as brain metastases are common in EGFR-mutant tumors and may be amenable to further radiation (43).

Our work shows the possible impact that mathematical modeling can have on clinical trials exploring the integration of new biological agents into current treatment approaches, a topic that has recently gathered increased attention (60). Mathematical models based on patient data have been used in the past to optimize radiotherapy (14) or to understand underlying disease dynamics (61), among a range of other applications (9, 62, 63). A recent example is adaptive drug therapy for metastatic castration-resistant prostate cancer (64, 65), which is currently being tested in a clinical trial (NCT02415621). The framework presented here could be translated into other disease sites or targetable mutations, however, this would require recalibration of model parameters, which involved several studies as noted in Table 1. To recalibrate the model, the following literature-obtained data are needed for the specific disease site: unrestricted tumor growth with time, failure rates after chemo/radiotherapy (both separately and together), failures rates after targeted therapy, and rates of noncancerous failures (from SEER; refs. 20, 21). With these data, along with literature estimates of preexisting resistance and mutation rates, one could independently tune each model parameters as was done in a previous study (18).

In conclusion, our mathematical framework provides an evolutionary argument against longer induction periods, as they could trigger the process of acquired drug resistance and may limit the efficacy of adjuvant therapy, possibly resulting in worse outcomes. According to our model, a shorter induction period of 1–2 weeks had a greater chance of controlling TKI-resistant cells with CRT, resulting in longer predicted PFS. Finally, the probability of observing a statistically significant increase in PFS due to a shorter induction period was stochastically derived as function of trial size, using randomly sampled heterogeneous patient populations. These model predictions are hypothesis generating and could have impact on clinical trial design. While this study has focused on optimizing TKI administration in combination with CRT for EGFR-mutant NSCLC, the generalized framework outlined in this article can be applied to oncogene-driven multimodal therapy designs in other cancers.

H. Willers reports grants from American Lung Association during the conduct of the study. A.N. Hata reports grants from Pfizer, Relay Therapeutics, Amgen, Eli Lilly, Roche/Genentech, Blueprint Medicines, and Novartis outside the submitted work. Z. Piotrowska reports personal fees from Blueprint Medicines, C4 Therapeutics, AstraZeneca, Spectrum, Takeda, Novartis, Immunogen, AbbVie, Guardant Health, Genentech, Eli Lilly, Incyte, and Medtronic and grants from Novartis, Takeda, Spectrum, AstraZeneca, Tesaro, and Cullinan Oncology outside the submitted work. L.V. Sequist reports grants and personal fees from AstraZeneca, Genentech, and Merrimack, grants from Novartis, Boehringer Ingelheim, Loxo Oncology, and Blueprint Medicines, and personal fees from Janssen outside the submitted work, as well as has a patent for treatment of EGFR-mutant lung cancer pending. C. Grassberger reports other compensation from Nanolive (scientific advisory board) outside the submitted work. No potential conflicts of interest were disclosed by the other authors.

D.M. McClatchy III: Data curation, software, formal analysis, validation, investigation, visualization, methodology, writing-original draft, writing-review and editing. H. Willers: Conceptualization, resources, supervision, funding acquisition, project administration. A.N. Hata: Conceptualization, resources, funding acquisition, methodology, project administration. Z. Piotrowska: Conceptualization, resources, methodology, project administration. L.V. Sequist: Conceptualization, supervision, funding acquisition, project administration. H. Paganetti: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, methodology, project administration, writing-review and editing. C. Grassberger: Conceptualization, resources, software, formal analysis, supervision, investigation, methodology, project administration, writing-review and editing.

This work was supported by grants from the NCI (CA21239 to H. Paganetti, CA197389 to A.N. Hata, and CA137008 to L.V. Sequist) and the American Lung Association (LCD-400286 to H. Willers).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Rosell
R
,
Carcereny
E
,
Gervais
R
,
Vergnenegre
A
,
Massuti
B
,
Felip
E
, et al
Erlotinib versus standard chemotherapy as first-line treatment for European patients with advanced EGFR mutation-positive non-small-cell lung cancer (EURTAC): a multicentre, open-label, randomised phase 3 trial
.
Lancet Oncol
2012
;
13
:
239
46
.
2.
Soria
JC
,
Ohe
Y
,
Vansteenkiste
J
,
Reungwetwattana
T
,
Chewaskulyong
B
,
Lee
KH
, et al
Osimertinib in untreated EGFR-mutated advanced non-small-cell lung cancer
.
N Engl J Med
2018
;
378
:
113
25
.
3.
Ortiz-Cuaran
S
,
Scheffler
M
,
Plenker
D
,
Dahmen
L
,
Scheel
AH
,
Fernandez-Cuesta
L
, et al
Heterogeneous mechanisms of primary and acquired resistance to third-generation EGFR inhibitors
.
Clin Cancer Res
2016
;
22
:
4837
47
.
4.
Ramirez
M
,
Rajaram
S
,
Steininger
RJ
,
Osipchuk
D
,
Roth
MA
,
Morinishi
LS
, et al
Diverse drug-resistance mechanisms can emerge from drug-tolerant cancer persister cells
.
Nat Commun
2016
;
7
:
10690
.
5.
Hata
AN
,
Niederst
MJ
,
Archibald
HL
,
Gomez-Caraballo
M
,
Siddiqui
FM
,
Mulvey
HE
, et al
Tumor cells can follow distinct evolutionary paths to become resistant to epidermal growth factor receptor inhibition
.
Nat Med
2016
;
22
:
262
9
.
6.
Sequist
LV
,
Waltman
BA
,
Dias-Santagata
D
,
Digumarthy
S
,
Turke
AB
,
Fidias
P
, et al
Genotypic and histological evolution of lung cancers acquiring resistance to EGFR inhibitors
.
Sci Transl Med
2011
;
3
:
75ra26
.
7.
Antonia
SJ
,
Villegas
A
,
Daniel
D
,
Vicente
D
,
Murakami
S
,
Hui
R
, et al
Durvalumab after chemoradiotherapy in Stage III non-small-cell lung cancer
.
N Engl J Med
2017
;
377
:
1919
29
.
8.
Lee
CK
,
Man
J
,
Lord
S
,
Links
M
,
Gebski
V
,
Mok
T
, et al
Checkpoint inhibitors in metastatic EGFR-mutated non-small cell lung cancer-a meta-analysis
.
J Thorac Oncol
2017
;
12
:
403
7
.
9.
Altrock
PM
,
Liu
LL
,
Michor
F
. 
The mathematics of cancer: integrating quantitative models
.
Nat Rev Cancer
2015
;
15
:
730
45
.
10.
Grassberger
C
,
Scott
JG
,
Paganetti
H
. 
Biomathematical optimization of radiation therapy in the era of targeted agents
.
Int J Radiat Oncol Biol Phys
2017
;
97
:
13
7
.
11.
Bozic
I
,
Reiter
JG
,
Allen
B
,
Antal
T
,
Chatterjee
K
,
Shah
P
, et al
Evolutionary dynamics of cancer in response to targeted combination therapy
.
Elife
2013
;
2
:
e00747
.
12.
Grassberger
C
,
McClatchy
DM
,
Geng
C
,
Kamran
SC
,
Fintelmann
F
,
Maruvka
YE
, et al
Patient-specific tumor growth trajectories determine persistent and resistant cancer cell populations during treatment with targeted therapies
.
Cancer Res
2019
;
79
:
3776
88
.
13.
Yoon
N
,
Vander Velde
R
,
Marusyk
A
,
Scott
JG
. 
Optimal therapy scheduling based on a pair of collaterally sensitive drugs
.
Bull Math Biol
2018
;
80
:
1776
809
.
14.
Leder
K
,
Pitter
K
,
LaPlant
Q
,
Hambardzumyan
D
,
Ross
BD
,
Chan
TA
, et al
Mathematical modeling of PDGF-driven glioblastoma reveals optimized radiation dosing schedules
.
Cell
2014
;
156
:
603
16
.
15.
Chmielecki
J
,
Foo
J
,
Oxnard
GR
,
Hutchinson
K
,
Ohashi
K
,
Somwar
R
, et al
Optimization of dosing for EGFR-mutant non-small cell lung cancer with evolutionary cancer modeling
.
Sci Transl Med
2011
;
3
:
90ra59
.
16.
Yu
HA
,
Sima
C
,
Feldman
D
,
Liu
LL
,
Vaitheesvaran
B
,
Cross
J
, et al
Phase 1 study of twice weekly pulse dose and daily low-dose erlotinib as initial treatment for patients with EGFR-mutant lung cancers
.
Ann Oncol
2017
;
28
:
278
84
.
17.
Gao
G
,
Ren
S
,
Li
A
,
Xu
J
,
Xu
Q
,
Su
C
, et al
Epidermal growth factor receptor-tyrosine kinase inhibitor therapy is effective as first-line treatment of advanced non-small-cell lung cancer with mutated EGFR: a meta-analysis from six phase III randomized controlled trials
.
Int J Cancer
2012
;
131
:
E822
9
.
18.
Geng
C
,
Paganetti
H
,
Grassberger
C
. 
Prediction of treatment response for combined chemo- and radiation therapy for non-small cell lung cancer patients using a bio-mathematical model
.
Sci Rep
2017
;
7
:
13542
.
19.
Yagishita
S
,
Horinouchi
H
,
Katsui Taniyama
T
,
Nakamichi
S
,
Kitazono
S
,
Mizugaki
H
, et al
Epidermal growth factor receptor mutation is associated with longer local control after definitive chemoradiotherapy in patients with stage III nonsquamous non-small-cell lung cancer
.
Int J Radiat Oncol Biol Phys
2015
;
91
:
140
8
.
20.
Edwards
BK
,
Noone
AM
,
Mariotto
AB
,
Simard
EP
,
Boscoe
FP
,
Henley
SJ
, et al
Annual report to the nation on the status of cancer, 1975–2010, featuring prevalence of comorbidity and impact on survival among persons with lung, colorectal, breast, or prostate cancer
.
Cancer
2014
;
120
:
1290
314
.
21.
Cho
H
,
Mariotto
AB
,
Mann
BS
,
Klabunde
CN
,
Feuer
EJ
. 
Assessing non-cancer-related health status of US cancer patients: other-cause survival and comorbidity prevalence
.
Am J Epidemiol
2013
;
178
:
339
49
.
22.
Enewold
L
,
Thomas
A
. 
Real-world patterns of EGFR testing and treatment with erlotinib for non-small cell lung cancer in the United States
.
PLoS One
2016
;
11
:
e0156728
.
23.
Lim
YJ
,
Chang
JH
,
Kim
HJ
,
Keam
B
,
Kim
TM
,
Kim
DW
, et al
Superior treatment response and in-field tumor control in epidermal growth factor receptor-mutant genotype of stage III nonsquamous non-small cell lung cancer undergoing definitive concurrent chemoradiotherapy
.
Clin Lung Cancer
2017
;
18
:
e169
78
.
24.
Mak
RH
,
Doran
E
,
Muzikansky
A
,
Kang
J
,
Neal
JW
,
Baldini
EH
, et al
Outcomes after combined modality therapy for EGFR-mutant and wild-type locally advanced NSCLC
.
Oncologist
2011
;
16
:
886
95
.
25.
Das
AK
,
Sato
M
,
Story
MD
,
Peyton
M
,
Graves
R
,
Redpath
S
, et al
Non-small-cell lung cancers with kinase domain mutations in the epidermal growth factor receptor are sensitive to ionizing radiation
.
Cancer Res
2006
;
66
:
9601
8
.
26.
Thames
HD
,
Bentzen
SM
,
Turesson
I
,
Overgaard
M
,
Van den Bogaert
W
. 
Time-dose factors in radiotherapy: a review of the human data
.
Radiother Oncol
1990
;
19
:
219
35
.
27.
Switzer
P
,
Gerstl
B
,
Greenspoon
J
. 
Karyometry in the estimation of nuclear population in pulmonary carcinomas
.
J Natl Cancer Inst
1974
;
52
:
1699
704
.
28.
Cuaron
J
,
Dunphy
M
,
Rimner
A
. 
Role of FDG-PET scans in staging, response assessment, and follow-up care for non-small cell lung cancer
.
Front Oncol
2012
;
2
:
208
.
29.
Williams
MJ
,
Werner
B
,
Barnes
CP
,
Graham
TA
,
Sottoriva
A
. 
Identification of neutral tumor evolution across cancer types
.
Nat Genet
2016
;
48
:
238
44
.
30.
Foo
J
,
Chmielecki
J
,
Pao
W
,
Michor
F
. 
Effects of pharmacokinetic processes and varied dosing schedules on the dynamics of acquired resistance to erlotinib in EGFR-mutant lung cancer
.
J Thorac Oncol
2012
;
7
:
1583
93
.
31.
Virtanen
P
,
Gommers
R
,
Oliphant
TE
,
Haberland
M
,
Reddy
T
,
Cournapeau
D
, et al
SciPy 1.0: fundamental algorithms for scientific computing in Python
.
Nat Methods
2020
;
17
:
261
72
.
32.
Oliphant
TE
.
A guide to NumPy
.
USA:
Trelgol Publishing
; 
2006
.
33.
Waskom
M
,
Botvinnik
O
,
O'Kane
D
,
Hobson
P
,
Ostblom
J
,
Lukauskas
S
, et al
mwaskom/seaborn: v0.9.0 (July 2018). v0.9.0 ed. seaborn.pydata.org. Zenodo
.
DOI:
.
34.
Davidson-Pilon
C
.
CamDavidsonPilon/lifelines: v0.24. v0.24 ed. Zenodo
.
DOI:
.
35.
Senan
S
,
Brade
A
,
Wang
LH
,
Vansteenkiste
J
,
Dakhil
S
,
Biesma
B
, et al
PROCLAIM: randomized Phase III trial of pemetrexed-cisplatin or etoposide-cisplatin plus thoracic radiation therapy followed by consolidation chemotherapy in locally advanced nonsquamous non-small-cell lung cancer
.
J Clin Oncol
2016
;
34
:
953
62
.
36.
Auperin
A
,
Le Pechoux
C
,
Rolland
E
,
Curran
WJ
,
Furuse
K
,
Fournel
P
, et al
Meta-analysis of concomitant versus sequential radiochemotherapy in locally advanced non-small-cell lung cancer
.
J Clin Oncol
2010
;
28
:
2181
90
.
37.
Komaki
R
,
Allen
PK
,
Wei
X
,
Blumenschein
GR
,
Tang
X
,
Lee
JJ
, et al
Adding erlotinib to chemoradiation improves overall survival but not progression-free survival in Stage III non-small cell lung cancer
.
Int J Radiat Oncol Biol Phys
2015
;
92
:
317
24
.
38.
Ready
N
,
Janne
PA
,
Bogart
J
,
Dipetrillo
T
,
Garst
J
,
Graziano
S
, et al
Chemoradiotherapy and gefitinib in stage III non-small cell lung cancer with epidermal growth factor receptor and KRAS mutation analysis: cancer and leukemia group B (CALEB) 30106, a CALGB-stratified phase II trial
.
J Thorac Oncol
2010
;
5
:
1382
90
.
39.
Solit
DB
,
She
Y
,
Lobo
J
,
Kris
MG
,
Scher
HI
,
Rosen
N
, et al
Pulsatile administration of the epidermal growth factor receptor inhibitor gefitinib is significantly more effective than continuous dosing for sensitizing tumors to paclitaxel
.
Clin Cancer Res
2005
;
11
:
1983
9
.
40.
Gandara
DR
,
Gumerlock
PH
. 
Epidermal growth factor receptor tyrosine kinase inhibitors plus chemotherapy: case closed or is the jury still out?
J Clin Oncol
2005
;
23
:
5856
8
.
41.
Morgensztern
D
,
Waqar
S
,
Subramanian
J
,
Gao
F
,
Trinkaus
K
,
Govindan
R
. 
Prognostic significance of tumor size in patients with stage III non-small-cell lung cancer: a surveillance, epidemiology, and end results (SEER) survey from 1998 to 2003
.
J Thorac Oncol
2012
;
7
:
1479
84
.
42.
Carmichael
J
,
Degraff
WG
,
Gamson
J
,
Russo
D
,
Gazdar
AF
,
Levitt
ML
, et al
Radiation sensitivity of human lung cancer cell lines
.
Eur J Cancer Clin Oncol
1989
;
25
:
527
34
.
43.
Magnuson
WJ
,
Yeung
JT
,
Guillod
PD
,
Gettinger
SN
,
Yu
JB
,
Chiang
VL
. 
Impact of deferring radiation therapy in patients with epidermal growth factor receptor-mutant non-small cell lung cancer who develop brain metastases
.
Int J Radiat Oncol Biol Phys
2016
;
95
:
673
9
.
44.
Wang
M
,
Morsbach
F
,
Sander
D
,
Gheorghiu
L
,
Nanda
A
,
Benes
C
, et al
EGF receptor inhibition radiosensitizes NSCLC cells by inducing senescence in cells sustaining DNA double-strand breaks
.
Cancer Res
2011
;
71
:
6261
9
.
45.
Gerlinger
M
,
Rowan
AJ
,
Horswell
S
,
Math
M
,
Larkin
J
,
Endesfelder
D
, et al
Intratumor heterogeneity and branched evolution revealed by multiregion sequencing
.
N Engl J Med
2012
;
366
:
883
92
.
46.
Clevers
H
. 
The cancer stem cell: premises, promises and challenges
.
Nat Med
2011
;
17
:
313
9
.
47.
Kreso
A
,
Dick
JE
. 
Evolution of the cancer stem cell model
.
Cell Stem Cell
2014
;
14
:
275
91
.
48.
Marusyk
A
,
Tabassum
DP
,
Altrock
PM
,
Almendro
V
,
Michor
F
,
Polyak
K
. 
Non-cell-autonomous driving of tumour growth supports sub-clonal heterogeneity
.
Nature
2014
;
514
:
54
8
.
49.
Kaznatcheev
A
,
Peacock
J
,
Basanta
D
,
Marusyk
A
,
Scott
JG
. 
Fibroblasts and alectinib switch the evolutionary games played by non-small cell lung cancer
.
Nat Ecol Evol
2019
;
3
:
450
6
.
50.
Seto
T
,
Kato
T
,
Nishio
M
,
Goto
K
,
Atagi
S
,
Hosomi
Y
, et al
Erlotinib alone or with bevacizumab as first-line therapy in patients with advanced non-squamous non-small-cell lung cancer harbouring EGFR mutations (JO25567): an open-label, randomised, multicentre, phase 2 study
.
Lancet Oncol
2014
;
15
:
1236
44
.
51.
Liccardi
G
,
Hartley
JA
,
Hochhauser
D
. 
Importance of EGFR/ERCC1 interaction following radiation-induced DNA damage
.
Clin Cancer Res
2014
;
20
:
3496
506
.
52.
Liccardi
G
,
Hartley
JA
,
Hochhauser
D
. 
EGFR nuclear translocation modulates DNA repair following cisplatin and ionizing radiation treatment
.
Cancer Res
2011
;
71
:
1103
14
.
53.
Wu
YL
,
Lee
JS
,
Thongprasert
S
,
Yu
CJ
,
Zhang
L
,
Ladrera
G
, et al
Intercalated combination of chemotherapy and erlotinib for patients with advanced stage non-small-cell lung cancer (FASTACT-2): a randomised, double-blind trial
.
Lancet Oncol
2013
;
14
:
777
86
.
54.
Guo
T
,
Ni
J
,
Yang
X
,
Li
Y
,
Li
Y
,
Zou
L
, et al
Pattern of recurrence analysis in metastatic EGFR-mutant NSCLC treated with osimertinib: implications for consolidative stereotactic body radiation therapy
.
Int J Radiat Oncol Biol Phys
2020
;
107
;
62
71
.
55.
Schiller
JH
,
Harrington
D
,
Belani
CP
,
Langer
C
,
Sandler
A
,
Krook
J
, et al
Comparison of four chemotherapy regimens for advanced non-small-cell lung cancer
.
N Engl J Med
2002
;
346
:
92
8
.
56.
Qiao
GB
,
Wu
YL
,
Yang
XN
,
Zhong
WZ
,
Xie
D
,
Guan
XY
, et al
High-level expression of Rad51 is an independent prognostic marker of survival in non-small-cell lung cancer patients
.
Br J Cancer
2005
;
93
:
137
43
.
57.
Torres-Roca
JF
,
Eschrich
S
,
Zhao
H
,
Bloom
G
,
Sung
J
,
McCarthy
S
, et al
Prediction of radiation sensitivity using a gene expression classifier
.
Cancer Res
2005
;
65
:
7169
76
.
58.
Scott
JG
,
Berglund
A
,
Schell
MJ
,
Mihaylov
I
,
Fulp
WJ
,
Yue
B
, et al
A genome-based model for adjusting radiotherapy dose (GARD): a retrospective, cohort-based study
.
Lancet Oncol
2017
;
18
:
202
11
.
59.
Chalkidou
A
,
Landau
DB
,
Odell
EW
,
Cornelius
VR
,
O'Doherty
MJ
,
Marsden
PK
. 
Correlation between Ki-67 immunohistochemistry and 18F-fluorothymidine uptake in patients with cancer: a systematic review and meta-analysis
.
Eur J Cancer
2012
;
48
:
3499
513
.
60.
Rockne
RC
,
Hawkins-Daarud
A
,
Swanson
KR
,
Sluka
JP
,
Glazier
JA
,
Macklin
P
, et al
The 2019 mathematical oncology roadmap
.
Phys Biol
2019
;
16
:
041005
.
61.
Werner
B
,
Scott
JG
,
Sottoriva
A
,
Anderson
AR
,
Traulsen
A
,
Altrock
PM
. 
The cancer stem cell fraction in hierarchically organized tumors can be estimated using mathematical modeling and patient-specific treatment trajectories
.
Cancer Res
2016
;
76
:
1705
13
.
62.
Grassberger
C
,
Ellsworth
SG
,
Wilks
MQ
,
Keane
FK
,
Loeffler
JS
. 
Assessing the interactions between radiotherapy and antitumour immunity
.
Nat Rev Clin Oncol
2019
;
16
:
729
45
.
63.
Walker
R
,
Enderling
H
. 
From concept to clinic: mathematically informed immunotherapy
.
Curr Probl Cancer
2016
;
40
:
68
83
.
64.
Zhang
J
,
Cunningham
JJ
,
Brown
JS
,
Gatenby
RA
. 
Integrating evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer
.
Nat Commun
2017
;
8
:
1816
.
65.
West
JB
,
Dinh
MN
,
Brown
JS
,
Zhang
J
,
Anderson
AR
,
Gatenby
RA
. 
Multidrug cancer therapy in metastatic castrate-resistant prostate cancer: an evolution-based strategy
.
Clin Cancer Res
2019
;
25
:
4413
21
.

Supplementary data