Preventing evolutionary rescue in cancer

First-line cancer treatment frequently fails due to initially rare therapeutic resistance. An important clinical question is then how to schedule subsequent treatments to maximize the probability of tumour eradication. Here, we provide a theoretical solution to this problem by using mathematical analysis and extensive stochastic simulations within the framework of evolutionary rescue theory to determine how best to exploit the vulnerability of small tumours to stochastic extinction. Whereas standard clinical practice is to wait for evidence of relapse, we confirm a recent hypothesis that the optimal time to switch to a second treatment is when the tumour is close to its minimum size before relapse, when it is likely undetectable. This optimum can lie slightly before or slightly after the nadir, depending on tumour parameters. Given that this exact time point may be difficult to determine in practice, we study windows of high extinction probability that lie around the optimal switching point, showing that switching after the relapse has begun is typically better than switching too early. We further reveal how treatment dose and tumour demographic and evolutionary parameters influence the predicted clinical outcome, and we determine how best to schedule drugs of unequal efficacy. Our work establishes a foundation for further experimental and clinical investigation of this evolutionarily-informed “extinction therapy” strategy.


Introduction
Just as species in an ecosystem interact, compete for resources, adapt to changing environmental conditions and undergo natural selection, so cancer clones rise and fall in a tumour ecosystem.Darwinian principles inevitably determine therapeutic responses [1] including the emergence of resistance, which, despite pharmaceutical advances, remains the greatest challenge in oncology.As cancer cells can use a variety of adaptive strategies to achieve resistance [2], targeting a single molecular mechanism often proves ineffective in the long term [3].Understanding intratumour evolutionary processes provides a rational foundation for developing treatment strategies that, by explicitly accounting for evolutionary dynamics, achieve better clinical outcomes [4,5,6].Mathematical modelling of clonal dynamics and the emergence of resistance is critical for optimising clinical treatment strategies based on evolutionary principles.Consequently, the historical development of evolutionary therapies has followed a trajectory that begins with a theoretical and mathematical exploration of associated eco-evolutionary models [7,8,9].
The clinical strategy we study here uses evolutionary rescue theory to inform the probability of tumour extinction under multiple treatment administrations or "strikes".Although it is more usual to consider evolutionary rescue in a conservation context, the same theory is applicable when extinction is the goal, such as Preventing evolutionary rescue in cancer in bacterial infections or cancer [10].Since an oncologist can influence the tumour environment, they can anticipate the evolutionary trajectories of cancer clones and, in theory, follow a strategy to avoid evolutionary rescue and so cure the patient [11].The key idea is that, even if a single strike fails to eradicate cancer cells due to resistant phenotypes, it can still render the population small and fragmented.Small populations are more vulnerable to stochastic extinction and less capable of adapting to environmental changes owing to loss of phenotypic heterogeneity [10].Cell proliferation may also slow due to Allee effects [12].Subsequent therapeutic strikes, if well timed, can exploit these weaknesses to drive the cancer cell population to extinction [13].
Combination cancer therapies are typically designed such that cells resistant to one treatment are likely to have collateral sensitivity to another [14].The main differences between multi-strike therapy and conventional combination therapy are in the timing of the strikes and the use of evolutionary principles to guide treatment.In combination or sequential therapy, the second or subsequent treatments are usually given during relapse, when the first treatment appears to have failed.Another conventional strategy is to simultaneously administer multiple drugs with collateral sensitivities from the beginning of treatment [15].In multi-strike therapy, the idea is instead to attack the cancer at its weakest point when it may well be clinically undetectable.It has been suggested in a prior study that the best time to give the second strike may be while the tumour is still shrinking in response to the first therapy [16].The success rate of multi-strike therapy is expected to be highly sensitive to the timing and severity of the second and any subsequent strikes.
Demonstrating its potential to improve cure rates across diverse cancer types, three clinical trials of multistrike therapy are already underway.A Phase 2 trial using conventional chemotherapy drugs in metastatic rhabdomyosarcoma started recruiting patients in 2020 and is expected to run until 2026 [17].A Phase 1 trial in metastatic prostate cancer (2022-27) involves agents that exploit the hormone sensitivity of cancer cells [18].A Phase 2 trial using targeted therapies in metastatic breast cancer began in 2024 [19].Further trials are in development.
Yet, despite this rapid progress to clinical evaluation, many critical questions regarding the timing of the subsequent strikes, the time until extinction, the effect of environmental and demographic factors, and most importantly, the conditions under which multi-strike therapy is a feasible alternative to other therapies remain unanswered.How effective is the first strike, and does it make the population vulnerable enough for further strikes to work?What is the probability that a population is rescued either by pre-existing mutants or those that arise during the treatment?How do outcomes vary with the cost of resistance, density dependence, and other factors that affect clonal growth rates?
We tackle these pressing questions in two ways.First, using ideas from evolutionary rescue theory, we develop and study the first analytical model of two-strike therapy.This simple, tractable mathematical model enables us to compute extinction probabilities and to identify the optimal time for the second strike.Second, we use extensive stochastic simulations to test the robustness of our analytical results and to study the effects of additional factors.We thus establish a necessary foundation for further theoretical, experimental, and clinical investigations of multi-strike therapy.

Methods
To obtain general, robust insights into the factors that determine a successful two-strike treatment strategy, we combine a deterministic analytical model and a stochastic simulation model.Both models involve two stressful environments E 1 and E 2 (corresponding to the two treatment strikes) and four cell types.Cells can be sensitive to both treatments (S), resistant to one treatment but sensitive to the other (R 1 and R 2 ), or resistant to both treatments (R 1,2 ).The time of switching to the second treatment is τ and the population size at this time is N (τ ).

Analytical methods
Our analytical modelling method is composed of two stages (Figure 1).First, we model the population dynamics during the first treatment as a set of numerically-solved differential equations.We then use those solutions to predict extinction probability using evolutionary rescue theory.
To calculate extinction probabilities due to standing genetic variation (pre-existing rescue mutants) and denovo rescue mutants, we must obtain the population composition at the beginning of the second strike.For this we use the system of differential equations given in Figure 1 (first box), describing logistic growth Preventing evolutionary rescue in cancer in environment E 1 of the four subpopulations S(t), R 1 (t), R 2 (t) and R 1,2 (t) that make up the tumour cell population N (t).The model includes mutations from less resistant to more resistant states while, for simplicity, ignoring negligible back mutations.For plausible parameter values, a tumour that grows from a single treatment-sensitive (S) cell is unlikely to harbour any doubly-resistant (R 1,2 ) cells at the time it is first treated (see Appendix A.1.1).If this were not so then a cure would be highly improbable.We therefore assume R 1,2 (0) = 0. Other default initial conditions and parameter values are listed in Table 1.By solving the differential equations numerically over the course of the first treatment (time 0 to time τ ), we determine the subpopulation sizes at the time of switching to the second treatment.2) is denoted by µ1(µ2).Growth rates gS for sensitive cells and gR for resistant cells depend on the intrinsic birth rate, intrinsic death rate and the cost of resistance (see Table 1).The treatment-induced death rate is denoted by δ1.Initial conditions are specified by the initial population sizes of S, R1, R2 and R1,2 cells.The total initial population N (0) is the sum of these four subpopulations.
Given the population composition at treatment switching time τ , we next compute the probabilities P SGV E (τ ) and P DN E (τ ) of no evolutionary rescue due to standing genetic variation and de-novo mutations, respectively (Part 2 in Figure 1; Appendix A.3).Because successful treatment requires the eradication of both preexisting and de-novo mutants during the second treatment period, E 2 , the tumour extinction probability P E (τ ) is the product of these two probabilities: where π e is the probability of establishment of a single resistant lineage, which depends on parameters b, d and c (see Appendix A.5 for the derivation).With Eq 2, we study the behaviour of extinction probability as a function of τ under different conditions.Per capita death rate due to treatment 1,2 2.0 S(0) Initial population of S cells

