Do mechanisms matter? Comparing cancer treatment strategies across mathematical models and outcome objectives

1.   Introduction

Treatment that delays or prevents the emergence of resistance can control cancers, potentially indefinitely, and provides a suitable strategy when eradication is impossible ^[1]. As with bacterial resistance to antibiotics or herbivore resistance to pesticides, high levels of treatment can lead to emergence of resistant strains that had been controlled by a combination of costs of resistance and competition with susceptible strains ^[2]. Adaptive therapies, where cessation of treatment precedes loss of efficacy, have been proposed as a way to delay emergence of resistance, an approach supported by mathematical modeling, laboratory experiments, and some preliminary results of clinical trials ^[3,4,5]. Modeling has played a key role in evaluating therapeutic timing, often providing evidence that reduced doses with treatment holidays can provide longer-term control ^[6].

The effectiveness of adaptive therapy depends on three assumptions: resistance is costly, resistant cells can be suppressed by competition with sensitive cells, and therapy reduces the population of sensitive cells. Under these assumptions, mathematical models have provided some support for adaptive therapy ^[5,7]. Our goal here is to investigate a wider range of mathematical models to better understand the role of these assumptions and of particular modeling choices in shaping the tradeoff between controlling cell populations and delaying the emergence of resistance. However, even for a relatively simple and well-understood cancer like prostate cancer, the mechanisms delaying resistance are not fully known. In particular, neither the role of costs nor that of competition with susceptible cell lineages has been clearly established.

For our range of models, we compare three strategies:

1. Maximum Tolerable Dose (MTD): a constant high dose determined by side effects.

2. Intermittent: a periodic scheduled dose with treatment holidays ^[8]. We here follow more common usage, rather than metronomic therapy, as recently reviewed ^[9].

3. Adaptive: initiation and termination of treatment based on the status of individual patient biomarkers, often with much earlier cessation of treatment than in intermittent strategies.

We begin by examining the original model by Zhang et al (2017) in more depth, and looking at a wider range of therapy strategies. Our central result is that all therapies follow the same tradeoff between total cancer cell burden and time to emergence of resistance.

To investigate a wider range of models, we first place the original model into a broader framework of game theoretic and competition models. We extend these in four ways that could alter responses to therapy.

● Allee effect: Control of resistance solely by competition with susceptible cell populations leads logically to the tradeoff between cancer cell burden and time to emergence of resistance. Inclusion of an Allee effect, where cell populations grow more slowly at low population sizes ^[10], could break that tradeoff.

● Healthy cells: Cancer cells compete both with each other and with unmutated cells ^[11] that respond differently to therapy and could alter the response to different therapeutic regimes.

● Immune response: Apparent competition, whereby species interact not through competition for space or resources, but as mediated by predators ^[12] or an immune response, could alter responses to therapy by introducing a delay and through their own responses to treatment ^[13,14].

● Resource competition: Like competition mediated through the immune system, explicit resource competition models introduce delays mediated through the dynamics of depletable factors. Models include consumer-resource dynamics ^[15,16,17] and more mechanistic models of androgen dynamics ^[18,19,20].

After summarizing the original model from ^[5], we present the alternative models, provide parameter values that scale dynamics to be comparable to the original, derive analytical results on the simplest of these to illustrate tradeoffs, and test the three treatment strategies. We hypothesize that intermittent and adaptive therapy will produce similar results in all cases, with the exception of the Allee effect where the rapid cessation of treatment could allow cancers to escape.

2.   Materials and methods

2.1.   The basic Zhang model (Zh)

We begin with the model published by Zhang et al (2017). This Lotka-Volterra model has three competing cell types, which we reletter for consistency with our later models: $S$ represents the population of androgen-dependent cells that are sensitive to treatment, $P$ the population of androgen-producing cells and $R$ the population of androgen-independent cells that are resistant to treatment ($x_1$, $x_2$ and $x_3$ respectively in the original).

Each cell type has an intrinsic per cell growth rate $r$ and carrying capacity $K$. The carrying capacity for $S$ cells is assumed to be proportional to $P$. The competition coefficients $a$ represent pairwise competitive effects between cells. Adaptive therapy is based on prostate specific antigen (PSA) dynamics described by

where $0.5$ represents the rate of PSA decay. This is much faster than the other time scales in the model, making PSA almost exactly proportional to total cell numbers.

