Authors: Srishti Patil, Armaan Ahmed, Yannick Viossat, Robert Noble
Categories: Article, cancer treatment, evolutionary rescue, evolutionary therapy, mathematical oncology, therapeutic resistance
Source: bioRxiv
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.
Keywords: mathematical oncology, evolutionary therapy, evolutionary rescue, therapeutic resistance, cancer treatment
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 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 multi-strike 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.
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 E1 and E2 (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 (R1 and R2), or resistant to both treatments (R1,2). The time of switching to the second treatment is τ and the population size at this time is N(τ).
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.
Figure 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 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.
To calculate extinction probabilities due to standing genetic variation (pre-existing rescue mutants) and de-novo 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 in environment E1 of the four subpopulations S(t), R1(t), R2(t) and R1,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 (R1,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 R1,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.
Given the population composition at treatment switching time τ, we next compute the probabilities PESGV(τ) and PEDN(τ) 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 pre-existing and de-novo mutants during the second treatment period, E2, the tumour extinction probability PE(τ) is the product of these two
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.
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 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.
To compare treatment outcomes for varied parameter values, we use a summary variable Nq 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 PE(τ) generally decreases as N(τ) increases, unless N(τ) is very close to the population nadir that would pertain in the absence of a second strike (Nmin). Therefore, for a given extinction probability q (with 0≤q≤max(PE(τ))), we can obtain a corresponding value Nq, 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(τ)≤Nq, 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. Nq is therefore defined for both before-nadir and after-nadir switching
where t(Nmin) is the time at which the nadir would be reached in the absence of a second strike.
Nq 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 N0.1(before) and N0.9(before) is slight, then the extinction probability must increase steeply between these two population sizes. A higher value of N0.9(before) implies a wider window of opportunity for implementing a successful second strike. More generally, we want the range of N(τ) values with low PE(τ) to be as small as possible. We plot Nq versus q to analyse the trend in extinction probabilities across the range of potential switching points. Several metrics for comparing parameter sets based on Nq are described in Appendix A.8.
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 strike) creates a stressful environment that we denote E1. After switching to the second treatment, the tumour enters the second stressful environment, E2.
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 (R1 and R2) and resistant to both treatments (R1,2). Consequently, R1 and R1,2 cells are resistant in E1, and R2 and R1,2 cells are resistant in E2. Even if they initially rescue the population, all R1 cells will eventually go extinct due to the second strike. Evolutionary rescue from the second strike must be due to surviving R2 or R1,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 Nmin, 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).
Figure (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 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) Nqvs 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.
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 Nq (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.
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 Nmin (Figure 2). According to the analytical model, the optimal switching point may, depending on parameter values, lie slightly before or after Nmin (Appendix A.1.2.4. Yet the difference between the optimal N(τ) and the Nmin is generally so small that it is not captured by our simulation results. The difference is significant only when R2 cells are initially abundant, and the cost of resistance is low (Figure 2, second column).
To explain why the optimal N(τ) is close to Nmin, we refer to Eq 2 and see that the maximum PE(τ) will be achieved by minimizing the sum of all the rates of generating rescue mutants. The last two integral terms are minimized slightly after Nmin (see Appendix A.1.2.2 for an approximate analytical expression). However, in the terms constituting PESGV (pre-existing rescue mutants at the beginning of the second strike), the decay in the R2 population relative to the generation of R1,2 mutants determines where the optimal N(τ) lies relative to Nmin. For further analysis, we focus on regions of high extinction probability (PE≥0.8) instead of the exact optima. Consistent with our analytical predictions, we observe that these high-PE regions lie around Nmin (Figures A.3 and A.2).
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 Nmin usually have higher extinction probabilities than those before Nmin (Figure 2(A), blue versus black).
We hypothesize that this result is due the pre-existing (or initially accumulated) R2 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 R1,2 mutants to accumulate. In most cases of interest, the generation of R1,2 mutants is slower than the decay of the R2 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.
Figure Effects of varying parameter values or initial conditions. (A) Nqversus q plots for several parameters and initial conditions. The x-axes show extinction probability threshold q, and y-axes are the (normalised) Nqbefore (black) and Nqafter (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.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.
Extinction probabilities depend on the number of R2 cells present at the time of switching, which in turn depends on both the cost of resistance and the number of R2 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.
A cost of resistance is expected to hasten the decay of R2 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 counter-intuitive 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 R1,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 Nmin, which makes it impossible to achieve high rates of extinction (Figure 2(A), second row, second column).
Figure 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.
Assuming that R2 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 R2 population breaks down (see Sections A.3 and A.4 for details).
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 (PE≥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-PE regions. This is what we observe in the case of no resistance cost and R2(0)=0 (Figure 3(A, C)). However, for our default parameter values, higher values of δ1 give smaller high-PE 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 Nq versus q plot for four δ1−δ2 combinations confirms that intermediate δ1 and high δ2 produce the best treatment outcome in terms of high-PE 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.
Figure (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) Nqvs 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.
The somewhat counter-intuitive result is explained by the interaction of the dose, the cost of resistance, and the R2 population. For lower δ1, the S cells decay more slowly, so the switching point N(τ) is reached later. This provides more time for the R2 population to decay but also more time for R1,2 mutants to arise. When the cost of resistance and the initial R2 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 R2 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-PE 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-PE 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 Nmin and the high-PE region, in this case, are relatively large.
Note that this result depends on the fact that we compare high-PE 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-PE 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-PE 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 R2 population, the optimal treatment combinations in the before-nadir and after-nadir regimes are similar (Figure 3(A,C)).
Using the analytical model, we compare different values of Nq (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 Nq 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.
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 E1 and E2 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 Nmin, when the R2 population has fallen to close to zero. Therefore, there is an upper bound on the extinction probability when we switch before the Nmin, but extinction probabilities for the after-nadir switching points are equal to 1 (blue curve in Figure 4(A), second row, first panel).
To compare treatment outcomes for different plausible combinations of birth and death rates, we plot heat maps of high-PE 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-PE regions within the “good” region remain more or less the same with changes in birth rate.
Third, we observe that the high-PE regions become larger as we increase the turnover (defined as the sum b+d) while keeping the intrinsic growth rate gS 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 gS constant, the growth rate gR 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 gS 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.
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 R1 cells. As explained in Section 3.3, the cost of resistance and the R2 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 PE(τ). 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 PE(τ).
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.
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.
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 Nmin, 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 Nmin 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-PE (≥ 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-PE 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-PE 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-PE region but a shorter high-PE 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.
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-PE (≥ 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 treatment-induced 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.
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.
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 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].
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.
SP and RN were supported by the National Cancer Institute of the National Institutes of Health under Award Number U54CA217376. YN benefited from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 955708. YN and RN were also supported by a London Mathematical Society Research in Pairs award (reference 42320). The opinions expressed in this document reflect only the author’s view and in no way reflect the European Commission’s opinions. The European Commission is not responsible for any use that may be made of the information it contains. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. All simulations were performed on City, University of London’s Hyperion cluster.
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 g1, g2, and g1,2 the baseline per-cell growth-rates of the resistant strains R1, R2 and R1,2, respectively. In the main text, for simplicity, these growth rates are assumed to be equal and denoted by gR, but this assumption is not needed here. Let γS=gS−μ1−μ2, γ1=g1−μ2 and γ2=g2−μ1 denote the per-cell growth-rates of subpopulations S, R1, and R2, net of outgoing per-cell mutations rates. Let ci=γ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≃gS, γi≃gi, and ci≃gS−gi, i=1,2.
During the administration of treatment 1, the model
This system of linear differential equations is easy to solve. Assuming γ2≠γS, we ^^1^^
where the approximation for R1 holds for t large enough. The quantity R1,2(t) is also easy to compute, but is not relevant. Indeed, it represents the expected number of R1,2 cells, while we are interested in the expected number of potential R1,2 rescue mutants, i.e, of R1,2 cells directly arising from mutations from R1 or R2 cells, as opposed to reproduction of previous R1,2 cells. The latter is given by
Here, R1,2∗(0)=∫−T0μ2R1(t)+μ1R2(t)dt where −T is the time of tumour initiation.
Explicitly,
where the most important terms are those of the first line.
Letting π1 and π1,2 denote the probability of establishment of R2 and R1,2 mutants, respectively, the total expected number of rescue mutants is given by
where the approximation neglects the impact of mutations from S to R1 on the dynamic of R1 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 R2 cells at time τ is approximately Poisson distributed, see Appendix A.4).
Letting R˜1(0)=R1(0)+μ1S0δ1−c1, R˜2(0)=R2(0)−μ2S0c2, and using (7) and (8), we obtain the following approximation of the expected number of rescue mutants when e−(δ1−γS)τ≪1:
where
and
is a constant that impacts the probability of extinction but not the comparison between switching times.^^2^^
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 R1(0) and R2(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 S˙=γSS and R˙i=γiRi+μiS,i=1,2. Using that eγit=[eγSt]γ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 Ri cells we assume.^^3^,^4^^ Integrating the analog of (11) between −T and 0 leads
Plugging the initial resistant populations obtained in (11) in the expressions for Ri(t) and R1,2∗(t) leads to (during treatment 1):
This leads
where αS and α1 were defined in (10). The key difference with (9) is that, since R˜2(0)≃0, the term in R˜2(0) disappeared, and the R2 population is fully hidden in a part of the term αSS(τ). 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.
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.
We first need a lemma.
Let A, B, α, β be positive real numbers. The quantity x(t)=Ae−αt+Beβt is minimal at tmin=1α+βln(αAβB) and its value is then
Moreover, if C>0, the function f(t)=Beβt/Ae−αt satisfies f(t)=Cf(tmin) when e(α+β)(t−tmin)=C, and we then have
Proof. Let y(t)=Ae−αt, z(t)=Beβt so that x(t)=y(t)+z(t). The function x is convex and is thus minimal when its derivative cancels. Since x′(t)=−αy(t)+βz(t), this occurs when z/y=α/β. This leads to e(α+β)tmin=αAβB, hence the formula for tmin. Moreover, z/y=α/β implies y/x=β/(α+β) and z/x=α/(α+β). Since ddtln(yz)=α+β, we get f(t)=f(tmin)e(α+β)Δt where Δt=t−tmin. Thus, f(t)=Cf(tmin) when C=e(α+β)Δt. We then have y(t)=y(tmin)e−αΔt=y(tmin)C−αα+β and similarly, z(t)=z(tmin)Cβα+β. Therefore,
Around the nadir, most tumour cells are either S cells or R1 cells,
The right-hand side may be minimized using Lemma A.1. Letting x˜1(0)=R˜1(0)/S0, this leads
where we used that δ1−γS+γ1=δ1−c1. With our default parameter values, we obtain Nmin≃2018. With the initial R1 population given by (11) and our default parameter values otherwise, we would get Nmin≃286. Adding competition, as in the main text, would somewhat decrease these numbers since this reduces all net growth rates.
This part of the analysis is also valid for the model with competition, as long as Nmin≪K. Let α=δ1−γS and β=γ1. At relevant switching times, N≃S+R1, S˙=−(δ1−γS)S=−αS and R˙1≃γ1R1=βR1, so N˙≃−αS+βR1. The nadir is thus reached at the time tmin such that R1/S=α/β.
Let a1=π1,2/(δ2−γ1) and as=π2/(δ2−γS). The expected number λDN of de novo rescue mutants satisfies
where we used that close to the nadir, N≃S+R1. In the absence of costs of resistance, that is, if γ1=γ2=γ1.2=γS (which also implies π2≃π1,2 and hence aS=a1), the quantity λDN(τ) is minimized at the nadir. If, however, we assume a cost of resistance for R1 cells and similar probabilities of establishments for R2 and R1,2 cells (π2≃π1,2), then aS>a1. Since at the nadir, N˙=0 and S˙<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 tDN such that aSS˙+a1R˙1=0, that is, −aSαS+a1βR1=0. It follows
Using Lemma A.1 with C=aSa1 thus implies that
With our standard parameter values, we get N(tDN)≃1.013Nmin. 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 1 when changing parameters (Figure A.8).
During treatment 1, the R2 population decreases, but the expected number of R1,2 rescue mutants increases. The expected number λSGV(τ)=π2R2(τ)+π1,2R1,2∗(τ) of pre-existing rescue mutants is minimal when these two forces exactly balance each other, i.e.
that is,
Assuming μ1≪δ1−γ2, this happens
If, as in the main text, γ1,2=γ2, hence π2=π1,2, this boils down
This may occur before or after the nadir, depending on the value of parameters, and in particular of the initial R2 the larger R2(0), the more important it is to wait that the R2 population decreases.
In the reference case Ri(0)=μtctS0, the expected number of preexisting rescue mutants is approximately
Assuming R˙1≃γ1R1, this is minimal when
Recalling that at the nadir, R1/S=δ1−γSγ1, and that the ratio R1/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 pre-existing rescue mutants reaches its minimum before the nadir. More precisely, by Lemma A.1 and using C=γ1π2c2π1,2, the expected number of pre-existing rescue mutants is minimized at the time tSGV such that tSGV≃0.983tmin, and the tumour size is then N(tSGV)≃1.005Nmin.
However, with a smaller cost of resistance c2, 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 R2 population than in our simulations. Thus, the fact that the optimal switching time is after the nadir does not require a large initial R2 population.
To see whether the optimal switching time is before or after the nadir, consider first the reference case where R2(0)=μ2c2S0. Since at relevant switching times, N≃S+R1, Eq. (12) may be
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,
This condition may be understood as to minimize the number of de novo rescue mutants, it is best to switch after the nadir if π1,2δ2−γ1<π2δ2−γS. To minimize the number of pre-existing rescue mutants, it is best to switch after the nadir if π1,2γ1<π2c2. 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 αSS˙+α1R˙1=0, which leads to
Denoting by topt the time at which the expected number of rescue mutants is minimized, we then get from Lemma A.1 with C=αSα1 that topt≃0.99994tmin and that N(topt)≃1.0005Nmin. 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.^^5^^
We 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.^^6^^ This makes switching after the nadir a safer move (see Figure 5).
For arbitrary initial resistant populations, replacing (12) by (9) leads
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 R2 population, hence R˜2(0) is quite positive. Since the term R˜2(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.
For a given switching time, and not taking into account treatment toxicity, a larger value of treatment 1 kill-rate δ1 is unambiguously 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.^^7^^
For a given switching size, the situation is different, at least when we switch before the nadir, and the initial R2 population is relatively large. In that case, as mentioned in the main text, lower values of δ1 provide an advantage by letting the R2 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 Ri and a higher growth rate of the resistant populations. It 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))/S0 is increasing with δ1 (lower for a low value of δ1).
But S(τ(δ1))/S0=e−(δ1−γS)τ(δ1), so
The lower δ1, the lower the quantity in the bracket, and the larger the exponent (assuming a cost of resistance). Thus a low value of δ1 leads to a low value of e−(δ1−γ2)τ(δ1), hence to a low value of the term R˜2(0)e−(δ1−γ2)τ when R˜2(0)>0.
Since a low value of δ1 increases the proportion of R1 cells when reaching size N¯, the value of δ1 also affects the term αSS+α1R1=N[αS+(α1−αS)R1/N]. However, as long as the ratio αS/α1 is close to 1, as with our default parameter values, this effect will be small.
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 R2 and R1,2 mutants are still they do not significantly affect the switching time, nor the proportion of R1 and S cells at that time, but reduce the expected number of rescue mutants, through a quicker decay of the R2 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 R1 cells when switching at a given size is more subtle. The coefficient α1 in (9) is equal
In our model, the number of mutations from R1 to R1,2 cells is proportional to the R1 population but not to its net growth rate. As a result, for a given resistant population R1(τ), R1 cells will generate more R1,2 mutants in phase 1 when incurring a large cost of resistance.^^8^^ This is reflected in the fact that the term 1/γ1 is a decreasing function of γ1, hence an increasing function of c1. They will, however, generate fewer R1,2 mutants in phase 2, as a higher cost of resistance means a quicker decay of the R1 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 c1, hence a lower growth-rate γ1, leads to more R1,2 mutants in phase 1.
When switching long after the nadir, the increase of α1 is the main effect of a cost of resistance c1. Indeed, pre-existing R2 mutants should be rare, and the fraction of R1 cells when switching is anyway close to 1, hence should vary little with c1. A larger cost of resistance is then paradoxically detrimental.
When switching before the nadir, other effects should be taken into 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 R2 population to decay (and in particular decreasing the term in R˜2(0) in (9) when R˜2(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 R1 population. A larger cost of resistance also leads to a lower proportion of resistant cells when switching, which affects the term αSS+α1R1, 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 c1 might counter-intuitively decrease the probability of survival when switching at a given size, also when switching before the nadir.
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 de-novo 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
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 text 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.
Given the populations S, R1, R2 and R1,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 R2 mutants. However, to estimate the rescue probability due to R1,2 mutants, we must consider the entire duration of the treatment because these cell types are resistant in both environments (E1 and E2). Therefore, in our model, the distribution of pre-existing R2 and R1,2 resistant cells is Poisson with a rate equal to λSGV=πe(R2(τ)+∫0τμ2R1(t)dt+∫0τμ1R2(t)dt), 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 E2 is equal to PESGV(τ)=exp[−λSGV].
To find the probability that no de-novo rescue mutants arise in E2, we assume that the generation of new mutants is a Poisson process, and the number of rescue mutants in E2 is Poisson distributed with a rate equal to λDN=πeμ2(∫τ∞S(t)dt+R1(t)dt). This rate is proportional to the rate of mutation from the sensitive cell fraction (S and R1) to the resistant cell fraction (R2 and R1,2) in E2.
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 R2 and R1,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 E1 (from t=0 to τ, given in Figure 1) and in E2 (from t=τ till extinction, given below). In (16) and (17), we ignore the mutations to R2 and R1,2 and the changes in these resistant populations because we use only the deterministic decay of the sensitive population to calculate extinction probabilities.
To compute extinction probabilities using our analytical model, we assume that the distribution of established rescue lineages (established lineages of R1,2 generated in either treatment phase and R2 in the second treatment) and R2 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 R2 cells in certain cases.
To investigate the distribution of R2 cells analytically, consider the stochastic master equation formulation of the Probability Mass Function (PMF) for the number of R2 cells (NR2) during treatment
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 R2 to R1,2 are ignored. Following Getz [34], (18) can be converted to a single partial differential equation involving the probability generating function (GNR2(z;t)=E[zR2]):
Equation (19) can be solved using the method of the characteristics [34]:
with k(z,t)=1p(t)(z−1)−∫0tbReff(τ)[p(τ)]−1dτ,p(t)=exp(∫0tbReff(τ)−dReff(τ)dτ) and with approximation (21) holding for k(z,t) large. There are two types of populations. Internal cells are R2 cells stemming from the initial R2(0) population and immigrant cells are R2 cells stemming from the mutant population. The random variables representing these populations are roughly independent at all times as the total R2 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 R2 cells that arise before switching to the second treatment, let NR2e denote the number of R2 cells existing at the time of switching that eventually get established. Since NR2e∣NR2∼Binomial(NR2,πE):
The probability of extinction from R2 cells arising before switching is just the mass of ℙ(NR2e=0)=GNR2e(0;t). In the limit where the term ∫0tbReff(τ)[p(t)]−1dτ is negligible, which occurs when bReff(t) is made small at fixed gReff(t)(b−cd+δ1≪1) and using approximation (21), we
Expression (22) is what we use in (15). However when bReff(t) is made large at fixed gReff(t)(orb−cd+δ1≈1), GNR2(1−πe;t) is strictly increased and can be made arbitrarily closer to 1^^9^^.
Consequently, the probability of extinction of the SGV R2 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 R2 population, we mean R2 cells descending from R2 cells initially present or from R2 cells arising from mutations of S cells during treatment 1.
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
where B is the standard Brownian motion, r∈ℝ is the intrinsic growth rate, and σ∈ℝ+ is the reproduction variance, defined as [31]:
We use a general result from diffusion theory that the probability u(z0) of the diffusion hitting an absorbing barrier zabs, where z0 is the initial condition (Z0=z0) is such that the function u solves the equation Gu(z0)=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 zabs=0, we will obtain the extinction probability of a resistant population starting with z0 cells by solving Gu(z0)=0, with boundary conditions u(zabs)=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.
With a cost of resistance, we get,
We use this result to compute the extinction probabilities with our deterministic analytical model in Section 2.1.
The Gillespie algorithm with birth rate b and death rate d is a numerical solution to the master stochastic
where δ is the Kronecker delta. Suppose we start with X0=1. This is a linear birth-death process, so the solution in terms of the probability-generating function G is [34]:
Note that G(0,t)=P(Xt=0). To obtain the probability of extinction, the relevant event is ∃t,Xt=0:
Now, we form a filtration {Xt1≠0}⊆{Xt2≠0},∀t1>t2, and
The limit limt→∞G(0,t)=db, so pext=1−(1−db)=db. Therefore,
With a cost of resistance, b in (30) should be replaced with b−c. This expression is consistent with the method of simulation we use. The difference between the extinction probabilities obtained using (29) and (30) is negligible. To see this, firstly note that the first order Taylor approximations of πe(b,d) and πeA(b,d) centered about (q0,q0) are b−dq0 and so behave quite similarly for small g=b−d. Their values are made most different as (b,d)→(∞,d) for any fixed d and in this limit become πe→1−e−2 and πeA→1. These values differ by ≈ 0.13, which is quite large, but the more relevant difference is between expressions of the form e−πez and e−πeAz. As (b,d)→(∞,d), the maximum difference is around 0.05 and is obtained at z≈1.07. Finally, any qualitative comparisons between probabilities of extinction remain the same independent of choice of πe or πeA.
The total population at time t is N(t)=S(t)+R1(t)+R2(t)+R1,2(t), as defined in Section 2.1. The code developed for the implementation of this model is designed to be highly versatile and easy to modify. It can be extended to analyse more complex systems with more treatments and corresponding cell types.
The stochastic simulation algorithm (SSA), commonly referred to as the Gillespie algorithm [37], is a Monte Carlo method initially devised in 1977 for modelling the temporal evolution of chemical reactions in a well-mixed system. We implement a version of the Gillespie algorithm to simulate the population dynamics of our system in the context of two-strike therapy. The idea is that given a set of rates corresponding to birth, death and mutation events, we can track the size of each subpopulation. These rates are specified for all cell types and can change with time or in response to the environment. The basic steps of our algorithm are as
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 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 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)
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→R1 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 R1 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 gM(t) is the growth rate of subpopulation M∈{S,R1,R2,R1,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 R1 cells in E1 will not include the treatment-induced death term, but in E2 will have a δ2 term.
We can write the intrinsic growth rates as the difference between intrinsic birth and death rates like
where b, d and c are the intrinsic birth rate, death rate and the cost of resistance. Now, we separate the positive and negative parts of the subpopulation growth rates to find the effective birth and death rates to use for our simulations.
We use Equations 36 and 37 in our simulations. Note that the separation of positive and negative terms in the last step is not unique. The only constraint is that the separation follows the condition required for the implementation of carrying capacity, i.e. when the total population is equal to K, we must have bS/Reff(t)=dS/Reff(t) (without treatment). We chose a simple way to implement this by taking a density-dependent birth rate and a density-independent death rate. Additionally, it is convenient and intuitive to group the death-rate terms together, like we do in the dS/Reff(t).
For mutation events, we consider the same mutation rate μEi for all cell types, which may change with the environment. Once the source population for an event is chosen (say, S cells), the target cell is chosen according to the rates of acquiring each type of mutation. These mutation probabilities from one cell type to the other can be modified.
For this study, we consider two environments E1 and E2, each corresponding to the two strikes or treatments. Now we ask the how does the environment change with time? In other words, what is the condition under which we switch treatments and apply the second strike? Two-strike therapy aims to minimize the probability of evolutionary rescue, and the expected number of rescue mutants depends on the population size and time at which we switch between treatments. The population size at the time of switching is N(τ).
It is necessary to understand the behaviour of this variable N(τ) in order to evaluate the efficacy and limitations of two-strike therapy. As reasoned by Gatenby et al. (2020) [16]and supported by our results, the hypothesis is that the optimal N(τ) would be close to Nmin. To test this, we run a set of simulations with the same parameter values but with different random seeds. For a single random seed, we run multiple simulations (using the algorithm in Appendix A.6.1), each with a different N(τ). To produce Figures 2(A) and A.2, we run simulations with 100 random seeds for a single set of parameters and then calculate extinction probabilities for various values of N(τ). The algorithm is as
The default parameter set used for the stochastic simulations and numerical solutions are listed in Table 1. The carrying capacity is set equal to the total initial population size for consistency with previous studies on extinction therapy [16]. Moreover, we show that changing the carrying capacity has negligible effects on our results (Section 3.8). The default cost of resistance is set to a high level (half of the intrinsic birth rate), but the case with zero cost of resistance is also studied in detail to elucidate the effect of this parameter (Section 3.3). The size of the initial resistant populations is also set to 10^−4^ times the initial sensitive population, which is higher than what we obtain using the growth model described in Appendix A.1. However, in our results, we consider the case with R2=0. Moreover, our default parameter value shows that two-strike therapy can work with a conservative estimate of the initial resistant population.
After devising a way to capture the drop in extinction probabilities as N(τ) increases, we wanted to compare these curves for different sets of parameters to determine the optimal conditions for two-strike therapy. After observing the (nearly) monotonously decreasing nature of the Nq vs q plot, we came up with three metrics to achieve this goal. The first metric is the normalised area under the Nq vs q curve. A larger area under the curve (AUC) can mean one of two things – first, that the drop in extinction probability occurs at a larger N(τ) or second, that the drop occurs gradually over a larger range of N(τ) values. The first scenario is favourable for ET, and the second one can be beneficial in some cases. If the drop is close to the initial population, then it is easier to implement ET, and if the drop is gradual, then there are chances to obtain high extinction probabilities at high N(τ) values. Clearly, this metric is not very accurate, but it is a reasonable choice considering the general trend of Nq vs q curves showing a sudden drop in extinction probabilities. Therefore, the set of parameters for which the area under the Nq curve is higher is expected to result in a better outcome in terms of ease of implementation, higher extinction probabilities, or both.
The second metric we propose allows us to find the population size at which the drop in extinction probability occurs. We define the critical second strike threshold as Nc=mean(N0.1,N0.9). The higher the Nc for a system, the easier it will be to implement two-strike therapy. This metric has the advantage of approximating the time of the drop, but it is accurate only if the drop is from a value as high as 0.9 to a value as low as 0.1. Other values for the definition of Nc may be considered, but it is hard to find appropriate values that work for a wide range of parameter sets. The third metric is similar to the second one and uses the N0.5/N(0) as a measure of the population size at which this drop occurs.
We compare the three metrics in Figure A.1 and see that they produce similar results. Furthermore, these results are consistent with the trends observed when comparing regions of high-PE (Figure 4(B)).