Stochastic simulations
To test the robustness of our analytical results, we separately obtain extinction probabilities using a stochastic simulation model with the same initial conditions as the ODE system and equivalent default parameter values (see Appendix A.6).The main difference between the models is that the analytical method uses evolutionary rescue theory to calculate extinction probabilities, whereas the computational approach uses the stochastic Gillespie algorithm to simulate birth, death, and mutation events.Each simulation ends with one of three outcomes: extinction, progression, or persistence (see Table A.1).The extinction probability is estimated as the proportion of extinction outcomes in a large number of simulations.

Comparing results across parameter values
To compare treatment outcomes for varied parameter values, we use a summary variable N q to describe how small the tumour must be at the time of switching treatment to achieve a given probability of extinction.This concept is based on our observation that the extinction probability P E (τ ) generally decreases as N (τ ) increases, unless N (τ ) is very close to the population nadir that would pertain in the absence of a second strike (N min ).Therefore, for a given extinction probability q (with 0 ≤ q ≤ max(P E (τ ))), we can obtain a corresponding value N q , which is the maximum population size threshold below which we achieve an extinction probability greater than or equal to q.In other words, if N (τ )≤ N q , then we will achieve an extinction probability of at least q.Any given switching point N (τ ) is reached twice in the trajectory of a population undergoing evolutionary rescue, once before and once after the start of relapse.N q is therefore defined for both before-nadir and after-nadir switching points: where t(N min ) is the time at which the nadir would be reached in the absence of a second strike.
N q values tell us when and how fast the extinction probability drops from a high value to a low one.For instance, if the difference between N (before) 0.1 and N (before) 0.9 is slight, then the extinction probability must increase steeply between these two population sizes.A higher value of N (before) 0.9 implies a wider window of opportunity for implementing a successful second strike.More generally, we want the range of N (τ ) values with low P E (τ ) to be as small as possible.We plot N q versus q to analyse the trend in extinction probabilities across the range of potential switching points.Several metrics for comparing parameter sets based on N q are described in Appendix A.8.

Results
To enable us to uncover general principles and determine the most important factors in a successful treatment strategy, we consider the simplest case of two strikes.Since further strikes can only increase extinction probabilities, we thus obtain conservative lower bounds on potential clinical benefits.The first treatment (or Preventing evolutionary rescue in cancer strike) creates a stressful environment that we denote E 1 .After switching to the second treatment, the tumour enters the second stressful environment, E 2 .
Corresponding to the two treatments, we consider four cell types -sensitive to both treatments (S cells), resistant to one of the treatments but sensitive to the other (R 1 and R 2 ) and resistant to both treatments (R 1,2 ).Consequently, R 1 and R 1,2 cells are resistant in E 1 , and R 2 and R 1,2 cells are resistant in E 2 .Even if they initially rescue the population, all R 1 cells will eventually go extinct due to the second strike.Evolutionary rescue from the second strike must be due to surviving R 2 or R 1,2 cells, called rescue mutants.
We obtain extinction probabilities using both a deterministic analytical model (Section 2.1) and a stochastic simulation model (Section 2.2), and we compare the two wherever possible.Unless mentioned otherwise, we use a default set of parameters and initial conditions (Table 1).Except in Section 3.4, we assume the two treatments induce identical death rates (that is, δ 1 =δ 2 =δ).For brevity, we will use "dose" and "treatment level" to refer to these treatment-induced death rates, which in reality also depend on pharmacodynamics and pharmacokinetics.
Our focus will be on the population size at the time of switching between the two treatments, denoted N (τ ).Since the optimal N (τ ) changes as we vary the parameters, the trend of extinction probabilities obtained at a fixed N (τ ) could differ from the trend obtained at the optimal N (τ ).The rationale for using a fixed N (τ ) for such comparisons is that it may, in practice, be impossible to determine the optimal N (τ ), which requires knowing the values of all the system parameters.
The fixed N (τ ) can be implemented on either side of the population nadir since a given population threshold is met twice in the characteristic U-shaped trajectory of a population undergoing evolutionary rescue (Figure 2(B)).Therefore we consider both before-nadir and after-nadir switching points.The population nadir reached in the absence of a second strike, which we denote N min , can be calculated by numerically solving the system of differential equations shown in Figure 1 (see Appendix A.1.2.1 for an analytical approximation).
To compare different parameter values and treatment conditions, we measure the regions of high extinction probability, defined as the range of N (τ ) values that give an extinction probability ≥ 0.8.Additionally, we use comparison metrics based on the quantity N q (Section 2.3 and Appendix A.8), which is defined as the maximum population size threshold that must be crossed to achieve an extinction probability greater than or equal to q.