This model simulates therapy by reducing the carrying capacity of $S$ and $P$. More specifically, $K_P$ is reduced by a factor of 100, and the carrying capacity of $S$ is changed from $1.5 P$ to $0.5 P$. Treatment thus results in an extremely rapid drop in these two populations. Treatment has no effect on $R$ cells which are thus released from the competitive suppression created by their smaller carrying capacity, and grow quickly until therapy is stopped. The full set of parameter values for our main simulations are given in Table 1 based on representative patient #1 ^[5].

This model includes the complexity of testosterone-producing ($P$) cells in addition to competition of sensitive and resistant cells. To test whether $P$ cells are essential to the result, we build a simpler version of the model (Zs) with two cell types, and implement therapy by directly reducing the carrying capacity of the sensitive cells, with equations

We use the same parameters as for the full model, except that $K_S = 1.5 \times 10^{4}$ in the absence of therapy and $50$ with therapy.

2.2.   General model framework

To capture the key assumptions of this model and examine the conditions that lead to success of adaptive therapy, we consider the following general framework of interaction between sensitive cells $S$ and resistant cells $R$, and an additional variable or variables $X$:

The additional dimension $X$ represents androgen-producing cells in the Zhang model, but could also be healthy cells, androgen, another resource or growth factor, or an immune response. The variable $u$ represents treatment, which will be a function of time for intermittent and adaptive therapy. The functions $r_R$ and $r_S$ describe growth as functions of population size to model competition, of $X$ to capture the tumor microenvironment via use of resources or immune attack, and of $u$ to represent the effects of treatment. The death terms $\delta_S$ and $\delta_R$ are included to separate births and deaths, and as a place for treatment effects $u$, here restricted to sensitive cells.

As outputs, we solve for the time of two types of treatment failure:

● The time $T_C$ when the total cancer cell population exceeds some threshold $C_{crit}$,

● The time $T_R$ when the population of resistant cells exceeds some threshold $R_{crit}$.

As the costs, we also track two outputs:

● Mean cancer cell numbers,

● The fraction of time being treated.

We use mean cancer cell numbers, which effectively assumes that costs are linear in the number of cells, both for simplicity and because risks of mutation and metastasis will be proportional to the number of cells. We do not attempt optimal control analysis ^[21], and present as results tradeoffs between the times, treatment, and cancer burden, looking for conditions where maximizing time to emergence reduces both treatment and cancer burden.

2.2.1.   Game theoretic model (GT)

The simplest version of this framework uses game theoretic models that focus on how strategy frequencies depend on frequency-dependent fitness. A basic model with density dependence is given by

where $C = S+R$. We place treatment costs in the growth rate of $S$ cells, and give both cell types the same carrying capacity and symmetric competitive effects.

2.2.2.   Lotka-Volterra model (LV)

Equation 2.5 is a the special case of a Lotka-Volterra model with equal competition coefficients. We here generalize to a model similar to equation 2.1, but with a more realistic approach to the effects of therapy and consequences of resistance.

Treatment increases the death rate $\delta_S$. We scale the competition coefficients describing the effect of each type on itself to 1, but can have asymmetric competitive effects and different carrying capacities for the two cell types.

2.2.3.   Allee Effect (AL)

The next model complements this framework with an Allee effect, whereby cancer cells grow more quickly when the population is above some threshold. We use the form

The parameter $b$ scales the strength of the effect, with $b = 1$ reducing to equation 2.6. Values of $0 < b < 1$ create a weak Allee effect where populations grow more slowly when rare, and $b < 0$ generate a strong Allee effect where populations decline when rare. The parameter $k_a$ scales the critical population size below which the Allee effect is strongest.

2.2.4.   Healthy Cells (HC)

As our first example of an additional dimension $X$, we consider interactions with healthy cells, denoted by $H$.

In this simple model, healthy cells are distinguished by their reduced growth and lack of sensitivity to treatment.

2.2.5.   Lottery model with cancer growth (LM)

Our other models include carrying capacities, which is unrealistic for cancers. To address this, we create an extended lottery model of competition for sites with healthy cells ^[22]. Assume that healthy tissue has $K$ sites, with $H$ occupied by healthy cells and the rest $E$ empty. The healthy cells turn over at rate $\delta_H$, and replicate at rate $r_H$ but only into empty locations. Then

