## Abstract

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.

## Introduction

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.

## Materials and Methods

### 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:

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:

Where

We can therefore rewrite equation (A) as

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 g_{1} is positive in the absence of treatment (g_{1off}), and negative in the presence of treatment (g_{1on}). We assume that g_{2} will be unaffected by therapy; therefore, g_{2} 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 (M_{off}), (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:

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:

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) n_{1}: cells sensitive to both *vA* and *vB*, n_{2}: resistant to only *vA*, n_{3}: resistant to only *vB*, or n_{4}: 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:

Where p_{R}i,j is the transition rate from cell type i to cell type j: p_{mi}*p_{r}i,j. Inputs into the model are: (1) matrix of parameters off-treatment (M_{off}), (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.

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 M_{off} 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.

## Results

### 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.

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(M_{off}*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(M_{off}*n[0]), and higher sum(M*$n[ 0 ]$) lead to a greater pulse advantage (Fig. 1,F–H). 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 *R*^{2} 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(M_{off}*n[0]) and sum(M*$n[ 0 ]$) as predictors, and the model still accurately predicted pulse advantage with an *R*^{2} 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(M_{off}*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.

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).

. | Better than continuous strategy . | Best option . | Better than continuous strategy . | Best 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 strategy . | Best option . | Better than continuous strategy . | Best 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(M_{off}*$\ 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.

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).

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:

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).

Patient . | Treatment . | Resistance mechanism . | Data type . | Estimated 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% |

Patient . | Treatment . | Resistance mechanism . | Data type . | Estimated 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 M_{off} 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.

## Discussion

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.

## Authors' Disclosures

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.

## Authors' Contributions

**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.

## Acknowledgments

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.