The optimal switching time is when the population size is close to its nadir
Our first aim is to find the optimal population size N (τ ), at which to switch from the first to the second treatment.Our analytical and stochastic models both show that the optimal N (τ ), in terms of maximizing extinction probability, is close to the population nadir N min (Figure 2).According to the analytical model, the optimal switching point may, depending on parameter values, lie slightly before or after N min (Appendix A.1.2.4.Yet the difference between the optimal N (τ ) and the N min is generally so small that it is not captured by our simulation results.The difference is significant only when R 2 cells are initially abundant, and the cost of resistance is low (Figure 2, second column).
To explain why the optimal N (τ ) is close to N min , we refer to Eq 2 and see that the maximum P E (τ ) will be achieved by minimizing the sum of all the rates of generating rescue mutants.The last two integral terms are minimized slightly after N min (see Appendix A.1.2.2 for an approximate analytical expression).However, in the terms constituting P SGV E (pre-existing rescue mutants at the beginning of the second strike), the decay in the R 2 population relative to the generation of R 1,2 mutants determines where the optimal N (τ ) lies relative to N min .For further analysis, we focus on regions of high extinction probability (P E ≥ 0.8) instead of the exact optima.Consistent with our analytical predictions, we observe that these high-P E regions lie around N min (Figures A.3 and A.2).

It is better to implement the second strike after the nadir than before
Given the practical impossibility of treating at the exact optimal time, we next compare outcomes for treating earlier or later (see Appendix A.6 for further details of the simulation algorithm).While the optimal switching point may lie slightly before or after the nadir, we observe that switching points after N min usually have higher extinction probabilities than those before N min (Figure 2(A), blue versus black).The default parameter set is given in Table 1.Extinction probabilities from simulations are computed using the outcomes of 100 independent runs, each with several switching points.Switching sizes smaller than Nmin are equivalent to the absence of a second treatment, and result in extinction probabilities equal to zero (black points to the left of red dashed lines).Error bars show 95% binomial proportion confidence intervals.In the first column, labels (1) and (2) correspond to curves in Figure 3(C).(B) An illustration of switching points before and after Nmin, implemented with the same random seed.See Appendix A.6 for a description of the algorithm for these simulations.(C) Nq vs q plots for four parameter sets and different treatment levels.Black curves show before-nadir switching points, and blue curves show after-nadir switching points.Labels indicate treatment levels δ1 = δ2.The same line style (solid, dashed, dotted, mixed) is used for a given treatment level for both before and after nadir curves.
We hypothesize that this result is due the pre-existing (or initially accumulated) R 2 population.If we delay switching, then this resistant subpopulation has longer to decay, which results in a smaller rescue population during the second treatment (Figure 4(A), first row, second panel).On the other hand, there is more time for doubly-resistant R 1,2 mutants to accumulate.In most cases of interest, the generation of R 1,2 mutants is slower than the decay of the R 2 population (see Section 3.3 for exceptions), and so the window of opportunity for effective treatment extends further to the right of the nadir than to the left.

The cost of resistance and the initial R 2 population size modulate extinction probabilities
Extinction probabilities depend on the number of R 2 cells present at the time of switching, which in turn depends on both the cost of resistance and the number of R 2 cells at the start of the first treatment.By removing one or both of the latter two factors, we can better understand how they affect treatment outcomes.

Preventing evolutionary rescue in cancer
A cost of resistance is expected to hasten the decay of R 2 mutations during the first treatment phase and so make the second treatment more effective.Accordingly, in most cases removing the cost of resistance reduces extinction probabilities (Figure 2(A,C), second column; Figure A.2; Figure 4(A), first row, first panel).
The exception is that in the case of high-dose treatment (δ = 2), extinction probabilities for late switching points (well after the nadir) can be slightly higher in the absence than in the presence of a resistance cost (but the optimal extinction probability is high even for severe resistance costs).The reason for this counterintuitive result is that in the absence of a cost of resistance, the first treatment phase is shorter, which gives less time for the generation of new R 1,2 mutants (see Appendix A.1.2.6 for a more detailed explanation).If we fix the switching time instead of the switching size, we see that a higher resistance cost is always beneficial (Figure 5, first and second columns).For intermediate-dose treatment (δ = 1.2), the main effect of removing the cost of resistance is to increase N min , which makes it impossible to achieve high rates of extinction (Figure 2(A), second row, second column).
Assuming that R 2 cells are absent at the start of the first treatment substantially increases extinction probabilities for before-nadir switching (Figure 2(A,C), third and fourth columns; Figure A.2).The result of switching before the nadir is then similar to -or, in the case of high-dose treatment, even slightly better than -the result of switching at the same tumour size after the nadir.
Although our analytical predictions are generally very close to our simulation results, they underestimate the probability of extinction in the case of intermediate treatment dose and no cost of resistance (Figure 2(A), second row, second and fourth columns).In this case, our modeling assumption of a Poisson-distributed R 2 population breaks down (see Sections A.3 and A.4 for details).

Higher doses do not necessarily maximize cure rates
We next relax our assumption of equal doses during the two treatment phases by examining alternative combinations of δ 1 (first treatment dose) and δ 2 (second treatment dose).In each case, we consider both the optimal switching point and the range of N (τ ) values that lead to a high extinction probability (P E ≥ 0.8).
We will refer to the default doses 2 and 1.2 as high and intermediate, respectively.
We first consider the before-nadir regime.The intuitive prediction is that higher doses should lead to larger high-P E regions.This is what we observe in the case of no resistance cost and R 2 (0) = 0 (Figure 3(A, C)).However, for our default parameter values, higher values of δ 1 give smaller high-P E regions (Figure 3(A), first row, second panel).Similarly, in the before-nadir regime, a lower dose results in higher extinction probabilities for all switching points not very close to the nadir (Figure 3(B,C) and 2(A)).A normalised N q versus q plot for four δ 1 -δ 2 combinations confirms that intermediate δ 1 and high δ 2 produce the best treatment outcome in terms of high-P E regions, because it gives a higher extinction probability at the same N (τ ) (Figure 3(C)).When the two doses are equal (along the black line in Figure 3(A)), we observe a similar trend (Figure A.7). Thus, for our default parameter values, an intermediate first dose paired with a high second dose gives the largest window of opportunity when switching before the nadir.
The somewhat counter-intuitive result is explained by the interaction of the dose, the cost of resistance, and the R 2 population.For lower δ 1 , the S cells decay more slowly, so the switching point N (τ ) is reached later.This provides more time for the R 2 population to decay but also more time for R 1,2 mutants to arise.When the cost of resistance and the initial R 2 population size are both large, the benefit of a lower δ 1 outweighs the disadvantages (see Appendix A.1.2.5 for a formal explanation).We do not observe this effect when the cost of resistance or the initial R 2 population size is set to zero (Figure A.4, Figure 3(A)).Note that if the dose is too low then the population size will never become small enough to permit stochastic extinction.
Although an intermediate δ 1 gives a larger high-P E region, the optimal extinction probability in all cases is obtained when both treatment levels are high (white contours in Figure 3(A)).Therefore, it is important to define what we need for a good treatment outcome.The best dose combination should not only lead to a high extinction probability at the optimal switching point, but it should also offer a large window of opportunity.An intermediate δ 1 allows a larger high-P E region than a high δ 1 at the cost of compromising on the optimal extinction probability (0.99 compared to 0.94).Moreover, the absolute values of N (τ ) giving a high extinction probability are also high with an intermediate δ 1 .This is because the N min and the high-P E region, in this case, are relatively large.
Note that this result depends on the fact that we compare high-P E regions based on switching sizes.The conclusion would be different if we were thinking in terms of switching times.A larger value of δ 1 is expected to lead to a larger high-P E time interval (Appendix A.1.2.5).Thus, the best dose in the before-nadir regime depends on how the therapy is implemented.
Outcomes for the after-nadir regime are best when both doses are high (Figure 3(A), first row, second panel).For our default parameter values, as expected, the high-P E region is also much larger for the after-nadir than for the before-nadir regime.When we eliminate the cost of resistance and the initial R 2 population, the optimal treatment combinations in the before-nadir and after-nadir regimes are similar (Figure 3(A,C)).

Two-strike therapy is feasible only in small tumours
Using the analytical model, we compare different values of N q (not normalised) for different initial population sizes N (0), bearing in mind that the resistant population size scales with N (0).We observe that the absolute values of N q for q close to 1 do not vary by more than an order of magnitude when N (0) ranges over three orders of magnitude, from 10 4 to 10 7 cells (Figure 4(A), second row, second panel).This implies that, within this range of initial tumour sizes, a high extinction probability can be achieved by applying the second strike at a sufficiently small population size (determined by the treatment dose, growth rates, and other parameters).Nevertheless, if N (0) is larger than 10 8 cells then the extinction probability never exceeds 0.4 (Figure 4(A), second row, second panel).There is, therefore, a limit on the size of tumours for which two-strike therapy is likely to succeed.(blue) values.The only plot which is not normalised by the initial population is the one showing variation in PE with N (0) (second row, second panel), assuming a constant proportion of initial resistant cells.The title of each plot indicates the parameter or initial condition that varies across curves.Solid curves correspond to the default parameter values (Table 1), and the same line style (solid, dashed, dotted, mixed) is used for both before and after nadir curves.In the first and last panels showing variation with changing cost of resistance and death rate, respectively, values beyond 0.9 and 0. Preventing evolutionary rescue in cancer

Mutation during treatment reduces the extinction probability
Next, we examine how ongoing mutation influences the treatment outcome.In our model, there are four types of mutation (Figure 1), and the total mutation rate is 2(µ 1 + µ 2 ).In Figure 4(A), second row, first panel, we see that increasing the total mutation rates in both E 1 and E 2 while keeping N (τ ) and the initial frequency of resistance unchanged, results in lower extinction probability.We observe the same trend if we change the mutation rate in only one environment (Figure A.5).This effect is due to higher mutation rates resulting in a larger rescue population size and hence, a higher probability of evolutionary rescue.For a total mutation rate as high as 10 −3 , the extinction probability never exceeds 0.2.On the other hand, the benefit of decreasing the mutation rate greatly diminishes after µ = 10 −5 in the before-nadir regime.In the extreme, unrealistic case of abundant pre-existing resistance and very low mutation rates, the optimal switching point would be long after N min , when the R 2 population has fallen to close to zero.Therefore, there is an upper bound on the extinction probability when we switch before the N min , but extinction probabilities for the after-nadir switching points are equal to 1 (blue curve in Figure 4(A), second row, first panel).

Extinction probability increases with death rate and turnover
To compare treatment outcomes for different plausible combinations of birth and death rates, we plot heat maps of high-P E regions in b-d space (Figure 4(B)).We observe that in the lower right region (high birth rates, low death rates), extinction probabilities are very low.This result also holds for alternative metrics for comparing parameters (Appendix A.8, Figure A.1).This leaves us with a diagonal band in the b-d space within which it is possible to attain high extinction probabilities.
Within this "good" region, we make three major observations.First, a higher death rate results in a higher extinction probability (Figure 4(A), second row, third panel).Second, as the birth rate increases, optimal extinction probability decreases and the optimal N (τ ) increases (Figure A.5).However, when the birth rate (of resistant cells) is close to the death rate (b − c ≈ d), extinction probabilities are close to 1.The high-P E regions within the "good" region remain more or less the same with changes in birth rate.
Third, we observe that the high-P E regions become larger as we increase the turnover (defined as the sum b + d) while keeping the intrinsic growth rate g S constant (dashed line in Figure 4(B)).Note that the cost of resistance is always a fixed fraction (0.5 by default) of the birth rate of S cells.It follows that when increasing turnover while keeping the growth rate g S constant, the growth rate g R of resistant cells decreases.This leads to a smaller rescue population, contributing to the increase in extinction probability.Another effect of turnover may relate to the establishment probability of resistant mutants.As noted in Appendix A.5, turnover appears in the expression for estimating the establishment probability π e .Higher turnover at a constant net growth rate g S leads to a lower π e .If it is harder for resistant lineages to establish, then there will be fewer rescue lineages, leading to better treatment outcomes.

Extinction probability is insensitive to carrying capacity
As the carrying capacity is increased from N (0) (default value), we see a slight increase in extinction probability, but this effect saturates before K = 10N (0).This is demonstrated in Figures 4(A), first row, third panel, using the analytical model and in Figure A.7 with stochastic simulation results.Systems with a lower K have an extra constraint on population growth since the initial population is closer to the carrying capacity.In our model, this results in a lower decay rate for S cells and a higher growth rate for R 1 cells.As explained in Section 3.3, the cost of resistance and the R 2 population size at a given N (τ ) affect the extinction probability.By default, if both factors are present, then an increase in K causes a slight increase in P E (τ ).The individual effects of both these factors are shown in Figure A.6.The same figure shows that, in the absence of both factors, an increase in the carrying capacity results in a small decrease in P E (τ ).

Discussion
We study a novel evolutionary therapy for cancer that aims to push tumours to extinction by exploiting stochasticity in small and vulnerable populations.This is done by applying multiple treatment strikes at appropriate times.A tumour that responds well to the primary therapy is primed for a second strike when it is small and susceptible to stochastic effects.The aim then is to "kick it while it's down" [11].This new strategy demands urgent theoretical investigation given that three clinical trials are already underway [17,18,19] and others are in development.

Preventing evolutionary rescue in cancer
Here, we have developed the first analytical model of a two-strike therapy derived from the principles of evolutionary rescue.This model is mathematically tractable and yields clearer explanations and more general results than previous approaches in modelling "extinction therapy" [16].We have also developed a complementary stochastic simulation model, which generally confirms the accuracy of our analytical predictions.We have sought to make both models as simple as possible, with minimal assumptions about parameter values and relations between different quantities.
We have used these new mathematical and computational models to investigate the optimal timing of the second strike and how the treatment outcome depends on crucial system parameters, including treatment levels and cost of resistance.The combination of analytical and computational analyses arms us with powerful tools to explore two-strike therapy in a wide range of scenarios, with a solid basis in eco-evolutionary theory.
When do we get the optimal extinction probabilities?The ability to analytically predict the optimal switching point for a large range of parameter values promises to aid the design of effective treatment schedules.By numerically solving our analytical model (Section 3.1), we have shown that the optimal N (τ ) is close to N min , the population nadir in the absence of a second strike.This result -which is supported by extensive simulations (Figure A.3) -is consistent with a hypothesis proposed in the previous investigation of two-strike extinction therapy [16].In contrast to the previous study, which suggested that striking before the N min is better, our results showed the optimal switching point may be slightly before or after the nadir (Section 3.1 and Appendix A.1.2.4).
Since it is unreasonable to expect switching to the second strike exactly at the optimal point, we examined high-P E (≥ 0.8) regions or windows of opportunity.We found that switching after the nadir generally results in a higher extinction rate than switching before (Section 3.2) because the window of opportunity extends further into the after-nadir regime.The difference in the extinction probabilities between before and after nadir switching points depends on the cost of resistance and the size of the initial population resistant to treatment 2 (Section 3.3).We conclude that it is generally better to wait slightly longer and risk missing the optimal N (τ ) than to apply the second strike too early.However, one should certainly not wait until the tumour becomes detectable again (as is the current practice) because that increases the probability that rescue mutants will emerge.
In practice, the tumour detection thresholds are much higher than the initial population sizes we consider in our model (∼ 1 million cells, see Section 3.5).Therefore, it is important to have an estimate of the population size as well as time windows in which switching to the second treatment gives a high extinction probability.In Figure 5, we compare the high-P E time window to the total duration in which the tumour is undetectable and see that the time window sizes correlate with the size of high-P E regions in most cases.A notable exception is when the cost of resistance is zero and the treatment level is high.In this case, we observed a larger high-P E region but a shorter high-P E time window compared to the case with a high cost of resistance (Figure 2(A), first row and Figure 5, fourth row).For clinical applications, the size of the window of opportunity, both in terms of time and population size, should be considered.
What are the optimal doses?The treatment levels during the two strikes (δ 1 and δ 2 ) are the easiest model parameters to control in practice.The higher the two doses, the higher the extinction probability at the optimal switching point (Figure 3.4).However, switching at the optimal size, which becomes smaller as the treatment dose is increased, may not be feasible.In this case, the treatment combination that gives a wide range of switching points with a high probability of extinction may be better.Surprisingly, at least with a high cost of resistance, we found that the largest high-P E (≥ 0.8) region in the before-nadir regime is obtained with an intermediate first dose paired with a high second dose (Section 3.4).This result emphasizes the importance of timing in two-strike therapy -a stronger treatment with a poorly chosen switching time can be worse than a weaker treatment given at the right time.An interesting implication of this result is that the two treatments need not both be very effective.An intermediate treatmentinduced cell death rate can give good treatment outcomes if the cost of resistance is high.Moreover, the optimal switching threshold is also relatively high for a less effective first treatment, which may be beneficial in practice.
What other tumour parameters determine the success of two-strike therapy?Our systematic exploration of the model parameter space reveals several noteworthy effects on treatment outcomes.First, although a high cost of resistance is predictably beneficial, we found that two-strike therapy can outperform conventional treatment even when this cost is small or non-existent (Section 3.3).Therefore, in common with adaptive therapy [20], two-strike therapy is not contingent on a cost of resistance.Moreover, we saw a variable response to a change in the cost of resistance depending on treatment level.For high treatment doses, we observed comparable extinction probabilities in the presence and absence of the cost of resistance, but for intermediate doses, a small cost of resistance gives much worse treatment outcomes than a high cost.
Second, mutation during either treatment is detrimental to treatment outcome (Section 3.6).This result suggests, for example, that mutagenic therapies may be less appropriate.Third, we find that higher death rates and higher turnover are beneficial, as has previously been shown for adaptive therapy [21].Fourth, although a higher carrying capacity allows more tumour growth, we found that for a given initial tumour size, changes in carrying capacities have little effect on treatment outcome (Section 3.8).
Although we have used simple models with minimal assumptions to ensure that our main findings are qualitatively robust, we have not explored all plausible functional forms.For example, the effect of changing the mutation rate might be different in a model in which mutations occur only at the time of cell division.Therefore, our results are subject to certain methodological assumptions and limitations.
When should two-strike therapy be used?Two-strike therapy holds most promise as an alternative to conventional therapy in cases where a very good initial response to treatment is typically followed by relapse.Our results suggest it is likely to succeed only in relatively small tumours (Section 3.5).Nevertheless, we expect that subsequent treatment strikes, following the same principle, would lead to higher extinction probabilities for larger tumours.Additionally, we have assumed a larger initial resistant population than expected for our default tumour size (see Appendix A.1.1),providing conservative estimates of extinction Preventing evolutionary rescue in cancer probabilities.Including an Allee effect [12] is also expected to increase extinction probabilities, making the therapy viable in a wider range of scenarios [16].Hence even two-strike therapy may be feasible in substantially larger tumours.Two-strike therapy may be a wise strategy when one of two available treatments is less effective than the other.Conversely, if resistant cells are abundant and have relatively high fitness, then this therapy is unlikely to succeed, and a long-term tumour control strategy such as adaptive therapy could be a better option [9,20].Even when it may be theoretically optimal, two-strike therapy crucially depends on the availability of effective treatments with low cross-resistance and methods for monitoring tumour burden over time [22].
Demonstrating the broad feasibility of two-strike therapy, the three clinical trials that are already underway in metastatic rhabdomyosarcoma [17], metastatic prostate cancer [18], and metastatic breast cancer [19] involve not only diverse cancer types but also very different classes of treatment, including chemotherapy, targeted therapies, and hormonal agents.Other proposed targets include locally advanced rectal adenocarcinoma [23] and paediatric sarcomas [22].

Conclusion and future directions:
We have shown that two-strike therapy is a theoretically sound concept that, in certain scenarios, could plausibly increase cancer cure rates.Our work provides a necessary foundation for further mathematical investigation and justification for experimental and clinical testing of this innovative strategy.
An important topic for further mathematical analysis is the prevention of evolutionary rescue with more than two strikes.Previous work on the optimal scheduling of multiple treatments [24,25,26] suggests that alternating two treatments is a theoretically sound approach.An alternative strategy, more in line with the original conception of extinction therapy, is to switch to a third treatment whenever possible.Other immediate directions for mathematical investigation include accounting for cross-resistance and considering alternative biological assumptions, such as modelling resistance as a continuous, plastic trait.

A.1 Analytic Model Without Competition
We study here the model without competition, obtained by letting the carrying capacity K go to infinity.This is a reasonable approximation of the true model (Figures 4A and A.7, bottom right subplot), and it may be solved explicitly.This allows us to better understand the impact of various parameters.
We denote by g 1 , g 2 , and g 1,2 the baseline per-cell growth-rates of the resistant strains R 1 , R 2 and R 1,2 , respectively.In the main text, for simplicity, these growth rates are assumed to be equal and denoted by g R , but this assumption is not needed here.Let γ S = g S − µ 1 − µ 2 , γ 1 = g 1 − µ 2 and γ 2 = g 2 − µ 1 denote the per-cell growth-rates of subpopulations S, R 1 , and R 2 , net of outgoing per-cell mutations rates.Let c i = γ S − γ i (i = 1, 2) denote the cost of resistance net of the mutation rates.Since the mutation rates are assumed to be much lower than the growth rates, γ S ≃ g S , γ i ≃ g i , and c i ≃ g S − g i , i = 1, 2.
During the administration of treatment 1, the model reads: This system of linear differential equations is easy to solve.Assuming γ 2 ̸ = γ S , we get:1 where the approximation for R 1 holds for t large enough.The quantity R 1,2 (t) is also easy to compute, but is not relevant.Indeed, it represents the expected number of R 1,2 cells, while we are interested in the expected number of potential R 1,2 rescue mutants, i.e, of R 1,2 cells directly arising from mutations from R 1 or R 2 cells, as opposed to reproduction of previous R 1,2 cells.The latter is given by Here where the most important terms are those of the first line.
Letting π 1 and π 1,2 denote the probability of establishment of R 2 and R 1,2 mutants, respectively, the total expected number of rescue mutants is given by where the approximation neglects the impact of mutations from S to R 1 on the dynamic of R 1 cells during the administration of treatment 2. The expected number of rescue mutants is the quantity that we seek to minimize (assuming that, in the stochastic analogue of the model, the number of R 2 cells at time τ is approximately Poisson distributed, see Appendix A.4).

A.1.1 A reference case for initial resistant populations
The initial resistant populations depend on the tumour and treatment history before the phase that we are interested in.For instance, if other treatments were given before and that resistance to these treatments is correlated to resistance to treatments 1 and 2, the initial resistant populations R 1 (0) and R 2 (0) could be large.Nevertheless, as a reference point, assume that the tumour starts at time −T with a single sensitive cell and no resistant cells, and that, until the administration of treatment 1, the tumour follows the dynamics we assumed, without treatment kill rate: Ṡ = γ S S and Ṙi = γ i R i +µ i S, i = 1, 2. Using that e γit = [e γ S t ] γi/γ S , it is then easy to see that This is equal to 5 with our default parameter values, much lower than the 100 initial R i cells we assume. 3,4 Integrating the analog of (11) between −T and 0 leads to: Plugging the initial resistant populations obtained in (11) in the expressions for R i (t) and R * 1,2 (t) leads to (during treatment 1): This leads to: where α S and α 1 were defined in (10).The key difference with ( 9) is that, since R2 (0) ≃ 0, the term in R2 (0) disappeared, and the R 2 population is fully hidden in a part of the term α S S(τ ).This allows us to study the worst-case scenario in which the treatment history leads to a large resistant population, giving us conservative bounds on the extinction probability.

A.1.2 Intuition about our findings, robustness
We now use previous formulas, with or without assuming (11), to compute some interesting quantities and get some intuition about our findings in the main text.

A.1.2.1 Computation of the nadir
Around the nadir, most tumour cells are either S cells or R 1 cells, The right-hand side may be minimized using Lemma A.1.Letting x1 (0) = R1 (0)/S 0 , this leads to: where we used that With our default parameter values, we obtain N min ≃ 2018.With the initial R 1 population given by ( 11) and our default parameter values otherwise, we would get N min ≃ 286.
Adding competition, as in the main text, would somewhat decrease these numbers since this reduces all net growth rates.

A.1.2.2 The expected number of de novo rescue mutants is minimized by switching slightly after the nadir.
This part of the analysis is also valid for the model with competition, as long as The nadir is thus reached at the time t min such that R 1 /S = α/β.
Let a 1 = π 1,2 /(δ 2 − γ 1 ) and a s = π 2 /(δ 2 − γ S ).The expected number λ DN of de novo rescue mutants satisfies where we used that close to the nadir, N ≃ S + R 1 .In the absence of costs of resistance, that is, if γ 1 = γ 2 = γ 1,2 = γ S (which also implies π 2 ≃ π 1,2 and hence a S = a 1 ), the quantity λ DN (τ ) is minimized at the nadir.If, however, we assume a cost of resistance for R 1 cells and similar probabilities of establishments for R 2 and R 1,2 cells (π 2 ≃ π 1,2 ), then a S > a 1 .Since at the nadir, Ṅ = 0 and Ṡ < 0, it follows that λ DN is still decreasing.To minimize it, it is thus best to switch after the nadir.
To be more precise, note that λ DN is minimal at the time t DN such that a S Ṡ + a 1 Ṙ1 = 0, that is, −a S αS + a 1 βR 1 = 0.It follows that: a1 thus implies that With our standard parameter values, we get N (t DN ) ≃ 1.013 N min .The expected number of de novo rescue mutants is minimized by switching just slightly after reaching the nadir.Testing other parameter values suggests that this conclusion is robust, in the sense that the ratio ( 13) varies relatively little from when changing parameters (Figure A.8).

A.1.2.3 The expected number of pre-existing rescue mutants may be minimized before or after the nadir
During treatment 1, the R 2 population decreases, but the expected number of R 1,2 rescue mutants increases.The expected number λ τ ) of pre-existing rescue mutants is minimal when these two forces exactly balance each other, i.e. when: Assuming µ 1 ≪ δ 1 − γ 2 , this happens when: If, as in the main text, γ 1,2 = γ 2 , hence π 2 = π 1,2 , this boils down to: This may occur before or after the nadir, depending on the value of parameters, and in particular of the R 2 population: the larger R 2 (0), the more important it is to wait that the R 2 population decreases.
In the reference case R i (0) = µi ci S 0 , the expected number of preexisting rescue mutants is approximately .
Recalling that at the nadir, R 1 /S = δ1−γ S γ1 , and that the ratio R 1 /S increases with time, it follows that the number of preexisting rescue mutant reaches its minimum before the nadir if and after the nadir otherwise.With our default parameters, the cost of resistance is very strong, and this inequality is satisfied (0.5 > 0.4).Thus, the number of preexisting rescue mutants reaches its minimum before the nadir.More precisely, by Lemma A.1 and using C = γ1π2 c2π1,2 , the expected number of pre-existing rescue mutants is minimized at the time t SGV such that t SGV ≃ 0.983 t min , and the tumour size is then However, with a smaller cost of resistance c 2 , the number of preexisting rescue mutants would reach its minimum after the nadir.The optimal switching time would then unambiguously be after the nadir, even though we assumed in this computation a much smaller initial R 2 population than in our simulations.Thus, the fact that the optimal switching time is after the nadir does not require a large initial R 2 population.

