Preventing evolutionary rescue in cancer

Abstract

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: cancer treatment; evolutionary rescue; evolutionary therapy; mathematical oncology; therapeutic resistance.

Figure 1:

A schematic describing the deterministic analytical model. Part 1 (left) shows the ODE population growth model during the first treatment (in $E_1$). 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 $\tau$. 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 $R_1\left(R_2\right)$. Cells resistant to both treatments are denoted by $R_{1,2}$. The per capita rate of acquiring resistance to treatment 1(2) is denoted by ${\mu }_1\left({\mu }_2\right)$. Growth rates $g_S$ for sensitive cells and $g_R$ 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 ${\delta }_1$. Initial conditions are specified by the initial population sizes of $S$, $R_1$, $R_2$ and $R_{1,2}$ cells. The total initial population $N\left(0\right)$ is the sum of these four subpopulations.

Figure 2:

(A) Comparing stochastic simulation results (dots) to analytical estimates (solid line) of extinction probabilities $P_E(\tau )$ for different values of $N(\tau )$ implemented before reaching $N_{\min }$ (black) and after crossing $N_{\min }$ (blue). Red dashed lines show the expected $N_{\min }$ (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 $R_2\left(0\right)$ (initial $R_2$ 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 $N_{\min }$ 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 $N_{\min }$, implemented with the same random seed. See Appendix A.6 for a description of the algorithm for these simulations. (C) $N_q$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 ${\delta }_1={\delta }_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:

(A) Heatmaps (obtained from the analytical model) showing high-$P_E$ regions (≥ 0.8) for different combinations of treatment levels ${\delta }_1$ and ${\delta }_2$. The default case is shown on the left, and the case with no cost of resistance and no initial $R_2$ 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 ${\delta }_1\ne {\delta }_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) $N_q$vs $q$ plots for the default case and the case with no cost of resistance and no initial $R_2$ 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:

Effects of varying parameter values or initial conditions. (A) $N_q$versus $q$ plots for several parameters and initial conditions. The $x$-axes show extinction probability threshold $q$, and $y$-axes are the (normalised) $N_q^{\text{before}}$ (black) and $N_q^{\text{after}}$ (blue) values. The only plot which is not normalised by the initial population is the one showing variation in $P_E$ with $N\left(0\right)$ (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-$P_E$ regions for different parameter values in $b\text{-}d$ space. Only non-negative growth rates (excluding the effects of treatment) are considered ($d\le b-c$, solid black line). The dashed black line indicates the set of birth and death rates corresponding to our default growth rate $\left(b-d=g_S=0.9\right)$. Solid white lines show optimal extinction probability contours.

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.