Therapeutic resistance is a fundamental obstacle in cancer treatment. Tumors that initially respond to treatment may have a preexisting resistant subclone or acquire resistance during treatment, making relapse theoretically inevitable. Here, we investigate treatment strategies that may delay relapse using mathematical modeling. We find that for a single-drug therapy, pulse treatment—short, elevated doses followed by a complete break from treatment—delays relapse compared with continuous treatment with the same total dose over a length of time. For tumors treated with more than one drug, continuous combination treatment is only sometimes better than sequential treatment, while pulsed combination treatment or simply alternating between the two therapies at defined intervals delays relapse the longest. These results are independent of the fitness cost or benefit of resistance, and are robust to noise. Machine-learning analysis of simulations shows that the initial tumor response and heterogeneity at the start of treatment suffice to determine the benefit of pulsed or alternating treatment strategies over continuous treatment. Analysis of eight tumor burden trajectories of breast cancer patients treated at Memorial Sloan Kettering Cancer Center shows the model can predict time to resistance using initial responses to treatment and estimated preexisting resistant populations. The model calculated that pulse treatment would delay relapse in all eight cases. Overall, our results support that pulsed treatments optimized by mathematical models could delay therapeutic resistance.

Cancers develop resistance to therapy in up to 70% of cases depending on the cancer type (1). In breast cancer, patients with HER2-amplified metastatic cancer initially respond to HER2-targeting drugs but the disease eventually progresses (2). This pattern of initial response followed by resistance occurs often in targeted therapies (3). Acquired and preexisting resistance to cytotoxic chemotherapies is also common in breast cancer (4). New therapies against resistance can take decades—despite extraordinary efforts (5)—before they are available to patients. The moving target of resistance to therapy is one of today's greatest challenges in cancer.

Resistance is due to tumor heterogeneity; some subpopulations of cancer cells have varying attenuated sensitivity to treatment (6–8). Even one subpopulation in addition to the dominant clone can alter tumor progression and treatment outcomes (9). Mathematical models can bridge the gap between known biological mechanisms and unknown outcomes. Models can predict synergistic drug combinations (10), elucidate how cancer cells evolve under chemotherapy (11), and calculate the likelihood of resistance (12). Mathematical models can also optimize treatments to delay resistance (13). These models consider preexisting tumor heterogeneity as well as de novo resistance mutations and their impact on the emergence of drug resistance.

Previous models suggest that it is possible to optimize treatment strategies with current therapy options to delay the time to resistance. Perhaps the most well-known of these strategies is “adaptive therapy,” in which the patient's response is monitored and treatment dose is consequently adjusted. Importantly, the dose is reduced if the tumor is responding to the therapy, to maintain a population of treatment-sensitive cells, which outcompete and inhibit the growth of resistant cells (14). A model of metastatic castrate-resistant prostate cancer predicted the efficacy of adaptive therapy in which tumors were treated until they dropped to 50% their original size, upon which treatment was ceased until the tumor returned to its original size. Adaptive therapy delayed the takeover by a resistant population, due to its fitness disadvantage compared with the sensitive population in the absence of therapy. A stable tumor size could thus be maintained for longer in the adaptive regimen than if the treatment was given continuously (15). In this model, it was important that the treatment caused a 50% drop in tumor volume or prostate-specific antigen. It was also assumed that resistance confers a fitness cost in the absence of the treatment.

The assumption that resistance has a fitness cost is supported by some data: Resistant non–small cell lung cancer cell lines grew slower than sensitive cells, and continuous low-dose inhibitors with occasional high-dose pulses delayed the onset of tumor resistance in this case (16). A fitness cost of resistance was assumed in mathematical models of breast cancer (17, 18) as well as general models of tumor progression (19). However, some mechanisms of resistance lack a fitness cost (20–22). Resistant subclones may exist at low frequency prior to treatment even if they are fitter than the dominant treatment-sensitive population, simply because they evolved recently (6). This key point was not accounted for in many previous models examining resistance but can alter model outcomes significantly.

Mathematical modeling by Norton and Simon predicted that maximizing the effective dose over as short amount of time as possible increases the chance of cure (23). Dose-dense scheduling for chemotherapy was therefore proposed, which delivers a higher total integrated dosage over a short period of time (23, 24) and limits the time cells have to regrow in between treatments (25). Dose-dense therapy was successful in clinical trials in primary breast cancer, where a higher density increased both disease-free and overall survival with no significant increases in toxicity (24). It is possible that combining the ideas of dose-dense therapy (increased dose) with adaptive therapy (well-timed breaks from treatment) could be an effective strategy.

Here, we develop a model of tumor resistance to study strategies that delay tumor progression. Our model shows that pulsed treatment—short, escalated doses followed by breaks—delays resistance compared with continuous treatment even in conditions where resistance has no fitness cost. Pulsed combination therapy or alternating between therapies were better strategies than sequential treatment when considering more than one therapy. A few measurements made at the start of treatment may be sufficient to parameterize a model and choose the best treatment strategy. We customized the model with parameters from patients treated at Memorial Sloan Kettering Cancer Center (MSKCC); the model predicted time to relapse under various treatment schedules, demonstrating how it could be applied in a clinical setting.

Model development

Model to choose the best schedule for a single therapy in a cancer with two subpopulations

We first model a tumor which has n cells of two kinds: sensitive or resistant to a therapeutic intervention v. Although there will be heterogeneity within each of these groups, we can model each subpopulation using their average properties. ${n_1}[ t ]$ is the number of sensitive tumor cells at time $t$⁠, and ${n_2}[ t ]$ is the number of resistant tumor cells at time $t$⁠, each population with its own average rate of mitosis (⁠${p_m}$⁠) and death (⁠${p_d}$⁠) which we assume are constant in time. Of the ${n_1}$ cells that divide, some daughter cells may acquire a resistance mutation at some low rate ${p_r}$⁠. We encompass these terms in ${g_1}$ and ${g_2}$ for ${n_1}$ and ${n_2}$⁠, respectively:

formula
formula

To simulate tumor growth, we first represent the model using a matrix differential equation, where net growth rates are on the diagonals and transition rates are on the off-diagonals:

formula

Where

formula

We can therefore rewrite equation (A) as

formula

See Supplementary methods for equation derivation and solutions.

Consider that this matrix of cell growths and transitions is different during or in the absence of treatment. We are interested in situations where the value of g1 is positive in the absence of treatment (g1off), and negative in the presence of treatment (g1on). We assume that g2 will be unaffected by therapy; therefore, g2 is always nonnegative even in conditions where resistance confers a relative fitness cost in the absence of therapy. The inputs to simulate the system are: (1) matrix of parameters off-treatment (Moff), (2) matrix of parameters on continuous treatment (M), and (3) vector of initial tumor composition at the start of treatment (⁠$n[ 0 ]$⁠, where t = 0 corresponds to the time of treatment initiation). See Supplementary Methods for detailed inputs and assumptions.

The assumptions of our model impose that if the rate of gaining resistance is greater than zero, or the number of preexisting resistant cells is greater than zero, the resistant cell population will increase over time, consistent with previous studies (26, 27). This outcome holds even when resistant cells have no fitness advantage in the absence of therapy, and even when resistance carries a relative fitness cost. Given the assumptions of our model, resistance and regrowth of the tumor is inevitable even if a treatment initially shrinks the tumor.

We determine the time to resistance (TTR) as the time it takes the tumor to return to its initial size $n\ = \ {n_1}[ 0 ] + \ {n_2}[ 0 ]$ following all possible treatment. We compare two options for scheduling. The first is continuous treatment, where therapy v is given continuously at dose $f$until resistance is reached. The second is pulsed treatment, where a sharp pulse of therapy v at a higher dose (⁠$f*x$⁠) is given for a defined period of time, followed by a break in treatment. For a fair comparison, we escalate the dose given during the on period by a factor $x$ such that the cumulative amount across the sum of on+off periods equals the dose given during that same period had we treated continuously with the lower dose:

formula

Thus, the total dose of the therapy across a period of time is equal between the two strategies, to avoid exceeding the maximum tolerated dose of a therapy in a given time frame. Note that $f$ is not an additional parameter but is already implicitly defined in our model inputs as the difference in sensitive cell growth rate off and on continuous treatment:

formula

See Supplementary Methods for further explanation.

For each simulated tumor with inputs in Eq. (C), the continuous therapy strategy is applied and TTR is calculated. Next, 50 iterations of the pulse treatment strategy are applied to the tumor, each one testing a different randomly generated combination of $x$ and $timeon$⁠. We then calculate the ratio of TTR for each pulse strategy over continuous strategy (TTRp/TTRc), where values greater than 1 indicate a benefit of the pulse treatment compared to continuous treatment. The above simulation is repeated for thousands of randomly generated tumors.

Model expanded to a two-therapy system in a four-subpopulation cancer

We next consider the case where a given tumor has two possible interventions, vA and vB. We can thus classify the tumor cells into four cell types: (i) n1: cells sensitive to both vA and vB, n2: resistant to only vA, n3: resistant to only vB, or n4: resistant to both vA and vB. Each cell type has its own growth rate and can transition from one type to the other by acquiring one level of resistance. Double-resistant cells cannot transition to any other cell type. Growth and transition rates can be conveniently described using matrices as before:

formula

Where pRi,j is the transition rate from cell type i to cell type j: pmi*pri,j. Inputs into the model are: (1) matrix of parameters off-treatment (Moff), (2) matrix of parameters on vA (MA), (3) matrix of parameters on vB (MB), (4) vector of initial tumor composition at the start of treatment (n[0]). See Supplementary Methods for detailed inputs and assumptions.

We consider “standard” treatment to be continuous treatment, in which a tumor is first treated continuously with vA at dose $fA$ until resistance to vA is reached, followed by vB at dose $fB$ until resistance to vB is reached (or vice versa). We then compare this sequential strategy to the following four treatment strategies: (i) Combination treatment, in which tumors are treated simultaneously with vA and vB, but the total dose does not change. To achieve this, the greater value between $fA$ and $fB$ is selected as the total dose, and $fA$ and $fB$ are scaled down at their respective ratios. (ii) Alternating treatment, in which $fA$ is given for a defined period of time less than the time it would take to reach resistance, followed by $fB$ for a defined period of time less than the time it would take to reach resistance. This is repeated, switching back and forth, until resistance to both vA and vB is reached. (iii) Pulsed-alternating treatment, in which a sharp dose of vA is given (⁠$x*fA$⁠) for a defined period of time (⁠$timeonA$⁠) followed by a break akin to Eq. (E). This is followed by a sharp dose of vB (⁠$x*fB$⁠) for a defined period of time (⁠$timeonB$⁠) followed by a corresponding break. These are repeated until resistance to both vA and vB is reached. (iv) Pulsed combination treatment, in which the simultaneous therapy created in option (i) is applied at a sharp pulse followed by a break as defined in Eq. (E).

For each simulated tumor, TTR under each strategy is simulated. For strategies that require random generation of $x$ and $timeon$ values, 50 iterations are tested as were done previously. The above simulation is repeated for thousands of randomly generated tumors.

The differential equations and assumptions in our model—including proportionally reducing the dose for each treatment in a combination cocktail, non-cross–resistant therapies, independent exponential growth of clones (no competition), dose-dependent reduction in growth rate, and no requirement for resistance bearing a fitness cost—have been adopted before in previous models (28–31).

Stochastic simulations

To simulate the effect of stochasticity in the model, we simulated the system using Euler's method. At each Euler step, we added noise centered around a preset amount and proportional to the number of cells in that cell type.

formula

Where $i$ is a cell type, $noisevariance$ is set to 0.1, and $randn$ is a random scalar drawn from the standard normal distribution. This leads to noise up to 100× the transition rate. This noise was added to the calculated value of ${n_i}[ t ]$⁠. All simulations and treatment strategies are defined as before and were implemented using this method.

Technical details of simulations

All simulations, analyses, and calculations were performed using custom MATLAB scripts and functions. See Supplementary Methods for details.

Patient data

Mutation and copy-number data, treatment data, and clinical metadata were obtained for patients with breast cancer at our institution from both published (32) and unpublished data. All data were in accordance with recognized ethical guidelines and approved by an institutional review board (IRB). Data were obtained from a repository at our institution from samples that had been collected under previous IRBs and previously published; therefore, a new IRB for this study was not required. The data usage in this study was approved by the 12–245 Data & Tissue Usage Committee and privacy review. A subset of patients was selected, and TTR for each treatment was determined by calculating the approximate number of days until the tumor marker value reached the value it had at the start of treatment. The proportion of ${n_1}[ 0 ]$ cells was approximated with pathology or sequencing data. Either a continuous or pulsed model was assumed for each patient's trajectory based on the treatment given, and simulated values of $n$ were compared with a patient's first two or three tumor marker data points to determine the likely M and Moff matrices for each. Predicted trajectories under various treatment strategies were modeled using these inputs. Because the optimal pulse strategy cannot be determined until the model is fit to the parameters determined from these data points, we assumed the actual dosing for the period of time during the first two data points for each patient and only began the optimal strategy subsequently. See Supplementary Methods for more details.

Data availability statement

All codes, simulation data generated in this study, and deidentified patient data are available in the following repository on github: https://github.com/dm2791/Optimal-strategy-and-benefit-of-non-continuous-therapy.

Two-subpopulation model showed pulse treatment preferred over continuous treatment

We first considered the case where only one therapy was applied, and the cells in the heterogenous tumor were therefore grouped into: (i) cells sensitive to this intervention and (ii) cells resistant to this intervention (Fig. 1A). A randomly simulated tumor showed the trajectory of the tumor under continuous treatment: the tumor initially shrank upon treatment but relapsed once resistant cells accumulated (Fig. 1B, black line). When the same tumor was given sharp escalated-dose pulses of the therapy at defined intervals, with no-therapy breaks in between, the TTR of the tumor was delayed (Fig. 1B, pink line). In this example, relapse was due to a high rate of transition to resistance, such that a quicker elimination of sensitive cells in the pulse therapy decreased the number of cells available to develop resistance (Supplementary Fig. S1A). In another example (Supplementary Fig. S1B and S1C), pulse therapy delayed TTR compared with continuous treatment although the transition rate to resistance was low. Unlike the previous example, this randomly simulated tumor had a preexisting resistant population, which rapidly increased and caused tumor relapse. Pulse therapy was effective in this case because it shrank the sensitive cell population more rapidly, so that even with accumulating resistant cells the overall tumor stayed smaller for a longer period of time.

Figure 1.

Single-treatment simulations show that pulsed treatment is always better than continuous treatment, and its benefit is higher for aggressive tumors. A, In the one-therapy model, sensitive and resistant populations each have an average rate of mitosis (m), death (d), and stagnating (s). Of the sensitive cells that divide, the daughter cells also have an average rate of acquiring resistance (r). B, Simulated tumor treated continuously with baseline dose or with pulse followed by an off period yielding about 1.5× delay in resistance. C, Histogram of 50,000 randomly simulated tumors. Time to resistance could be increased in 100% of them by using a pulse treatment strategy. D, PCA of data inputs from 50,000 randomly simulated tumors, shaded by log(pulse advantage). E, The best predictors of pulse advantage. F, Correlation between ${g_1}off/{g_2}off$ and log(pulse advantage), color indicating local density of data. G, Correlation between log(sum(Moff*$n[ 0 ]$⁠)) and log(pulse advantage), color indicating local density of data. H, Correlation between sum(M*$n[ 0 ]$⁠) and log(pulse advantage), color indicating local density of data.

Figure 1.

Single-treatment simulations show that pulsed treatment is always better than continuous treatment, and its benefit is higher for aggressive tumors. A, In the one-therapy model, sensitive and resistant populations each have an average rate of mitosis (m), death (d), and stagnating (s). Of the sensitive cells that divide, the daughter cells also have an average rate of acquiring resistance (r). B, Simulated tumor treated continuously with baseline dose or with pulse followed by an off period yielding about 1.5× delay in resistance. C, Histogram of 50,000 randomly simulated tumors. Time to resistance could be increased in 100% of them by using a pulse treatment strategy. D, PCA of data inputs from 50,000 randomly simulated tumors, shaded by log(pulse advantage). E, The best predictors of pulse advantage. F, Correlation between ${g_1}off/{g_2}off$ and log(pulse advantage), color indicating local density of data. G, Correlation between log(sum(Moff*$n[ 0 ]$⁠)) and log(pulse advantage), color indicating local density of data. H, Correlation between sum(M*$n[ 0 ]$⁠) and log(pulse advantage), color indicating local density of data.

Close modal

We compared continuous and pulsed treatment in 50,000 tumor models with randomly generated input parameters. The TTR could be extended by using a pulse strategy rather than continuous in 100% of cases. In the worst cases, the improvement was marginal, but in the best cases, the pulse treatment delayed TTR by >140× compared to the continuous treatment (Fig. 1C). The extent to which the pulse strategy wins we call the “pulse advantage.” Importantly, in a given time frame, we do not allow the total dose administered in the pulse strategy to exceed the total dose administered in the continuous strategy.

Machine learning shows pulsed strategy more advantageous on aggressive tumors

While the pulse strategy is always better than continuous, we find that there is a fairly wide range of the pulse advantage (Fig. 1D; Supplementary Fig. S1D). To ascertain the properties of a tumor that predispose to a high pulse advantage, we first calculated the following values for each simulated tumor: (i) sum(Moff*n[0]): the net growth rate of the whole tumor in the absence of treatment assuming its t = 0 composition; (ii) ${g_1}$off/${g_2}$off: the ratio of growth rates of sensitive and resistant cells in the absence of treatment; (iii) sum(M*$n[ 0 ]$⁠): the net growth rate of the whole tumor on continuous treatment assuming its t = 0 composition; (iv) ${g_1}$on/${g_2}$on: the ratio of growth rates of sensitive and resistant cells on continuous treatment; (v) $f$⁠: ${g_1}$off-${g_1}$on, the effect of the continuous treatment dose on sensitive cells; (vi) ${n_2}[ 0 ]/{n_1}[ 0 ]$⁠: the ratio of resistant to sensitive cells at t = 0 (time of therapy initiation); (vii) ${p_R}on/{p_R}off$⁠: the ratio of transition rates while on therapy versus off therapy. These calculations allow us to characterize the properties of each simulated tumor as relative values, and some of these quantities may be measurable in a clinical setting.