A.1.2.4 When to switch in total: before or after the nadir?
To see whether the optimal switching time is before or after the nadir, consider first the reference case where R 2 (0) = µ2 c2 S 0 .Since at relevant switching times, N ≃ S + R 1 , Eq. ( 12) may be rewritten: Since at the nadir dN/dτ = 0, and R is increasing, it follows that dRM/dτ has then the sign of α 1 − α S .Thus, it is best to switch after the nadir if α 1 < α S , that is, if: This condition may be understood as follows: to minimize the number of de novo rescue mutants, it is best to switch after the nadir if To minimize the number of pre-existing rescue mutants, it is best to switch after the nadir if π1,2 γ1 < π2 c2 .If both conditions hold, or if they hold in an appropriate average sense, then it is globally best to switch after the nadir.
With our default parameter values except for the initial resistant populations, α S /α 1 ≃ 0.93 < 1, so it is best to switch a bit before the nadir.More precisely, from (12), we get that the expected number of rescue mutants is minimized when α S Ṡ + α 1 Ṙ1 = 0, which leads to Denoting by t opt the time at which the expected number of rescue mutants is minimized, we then get from Lemma A.1 with C = α S α1 that t opt ≃ 0.99994 t min and that N (t opt ) ≃ 1.0005 N min .Thus, the optimal switching time is very slightly before the nadir, but so slightly that this could be due to approximation errors in our formulas, and could be easily reversed by slight changes in parameters. 5e conclude that whether the optimal switching time is before or after the nadir is not a robust finding, and depends on parameters.What seems robust, however, is that it is best to switch close to the nadir.Since the time at which the nadir is reached may not be observed in practice, it is important to note that due to the assumed cost of resistance, tumour size and the number of expected rescue mutants tends to decrease relatively quickly before the nadir, and then to rebound more slowly. 6This makes switching after the nadir a safer move (see Figure 5).
For arbitrary initial resistant populations, replacing (12) by ( 9) leads to: The previous argument shows that, with our default parameter values, the bracket is minimized very slightly before the nadir.However, in our simulations, we chose a relatively large initial R 2 population, hence R2 (0) is quite positive.Since the term R2 (0)e −(δ1−γ2)t decreases with time, switching after the nadir should give better results than switching before.This is indeed what we observed.