This model maintains a constant number of sites $E + H = K$, and we think of the equilibrium healthy cell population as corresponding to a physiological optimum.

Cancer cells can differ from healthy cells in several ways: they may replicate more quickly, reproduce into sites occupied by healthy cells, and reproduce into sites occupied by other cancer cells and increase the total cell population. With two cancer cell types $S$ and $R$, we have a total cell population of $N = H+S+R$ and assume that $E = K-N$ if $N < K$ and $E = 0$ otherwise. The cells follow the equations

The competition coefficients $a$ describe the probability that a cancer cell reproduces into an occupied site. If that site is occupied by a cancer cell, we assume that the total cell population will increase, and the parameter $\eta$ represents the probability that a cancer cell kills a healthy cell that it overgrows.

2.2.6.   Immune response (IC)

A simple model of the immune system obeys

We assume that immune cells are induced to replicate by presence of cancer cells, and with growth limited by a carrying capacity. Immune cells directly kill cancer cells at rates $\eta_S$ and $\eta_R$. The replication and death rates of immune cells could be altered by treatment or by the effect of treatment in priming the immune response, although we do not address these factors here.

2.3.   Models with androgen dynamics

We consider a range of models that include androgen-dependent growth by sensitive cells, and study androgen deprivation therapies that reduce the supply of this resource.

2.3.1.   Resource competition (RC)

We adapt the basic resource competition model from population biology ^[23] by treating androgen as a resource. Upon activation, androgen receptor is translocated to the nucleus, and we take the subsequent chemical transformations that occur within the cell to mean that any androgen that is used is destroyed in the process ^[24,25]. The model tracks the two cell types and the androgen level $A$ with the equations

where we model competition as in the basic game theory model (GT). Growth rates depend on androgen levels, and androgen is supplied externally at rate $\sigma_A$ that is reduced by treatment, and used by other cells at rate $\delta_A$. For simplicity, we assume non-saturating per cell absorption rates, but saturating growth

The parameters $\epsilon_S$ and $\epsilon_R$ describe growth efficiency of the two types, and $k_S$ and $k_R$ the half-saturation constants of growth as a function of androgen uptake. Androgen deprivation treatment reduces the supply rate of androgen, and has a larger effect on susceptible cells if $k_S > k_R$ or $c_S < c_R$.

2.3.2.   Androgen Dynamics (A3)

These more detailed models include androgen-dependent ($S$), androgen-producing ($P$), and androgen-independent cells ($R$) along with explicit dynamics of production and use. The dynamics follow

where the division rates of $S$ and $P$ cells depend on their intracellular androgen concentrations $A_S$ and $A_P$ respectively.

The androgen concentrations derive from an accounting of androgen production and diffusion

$A_E$ is the external androgen concentration, $\eta$ is diffusion into and out of cells, $\mu$ is androgen use by cells, $\rho$ is production by $P$, and $\sigma$ is residual production outside the prostate. The equilibrium of the androgen system is

We assume that the dynamics of equation 2.15 are sufficiently fast to use these equilibrium values in equation 2.14.

2.3.3.   Androgen model without androgen-producing cells (A2)

For comparison with the two-dimensional models, we simplify the system to exclude testosterone-producing cells by setting $P = 0$ in equation 2.14.

2.4.   Parameter values

We choose parameter values to hit three targets.

1. Resistant cells do not invade without therapy

2. Sensitive cells grow to about 10,000 without therapy

3. Resistant cells invade with MTD therapy

We set initial conditions to $S(0) = 1000$ and $R(0) = 1.0 \times 10^{-10}$. The threshold for $R$ cell invasion is set to $R_{crit} = 1.0 \times 10^{-4}$ and for total cell numbers to $C_{crit} = 0.5 M_{max}$ where $M_{max}$ is the maximum PSA level (equal to the total cancer cell population) that occurs in 10,000 days in the absence of treatment. We assume that treatment increases the death rate of sensitive cells by a factor of 50 for the models without resources to match the strong effect of treatment in the original model ^[5] except for the Allee effect model where we use a factor of 10 to avoid driving the population below the threshold too quickly. For the models with explicit resources (RC, A2 and A3) we reduce resource availability by 90% as an upper bound of observed effects ^[26,27]. For treatment, we run a range of intermittent therapies with periods $t_P$ ranging from 100 to 1000 days and treatment duration $t_D$ running from 10 to 400 days, constrained to $t_D < t_P$. Treatment begins at time $t_P - t_D$. To implement adaptive therapy, we compare the levels of PSA to two critical levels. Therapy turns on when the total cell population increases above a fraction $M_{hi}$ of $M_{max}$ ranging from 0.2 to 0.9, and turns off when the total cell population decreases below a fraction $M_{lo}$ of $M_{max}$ ranging from 0.1 to 0.8.