We then used a machine learning approach to identify which of these tumor features were most informative to determine the success of a treatment strategy. The algorithm identified that parameters (1), (2), and (3) provide the greatest predictive power in determining the degree of pulse advantage (Fig. 1E). In general, a smaller ${g_1}$off/${g_2}$off ratio, higher sum(Moff*n[0]), and higher sum(M*$n[ 0 ]$⁠) lead to a greater pulse advantage (Fig. 1,FH). Hence, tumors in which resistance confers a fitness advantage in the absence of treatment are more likely to have a large pulse advantage, as are tumors that are more aggressively growing in the absence of therapy and/or are poor responders to the given therapy. This is likely because a stronger dose can more effectively shrink the population that generates resistance (23, 33). A model built with these three predictors alone accurately predicted pulse advantage in a test dataset of an additional 10,000 randomly simulated tumors with an R2 of 0.9 (Supplementary Fig. S1E).

The tumors for which the model failed to accurately predict pulse advantage (Supplementary Fig. S1F) tended to have weaker responses to continuous treatment (Supplementary Fig. S1G). This is interesting, because this is the metric that the machine learning analysis had deemed most predictive of pulse success, but also shows that tumors at one extreme end of this metric tend to be harder to predict.

Determining the fitness cost or benefit of resistance (the quantity ${g_1}$off/${g_2}$off) may be difficult to clinically ascertain unless very well-characterized mechanisms of resistance dominate the tumor. However, it is easier to determine the growth of the bulk tumor on and off treatment, for example, by measuring tumor size at multiple timepoints. Therefore, we trained a model with only sum(Moff*n[0]) and sum(M*$n[ 0 ]$⁠) as predictors, and the model still accurately predicted pulse advantage with an R2 of 0.8 (Supplementary Fig. S2A). As before, the more poorly responding tumors tended to require more data (i.e., ${g_1}$off/${g_2}$off) to generate good predictive power, while well-responding tumors are equally well-predicted in both models (Supplementary Fig. S2B and S2C).

In Supplementary Fig. S1D, there appears to be an inflection point at y = 0.79 (in log scale, corresponding to a 2.2-fold increase in time to relapse using a pulse strategy). We used this inflection point to binarize pulse advantage into “high” and “low” categories and trained a binary classifier. When tested on an independent set of 10,000 tumors, the model predicted tumors with high/low pulse advantage with an accuracy of 0.97 (AUC = 0.99) when using all the interaction terms calculated above, and 0.96 (AUC = 0.98) when using only sum(Moff*n[0]) and sum(M*$n[ 0 ]$⁠) as predictors (Supplementary Fig. S2D).