Preventing evolutionary rescue in cancer
A.1.2.5 A lower value of δ 1 may reduce the expected number of rescue mutants when switching before the nadir and at a given size.
For a given switching time, and not taking into account treatment toxicity, a larger value of treatment 1 kill-rate δ 1 is unambiguously beneficial: it leads to a lower sensitive population, hence also to fewer mutations from sensitive to resistant cells, and lower resistant populations.It thus leads to a lower probability of evolutionary rescue. 7or a given switching size, the situation is different, at least when we switch before the nadir, and the initial R 2 population is relatively large.In that case, as mentioned in the main text, lower values of δ 1 provide an advantage by letting the R 2 population decay more between the beginning of treatment 1 and the time when the switching tumour size is reached.We now explain this formally.
Assume that the protocol is to switch treatment the first time that tumour reaches size N , and consider various values of δ 1 such that this happens.As explained above, a lower value of δ 1 leads to a lower rate of decay of the sensitive population, hence also to more mutations from S to R i and a higher growth rate of the resistant populations.thus leads to a larger resistant population size and a larger tumour size for a given sensitive population size.It follows that to reduce tumour size to size N , the sensitive population must be reduced by a larger factor when δ 1 is low.Thus, letting τ (δ 1 ) denote the time at which the switching size N is reached, the ratio S(τ (δ 1 ))/S 0 is increasing with δ 1 (lower for a low value of δ 1 ).
Since a low value of δ 1 increases the proportion of R 1 cells when reaching size N , the value of δ 1 also affects the term However, as long as the ratio α S /α 1 is close to 1, as with our default parameter values, this effect will be small.