2.5.   Derivation of the tradeoff curve

Consider the basic model with two competing types,

where $u$ represents the level of drug. If $R$ is rare, then it will have a negligible effect on $S$ which will follow its own dynamics with solution $S(t)$. $R$ will obey the linear equation

Integrating,

where $\bar{S}$ is the mean of $S$ from time 0 to $t$. We solve for $T_R$ with $R(T_R) = R_{crit}$, or

Because the models differ in their details, we include this relationship by fitting the parameters $a$ and $b$ to the predicted linear relationship of $\bar{S}$ with $1/T_R$,

2.6.   Computational methods

We solve models with the package deSolve in R ^[28] as written with one exception. In equation 2.14, there are two regulation terms based on androgen and carrying capacity, and the models behave pathologically when both are negative. In this case, we use only the androgen-based growth term. To simulate adaptive therapy, we include the auxiliary equation

where $M$ represents the level of the marker, like PSA, which is set equal to the total number of cancer cells. Therapy is on when $U > \epsilon$ and off when $U < \epsilon$. $U$ increases when the marker $M$ is above the threshold to turn on (the fraction $M_{hi}$ of the maximum value $M_{max}$) and decreases otherwise. We adjusted the parameter values to $r_U = 20.0$, $\epsilon = 1.0 \times 10^{-8}$ and $U_{crit} = 10 \epsilon$ in order to have therapy remain on until the value of $U$ decreases below a fraction $M_{lo}$ of $M_{max}$.

3.   Results

We simulate adaptive and intermittent therapy, including as special cases MTD and No Therapy, and record $T_R$, the time when R cells emerge ($R(t) > R_{crit}$) and $T_C$, the time when the total cancer cell population exceeds its threshold of $C_{crit}$. At each of these times, we record the average cell population until that time, and the fraction of time under treatment.

Intermittent therapy and adaptive therapy each have two hyperparameters that control the timing of treatment. For intermittent therapy, they are treatment period and duration, while for adaptive therapy they are $M_{hi}$ and $M_{lo}$, the fractions of the maximum PSA where treatment turns on and off respectively. By varying these hyperparameters, we obtain the tradeoff curves relating time to escape, cell burden, and treatment burden.

To illustrate the dynamics, we show solutions of the equations for three of the models with representative parameters (Figure 1). In both the full and simplified Zhang models (Zh and Zs), MTD leads to the most rapid decline of sensitive cells and escape of resistant cells (blue curves), No Therapy leads to rapid increase of sensitive cells with no escape of resistant cells (red curves), and the greatest delay in emergence of resistance depending on the parameter choices. For the full range of hyperparameters, both intermittent and adaptive therapy generate repeated oscillations of sensitive cells that are eventually invaded and replaced by resistant cells. However, with Allee effect, the results are nearly the opposite. MTD drives the cancer cell population below the threshold and eliminates both sensitive and resistant cells. Resistance emerges most quickly with this choice of adaptive therapy, and intermittent therapy maintains both cell populations in a long-term oscillation.

To summarize across the hyperparameters, we illustrate the relationships of time to escape of resistant and total cancer cells with mean cancer cell burden and mean tumor burden until that time beginning with the Zh model (Figure 2). Due to the simple linear dynamics of resistant cells when rare, intermittent and adaptive therapy follow the same tradeoff curves in relation to the time of escape of resistant cells (equation 2.19). MTD gives the endpoint of this curve, with the most rapid escape time, minimum cell burden and maximum treatment. The time to escape of the total cell population follows a more complex relationship, with some intermittent therapies acting much like No Therapy, with rapid cell population escape and a relatively high tumor burden. Adaptive therapy deviates only slightly from intermittent therapies, with a slight benefit of delaying escape at the cost of higher cell populations.