For the majority of tumors, every combination of strength (dose escalation, which affects ${g_1}on$⁠) and duration of the pulse yielded a pulse advantage greater than 1, but for some tumors, certain combinations did not improve TTR compared to continuous treatment (Supplementary Fig. S2E). Time spent off treatment had the highest frequency of being the strongest correlate with pulse advantage, with higher values of $timeoff$ (equal to $timeon*( {x - 1})$ correlating with higher pulse advantage (Supplementary Fig. S2F).

Sequential treatment leads to resistance more quickly than other strategies

We next expanded our model to consider cases with two therapies (Fig. 2A). We now have at least five potential treatment strategies, and consider the following: (i) continuous, in which therapy vA is given until resistance is reached, sequentially followed by therapy vB until resistance is reached; (ii) combination treatment, in which vA and vB are given simultaneously and continuously until resistance is reached, but the dose of the combined cocktail does not exceed the dose of either therapy alone; (iii) alternating treatment, in which defined intervals of vA and vB are oscillated back and forth with no dose escalation; (iv) combination pulse, in which the combined cocktail in (ii) is given at an escalated dose pulse followed by a period of no therapy and this sequence is repeated until full resistance is reached; and (v) alternating pulse, in which the defined vA interval is given at an escalated dose pulse, followed by a period of no therapy, followed by a vB pulse, followed by no therapy, etc.

Figure 2.

Simulations with two-treatment model show that sequential continuous treatments lead to resistance sooner than alternating or pulsed strategies. A, In the two-therapy model, double sensitive, resistant to vA, resistant to vB, and double resistant populations each have an average rate of mitosis (m), death (d), and stagnating (s). Of the sensitive cells that divide, the daughter cells also have an average rate of acquiring resistance (r). B, PCA of tumor properties from 100,000 randomly simulated tumors, colored by the treatment strategy that led to the longest time to resistance. Points closer together indicate tumors that have similar properties. Inset: pie chart depicting the percent that each strategy led to the longest time to resistance. C, Simulated tumor treated continuously with vA until no longer effective, followed by vB until no longer effective (black line), or with pulse of vA+vB followed by off period (blue line). This yields about a 3.5× delay in resistance. D, Pie chart of the data from (B) recalculated to include only strategies without a pulse, as percent of simulated tumors. E, Simulated tumor treated continuously with vA until no longer effective, followed by vB until no longer effective (black line), or with alternating vA and vB at defined short intervals. This yields about a 3.7× delay in resistance. F, Biplot of the advantage of each strategy relative to continuous. Vectors with similar direction and magnitude correspond to variables that have similar response profiles. G, The best predictors of combination pulse advantage are (1) relative fitness of double resistant cells in the absence of treatment (⁠${g_4}off/( {{g_1}off + {g_2}off + {g_3}off} )$⁠), (2) relative fitness of cells resistant to one therapy in the absence of treatment(⁠$({g_1}off + {g_3}off)/({g_2}off + {g_4}off)$⁠), (3) growth rate of the bulk tumor on continuous treatment vA (sum(MA*$\ n[ 0 ]$⁠), a negative value). H, The best predictors of alternating advantage are (1) relative fitness of double resistant cells in the absence of treatment (⁠${g_4}off/( {{g_1}off + {g_2}off + {g_3}off} )$⁠), (2) proportion of double-sensitive cells at the start of treatment (⁠${n_1}[ 0 ]$⁠), (3) proportion of double-resistant cells at the start of treatment (⁠${n_4}[ 0 ]$⁠).

Figure 2.

Simulations with two-treatment model show that sequential continuous treatments lead to resistance sooner than alternating or pulsed strategies. A, In the two-therapy model, double sensitive, resistant to vA, resistant to vB, and double resistant populations each have an average rate of mitosis (m), death (d), and stagnating (s). Of the sensitive cells that divide, the daughter cells also have an average rate of acquiring resistance (r). B, PCA of tumor properties from 100,000 randomly simulated tumors, colored by the treatment strategy that led to the longest time to resistance. Points closer together indicate tumors that have similar properties. Inset: pie chart depicting the percent that each strategy led to the longest time to resistance. C, Simulated tumor treated continuously with vA until no longer effective, followed by vB until no longer effective (black line), or with pulse of vA+vB followed by off period (blue line). This yields about a 3.5× delay in resistance. D, Pie chart of the data from (B) recalculated to include only strategies without a pulse, as percent of simulated tumors. E, Simulated tumor treated continuously with vA until no longer effective, followed by vB until no longer effective (black line), or with alternating vA and vB at defined short intervals. This yields about a 3.7× delay in resistance. F, Biplot of the advantage of each strategy relative to continuous. Vectors with similar direction and magnitude correspond to variables that have similar response profiles. G, The best predictors of combination pulse advantage are (1) relative fitness of double resistant cells in the absence of treatment (⁠${g_4}off/( {{g_1}off + {g_2}off + {g_3}off} )$⁠), (2) relative fitness of cells resistant to one therapy in the absence of treatment(⁠$({g_1}off + {g_3}off)/({g_2}off + {g_4}off)$⁠), (3) growth rate of the bulk tumor on continuous treatment vA (sum(MA*$\ n[ 0 ]$⁠), a negative value). H, The best predictors of alternating advantage are (1) relative fitness of double resistant cells in the absence of treatment (⁠${g_4}off/( {{g_1}off + {g_2}off + {g_3}off} )$⁠), (2) proportion of double-sensitive cells at the start of treatment (⁠${n_1}[ 0 ]$⁠), (3) proportion of double-resistant cells at the start of treatment (⁠${n_4}[ 0 ]$⁠).

Close modal

We simulated 100,000 randomly generated tumors and applied each of the treatment strategies described above. We found that continuous treatment was never the best strategy for a particular tumor, and combination pulse was most often the best strategy (Fig. 2B and C; Supplementary Fig. S3A and S3B). Alternating pulse was the second most likely to be the best, and—interestingly—was the only strategy that was always superior to continuous treatment (Fig. 2B; Supplementary Fig. S3C–S3E; Table 1 columns 1–2). We can consider the drug cocktail in combination pulse therapy to behave the same as a single therapy. The same intuition therefore applies: pulses eliminate populations that can generate resistance faster than sequential treatment does. Alternating therapies likewise delays resistance because cells resistant to the first therapy can be killed during the second therapy and vice versa. Alternating at short intervals delays the chance for a double-resistant clone to dominate the population. This is similar to results from previous models but has personalized schedules for each patient (29).

Table 1.

Best treatment strategies as revealed by simulations with the two-therapy model.

Better than continuous strategyBest optionBetter than continuous strategyBest option
Continuous — 0% — 0% 
Combination 44% 0% 80% 3% 
Alternating 99% 18% 100% 4% 
Combination pulse 95% 62% 95% 68% 
Alternating pulse 100% 20% 100% 25% 
Better than continuous strategyBest optionBetter than continuous strategyBest option
Continuous — 0% — 0% 
Combination 44% 0% 80% 3% 
Alternating 99% 18% 100% 4% 
Combination pulse 95% 62% 95% 68% 
Alternating pulse 100% 20% 100% 25% 

Note: Columns 1 and 2, deterministic version; columns 3 and 4, stochastic version.

If we consider only the treatment strategies that do not include a pulse (strategies i, ii, and iii) alternating between vA and vB (strategy iii) was the best option in the majority of cases (Fig. 2D and E; Supplementary Fig. S3F–S3G). This is important, given that there may be cases where even a short-term dose escalation is not feasible, although the overall dose does not change (23). While combination pulse had the highest mean advantage and mean success rate, alternating pulse had the highest maximum advantage (Supplementary Fig. S3H–S3J).

We noted that most tumors that were predicted to do well on a combination pulse treatment were also predicted to do well on an alternating pulse treatment (Fig. 2F), and were curious whether similar tumor properties lead to a high advantage in both strategies. We therefore calculated the following relative parameters for each simulated tumor: (1) sum(Moff*$\ n[ 0 ]$⁠), (2) sum(MA*$\ n[ 0 ]$⁠), (3) sum(MB*$\ n[ 0 ]$⁠), (4–7) the percentage of each cell type at t = 0, (8–11) the growth rate of each cell type relative to the others when off treatment, (12–15) the growth rate of each cell type relative to the others when on vA (continuous dose), (16–19) the growth rate of each cell type relative to the others when on vB (continuous dose), (20–21) $fA$ and $fB$⁠, and (22–23) the increase in transition rate when on vA or vB relative to off-treatment. We then used a machine learning algorithm to determine the parameters that best predicted the extent of the advantage of either of these strategies.

As hypothesized, we found that the top predictors were similar in both cases (Fig. 2G; Supplementary Fig. S3K). The relative fitness cost of resistance to both therapies when off treatment, as well as the tumor's initial response to baseline doses of vB and/or vA were important predictors in both strategies to determine the extent of their advantage over continuous treatment. As with the single-treatment model, a poorer response to treatment increases the benefit of either of these pulse strategies. In addition, knowing the relative fitness cost or benefit of resistance to just one of the therapies off-treatment added predictive power to the advantage of combination pulse treatment. On the other hand, the top predictors for determining the extent of the advantage of alternating treatment included the proportion of double sensitive and double resistant cells before starting therapy (Fig. 2H). A lower proportion of double resistant cells at the start of therapy increased the benefit of an alternating treatment strategy over continuous.

Due to the similarities between combination pulse and alternating pulse, we grouped these together as a “pulse” strategy and trained a classifier to predict the optimal strategy (pulse, alternating, combination, or continuous) for 60,000 simulated tumors. The best predictors were the tumor's initial response to baseline doses of vA, the proportion of double sensitive cells, and the proportion of cells only sensitive cells to vB (Supplementary Fig. S3L). When tested on an additional data set of 40,000 tumors, these three parameters predicted the best treatment strategy with an accuracy of 82%. However, this accuracy rate is equal to the percent of cases in which either combination pulse or alternating pulse was the best strategy.

Model results are robust to stochasticity

While the ODE models allowed us to gain insight into the various treatment strategies, biological systems are stochastic and noisy, making predictions even more difficult. To determine how stochasticity might affect our results we conducted simulations that introduced random noise at each time step.

We first implement this method in the 2-cell-type system. As with the ODE model, the pulse treatment delayed time to relapse (Fig. 3A; Supplementary Fig. S4A). However, unlike the ODE model, the fluctuations allowed in the Euler method could now lead to tumor extinction in some cases (Fig. 3B; Supplementary Fig. S4B). The pulse treatment strategy was always better than the continuous treatment strategy for tumors that would relapse, and we could now achieve tumor extinction in nearly 20% of cases (Fig. 3C). Moreover, when we added noise into the model, the maximum pulse advantage increased (Fig. 3C). Interestingly, tumors that were “cured” did not cluster with tumors that had the highest pulse advantage, suggesting that the properties that lead to extinction are independent of the properties previously identified that lead to a high pulse advantage. To identify these properties, we applied a binary machine learning classifier and found that the greatest predictor of tumor extinction is the ratio of sensitive/resistant cell growth when on continuous treatment (Fig. 3D). The more negative values (therefore the sharpest decline in sensitive cell growth) lead to a greater chance of tumor extinction, due to the ability of random fluctuation to eliminate tumor cells when they are present at smaller numbers.

Figure 3.

Stochastic model shows simulation results are robust to noise. A, Simulated tumor treated continuously with baseline dose or pulse treatment, allowing random fluctuation at each time step. This yields about a 3× delay in resistance. B, The same simulated tumor as (A), treated continuously with baseline dose or pulse treatment and allowing random fluctuation at each time step. This time, the pulse treatment led to extinction of the tumor due to fluctuation that eliminated remaining tumor cells. C, PCA of data inputs from 15,000 randomly simulated tumors in single-therapy model, shaded by degree of pulse advantage or whether tumor extinction was achieved. Inset: proportion of events where relapse delay or tumor extinction was achieved. Time to resistance could be increased or tumor extinction could be achieved in 100% of them by using a pulse treatment strategy compared to a continuous strategy. D, The best predictors of whether a tumor will reach extinction are (1) sensitive cell response to therapy relative to resistant cell growth (⁠${g_1}on/{g_2}on$⁠), (2) relative growth rates of each cell type in the absence of treatment (⁠${g_1}off/{g_2}off$⁠), (3) effective treatment dose (⁠$f$⁠). E, PCA of data inputs from 15,000 randomly simulated tumors in two-therapy model, colored by strategy that had the longest time to relapse or led to tumor extinction (“cured”). Inset: top pie chart: the percent that each strategy was the most successful; bottom pie chart: of the 1% of overall events in which tumor extinction was achieved, the percent that each treatment strategy contributed to the number of extinction events.

Figure 3.

Stochastic model shows simulation results are robust to noise. A, Simulated tumor treated continuously with baseline dose or pulse treatment, allowing random fluctuation at each time step. This yields about a 3× delay in resistance. B, The same simulated tumor as (A), treated continuously with baseline dose or pulse treatment and allowing random fluctuation at each time step. This time, the pulse treatment led to extinction of the tumor due to fluctuation that eliminated remaining tumor cells. C, PCA of data inputs from 15,000 randomly simulated tumors in single-therapy model, shaded by degree of pulse advantage or whether tumor extinction was achieved. Inset: proportion of events where relapse delay or tumor extinction was achieved. Time to resistance could be increased or tumor extinction could be achieved in 100% of them by using a pulse treatment strategy compared to a continuous strategy. D, The best predictors of whether a tumor will reach extinction are (1) sensitive cell response to therapy relative to resistant cell growth (⁠${g_1}on/{g_2}on$⁠), (2) relative growth rates of each cell type in the absence of treatment (⁠${g_1}off/{g_2}off$⁠), (3) effective treatment dose (⁠$f$⁠). E, PCA of data inputs from 15,000 randomly simulated tumors in two-therapy model, colored by strategy that had the longest time to relapse or led to tumor extinction (“cured”). Inset: top pie chart: the percent that each strategy was the most successful; bottom pie chart: of the 1% of overall events in which tumor extinction was achieved, the percent that each treatment strategy contributed to the number of extinction events.

Close modal

We further implemented stochasticity to the four-cell-type model and found a similar distribution in winning strategies as in Fig. 2C, although with a stronger propensity to favor pulse treatments (Fig. 3E, top pie chart; Table 1 columns 3–4). This was likely due to pulse treatments reducing the sensitive cell population to smaller minima than other treatment strategies, allowing random fluctuation to eliminate these small populations. Interestingly, the combination pulse strategy clusters separately from the other strategies (Supplementary Fig. S4C), unlike in our ODE model, which may in part be due to the strongest likelihood of combination pulse leading to tumor extinction (Fig. 3E, bottom pie chart).

The stochastic simulations showed that our main results are robust to random noise: pulse strategies win over continuous strategies, and while combination pulse is the most likely to lead to the longest relapse delay, alternating strategies are always better than continuous strategies in terms of TTR.

Patient data illustrate clinical significance of models

Our models have thus far shown the theoretical benefit of noncanonical treatment strategies in order to delay therapeutic resistance in tumors. It is critical to demonstrate how these models can be applied to currently available clinical data. We turned to a breast cancer cohort from our institution (32), consisting of more than 1500 patients of whom 60% developed resistance to at least one therapy and 46% developed resistance to two or more therapies (Supplementary Fig. S5A). In fact, therapeutic resistance was the most common reason for stopping a treatment (Supplementary Fig. S5B). However, the time to reach resistance was highly variable, highlighting that predicting the emergence of resistance is important and nontrivial (Supplementary Fig. S5C).

From four patients in this cohort, we obtained cancer antigen 15.3 (CA.15.3) data, a tumor marker used to monitor metastatic breast cancer progression (34), giving us a longitudinal proxy for tumor growth or decline (ref. 35; Supplementary Fig. S5D–S5G).

Two of these patients were given treatments that are generally administered continuously (36, 37). From the first two CA.15.3 measurements for each therapy, we calculated the initial slope of tumor regression. Our model had found that higher (less negative) values of sum(M*$\ n[ 0 ]$⁠)—the initial slope of the bulk tumor in the presence of continuous therapy—led to a greater pulse advantage (Fig. 1H). We would therefore predict that patient 288 treated with an everolimus-exemestane combination would have the least benefit of a pulse strategy, and that patient 468 treated with letrozole would have the greatest benefit of a pulse strategy (Fig. 4A).

Figure 4.

Patient data can be used to select model and predict benefit of pulsed treatment. A, Slope of tumor markers between first two data points for each treatment that is dosed daily. B, Actual (tumor marker) and predicted (modeled) tumor burden over time for patient 468 treated with letrozole, as well as predicted trajectory on pulse treatment. Vertical lines mark time to relapse for each scenario (TTR_actual, TTR_pred, and TRR_pulse, respectively). Input parameters for trajectory predictions were: 10% preexisting resistance (see Table 2), the first 2 CA.15.3 values in this treatment window, and the first 2 CA.15.3 values in an off-treatment window for this patient (see Supplementary Fig. S5F). Transition rate was assumed to be 10–4. C, Bar plot of the ratio of predicted over actual time to resistance for each treatment (teal bars), overlaid with bar plot of the pulse advantages for each treatment (navy blue bars). Mean values for each depicted by dashed lines. D, Actual and predicted tumor burden over time for patient 422 treated with eribulin followed by paclitaxel, as well as predicted trajectory on continuous, alternating, and combination pulse treatments. Vertical lines mark time to relapse for each scenario (TTR_actual, TTR_pred, TRR_cont, TTR_alternating, and TTR_combopulse, respectively). Input parameters for trajectory predictions were: 11% preexisting double resistance, 28% pre-existing resistance due to PIK3CA mutation, and 19% preexisting resistance due to PTEN mutation (see Table 2), the first 2 CA.15.3 values in both eribulin and paclitaxel treatment windows, and the first 2 CA.15.3 values in an off-treatment window for this patient (see Supplementary Fig. 5E). Transition rates were assumed to be 10–4.

Figure 4.

Patient data can be used to select model and predict benefit of pulsed treatment. A, Slope of tumor markers between first two data points for each treatment that is dosed daily. B, Actual (tumor marker) and predicted (modeled) tumor burden over time for patient 468 treated with letrozole, as well as predicted trajectory on pulse treatment. Vertical lines mark time to relapse for each scenario (TTR_actual, TTR_pred, and TRR_pulse, respectively). Input parameters for trajectory predictions were: 10% preexisting resistance (see Table 2), the first 2 CA.15.3 values in this treatment window, and the first 2 CA.15.3 values in an off-treatment window for this patient (see Supplementary Fig. S5F). Transition rate was assumed to be 10–4. C, Bar plot of the ratio of predicted over actual time to resistance for each treatment (teal bars), overlaid with bar plot of the pulse advantages for each treatment (navy blue bars). Mean values for each depicted by dashed lines. D, Actual and predicted tumor burden over time for patient 422 treated with eribulin followed by paclitaxel, as well as predicted trajectory on continuous, alternating, and combination pulse treatments. Vertical lines mark time to relapse for each scenario (TTR_actual, TTR_pred, TRR_cont, TTR_alternating, and TTR_combopulse, respectively). Input parameters for trajectory predictions were: 11% preexisting double resistance, 28% pre-existing resistance due to PIK3CA mutation, and 19% preexisting resistance due to PTEN mutation (see Table 2), the first 2 CA.15.3 values in both eribulin and paclitaxel treatment windows, and the first 2 CA.15.3 values in an off-treatment window for this patient (see Supplementary Fig. 5E). Transition rates were assumed to be 10–4.

Close modal

To test our predictions, we first aimed to accurately extrapolate time to relapse for each patient's treatments. While we have the benefit of hindsight in these examples, in actual clinical scenarios decisions regarding treatment schedule must be made with limited early information. We therefore returned to our model, starting first with the simple case of a single therapy:

formula

We now have information for the left side of the equation since we know ${n_1}$[t] + ${n_2}[ t ]$ at several timepoints. We thus want to reverse engineer the growth matrix. To do this, we also need to find ${n_1}$[0] and ${n_2}[ 0 ]$⁠. We approximated these values using mutational, copy number, and histological data to estimate the proportion of pre-existing resistance (Table 2): loss of ER expression can cause resistance to letrozole (38), oncogenic myc can cause resistance to nab-paclitaxel (39, 40), ESR1 mutations can cause resistance to exemestane (41), and mutations in the PI3K pathway can cause resistance to gemcitabine–carboplatin, eribulin, and paclitaxel (42–45). Although there may be other mechanisms of resistance for each of these therapies, the aforementioned genetic alterations were present in the patient samples treated with each respective therapy. Other alterations can also cause resistance to aromatase inhibitors, but in some cases sequencing was done after the aromatase inhibitor treatment window so we could not use this information in our predictive model (32).

Table 2.

Four patients with breast cancer treated at MSKCC.

PatientTreatmentResistance mechanismData typeEstimated proportion preexisting resistance
468 Letrozole ER loss IHC 10% 
468 Nab-paclitaxel Myc amplification Copy-number analysis 13.4% 
15 Gemcitabine-carboplatin PIK3CA gain-of-function missense mutation Sequencing 41% 
288 Letrozole ER loss IHC 14% 
288 Everolimus-exemestane ESR1 missense mutation (exemestane resistance) Sequencing 22.6% 
422 Eribulin PIK3CA gain-of-function missense mutation Sequencing 35% 
422 Paclitaxel PTEN loss-of-function missense mutation Sequencing 24% 
PatientTreatmentResistance mechanismData typeEstimated proportion preexisting resistance
468 Letrozole ER loss IHC 10% 
468 Nab-paclitaxel Myc amplification Copy-number analysis 13.4% 
15 Gemcitabine-carboplatin PIK3CA gain-of-function missense mutation Sequencing 41% 
288 Letrozole ER loss IHC 14% 
288 Everolimus-exemestane ESR1 missense mutation (exemestane resistance) Sequencing 22.6% 
422 Eribulin PIK3CA gain-of-function missense mutation Sequencing 35% 
422 Paclitaxel PTEN loss-of-function missense mutation Sequencing 24% 

We then ran several simulations with combinations of ${g_1}on$ and ${g_2}on$ selected randomly, and compared simulated tumor size to the first two CA.15.3 data points. We selected the best match and then continued the simulation to predict what the time to relapse would be given this growth matrix.

For the treatments that are administered daily (letrozole and everolimus–exemestane), we fit a continuous model to the data. For treatments that are administered once every several days (nab-paclitaxel, gemcitabine–carboplatin, eribulin, paclitaxel), we fit a pulsed model to the data. Using the most likely growth matrices for each treatment and estimated initial vectors as described above, we predicted actual time to relapse for each of these treatment courses (Fig. 4B, black line, Fig. 4C, teal bars; Supplementary Fig. S6A and S6B, black lines, and S6C–S6F, solid pink lines).

Likewise, we estimated Moff matrices for each of the patients. We then used these inputs for the single-therapy model described in Fig. 1. For the treatments that are administered continuously, we predicted the theoretical pulse advantage for each scenario (Fig. 4B; Supplementary Fig. S6A and S6B, pink lines). We implemented the constraint that the pulse strategy should not increase the baseline dose by more than three times, to help keep pulsed doses in the similar range as standard doses and not deviate too far on the pharmacokinetic curve. Pulses only start after the model has been parameterized by the patient's data. As expected, patient 468 treated with letrozole had the highest predicted pulse advantage. However, while patient 288 treated with everolimus–exemestane had the lowest predicted pulse advantage based on initial tumor regression alone, patient 288 treated with letrozole had the actual lowest pulse advantage, perhaps due to the very long time to relapse under continuous treatment or incorrect estimations of $n$[0]. Still, by using only the first few tumor marker data points and single-timepoint histological or mutation data, we (i) predicted actual time to relapse with reasonable accuracy for single therapies treated continuously, and (ii) determined that pulse strategies would delay relapse in these clinical cases.

For the treatments that are already administered as “pulses,” we used our model to determine how the TTR would change had the therapies been—theoretically—given continuously [at a correspondingly lower dose as calculated by Eq. (E)]. As expected, we found that the theoretical continuous treatment strategy would have led to earlier time to relapse than the current pulse schedule (Supplementary Fig. S6C–S6F, black lines). We then used our model to predict what the optimal pulse strategy would be for each of these scenarios, and found that we could theoretically extend relapse time even further by optimizing the schedule, even if this optimized schedule only starts after full normal chemotherapy cycle(s) (Supplementary Fig. S6C–S6F, dotted pink lines). Overall, across all patient trajectories we predicted that optimal pulse treatments would delay resistance compared to continuous treatment 1.25× on average, corresponding to about 8 months (Fig. 4C, blue bars).

Our pulse schedules implicitly assume a negligible half-life, with off-treatment periods reverting to the growth dynamics of the untreated tumor. However, the half-life of some of these drugs could affect the optimal pulse strategy employed. Letrozole has a half-life often in the range of 2 to 4 days, and we therefore tested different possible half-lives surrounding this range for patient 468 treated with letrozole (46). As we increase the half-life, the “pulsed” curve begins to look smoother and more closely resembles a continuous-dose curve, due to longer clearance times, which leads to extended anti-tumor activity (Supplementary Fig. S7A–S7F). Therefore, even in the “off” periods of the pulse schedule, some drug remains and contributes to tumor regression. Specific pharmacokinetics depend on the physiology of the individual; therefore, accurate modeling of pharmacokinetics requires additional clinical data. These examples nonetheless provide a proof of concept of how the same drug may cause different responses due to half-life variability. In this case, the differences in half-life did not affect the overall time to relapse (Supplementary Fig. S7A–S7F) because relapse is driven by growth of resistant cells, which are not responsive to letrozole in any of the scenarios. A caveat here is that when the half-life was longer, the effective dose accumulated to higher levels because clearance took longer. This may be acceptable for many patients, akin to how dose-dense scheduling increases integrated dose yet is well-tolerated and an accepted treatment schedule for chemotherapy (23); however, adjusting pulse schedules to account for the “extra” lingering dose during off-periods may be required for other patients.

To apply our multiple-therapy model to clinical data, we used patient 422 as an example. This patient was consecutively treated with eribulin followed by paclitaxel; each treatment was administered until resistance was reached, akin to the sequential treatment strategy described in Fig. 2. By using predicted growth matrix information from the single-treatment predictions for both eribulin and paclitaxel, we estimated the growth matrices for a four-cell type system (Eq. (H)). We predicted the growth trajectory of the tumor under sequential treatment, where each therapy was pulsed according to standard on/off schedules as in Supplementary Fig. S6E–S6F. We found that the predicted trajectory closely follows the actual time to resistance for this patient (Fig. 4D, pink line), and that continuous sequential treatment would have shortened time to relapse (Fig. 4D, black line). We also predict that alternating between eribulin and paclitaxel at fixed intervals without any breaks in treatment (Fig. 4D, red line) or combination pulse treatment (Fig. 4D, blue line) could have extended overall time to relapse. Despite the many assumptions in this example, we were able to apply our model to predict tumor trajectory over the course of more than one treatment and predicted the best way to schedule the two therapies for relapse delay.

Here, we presented a mathematical model to study treatment strategies and showed how they may be applied to determine the delay in tumor resistance. By generating large numbers of simulation results and using machine learning to select the most relevant input features we can draw some general conclusions.

Pulsed treatment is always better than continuous dosing when considering monotherapy, and aggressive cancers are more likely to have a meaningful benefit from this strategy. Simply knowing the growth rate of the tumor in the absence of treatment and how well it initially responds to treatment could be sufficient to predict the benefit of pulse therapy. This could help stratify patients according to whether pulse strategy would likely carry a meaningful benefit, meriting more in-depth tumor characterization to identify the ideal pulse strength and duration for the patient. We show how these quantities can be calculated from clinical data, and how our model can predict the optimal pulse strategy using limited early information from the patient. Our model offers a priori proactive-based strategies, contrasting with reactive-based techniques requiring frequent monitoring and treatment adjustment in real-time (18).

Pulse therapy and adaptive therapy are complementary treatment strategies: while both have been shown to outperform canonical regimens, they could each fulfill a different treatment requirement. Adaptive therapy works well when sensitive clones outcompete resistant clones in the absence of treatment, favoring dose reduction. Pulse therapy works well when resistant clones cannot be outcompeted by other clones, and temporary dose escalation can limit the population that can acquire resistance.

Our model compares pulse and continuous strategies assuming the cumulative dose is equal between the two in the same period of time. In our stochastic model (Fig. 3), there are some cases where the total dose over the lifetime of a patient is actually less in pulse treatment, since the disease is cured. Additionally, the standard dose selected for a therapy is the average maximum tolerated dose among patients in a trial, which means that many patients are underdosed at the standard dosing level (47). Because our application of the model is at a personalized medicine level, pulses can be adjusted for each individual. Many targeted therapies such as letrozole are given as a daily oral medication and could theoretically be taken at higher doses followed by a scheduled break. In fact, a phase I trial for patients with HER2-positive breast cancer with brain metastases found that a dosing schedule of lapatinib similar to our pulse strategy was well tolerated and promising in efficacy compared to continuous treatment (48, 49). We speculate that any therapy which includes a loading dose—a short-term higher dose at the start of treatment—is likely to tolerate short-term dose escalations, since the initial high dose was deemed safe. Fulvestrant, a breast cancer drug, is one such example (50).

An important caveat is that not all therapies may realistically change their scheduling or allow even short-term escalations. Chemotherapies such as paclitaxel are given as intravenous injections every few weeks and might not be reasonably scheduled in different specific pulses due to toxicity and impracticality. It is also possible that increasing the dose intensity of some drugs will have diminishing returns (23). However, a study of metastatic prostate cancer found that higher-dose docetaxel administered every 3 weeks was slightly superior to lower-dose docetaxel administered weekly, supporting the feasibility and efficacy of pulse treatments in the future (51).

Our simulations indicated that when multiple therapies are available, alternating between them at fixed intervals is always better than sequential treatment, and pulsed-simultaneous treatment with the therapies is most likely to give the largest advantage over continuous sequential treatment. In line with previous studies, we define our “combination therapy” to be a weighted average of the two therapies such that the total dose does not exceed the dose of the more potent therapy in order to not increase systemic toxicity (3, 28, 30, 52). Consistent with prior reports (26), we find that combination therapy is usually better than sequential therapy; however, we found it is least likely to have a benefit over continuous treatment compared to other strategies and was never the best option. The alternating strategy did not employ temporary dose increases yet also proved better than the combination strategy. This suggests that targeting multiple different resistance pathways simultaneously is less favorable than targeting one resistant population at a time if doses are held constant. Our results are consistent with other findings that demonstrate the benefit of switching between therapies and that alternating therapies can be better than combination treatment (29, 53, 54); our model additionally does not require frequent monitoring of tumor growth or heterogeneity over time, and allows personalized schedules for each patient.

Lack of clustering on our PCA plot (Fig. 2B) showed that tumors with similar properties can have different optimal treatment strategies because the system is very sensitive to initial conditions. However, given that most strategies were superior to continuous treatment >95% of the time, we can assume that even if we apply the second-best strategy to a tumor it will likely still be a successful option.

We chose to use a model of exponential tumor growth for its simplicity; this may be unrealistic in situations where tumors slow their growth when they approach a carrying capacity. Nonetheless, the exponential model has previously been shown to align well with data from clinical trials (55) and general mathematical models of tumor progression (53). As justified previously (31), if we consider the modeled system to be the entire patient, the idea of a carrying capacity seems less relevant due to metastatic expansion; in our clinical examples, the patients are being treated for metastatic disease and total tumor burden is measured. Moreover, we model tumor size as it grows and declines around detection level, well below carrying capacity, and previous models of resistance found that the incorporation of carrying capacity did not alter the main results (15). Logistic or Gompertz models could be more realistic for modeling primary tumors or large metastatic lesions as they reach their own local carrying capacity (22).

Previous work used mathematical approximations in addition to simulations in order to determine endpoints such as tumor size and survival (56, 57). While effective, mathematical analysis is increasingly complex for scenarios such as pulse therapy, which is a piece-wise function with unknown number of cycles. The machine learning approaches we use are not only time-tested, but allow us to retain intuition about the biological relevance of the predictor variables. Furthermore, the features selected by the machine learning algorithm were clinically relevant parameters (both sum(Moff*n[0]) and sum(M*n[0]) can be determined from CA.15.3 data from patients, allowing us to predict with just two values which patients are most likely to benefit from pulse therapy). Therefore, unlike purely mathematical approaches, we do not necessarily require each clone's individual growth rate measurement in order to fit the model to clinical scenarios. Machine learning can often be applied to real-world data; we can envision a future where our approach can be parameterized to patient data when we have an abundance of clinical data. Previous work used decision trees to simulate treatment schedules multiple steps in advance (30). Incorporating both high-dose monotherapies and combination therapies within complex sequences of treatment schedules led to the highest cure rates in simulated tumors. Their work further exemplifies the benefit of a priori predictions and that more detailed data would allow increasingly advanced predictive models.

Interestingly, a prior report found continuous treatment with lapatinib to be more effective than pulse treatment to lower tumor burden in glioblastoma (58). However, it is unclear whether the optimal treatment strategy to reduce tumor burden in the time frame studied would also be the best strategy for delaying resistance if a durable remission is not achieved. It is also possible that the one-size-fits-all method of pulsed treatments is unlikely to be effective. Likewise, a study of cancer evolution that modeled strategies very similar to ours found that the optimal strategy could vary depending on the pharmacokinetics of the drugs used (59), implying that tumor properties alone were not sufficient to determine to the best outcome. This study also used the least number of cancer cells at given time as the metric for success, which does not always correlate with TTR. The inclusion of pharmacokinetics is still, however, an important aspect of treatment strategies and would improve our model and clinical application. A switch frequency for pulse treatments that is shorter than the half-life for the drug will not truly “pulse” the tumor, and designing the optimal pulse strategy for a patient should take pharmacokinetics into account.

The mathematical model used to simulate tumor growth and drug treatment could also be further improved by allowing interactions between cell types that go beyond acquisition of resistance (60, 61). Allowing reversible transitions in our model, such as transcriptional changes that can alter the sensitivity of a cell without the need for a permanent mutation, will expand the applicability to multiple modes of resistance. This could also further increase the benefit of breaks between therapy doses, by allowing periods of resensitization to the therapy via phenotype switching (62). Note, however, that the more detail a model has the harder it is to reverse-engineer the model's input parameters from the clinical data available. This greatly limits our ability to make useful predictions and assist in cancer therapy.

Here, we showed how we could use early tumor information from patients to predict time to relapse and how pulse strategies may extend this time to relapse. However, there are several caveats: first, small differences in starting conditions vastly alter the prediction timeline. Therefore, data that provide estimates of $n[ 0 ]$ must be very accurate to correctly reverse-engineer the growth dynamics matrix. This not only requires precise genomic and histological information, but also better knowledge of the penetrance of the presumed resistance mutations. More research on resistance mutations and their effects on the population is needed. Second, cancer antigen 15.3 (CA.15.3) is a proxy for overall tumor burden but does not necessarily give us information about the size of individual tumors/metastases. Detailed longitudinal information—from imaging, for example—will better inform us of tumor size which may follow different growth patterns than those estimated by blood markers. Third, each prediction simulation had a wide standard deviation, suggesting that a higher density of data is required to make precise predictions. Despite these significant limitations, our study indicates that it is possible to predict resistance a priori using currently available clinical data. We therefore propose that accurate predictions may soon be within our reach with a concerted effort: improved models as described above, more frequent clinical sampling, and better biological knowledge of the effects of resistance mutations will likely enable us to make the best proactive choices for patients with cancer.

H.I. Scher reports personal fees from Ambry Genetics Corporation, Bayer, Konica, Minolta, Inc., Pfizer Inc., Sun Pharmaceuticals Industries, Inc., Amgen, WCG Oncology, Elsevier, Ltd., Sidney Kimmel Cancer Center, Jefferson Health, and Asterias Biotherapeutics; grants from AIQ Pharma, Epic Sciences, Illumina, Inc., Janssen, Menarini Silicon Biosystems, ThermoFisher, Prostate Cancer Foundation, National Cancer Institute, and from the Department of Defense; nonfinancial support from ESSA Pharma Inc., Amgen, Mabvax, BioNTech, WCG Oncology, Phosplatin, and nonfinancial support from Epic Sciences outside the submitted work. P. Razavi reports grants from Grail Inc., Illumina, grants and personal fees from Novartis, grants and personal fees from Archer Dx, grants and personal fees from AstraZeneca, personal fees from Foundation Medicine and Natera, grants and personal fees from Epic Sciences, grants and personal fees from Tempus Labs, and grants from Guardant outside the submitted work. J. Xavier reports grants from the NIH during the conduct of the study. No disclosures were reported by the other authors.

D. Mathur: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. B.P. Taylor: Software, methodology, writing–review and editing. W.K. Chatila: Resources, data curation. H.I. Scher: Supervision, project administration, writing–review and editing. N. Schultz: Supervision. P. Razavi: Resources, data curation, supervision, writing–review and editing. J.B. Xavier: Conceptualization, resources, supervision, funding acquisition, methodology, project administration, writing–review and editing.

This work was supported by NIH Research Program Grants under award numbers R01CA229215 and U54 CA209975 (to J. Xavier). D. Mathur is also supported by the Alan and Sandra Gerry Metastasis and Tumor Ecosystems Center (GMTEC) Gerry Fellowship and the 1F32CA260841-01 NRSA award from the NCI. This work was funded in part by the Marie-Josée and Henry R. Kravis Center for Molecular Oncology and the NCI Cancer Center Core grant No. P30-CA008748. We gratefully acknowledge the members of the Molecular Diagnostics Service in the Department of Pathology.

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.
Zahreddine
H
,
Borden
KLB
.
Mechanisms and insights into drug resistance in cancer
.
Front Pharmacol
2013
;
4
:
28
.
2.
Pernas
S
,
Tolaney
SM
.
HER2-positive breast cancer: new therapeutic frontiers and overcoming resistance
.
Ther Adv Med Oncol
2019
;
11
:
1758835919833519
.
3.
Wang
Z
,
Deisboeck
TS
.
Dynamic targeting in cancer treatment
.
Front Physiol
2019
;
10
:
96
.
4.
Nedeljković
M
,
Damjanović
A
.
Mechanisms of chemotherapy resistance in triple-negative breast cancer-how we can rise to the challenge
.
Cells
2019
;
8
:
957
.
5.
Van Norman
GA
.
Drugs, devices, and the FDA: part 1: an overview of approval processes for drugs
.
JACC Basic Transl Sci
2016
;
1
:
170
9
.
6.
Chowell
D
,
Napier
J
,
Gupta
R
,
Anderson
KS
,
Maley
CC
,
Sayres
MAW
.
Modeling the subclonal evolution of cancer cell populations
.
Cancer Res
2018
;
78
:
830
9
.
7.
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
.
8.
Inukai
M
,
Toyooka
S
,
Ito
S
,
Asano
H
,
Ichihara
S
,
Soh
J
, et al
.
Presence of epidermal growth factor receptor gene T790M mutation as a minor clone in non-small cell lung cancer
.
Cancer Res
2006
;
66
:
7854
8
.
9.
Kansal
AR
,
Torquato
S
,
Chiocca
EA
,
Deisboeck
TS
.
Emergence of a subpopulation in a computational model of tumor growth
.
J Theor Biol
2000
;
207
:
431
41
.
10.
Sun
X
,
Bao
J
,
Shao
Y
.
Mathematical modeling of therapy-induced cancer drug resistance: connecting cancer mechanisms to population survival rates
.
Sci Rep
2016
;
6
:
22498
.
11.
Almendro
V
,
Cheng
Y-K
,
Randles
A
,
Itzkovitz
S
,
Marusyk
A
,
Ametller
E
, et al
.
Inference of tumor evolution during chemotherapy by computational modeling and in situ analysis of genetic and phenotypic cellular diversity
.
Cell Rep
2014
;
6
:
514
27
.
12.
Iwasa
Y
,
Nowak
MA
,
Michor
F
.
Evolution of resistance during clonal expansion
.
Genetics
2006
;
172
:
2557
66
.
13.
Michor
F
,
Beal
K
.
Improving cancer treatment via mathematical modeling: surmounting the challenges is worth the effort
.
Cell
2015
;
163
:
1059
63
.
14.
Gatenby
RA
,
Silva
AS
,
Gillies
RJ
,
Frieden
BR
.
Adaptive therapy
.
Cancer Res
2009
;
69
:
4894
903
.
15.
Zhang
J
,
Cunningham
JJ
,
Brown
JS
,
Gatenby
RA
.
Integrating evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer
.
Nat Commun
2017
;
8
:
1816
.
16.
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
.
17.
Silva
AS
,
Kam
Y
,
Khin
ZP
,
Minton
SE
,
Gillies
RJ
,
Gatenby
RA
.
Evolutionary approaches to prolong progression-free survival in breast cancer
.
Cancer Res
2012
;
72
:
6362
70
.
18.
Gallaher
JA
,
Enriquez-Navas
PM
,
Luddy
KA
,
Gatenby
RA
,
Anderson
ARA
.
Spatial heterogeneity and evolutionary dynamics modulate time to recurrence in continuous and adaptive cancer therapies
.
Cancer Res
2018
;
78
:
2127
39
.
19.
Greene
JM
,
Gevertz
JL
,
Sontag
ED
.
Mathematical approach to differentiate spontaneous and induced evolution to drug resistance during cancer treatment
.
JCO Clin Cancer Inform
2019
;
3
:
1
20
.
20.
Bozic
I
,
Nowak
MA
.
Resisting resistance
.
Annu Rev Cancer Biol
2017
;
1
:
203
21
.
21.
Le
X
,
Antony
R
,
Razavi
P
,
Treacy
DJ
,
Luo
F
,
Ghandi
M
, et al
.
Systematic functional characterization of resistance to PI3K inhibition in breast cancer
.
Cancer Discov
2016
;
6
:
1134
47
.
22.
Grassberger
C
,
McClatchy
D
,
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
.
23.
Simon
R
,
Norton
L
.
The Norton-Simon hypothesis: designing more effective and less toxic chemotherapeutic regimens
.
Nat Clin Pract Oncol
2006
;
3
:
406
7
.
24.
Citron
ML
,
Berry
DA
,
Cirrincione
C
,
Hudis
C
,
Winer
EP
,
Gradishar
WJ
, et al
.
Randomized trial of dose-dense versus conventionally scheduled and sequential versus concurrent combination chemotherapy as postoperative adjuvant treatment of node-positive primary breast cancer: first report of Intergroup Trial C9741/Cancer and Leukemia Group B Trial 9741
.
J Clin Oncol
2003
;
21
:
1431
9
.
25.
Norton
L
.
Theoretical concepts and the emerging role of taxanes in adjuvant therapy
.
Oncologist
2001
;
6
:
30
5
.
26.
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
.
27.
Tomasetti
C
,
Levy
D
.
An elementary approach to modeling drug resistance in cancer
.
Math Biosci Eng
2010
;
7
:
905
18
.
28.
Beckman
RA
,
Schemmann
GS
,
Yeang
C-H
.
Impact of genetic dynamics and single-cell heterogeneity on development of nonstandard personalized medicine strategies for cancer
.
Proc Natl Acad Sci U S A
2012
;
109
:
14586
91
.
29.
Goldie
JH
,
Coldman
AJ
,
Gudauskas
GA
.
Rationale for the use of alternating non-cross-resistant chemotherapy
.
Cancer Treat Rep
1982
;
66
:
439
49
.
30.
Yeang
C-H
,
Beckman
RA
.
Long range personalized cancer treatment strategies incorporating evolutionary dynamics
.
Biol Direct
2016
;
11
:
56
.
31.
Beckman
RA
,
Kareva
I
,
Adler
FR
.
How should cancer models be constructed?
Cancer Control
2020
;
27
:
1073274820962008
.
32.
Razavi
P
,
Chang
MT
,
Xu
G
,
Bandlamudi
C
,
Ross
DS
,
Vasan
N
, et al
.
The genomic landscape of endocrine-resistant advanced breast cancers
.
Cancer Cell
2018
;
34
:
427
38
.
33.
Citron
ML
.
Dose-dense chemotherapy: principles, clinical results and future perspectives
.
Breast Care
2008
;
3
:
251
5
.
34.
Harris
L
,
Fritsche
H
,
Mennel
R
,
Norton
L
,
Ravdin
P
,
Taube
S
, et al
.
American Society of Clinical Oncology 2007 update of recommendations for the use of tumor markers in breast cancer
.
J Clin Oncol
2007
;
25
:
5287
312
.
35.
Cervino
AR
,
Saibene
T
,
Michieletto
S
,
Ghiotto
C
,
Bozza
F
,
Saladini
G
, et al
.
Correlation between cancer antigen 15.3 value and qualitative and semiquantitative parameters of positron emission tomography/computed tomography in breast cancer patients
.
Curr Radiopharm
2014
;
7
:
20
8
.
36.
López
AM
,
Pruthi
S
,
Boughey
JC
,
Perloff
M
,
Hsu
C-H
,
Lang
JE
, et al
.
Double-blind, randomized trial of alternative letrozole dosing regimens in postmenopausal women with increased breast cancer risk
.
Cancer Prev Res
2016
;
9
:
142
8
.
37.
Hortobagyi
GN
.
Everolimus plus exemestane for the treatment of advanced breast cancer: a review of subanalyses from BOLERO-2
.
Neoplasia
2015
;
17
:
279
88
.
38.
Lamb
HM
,
Adkins
JC
.
Letrozole. A review of its use in postmenopausal women with advanced breast cancer
.
Drugs
1998
;
56
:
1125
40
.
39.
Parasido
E
,
Avetian
GS
,
Naeem
A
,
Graham
G
,
Pishvaian
M
,
Glasgow
E
, et al
.
The sustained induction of c-MYC drives nab-paclitaxel resistance in primary pancreatic ductal carcinoma cells
.
Mol Cancer Res
2019
;
17
:
1815
27
.
40.
Gupta
N
,
Gupta
P
,
Srivastava
SK
.
Penfluridol overcomes paclitaxel resistance in metastatic breast cancer
.
Sci Rep
2019
;
9
:
5066
.
41.
Reinert
T
,
Saad
ED
,
Barrios
CH
,
Bines
J
.
Clinical implications of ESR1 mutations in hormone receptor-positive advanced breast cancer
.
Front Oncol
2017
;
7
:
26
.
42.
Yang
XL
,
Lin
FJ
,
Guo
YJ
,
Shao
ZM
,
Ou
ZL
.
Gemcitabine resistance in breast cancer cells regulated by PI3K/AKT-mediated cellular proliferation exerts negative feedback via the MEK/MAPK and mTOR pathways
.
Onco Targets Ther
2014
;
7
:
1033
42
.
43.
Eckstein
N
.
Platinum resistance in breast and ovarian cancer cell lines
.
J Exp Clin Cancer Res
2011
;
30
:
91
.
44.
Hayasaka
N
,
Takada
K
,
Nakamura
H
,
Arihara
Y
,
Kawano
Y
,
Osuga
T
, et al
.
Combination of eribulin plus AKT inhibitor evokes synergistic cytotoxicity in soft tissue sarcoma cells
.
Sci Rep
2019
;
9
:
5759
.
45.
Liu
L
,
Meng
T
,
Zheng
X
,
Liu
Y
,
Hao
R
,
Yan
Y
, et al
.
Transgelin 2 promotes paclitaxel resistance, migration, and invasion of breast cancer by directly interacting with PTEN and activating PI3K/Akt/GSK-3β pathway
.
Mol Cancer Ther
2019
;
18
:
2457
68
.
46.
Buzdar
AU
,
Robertson
JFR
,
Eiermann
W
,
Nabholtz
J-M
.
An overview of the pharmacology and pharmacokinetics of the newer generation aromatase inhibitors anastrozole, letrozole, and exemestane
.
Cancer
2002
;
95
:
2006
16
.
47.
Gurney
H
.
How to calculate the dose of chemotherapy
.
Br J Cancer
2002
;
86
:
1297
302
.
48.
Morikawa
A
,
de Stanchina
E
,
Pentsova
E
,
Kemeny
MM
,
Li
BT
,
Tang
K
, et al
.
Phase I study of intermittent high-dose lapatinib alternating with capecitabine for HER2-positive breast cancer patients with central nervous system metastases
.
Clin Cancer Res
2019
;
25
:
3784
92
.
49.
Morikawa
A
,
De Stanchina
E
,
Patil
S
,
Chandarlapaty
S
,
Li
BT
,
Norton
L
, et al
.
Abstract P4-14-24: Optimization of intermittent high dose lapatinib administration with or without capecitabine: a rational approach to drug dosing and scheduling using Norton-Simon modeling
.
Cancer Res
2016
;
76
:
P4–14–24–P4–14–24
.
50.
McCormack
P
,
Sapunar
F
.
Pharmacokinetic profile of the fulvestrant loading dose regimen in postmenopausal women with hormone receptor-positive advanced breast cancer
.
Clin Breast Cancer
2008
;
8
:
347
51
.
51.
Berthold
DR
,
Pond
GR
,
Soban
F
,
de Wit
R
,
Eisenberger
M
,
Tannock
IF
.
Docetaxel plus prednisone or mitoxantrone plus prednisone for advanced prostate cancer: updated survival in the TAX 327 study
.
J Clin Oncol
2008
;
26
:
242
5
.
52.
Zhao
B
,
Pritchard
JR
,
Lauffenburger
DA
,
Hemann
MT
.
Addressing genetic tumor heterogeneity through computationally predictive combination therapy
.
Cancer Discov
2014
;
4
:
166
74
.
53.
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
.
54.
Jonsson
VD
,
Blakely
CM
,
Lin
L
,
Asthana
S
,
Matni
N
,
Olivas
V
, et al
.
Novel computational method for predicting polytherapy switching strategies to overcome tumor heterogeneity and evolution
.
Sci Rep
2017
;
7
:
44206
.
55.
Yin
A
,
Moes
DJAR
,
van Hasselt
JGC
,
Swen
JJ
,
Guchelaar
H-J
.
A review of mathematical models for tumor dynamics and treatment resistance evolution of solid tumors
.
CPT Pharmacometrics Syst Pharmacol
2019
;
8
:
720
37
.
56.
Foo
J
,
Leder
K
.
Dynamics of cancer recurrence
.
Ann Appl Probab
2013
;
23
.
57.
Yamamoto
KN
,
Liu
LL
,
Nakamura
A
,
Haeno
H
,
Michor
F
.
Stochastic evolution of pancreatic cancer metastases during logistic clonal expansion
.
JCO Clin Cancer Inform
2019
;
3
:
1
11
.
58.
Stein
S
,
Zhao
R
,
Haeno
H
,
Vivanco
I
,
Michor
F
.
Mathematical modeling identifies optimum lapatinib dosing schedules for the treatment of glioblastoma patients
.
PLoS Comput Biol
2018
;
14
:
e1005924
.
59.
Chakrabarti
S
,
Michor
F
.
Pharmacokinetics and drug interactions determine optimum combination strategies in computational models of cancer evolution
.
Cancer Res
2017
;
77
:
3908
21
.
60.
Nowak
MA
,
Sigmund
K
.
Evolutionary dynamics of biological games
.
Science
2004
;
303
:
793
9
.
61.
Park
DS
,
Robertson-Tessi
M
,
Luddy
KA
,
Maini
PK
,
Bonsall
MB
,
Gatenby
RA
, et al
.
The goldilocks window of personalized chemotherapy: getting the immune response just right
.
Cancer Res
2019
;
79
:
5302
15
.
62.
Akhmetzhanov
AR
,
Kim
JW
,
Sullivan
R
,
Beckman
RA
,
Tamayo
P
,
Yeang
C-H
.
Modelling bistable tumour population dynamics to design effective treatment strategies
.
J Theor Biol
2019
;
474
:
88
102
.
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs International 4.0 License.

Supplementary data