A.1.2.6 Impact of a cost of resistance
A key benefit of costs of resistance, if they are incurred even before treatment begins, is to decrease the initial proportion of resistant cells.This was addressed when deriving (11).By contrast, we focus here on the impact of costs of resistance for fixed initial resistant populations.
For a given switching time, all costs of resistance are beneficial.They do not impact the sensitive population, and reduce resistant populations, hence the expected number of rescue mutants.It follows that any cost of resistance increases the range of switching times leading to a large probability of extinction.
For a given switching size, resistance costs for R 2 and R 1,2 mutants are still beneficial: they do not significantly affect the switching time, nor the proportion of R 1 and S cells at that time, but reduce the expected number of rescue mutants, through a quicker decay of the R 2 population in the first phase, and the reduction of the probability of establishment of potential rescue mutants in both phases.
The effect of a cost of resistance for R 1 cells when switching at a given size is more subtle.The coefficient α 1 in ( 9) is equal to: In our model, the number of mutations from 1 to R 1,2 cells is proportional to the R 1 population but not to its net growth rate.As a result, for a given resistant population R 1 (τ ), R 1 cells will generate more R 1,2 mutants in phase 1 when incurring a large cost of resistance. 8This is reflected in the fact that the term 1/γ 1 is a decreasing function of γ 1 , hence an increasing function of c 1 .They will, however, generate fewer R 1,2 mutants in phase 2, as a higher cost of resistance means a quicker decay of the R 1 population.This is reflected in the fact that the term 1 δ2−γ1 increases with γ 1 .The net effect depends on whether γ 1 is smaller or larger than δ 2 /2 (indeed, the denominator γ 1 (δ 2 − γ 1 ) is maximal when γ 1 = δ 2 /2).With our default parameter values, γ 1 < δ 2 /2, and so the dominant effect is that a larger cost of resistance c 1 , hence a lower growth-rate γ 1 , leads to more R 1,2 mutants in phase 1.
When switching long after the nadir, the increase of α 1 is the main effect of a cost of resistance c 1 .Indeed, pre-existing R 2 mutants should be rare, and the fraction of R 1 cells when switching is anyway close to 1, hence should vary little with c 1 .A larger cost of resistance is then paradoxically detrimental.
When switching before the nadir, other effects should be taken into account: a larger cost of resistance means that it takes a bit more time to reach the switching size, giving a bit more time for the R 2 population to decay (and in particular decreasing the term in R2 (0) in ( 9) when R2 (0) > 0).This effect will be smaller long before the nadir, as the time at which the switching size is reached then depends little on the R 1 population.A larger cost of resistance also leads to a lower proportion of resistant cells when switching, which affects the term α S S + α 1 R 1 , though, as said before, this effect should remain small when the ratio α S /α 1 is close to 1.We conclude that, though the net effect depends on parameter values, a larger cost of resistance c 1 might counter-intuitively decrease the probability of survival when switching at a given size, also when switching before the nadir.