The two-dimensional extensions follow similar dynamics with the exception of the Allee effect model AL (Figure 3). In the Zs, GT and LV models, the mean cancer cell burden follows the tradeoff curve (equation 2.19) for both intermittent and adaptive therapy. Overall, the two therapies behave similarly, although adaptive therapy can delay resistance at the cost of a slightly higher tumor burden. With the Allee effect, the results are quite different. MTD drives cells below the threshold, and prevents both resistant cells and total cells from reaching their thresholds. Adaptive therapy can behave quite poorly, leading to escape times nearly as short as those with No Therapy and with a high total cell burden.

In Figure 4, we illustrate the results of the HC, LM, and IC models for which total cells are comprised of an additional cell type: healthy (HC, LM) or immune (IC). Despite their greater complexity, the results from these models also closely follow the predicted tradeoff (equation 2.19), although with a greater deviation for the lottery model (LM) that lacks a carrying capacity. As before, adaptive therapy largely overlaps with intermittent therapy, deviating in producing longer times to total cell escape at the expense of greater tumor cell burden, which is quite large for LM and IC models.

The models with androgen dynamics, RC, A2, and A3, follow broadly similar patterns (Figure 5). The relationship of mean cancer burden with time to escape of resistant cells persists robustly. As before, adaptive therapy can deviate from intermittent therapy when considering time to escape of the total cancer cell number, with a more complex structure with two distinct branches evident with intermittent therapies that use a low level of treatment but not for any value of adaptive therapy.

With the exception of the success of MTD with a strong Allee effect, no universal therapy can achieve all three objectives of lowering mean treatment, delaying time to emergence of resistant cells, and reducing total tumor burden. All strategies follow the relationship of mean cancer burden with time to emergence of resistant cells, and the deviations between intermittent and adaptive therapies are relatively minor for time to escape of total cancer cells. The inclusion of androgen dynamics has the largest effect on whether intermittent, adaptive or MTD therapy has higher tumor cell burden with an earlier time to emergence of resistant cells.

4.   Discussion

Adaptive therapy is based on three key assumptions: resistance is costly, resistant cells can be suppressed by competition with sensitive cells, and therapy is effective in reducing the population of sensitive cells. Without additional interactions, the only factor impeding resistant cells is competition with sensitive cells, leading to a tradeoff between reducing overall cancer burden and delaying emergence of resistance. We use a suite of simple models to find whether any of three therapy strategies can break this tradeoff: 1) constant therapy with maximum tolerable dose (MTD), 2) intermittent therapy on a fixed schedule, and 3) adaptive therapy on a patient specific schedule. We seek to test whether adaptive therapy is a feasible way to reduce overtreatment.

We have three key results. First, with the exception of models that include a strong Allee effect, all models closely follow a tradeoff curve between cancer cell burden and time to emergence of resistant cells, due to their suppression by sensitive cells. The Allee effect, meaning that cancer cell populations decline if they fall below some threshold, adds an additional type of control of resistant cells, and breaks the relationship. Second, again with the exception of models with an Allee effect, the tradeoff among time to total cancer cell number escape, average cancer cell burden, and total treatment, is similar but not identical over a range of intermittent and adaptive therapies. In most cases, some adaptive therapies do delay tumor growth, but at the cost of higher cell populations. With the Allee effect, some adaptive therapies perform quite poorly because the threshold for stopping therapy (a fraction $M_{lo}$ of $M_{max}$) is above the Allee threshold. Third, and most importantly, no therapeutic choice robustly breaks the three-way tradeoff among delaying emergence of resistance, delaying cancer growth, and minimizing treatment.

These simple models leave out many factors known to shape response to therapy. Most clearly, these models use ordinary differential equations that neglect spatial interactions known to be critical in shaping cell interactions and response to therapy ^[29,30,31]. The stochasticity neglected by differential equations would have the strongest effects when cell numbers are low, such as during the initial invasion of resistant cells or when populations are close to the threshold created by the Allee effect.

In addition, if resistance is induced by therapy rather than arising from mutations, resistance may be much more difficult to suppress ^[32]. Reversibility of these responses can create complex responses to therapeutic timing ^[33]. Much cancer resistance derives from phenotypic plasticity ^[34] that leads to more rapid emergence of resistance than the dynamics used here, and would require a different set of objective functions to evaluate.

Heterogeneity of cells with more than two states ^[35] can alter responses ^[30], and these states can be induced by a variety of intracellular changes including ABC transporter upregulation ^[36] and the evolution of mutation and genetic instability ^[37]. Androgen dynamic models begin to link intracellular with tissue-level models, but do not include specifics of pharmacokinetics and pharmacodynamics that generate differences among cells ^[36].

Our model of the immune system is highly simplified, and more realistic interactions could include the positive effects of the immune system on tumors ^[38], and additional thresholds, such as a lack of immune response to tumors below a particular size or an inability of the immune system to suppress tumors larger than some upper bound for control ^[39].

Use of multiple therapies, such as cytotoxic and cytostatic ^[40] can help create a double bind, such as between chemotherapeutic agents and those that attack glycolysis in hypoxic tumors ^[41] or force cancer into cycles of futile evolution ^[42]. The original paper ^[5] has been extended to include cells that are testosterone-independent and resistant to the chemotherapy docetaxel ^[43]. Oncolytic viruses interact with cells and the immune system with feedbacks that have only begun to be modeled ^[44]. Anti-angiogenic drugs have complex interaction with the dynamics of vasculature ^[14]. A modeling approach based on non-small cell lung cancer ^[45] that includes glycolytic cells and vascular overproducers, shows that adaptive therapy targeting glycolytic cells has potential to be effective ^[46].

We here examine only a few preset strategies, and optimization approaches could greatly refine these. The original model ^[5] itself has been studied this way ^[21]. Alternative models of prostate cancer ^[33] have been studied to retrospectively compute optimal treatments for 150 patients ^[47]. Using an alternative competitive framework with Gompertz growth, optimization of the dose and timing of treatment can maintain tumors below a tolerable size ^[48].

Optimization has been applied in many more complicated models, such as those including vasculature and the immune system and anti-angiogenic drugs and immunotherapy ^[49], and in models that include a continuum of internal resistance states and a population of healthy cells ^[40]. Building on this framework, a simplified two-state model identified strategies that give the full dose, then a smaller dose, and then a zero dose ^[50]. In a comparison of different treatment goals, a model of sensitive, damaged and resistant cells proposed to treat early to minimize integrated cell numbers, but early and late to minimize total cancer cell numbers at a fixed terminal time ^[51]. In a detailed model of colon cancer, ^[52] consider various dosing schedules to address the case where cancer cells can be in a quiescent state that is released by treatment of normal cells.

Any optimization requires choosing an objective function in a single ordered currency. Combining the measures used here (cancer cell burden, emergence of resistance, and costs of therapy) into a single currency would require sufficient clinical information to combine these into survivorship or quality-of-life adjusted years ^[53].

Applying even the simpler adaptive therapies here requires fitting to data on individual patients. A comparison of a simple model ^[33], a more complex model with basic androgen dynamics ^[20], and a detailed model of androgen dynamics ^[19] found that all fit data reasonably well, although with some exceptions ^[54]. Whether PSA dynamics alone are sufficient to resolve differences across patients, in particular patients who have different types of emerging resistance, seems unlikely, and methods almost certainly will need to be complemented with sequencing data ^[55]. If the models can be fit to the dynamics, adaptive therapies may be more robust to patient variability than prescribed timing of intermittent therapy.

Our evaluation shows that different objectives, delaying emergence of resistance, limiting total cancer cell numbers, and minimizing treatment, are likely to be related, and that no treatment can achieve all three. To address this, the choice of treatment strategy must be based on an objective function that weights different outcomes and patient goals, some of which are typically not included explicitly in mechanistic models, such as side effects ^[36].

Future therapies will need more sophisticated approaches that take into account multiple drug effects, differences among patients including increased clearance under low doses ^[6], and thinking that anticipates cancer responses ^[56,57]. However, given the limited and noisy data we have on patients, we argue that the simple models and principles presented here will remain useful as we move toward the next generation of cancer therapy.

Acknowledgments