A.2 Using results from evolutionary rescue theory
While the conventional goal of evolutionary rescue is to maximize the probability of rescue, the same mathematical formulations can be applied with the aim of extinction.For our deterministic analytical model, we use the theory of evolutionary rescue to find the probability of extinction in cancer populations undergoing two-strike therapy.A population is rescued from extinction when one or more resistant lineages are established in the population.
We assume that the generation of resistant mutants in the population can be approximated by a Poisson process [10] (see Appendix A.4 for further investigation into the validity of this assumption in the context of our model).The resistant mutants that lead to evolutionary rescue are called rescue mutants.The probability of establishment of a resistant mutant can be approximated by Haldane's 2s [28], where s is the selective advantage of the mutant over the wild-type, which depends on the degree of stress.For a more general, continuous-time estimate, we use diffusion approximation to get the establishment probabilities of resistant lineages starting with a single cell (derived in Appendix A.5).The establishment probability is multiplied by the rate of generation of resistant mutants to obtain the rate of generation of rescue mutants.Therefore, once the rate of the Poisson process is calculated, it is easy to obtain the probability of extinction, which will be equal to the probability that no rescue mutants are generated.
We assume density independence of resistant lineages and use branching processes to model the emergence of new resistant (and rescue) mutants [29].The growth model we use is described in Appendix A.5. Density independence is a reasonable assumption because rescue lineages become important only when the population is much smaller than the carrying capacity.Additionally, we make the assumption that the probability of establishment of a resistant mutant is constant throughout.There is extensive literature on what happens when these assumptions do not hold.Analytical results can be derived in all such situations using stochastic methods and taking continuous time approximations [30,31,32].
An important distinction is that the probability of rescue by pre-existing mutants (from before the onset of stress, called standing genetic variation) is different from that of new mutants (via de-novo mutations) [29].This is because the pre-existing mutant (SGV) lineages typically get more time to grow than the denovo mutants (DN).There are several expressions by different authors [31,30,29] for the probability of evolutionary rescue following a single change in environment (treatment strike in our case) by both these classes of mutants, but considering the common conceptual basis, they all reduce to the following form: where N (0) is the population size at the onset of stress (t = 0), µ is the total per capita, per unit time mutation rate after the onset of stress, f (0) is the frequency of resistant variants at t = 0 and t ext is the time at which the population would go extinct if it is not rescued.Thus, the first and second terms, respectively, represent evolutionary rescue by a population that only becomes resistant to the stress induced at t = 0 and evolutionary rescue by a rescue population generated by mutation from N .
In Section 2.1, we use these results extensively in the context of two-strike therapy.The principles of evolutionary rescue provide a foundation for building the theoretical formulation of two-strike therapy.There are, of course, significant differences to account for -firstly, most evolutionary rescue models consider either one abrupt change in the environment or a continuous, gradual change [33].However, our approach is based on using two or more subsequent strikes, all of which are abrupt changes in the environment.Second, most existing models consider a single resistant variant (an exception is G. Martin et al. ( 2013) [31]), while the existing model for extinction therapy [16] works with a continuum of resistance effects.We choose to develop the simpler case in which we have discrete genetic variation.We, therefore, theoretically understand two-strike therapy as the prevention of evolutionary rescue.

A.3 Derivation of extinction probabilities
Given the populations S, R , where π e is the probability of establishment of a single resistant lineage and µ 1 and µ 2 are the rates of acquiring resistance to treatments 1 and 2. The probability of establishment depends on b, d and c (see Appendix A.5 for the derivation).Following the Poisson assumption, the probability that all pre-existing mutants go extinct in E 2 is equal to To find the probability that no de-novo rescue mutants arise in E 2 , we assume that the generation of new mutants is a Poisson process, and the number of rescue mutants in E 2 is Poisson distributed with a rate equal to λ DN = π e µ 2 ( ∞ τ S(t)dt + R 1 (t)dt).This rate is proportional to the rate of mutation from the sensitive cell fraction (S and R 1 ) to the resistant cell fraction (R 2 and R 1,2 ) in E 2 .
The total extinction probability as a function of τ is given by the product of the probability of extinction of pre-existing and de-novo mutants, The first two terms inside the exponent account for the expected number of pre-existing R 2 and R 1,2 rescue mutants, and the last term accounts for the expected number of de-novo rescue mutants.To make the computation easier, we calculate the approximate values of the last two integral terms by integrating them until the population size is equal to one.
To compute the expression in (15) numerically, we evolve the deterministic logistic growth equations for dynamics in E 1 (from t = 0 to τ , given in Figure 1) and in E 2 (from t = τ extinction, given below).In ( 16) and ( 17), we ignore the mutations to R 2 and R 1,2 and the changes in these resistant populations because we use only the deterministic decay of the sensitive population to calculate extinction probabilities. )

A.4 Validity of the Poisson assumption
To compute extinction probabilities using our analytical model, we assume that the distribution of established rescue lineages (established lineages of R 1,2 generated in either treatment phase and R 2 in the second treatment) and R 2 cells at the switching time is Poisson (rates given in Section A.3).This assumption is generally valid for the established rescue lineages, as they closely follow a pure Poisson point process, but it can break down for the pre-existing R 2 cells in certain cases.
To investigate the distribution of R 2 cells analytically, consider the stochastic master equation formulation of the Probability Mass Function (PMF) for the number of R 2 cells (N R2 ) during treatment 1: where δ 0,n denotes the Kronecker Delta function evaluated at (0, n) (δ 0,n = 1 if n = 0, and δ 0,n = 0 otherwise) and the S population is approximated as a deterministic function of time.This represents a birth-death-mutation process with time-inhomogenous event rates.Mutations from R 2 to R 1,2 are ignored.
Following Getz [34], (18) can be converted to a single partial differential equation involving the probability generating function Equation ( 19) can be solved using the method of the characteristics [34]: with k(z, t) = τ ) dτ and with approximation (21) holding for k(z, t) large.There are two types of populations.Internal cells are 2 cells stemming from the initial R 2 (0) population and immigrant cells are R 2 cells stemming from the mutant population.The random variables representing these populations are roughly independent at all times as the total R 2 population contributes little to competition, so in (20) the first term represents the PGF of immigrant cells and the second term internal cells.
To obtain the probability of extinction of all R 2 cells that arise before switching to the second treatment, let N e R2 denote the number of R 2 cells existing at the time of switching that eventually get established.Since The probability of extinction from R 2 cells arising before switching is just the mass of (0; t).In the limit where the term 1) and using approximation (21), we get: Expression ( 22) is what we use in (15).However when b eff R (t) is made large at fixed g eff R (t) (or b−c d+δ1 ≈ 1), G N R 2 (1 − π e ; t) is strictly increased and can be made arbitrarily closer to 1 9 .
Preventing evolutionary rescue in cancer Consequently, the probability of extinction of the SGV R 2 population calculated by our analytical model (Eq 2) is expected to be an underestimate of the probability of extinction of this population observed in our stochastic simulations especially when birth and death rates are comparable.For example, we observe this discrepancy in simulations with a zero cost of resistance and intermediate treatment dosages (second column of Figure 2 and Figure A.2).By SGV R 2 population, we mean R 2 cells descending from R 2 cells initially present or from R 2 cells arising from mutations of S cells during treatment 1.

A.5 Derivation of the establishment probability
In this section, we provide the derivation of the establishment probability of a single lineage of resistant cells.This result is taken from Lambert, 2006 [35].As described in Appendix A.2, the establishment probability is the probability with which resistant mutants establish themselves in the population and lead to evolutionary rescue.These mutants are then called rescue mutants.
We use the Continuous-time, real-valued Branching (CB) process to model the population dynamics of the resistant lineages.We chose this process because CB processes are the only diffusion processes that satisfy the additive property of the well-known BGW processes, which are commonly used to model stochastic population dynamics.Therefore, a CB process can be used to model the total population size summed over several lineages evolving independently.In our case, we use a CB process which is a continuous function of time, also called branching diffusion.As described in Lamperti, 1967 [36], branching diffusions are strong solutions of SDEs of the form: where B is the standard Brownian motion, r ∈ R is the intrinsic growth rate, and σ ∈ R + is the reproduction variance, defined as [31]: We use a general result from diffusion theory that the probability u(z 0 ) of the diffusion hitting an absorbing barrier z abs , where z 0 is the initial condition (Z 0 = z 0 ) is such that the function u solves the equation Gu(z) = 0 subject to appropriate boundary conditions.Here, G is the infinitesimal generator of the diffusion and characterizes the behaviour of the diffusion at small time intervals.For our one-dimensional diffusion (Equation 23), G is of the form, For the absorbing barrier z abs = 0, we will obtain the extinction probability of a resistant population starting with z 0 cells by solving Gu(z 0 ) = 0, with boundary conditions u(z abs ) = 1 and u(∞) = 0.The solution of this differential equation is, As explained in Martin et al., 2013 [31], the reproduction variance for a simple birth-death process with the birth rate b and death rate d can be approximated by b + d (turnover).Consequently, we obtain the following expression for the establishment probability, defined as one minus the extinction probability, for a resistant lineage starting from a single cell.
We use this result to compute the extinction probabilities with our deterministic analytical model in Section 2.1.