The authors thank the members of sLaM, Jason Griffiths, Andrea Bild, and Sandy Anderson for useful discussion. FRA acknowledges support from NIH-CSBC: U54: Combating subclonal evolution of resistant cancer phenotypes (U54 CA209978). CKB received support from the Modeling the Dynamics of Life Fund at the University of Utah and Research Experiences for Undergraduates from the Department of Mathematics at the University of Utah. RST received support from the Department of Mathematics at the University of Utah. We thank two anonymous reviewers for perceptive comments that greatly improved the manuscript.

Conflict of interest

The authors have no conflicts of interest to disclose.

Appendix A.   Parameter values and justification

Model Zh: These parameters are chosen to match those in Zhang et al ^[5], but with time in units of days. $r_S = 0.0278$, $r_P = 0.0355$, $r_R = 0.0665$, $K_S = 1.5*P(0)$, $K_P = 10000$, $K_R = 10000$, $a_{SS} = 1$, $a_{SP} = 0.7$, $a_{SR} = 0.8$, $a_{PS} = 0.4$, $a_{PP} = 1$, $a_{PR} = 0.5$, $a_{RS} = 0.6$, $a_{RP} = 0.9$, $a_{RR} = 1$, $P(0) = 300$. Effects of therapy: Reduce $K_P$ by 99% and reduce $K_S$ from $1.5 K_P$ to $0.5 K_P$.

Model Zs: These match those of Model Zh without the P cells. $r_S = 0.0278$, $r_R = 0.0665$, $K_S = 15000$, $K_R = 10000$, $a_{SS} = 1$, $a_{SR} = 0.8$, $a_{RS} = 0.6$, $a_{RR} = 1$. Effects of therapy: Reduce $K_S$ by 99%.

Model GT: $r_S = 0.05$, $r_R = 0.04$, $\delta_S = 0.01$, $\delta_R = 0.01$, $K = 10000$. Effects of therapy: Reduce growth of $S$ cells to 0, with $R$ cells unaffected.

Model LV: $r_S = 0.0278$, $r_R = 0.0278$, $K_S = 12000$, $K_R = 10000$, $a_{SS} = 1$, $a_{SR} = 0.8$, $a_{RS} = 0.6$, $a_{RR} = 1$, $\delta_S = 0.001$, $\delta_R = 0.001$. Effects of therapy: Increase death rate of $S$ cells by a factor of 50.

Model AL: $b = -0.1$, $k_a = 5000$, $r_S = 0.015$, $r_R = 0.015$, $K_S = 12000$, $K_R = 10000$, $a_{SS} = 1$, $a_{SR} = 0.8$, $a_{RS} = 0.6$, $a_{RR} = 1$, $\delta_S = 0.001$, $\delta_R = 0.001$. Effects of therapy: Increase death rate of $S$ cells by a factor of 10.

Model HC: $r_S = 0.0278$, $r_H = 0.0139$, $r_R = 0.0555$, $K_S = 12000$, $k_H = 10000$, $K_R = 10000$, $a_{SS} = 1.0$, $a_{SH} = 0.6$, $a_{SR} = 1.0$, $a_{HS} = 0.6$, $a_{HH} = 1.0$, $a_{HR} = 0.6$, $a_{RS} = 1.0$, $a_{RH} = 0.6$, $a_{RR} = 1.0$, $\delta_S = 0.001$, $\delta_H = 0.001$, $\delta_R = 0.001$, $H(0) = 10000$. Effects of therapy: Increase death rate of $S$ cells by a factor of 50.

Model LM: $r_H = 0.0139$, $r_S = 0.0278$, $r_R = 0.05$, $K = 10000$, $H(0) = 10000$, $a_{SS} = 0.008$, $a_{SH} = 0.008$, $a_{SR} = 0.008$, $a_{RS} = 0.0$, $a_{RH} = 0.008$, $a_{RR} = 0.008 \; \delta_H = 0.001$, $\delta_S = 0.001$, $\delta_R = 0.001$, $\eta = 0.3$. Effects of therapy: Increase death rate of $S$ cells by a factor of 50.

Model IC: $c = 0.4$, $r_S = 0.0278$, $r_I = r_S$, $r_R = r_S$, $K_S = 15000$, $K_I = 1000$, $K_R = 10000$, $a_{SS} = 1$, $a_{SR} = 0.8$, $a_{RS} = 0.6$, $a_{RR} = 1$, $\delta_S = 0.001$, $\delta_I = 0.001$, $\delta_R = 0.001$, $\beta_S = 2/K_S$, $\beta_R = 2/K_S$, $\eta_S = c*((r_S-\delta_S)/K_I)$, $\eta_R = \eta_S$, $I(0) = 1000$. We give immune cells a carrying capacity of $K_I = 1000$ for convenience. We set $\beta_S = \beta_R = 2/K_S$ so that the immune system is strongly primed when cancer cell populations are large. We pick