Preventing evolutionary rescue in cancer
Stopping condition Outcome condition Outcome N (t) = 0 N (t) = 0 Extinction t ≥ T N (t) ≥ N (0) Progression N (t) < N (0) Persistence N (t) ≥ min(0.99K,N max ) N (t) ≥ N (0) Progression N (t) < N (0) Persistence Table A.1: Stopping conditions and corresponding outcomes of a simulation.In the second condition, T is the maximum simulation time defined at the beginning.If we see a significant number of outcomes with persistence (more than 10%), it means that T is not large enough and the simulation is run again with a higher T value.In the third condition, the threshold 0.99K is arbitrary.The outcome remains the same as long as the threshold is close to K. We take the minimum of two quantities for cases where K is much larger than the initial population size, and a threshold of Nmax is enough to declare the outcome.In all cases, K ≥ N (0), and Nmax > N (0), so persistence in the last condition is not observed.However, the condition is mentioned for the sake of completion.
to the total rate i ω i (t).Alternatively, one can generate a random number z 1 from the uniform between 0 and 1 and use the following formula to determine the time interval for the next event: t int = − ln(z 1 )/ i ω i (t) 3. Calculate event probabilities for the next event by dividing the individual event rates by the total rate: p i (t) = ω i (t)/ i ω i (t).Use these probabilities to select the next event by generating a random number z 2 between 0 and 1.The chosen event would be the largest j such that z 2 − j i=0 p i (t) > 0. 4. Implement the chosen event by updating the population of the corresponding cell types.Increment time t = t + t int .For example, if the chosen event is the birth of an R 2 cell, then R 2 (t + t int ) = R 2 (t) + 1. 5. Repeat steps 2-4 till a stopping condition is reached.
Note that the effects of carrying capacity and treatment are included while specifying birth and death rates in Section A.6.2.Simulations are stopped under one of three conditions: if the population goes extinct, if it exceeds the maximum simulation time, or if it exceeds some maximum population size.Similarly, the outcome of one run of a simulation can be one of three possibilities: extinction (N (t) = 0), progression (N (t) ≥ N (0)) or persistence (0 < N (t) < N (0)).Note that there is not a one-to-one correspondence between stopping conditions and simulation outcome (see Table A.1)

A.6.2 Determining rates of demographic event
Following our variant of the Gillespie algorithm, one must define all possible demographic events at the beginning of the simulation and define rates corresponding to those events at each time step.Note that an individual event includes the type of event and the type of cell.For example, a mutation event S− →R 1 is one individual event and has a rate specified for it.Similarly, the birth of an S cell is a different event than the birth of an R 1 cell.All simulation parameters used to compute the individual event rates are listed in Table 1.
We derive the birth and death rates for all cell types from the deterministic Logistic model for population growth.
where g M (t) is the growth rate of subpopulation M ∈ {S, R 1 , R 2 , R 1,2 }.The growth rates of sensitive and resistant subpopulations differ by the cost of resistance, where the total treatment-induced death rate in environment i is δ i , and the presence of this term depends on the sensitivity of different cell types in both environments.For example, the growth rate of R 1 cells in E 1 will not include the treatment-induced death term, but in E 2 will have a δ 2 term.

2 :Figure 1 :
Figure 1: A schematic describing the deterministic analytical model.Part 1 (left) shows the ODE population growth model during the first treatment (in E1).Part 2 (right) uses input from the ODE model and evolutionary rescue theory to calculate extinction probabilities if the second treatment is given at time τ .The solid curve shows the population trajectory when only treatment 1 is applied.Sensitive cells are denoted by S. Cells resistant to treatment 1(2) and sensitive to treatment 2(1) are denoted by R1(R2).Cells resistant to both treatments are denoted by R1,2.The per capita rate of acquiring resistance to treatment 1(2) is denoted by µ1(µ2).Growth rates gS for sensitive cells and gR for resistant cells depend on the intrinsic birth rate, intrinsic death rate and the cost of resistance (see Table1).The treatment-induced death rate is denoted by δ1.Initial conditions are specified by the initial population sizes of S, R1, R2 and R1,2 cells.The total initial population N (0) is the sum of these four subpopulations.
parameters and initial conditions used in the analytical and stochastic simulation models, along with their default values.Note that for the analytical model, we use the values of growth rates for sensitive and resistant cells, gS = b − d and gR = b − c − d, respectively.

Figure 2 :
Figure 2: (A) Comparing stochastic simulation results (dots) to analytical estimates (solid line) of extinction probabilities PE(τ ) for different values of N (τ ) implemented before reaching Nmin (black) and after crossing Nmin (blue).Red dashed lines show the expected Nmin (calculated with the analytical model).Results for two treatment levels are shown (rows).Columns show results with different values for the parameter c (cost of resistance) and R2(0) (initial R2 population size).The default parameter set is given in Table1.Extinction probabilities from simulations are computed using the outcomes of 100 independent runs, each with several switching points.Switching sizes smaller than Nmin are equivalent to the absence of a second treatment, and result in extinction probabilities equal to zero (black points to the left of red dashed lines).Error bars show 95% binomial proportion confidence intervals.In the first column, labels (1) and (2) correspond to curves in Figure3(C).(B) An illustration of switching points before and after Nmin, implemented with the same random seed.See Appendix A.6 for a description of the algorithm for these simulations.(C) Nq vs q plots for four parameter sets and different treatment levels.Black curves show before-nadir switching points, and blue curves show after-nadir switching points.Labels indicate treatment levels δ1 = δ2.The same line style (solid, dashed, dotted, mixed) is used for a given treatment level for both before and after nadir curves.

Figure 3 :
Figure 3: (A) Heatmaps (obtained from the analytical model) showing high-PE regions (≥ 0.8) for different combinations of treatment levels δ1 and δ2.The default case is shown on the left, and the case with no cost of resistance and no initial R2 population is on the right.For each case, both before-nadir and after-nadir switching points are considered.White lines indicate optimal extinction probability contours.(B) Extinction probabilities for two combinations of treatment levels where δ1 ̸ = δ2.Dots show simulation results and solid lines indicate analytical model predictions.Extinction probabilities are obtained from 100 paired simulations with different random seeds.Error bars show 95% binomial proportion confidence intervals.(C) Nq vs q plots for the default case and the case with no cost of resistance and no initial R2 population.In both cases, four treatment dose combinations are shown (labels refer to treatment combinations shown in (B) and Fig 2(A)).Black(blue) lines show before(after)-nadir switching points.The same line style (solid, dashed, dotted, mixed) is used for a given treatment level for both before and after nadir curves.

Figure 4 :
Figure 4: Effects of varying parameter values or initial conditions.(A) Nq versus q plots for several parameters and initial conditions.The x-axes show extinction probability threshold q, and y-axes are the (normalised) N before q 5 are not considered since the resistant cell growth rate becomes negative.(B) Heatmap of high-PE regions for different parameter values in b-d space.Only non-negative growth rates (excluding the effects of treatment) are considered (d ≤ b − c, solid black line).The dashed black line indicates the set of birth and death rates corresponding to our default growth rate (b − d = gS = 0.9).Solid white lines show optimal extinction probability contours.

Figure 5 :
Figure 5: Time windows of high extinction probability (solid lines obtained from the analytical model and dots with 95% confidence intervals from 100-replicate Gillespie runs).Each subplot shows extinction probability trajectories under different modifications of the default parameter values.Population sizes are shown in the coloured bar, and the dotted red line shows when the nadir is achieved in the absence of a second strike.

Figure A. 3 :Figure A. 4 :
Figure A.3: Extinction probabilities for several switching points relative to the Nmin (different for each of 100 independent runs).Black(blue) dots show simulation results for before(after)-nadir switching points.
Mutation rate for acquiring resistance to treatment 1, 2 2.5 × 10 −6 δ 1 , δ 2 1, R 2 and R 1,2 at time τ , we first compute the probability of no evolutionary rescue due to standing genetic variation.From evolutionary rescue theory (Appendix A.2), we know that the distribution of pre-existing rescue variants can be reasonably approximated by a Poisson distribution with a rate proportional to the number of rescue mutants at the beginning of the second strike (see Appendix A.4 for exceptions).This gives us a good estimate of the rescue probability due to R 2 mutants.However, to estimate the rescue probability due to R 1,2 mutants, we must consider the entire duration of the treatment because these cell types are resistant in both environments (E 1 and E 2 ).Therefore, in our model, the distribution of pre-existing R 2 and R 1,2 resistant cells is Poisson with a rate equal to λ SGV