where $c$ determines the strength of the immune response. Effects of therapy: Increase death rate of $S$ cells by a factor of 50.

Model RC: $c_S = 5e-5$, $c_R = c_S$, $\sigma_A = 50$, $\epsilon_S = 0.01$, $\epsilon_R = 0.015$, $K_S = 15000$, $K_R = 10000$, $k_S = 0.005$, $k_R = 0.1*k_S$, $\delta_S = 0.001$, $\delta_R = 0.001$, $\delta_A = 0.5$, $A(0) = 100$. Effects of therapy: Reduce the resource supply rate $\sigma_A$ by 90%.

Model A3: $r_0 = 0.00208$, $s_S = 0.0313$, $s_P = 0.00581$, $r_R = 0.00565$, $K_S = 15000$, $K_P = 10000$, $K_R = 10000$, $\delta_A = 6000$, $\rho = 10$, $\mu = 9$, $\sigma = 10000$, $\eta = 1$, $P(0) = 1.0$. Effects of therapy: Reduce the resource supply rate $\sigma$ and resource production rate $\rho$ by 90%. We present the details of the derivation of these values in Appendix B.

Model A2: $r0 = 0.00208$, $sS = 0.0313$, $r_R = 0.00565$, $K_S = 15000$, $K_R = 10000$, $\delta_A = 6000$, $\rho = 10$, $\mu = 9$, $\sigma = 10000$, $\eta = 1$. Effects of therapy: Reduce the androgen supply rate $\sigma$ by 90%.

Appendix B.   Calculation of parameter values for Model A3

To find the parameter values that roughly match our targets, we aim to have

● $S$ and $P$ can coexist with total population of $1.0 \times 10^{4}$,

● About 90% of androgen comes from $P$ cells, and is used internally and fairly quickly,

● With treatment that reduces $\sigma$ and $\rho$ by 90%, populations decline quickly.

The key idea is that the growth functions take the form

and must satisfy the conditions that $r_S(A_S^*) = \tilde{r_{S}}$ and $r_P(A_P^*) = \tilde{r_{P}}$ taken as the values from Model Zh. To find $r_0$, we choose a value of $A_{S, \rm{crit}} < A_S^*$ which is point of zero growth. We then find

We can then find $A_{P, \rm{crit}}$ by assuming that $S$ and $P$ cells have the same value of $r_0$ as

The equilibrium equations for $S$ and $P$ take the form

To coexist, we must have one type regulated by the carrying capacity and the other by androgen. Because $K_P < K_S$, $P$ must be regulated by carrying capacity, so $S+P = K_P$. Then $S$ will be regulated by androgen, meaning that $A_S = A_{S, \rm{crit}}$ at equilibrium. We have that

by setting $\delta_A = 0$ and simplifying by using $S+P = K_P$. This will have a positive root for $P^*$ as long as $A_{S, \rm{crit}}$ is sufficiently large. Otherwise, $P$ will be excluded from the system because external production alone can maintain $S$ at a high enough level to exclude it. The root is

We suppose $\eta = 1.0$ and $\mu = 9.0$, meaning that most androgen is used internally and fairly quickly. Results are not sensitive to $\delta_A$ and we pick $\delta_A = 6000.0$ for convenience. The value of $\rho$ can be scaled out and we set $\rho = 10.0$ for convenience. For 90% of androgen to come from $P$ cells, $\sigma = 0.1 \cdot 0.1 \rho P$. To find equilibrium androgen values, we need equilibrium values of $S$ and $P$. In the Zhang model with $a_{SP} = a_{PS} = 1$, we have $P^* = 2 K_P/3 = 20000/3$ and $S^* = K_P/3 = 10000/3$. By picking $\sigma = 10000$, and $A_E^* = 0.511$, $A_S^* = 0.0511$, $A_P^* = 1.0511$. We assume that $K_S = 15000$, and $K_P = K_R = 10000$.