Towards Multidrug Adaptive Therapy

Dobzhansky's now-famous quote that “nothing in biology makes sense except in the light of evolution” succinctly explains a worldview that has been widely adopted by the cancer biology community (1). Taken one step further, others have claimed that “nothing in evolution makes sense except in the light of ecology,” which provided the basis for designing adaptive cancer therapies centered on principles from evolution and ecology (2, 3, 4).

Cancer is an evolutionary and ecological process (5, 6) driven by random mutations (7, 8) responsible for the genetic diversity and heterogeneity that typically arises via waves of clonal and subclonal expansions (9, 10). Clones and subclones compete and Darwinian selection favors highly proliferative cell phenotypes, which in turn drive rapid tumor growth (5, 6).

Recent emphasis on personalized medicine has largely focused on the development of therapies that target specific mutations. These targeted therapies do extend patient lives but cancer cells tend to evolve resistance within months or years (11, 12). Prior to therapy, preexisting resistant cell types are suppressed and kept in check by competitively superior, therapy-sensitive cell types. There is some evidence of “cost” to incurring resistant mutations. In one study, cells sensitive (MCF7) and resistant (MCF7Dox) to doxorubicin cocultured in vitro showed that sensitive MCF7 cells rapidly outcompeted the resistant MCF7Dox line after only a few generations, illustrating the “cost” of resistance cell lines cocultured with drug-sensitive cell lines (13).

With a targeted therapy suppressing sensitive cells, these resistant cell types may experience release from competition (14). If total eradication of all cancer cells is not accomplished, the tumor will relapse derived from resistant cells that survived initial therapy (15, 16). Upon relapse, a second drug may be administered. Yet continuous use of this subsequent targeted therapy may inevitably result in the emergence of the corresponding resistant clones (Fig. 1A). This approach ignores considerations of heterogeneity and therapy as a selection event in somatic evolution (17, 18, 9).

Historically, the problem of drug resistance has often enlisted the help of mathematical modeling in designing optimal therapy schedules. For example, Goldie and Coldman were the first to propose a mathematical model relating drug sensitivity of tumors to their mutation rates (19). The model had two clinical implications: first, smaller tumors are more likely to be curable (without preexisting resistance) and second, as many effective drugs as possible should be applied as soon as possible. In this article, we will explore sequential and concomitant therapy, with the goal of designing patient-specific therapeutic schedules, which drive tumor evolution into cycles, explained below.

Eradicating most disseminated cancers may be impossible, undermining the typical treatment goal of killing as many tumor cells as possible (20). Previous schools of thought saw maximum cell-kill as either leading to a cure or, at worst, maximizing life extension. Attempting to kill as many cancer cells as quickly as possible may be evolutionarily unsound and facilitate resistance evolution and loss of therapy efficacy. Evolution matters. A recent (2012) systematic literature analysis of cancer relapse and therapeutic research showed that while evolutionary terms rarely appeared in articles studying therapeutic relapse before 1980 (< 1%), the use of evolutionary terms has steadily increased more recently, due to the potential benefits of studying therapeutic relapse from an evolutionary perspective (21).

A new paradigm in the war on cancer replaces the “treatment-for-cure” strategy with “treatment-for-contain”—receiving cues from agriculturists who have similarly abandoned the goal of complete eradication of pests in favor of more limited and strategic application of insecticides for control (22). This ecoevolutionary inspired paradigm for cancer treatment known as “adaptive therapy” capitalizes on subclonal competitive interactions. Resistance may confer some fitness costs due to increased rates of DNA repair or other costly activities required to pump toxic drugs across cell membranes. Cancer cell resistance mechanisms, whether mitigation, detoxification, or re-routing metabolic pathways, divert finite resources that would otherwise be available for cell proliferation or other avenues for cell survival (20, 23, 14).

The goal of adaptive therapy is to maintain a controllable stable tumor burden by allowing a significant population of treatment-sensitive cells to survive (see Fig. 1B, blue). These readily treatable sensitive cells serve to suppress the proliferation of the less-fit resistant populations (see Fig. 1B, red). Adaptive therapies have been tested experimentally (24, 24, 25) and are currently being applied across multiple clinical trials (NCT02415621; NCT03511196; NCT03543969; NCT03630120) at the Moffitt Cancer Center (Tampa, FL; ref. 26). These adaptive therapies capitalize on competition for space and resources between drug-sensitive and slow-growing drug-resistant populations (27, 13, 28). These trials have shown initial promise in prostate cancer compared with contemporaneous patient cohorts (26), but it should be noted that evolutionary modeling-based dosing schedule did not improve progression-free survival in anticancer TKI regimens for patients with EGFR-mutant lung cancers (29).

This adaptive approach means that each patient's treatment is truly personalized on the basis of the tumor's state and response rather than a one-size-fits-all fixed treatment regime (30). There remain several open challenges in designing adaptive therapies. First, treatments can aim to steer and control the ecoevolutionary dynamics to where the tumor finds itself in an evolutionary dead end (31), an evolutionary double-bind (32, 33, 34), or evolutionarily stable control (13, 35). Yet, for clinical practice, how to design such therapies remains difficult (36). Second, it is not yet clear how to extend these evolutionarily enlightened therapeutic concepts to multiple treatments. Two schematic examples are illustrated in Fig. 1C and D. Is it evolutionarily optimal to reproduce the single-drug adaptive therapy in a sequential (Fig. 1C) setting or in a concomitant setting (Fig. 1D)? Synergizing treatments such that evolving resistance to one drug makes cells more susceptible to another requires mathematical modeling as well as improved monitoring methods (37, 38).

Figure 2 provides a schematic of steering the tumor into “cycles” of tumor evolution. A cycle is defined as a treatment regimen that steers the tumor into periodic and repeatable temporal dynamics of tumor composition. As seen in Fig. 2, a weekly biopsy shows the frequencies (pie charts) of four phenotypes (blue, green, red, and yellow) in the evolving tumor over time. At each time step, the state of the tumor is given by a vector, x, indicating the frequency of each cell type composing the tumor. In this example, the fifth week produces a tumor composition, which is approximately equivalent to the tumor composition at the start of therapy (x₅ = x₁). Theoretically, such cycles could be repeated ad infinitum to steer and trap tumor evolution in a repeatable (and controllable) cycle.

A trial may use sequences of n available drugs either alone or in combination (shown for n = 2 in Fig. 2). We employ the terminology “treatment” to indicate the 2ⁿ possible combinations: no drug, single drug, or combination therapy. These treatments can be administered in any arbitrary sequence with the goal of controlling m cell types. An adaptive trial for metastatic castrate-resistant prostate cancer (NCT02415621) uses only two treatments: (i) Lupron and (ii) Lupron + abiraterone. Likewise, an adaptive trial for advanced BRAF-mutant melanoma (NCT03543969) administers vemurafenib and cobimetinib in combination, followed by no treatment. Both trials use only two treatments out of the four (22) combinations possible with two drugs (no treatment; first drug only; second drug only; first and second drug in combination therapy). Opening up trial design to include the full range of complexity (i.e., 2ⁿ) may allow for greater tumor control, but the treatment administered at each clinical decision time point must be chosen with care and forethought, to steer the tumor into a desirable evolutionary state. To that end, the purpose of this study is to introduce the following three concepts to consider in the pursuit of designing multidrug adaptive therapies:

1. Frequency-dependent “cycles” of tumor evolution

2. Treatment-dependent evolutionary “absorbing region”

3. Frequency-dependent evolutionary velocities

In the next sections, we will illustrate these three concepts through a case study of 4 patients with metastatic castrate-resistant prostate cancer. This case study develops these generalizable concepts for designing multidrug adaptive therapy treatment schedules for three cell types (m = 3 phenotypes; Fig. 2) under two drug treatments (n = 2 drugs; Fig. 2).

Evolutionary game theory is a mathematical framework that models frequency-dependent selection for strategies (phenotypes) among competing individuals. Competition between individuals is typically characterized by a “payoff matrix,” which defines the fitness of an individual based upon interactions with another individual or the population at large (39, 40). As a game, the payoff to an individual depends both on its strategy and the strategies of others in the population. As an evolutionary game, the payoffs to individuals possessing a particular strategy influences the changes in that strategy's frequency. A strategy that receives a higher than population-wide average payoff will increase in frequency at the expense of strategies with lower than average payoffs. Such frequency-dependent mathematical models have shown success in modeling competitive release in cancer treatment (14), designing optimal cancer treatment (41, 42), evolutionary double binds (33), glioblastoma progression (43), the emergence of invasiveness in cancer (44) as well as in cocultures of alectinib-sensitive and alectinib-resistant non–small cell lung cancer (45).

In principle, there exist many frequency-dependent models of cell–cell competition, which could adequately characterize this system: replicator dynamics models, stochastic Moran process models, spatially explicit game theoretic representations, even normalized population dynamics models. Here, we simply require a model that analyzes trajectories of relative population sizes rather than absolute population sizes of m cell types under treatment from combinations of n drugs. For details on the specific implementation and parameterization of the model used in Fig. 3A–D, see Supplementary Information and refs. 26, 46.

The model presented below is a simplified frequency-dependent dynamics mathematical model (a qualitative extension of the model behind the first adaptive therapy clinical trial in metastatic prostate cancer; see ref 26). Using the simplifying assumption of a (relatively) constant tumor volume allows us to focus on the frequency dynamics within the tumor, ignoring population dynamics. Tracking frequency dynamics that are themselves frequency dependent allows us to use a game theoretic modeling framework.

We can use the 3 by 3 payoff matrix that describes the outcomes of interactions between the different cell types. The expected payoff to an individual of a given cell type is influenced by the frequency of cell types in the population (Eq. A). This can be thought of as the “inner game” (47). The replicator equation from game theory, then translates these strategy-specific payoffs into the evolutionary dynamics described by changes in the frequencies of each cell type within the tumor (Eq. B). This is the “outer game” where payoffs become translated into fitness.

where B = 1 - A, noted in Supplementary Information (26). The variables x₁, x₂, x₃ are the corresponding frequency of dependent (T⁺), producers (TP) and independent (T⁻) cells, respectively, such that ∑_i x_i = 1. The prevalence of each cell type changes over time according to the changing payoffs, f_i, compared with the average payoff of all three populations <f> = ∑_if_ix_i. Unless otherwise noted parameters for untreated dynamics are: b_1,2 = 0.2; b_1,3 = 0.6; b_2,1 = 0.3; b_2,3 = 0.5; b_3,1 = 0.4; b_3,2 = 0.1, for Lupron only are: b_1,2 = 0.4; b_1,3 = 0.3; b_2,1 = 0.6; b_2,3 = 0.5; b_3,1 = 0.2; b_3,2 = 0.1, and for Lupron + abiraterone are: b_1,2 = 0.5; b_1,3 = 0.1; b_2,1 = 0.6; b_2,3 = 0.2; b_3,1 = 0.4; b_3,2 = 0.3.

The expected payoff to a cell type is calculated as the product of the payoff entries and prevalence of each population (Eq. B). Treatments are assumed to alter the carrying capacity for each cell type. Under Lupron, each producer cell is assumed to support 1.5 dependent cells, limiting the T⁺ to rely on producers in the absence of systemic testosterone. During Lupron + abiraterone treatment, each producer is capable of supporting 0.5 T⁺ cells and the carrying capacity of TP cells significantly drops due to local antiandrogen effects.

To simulate therapy, we introduced a weighting term, w_i. The weighting term adjusts payoffs to a cell type by its carrying capacity (capacity of the tumor environment to support a given cell type). each cell type's weighting term equals a given cell type's carrying capacity, K_i, normalized by a maximum carrying capacity (w_i = K_i/K_max). The following parameters are used: no treatment: K₁ = 1.5 × 10⁴, K₂ = 10⁴, K₃ = 10²; Lupron: K₁ = 1.5 × 10⁴; K₂ = 10⁴; K₃ = 10⁴; Lupron + abiraterone: K₁ = 0.5 × 10⁴; K₂ = 10²; K₃ = 10⁴; K_max = 1.5 ×10⁴ cells (26).

Each subpopulation is assumed to contribute to changes in prostate-specific-antigen (PSA) levels, as follows:

The PSA is normalized over time to show changes relative to treatment initialization (see Fig. 3). Data was obtained via the clinical trial protocol (NCT02415621), which was approved by central Institutional Review Board and monitored by Moffitt Cancer Center's protocol monitoring committee. Written informed consent was obtained from all patients prior to enrollment in the trial. Best fits were performed using lsqcurvefit function in MATLAB 2018a. All parameters were held constant between patients (see Supplementary information) except for the initial conditions.

In progressive prostate cancer, continuous or intermittent androgen deprivation are common treatments that generally lead to a substantial decline in tumor burden. Eventually, tumor burden rebounds as a result of the rise of castrate-resistant cells (48). Upon castrate resistance, adaptive therapy using abiraterone has shown considerable treatment success.

In Fig. 3, we present data of castrate-resistant patients from NCT02415621. The 4 patients (Fig. 3, panels A–D) were selected who have undergone at least five adaptive doses of treatment at the time data was acquired (26, 36). The adaptive protocol is as follows: Lupron is administered continuously (limiting systemic testosterone production) while abiraterone (antiandrogen) is administered only until the PSA blood biomarker declines to 50% of the pretreatment value. abiraterone treatment is resumed after the PSA climbs to the original pretreatment value.

Previous studies on controlling resistance in cancer therapy have emphasized the value of continuous monitoring of evolving populations (37) and the importance of developing a resistance management plan (49, 50, 51). Here, we perform an “after-action analysis” on all patients who progress due to resistance (52, 53). This after-action reporting aids understanding of the mechanisms of treatment failure, as well as identifies improved resistance management strategies.

To facilitate this after-action analysis, we parameterize a previously published Lotka–Volterra model of treatment dynamics by fitting to PSA data from all patients (Fig. 3, dashed red lines). This population dynamics model has three cell types: those that require testosterone for growth, T⁺, cells that require testosterone but produce their own, TP, or cells that are independent of testosterone (and thus are resistant to abiraterone), T⁻ (26, 46, 36). Mathematical modeling allows for continuous prediction of tumor subpopulations (Fig. 3, pie charts).

While the ongoing adaptive trial (NCT02415621) has shown promise in extending median time to progression of patients with metastatic prostate cancer, resistance is still expected to inevitably occur. The retrospective, after-action analysis presented in Fig. 3 indicates the emergence of resistance for all 4 patients shown (see Fig. 3 embedded pie charts; see also detailed analysis shown for each patient in Supplementary Fig. S1, panels A, B, C, D).

On the basis of the model predictions developed in the after-action analysis, all 4 patients eventually fail due to resistance. This indicates that a secondary resistance management plan is still necessary to mitigate this emergence of resistance. The resistance management plan introduced here is additional periods of no treatment to mitigate the emergence of the abiraterone-resistant testosterone-independent (T⁻) subpopulation. The timing of these off-treatment periods are constructed such that the tumor subpopulations return to the original proportion and thus complete a cycle, as proposed in the above Fig. 2.

This goal of this after-action analysis is to utilize retrospective patient data to parameterize a model to test the efficacy of this evolutionary cycle design to mitigate resistance. Importantly, our after-action analysis predicts that evolutionary cycling prolongs the relapse of the resistant T⁻ subpopulation while maintaining a controllable tumor volume, for each patient (Fig. 3, blue lines; see also Supplementary Fig. S1). The three-pronged approach of continuous monitoring, after-action reporting, and resistance management planning provide a key step forward to implementing this evolutionary cycling concept in the clinic.

To help elucidate the nature of such cycles, we simplify our model system by ignoring (for the moment) population dynamics to leverage several mathematical conveniences of a frequency-dependent game theoretic modelling framework, explained below. Frequency-dependent models of tumor evolution are particularly suited for studying tumor control (i.e., “cycles” of tumor evolution), determining the set of possible evolutionary dynamics (i.e., evolutionary absorbing region), and timing of evolution (i.e., evolutionary velocity).

To track the state vector of the tumor, x, we set the variables x₁, x₂, x₃, to the frequency of dependent (T⁺), producers (TP), and independent (T⁻) cells, respectively. The temporal dynamics of these different populations under treatment can be characterized inside a trilinear coordinate simplex (Fig. 4A), which gives a representation of every possible value of x, on a triangle. The corners represent a tumor consisting of 100% of either T⁺ (top corner), TP (left corner), and T⁻ (right corner). Figure 4 shows the temporal dynamics under each treatment.

Every treatment scenario has an associated long-time equilibrium state: the tumor composition where all temporal dynamics eventually converge. For example, continuous treatment of Lupron leads to roughly equal fraction of T⁺ and TP (blue circle, Fig. 4B). Continuous treatment of Lupron + abiraterone leads to full saturation of T⁻ (blue circle, Fig. 4C). Similarly, administering no treatment leads to an equilibrium of mostly T⁺ and roughly equal (but small) proportions of T⁻ and TP (purple circle Fig. 4D).

This leads us to a key insight. We are able to ignore the tumor volume information (for the moment) because of one convenient fact: an optimal treatment schedule will necessarily be a schedule, which successfully avoids all equilibria. Equilibria associated with “no treatment” will lead to tumor saturation and presumably, death. Equilibria associated with any treatment will lead to a fully resistant tumor (and also presumably, death).

The goal of adaptive therapy in prostate cancer is to delay the onset of the resistant T⁻ population by well-timed switching between each treatment. This is equivalent to switching between each triangular phase portrait in Fig. 4 before reaching the resistant equilibrium state of any given treatment. In theory, a periodic (closed) cycle can be constructed by switching between treatments at carefully chosen times to design schedules that are superior to continuous MTD (54).

Figure 5A details the process by which an evolutionary cycle is constructed. Each patient in Fig. 3 has an associated initial condition within the triangle. A sequence of treatments is constructed for each patient, which arrives back at the same initial condition: an evolutionary cycle. This treatment sequence can be repeated, controlling the tumor. The dashed lines on Fig. 5A show dynamics under continuous treatment, which lead to respective equilibrium states, spiraling out of the evolutionary cycle. The identical evolutionary cycle paradigm is shown with cell fraction over time in Fig. 5B. By appropriate treatment switching, the three cell types remain in competition with each other, and no cell type is able to dominate but instead are balanced indefinitely in closed periodic cycles thereby avoiding the emergence of resistance. The tumor undergoes three such evolutionary cycles, where the tumor “resets” back the initial state before treatment (pie charts, Fig. 5).

Many such cycles may exist in the state space shown in Fig. 4D: for example, traveling down any Lupron trajectory (blue), switching to a no treatment trajectory (purple) and so on, repeated ad infinitum. Which cycles are preferable? To answer this question, the next section introduces concepts of evolutionary absorbing region and evolutionary velocity.

Phase portraits can be drawn for each pairwise treatment combination (Fig. 6). Each treatment has an associated equilibrium (solid circles), which can be connected by a single evolutionary trajectory. These two connecting trajectories represent a bounding region of state space where a periodic cycle resulting from sequential administration of the two treatments is guaranteed on the outer rim of this region (directionality shown with black arrows). As can be seen in Fig. 6, all external trajectories tend toward either treatment's equilibrium, or toward the inside of the bounding domain. The implications are clear: sequential treatment between any two drugs for a sufficiently long time results in limiting the “evolutionary absorbing region” of the tumor composition.

Herein lies a second key advantage of frequency-dependent game theoretic models: parameterization robustness. Shown in the inset of each subfigure (Fig. 6), are 36 random parameterizations (subject to the inequalities in Supplementary Tables S1 and S2), overlaid in transparent shading. Absorbing regions overlap significantly due to the fact that model dynamics are more sensitive to relative parameter values (i.e., the inequalities in Supplementary Tables S1 and S2) than absolute values chosen.

At least one cycle exists for each pairwise treatment scenario: the outer rim of the absorbing state. There may be other cycles, existing only in the interior of the region. For some pairwise combinations, the bounding region is quite large (see Fig. 6A: Lupron, Lupron + abiraterone); for others quite small (see Fig. 6B: no treatment, Lupron or C: no treatment, Lupron + abiraterone). Once treatment drives the tumor composition into one of these bounding regions, it is impossible to traverse outside without the addition of a third treatment.

Considering the three treatment scenarios together (see Fig. 6D: no treatment, Lupron, Lupron + abiraterone) expands the evolutionary absorbing region. Arbitrary cycles can be constructed by clever sequences of treatment inside this expanded evolutionary absorbing region. To illustrate which regions in which cycles may be preferable, the next section will describe evolutionary velocity.

The evolutionary “velocity” (and each of its velocity components) of each treatment scenario can be viewed on the same trilinear simplex, shown for total velocity magnitude (Fig. 7A–C), T⁺ velocity (Fig. 7D–F), TP velocity (Fig. 7G–I), and T⁻ velocity (Fig. 7J–L). It is useful to compare no treatment (7, left column), Lupron (7, middle column), and Lupron + abiraterone (7, right column). A low velocity indicates a slow change to the tumor composition (calculated as the magnitude of the resultant vector of all components in Eq. A (see Supplementary Information). Several recent studies have indicated the importance of studying temporal effects of drug administration. Using human breast cancer explants, in vitro cells, mouse in vivo studies combined with mathematical modeling, one study showed that a temporary exposure to a taxane induces phenotypic cell state transition toward a favored transient chemotherapy-tolerant state, which disappears when the drugs are administered concomitantly (55). A separate study noted the appearance of weakly proliferative, resistant cells under high-dose therapy is transient and controllable by nongenetic, stem-like (reversible) characteristics that depend on timing and length of drug administration (56).

In short: timing matters. Introduction of concomitant therapy or altering the timing of sequential therapy changes the evolutionary velocity of the underlying tumor composition. One can imagine scenarios where a fast or slow evolutionary velocity is desired. A tumor with a high composition of resistant cells may need to navigate to a fast dynamics region to rapidly decrease the resistant subpopulation. Alternatively, an adaptive regime may capitalize on slow velocities regions that ensure slow evolutionary dynamics on treatment holidays.

Frequency dynamics models allow for monitoring and control of the velocity of a single-cell type (see Fig. 7D–L). For example, under continuous Lupron + abiraterone treatment, the velocity of the T⁻ population is positive for most of the state space (see Fig. 7L). To control the T⁻ population, it is necessary to switch to a new treatment: Fig. 7J or Fig. 7K. Depending on the current tumor composition, no treatment or Lupron may be desirable for negative T⁻ velocity (blue). For example, while certain adaptive evolutionary cycles may be technically feasible, clinical considerations of fast evolutionary velocities may preclude certain schedules requiring high frequency hospital visits or impractically slow wash-out times of particular treatments.

In the design of multidrug adaptive therapy, the optimal method of combining therapies in an additive or sequential manner is unclear. A key observation is that the goal of adaptive therapy is to maintain a stable tumor volume in favor of designing therapies that alter tumor composition. We advocate for the use of frequency-dependent competition models (and in particular evolutionary game theory) to design novel multidrug adaptive therapy regimens. This is a promising modeling approach, and more work must be done combining population models of changing tumor volume with frequency-dependent models of changing tumor fraction to eliminate clinically unfeasible treatment schedules from consideration.

While it is clear that tumors evolve in response to treatment, it has proved difficult to exploit evolutionary principles to steer tumor evolution in an ideal direction. One limitation of current adaptive therapy clinical trials is the limited selection of drugs and limited monitoring methods. Even despite the lack of monitoring the exact state of the tumor composition, the technique of maintaining a substantial drug-sensitive population for extended tumor control has shown promise in mathematical models as well as preclinical and clinical trials. However, each of the on-going or planned adaptive therapy clinical trial utilizes less than the full range of treatment options available. For example, prostate cancer adaptive therapy sequentially administers Lupron and both Lupron + abiraterone, while ignoring no treatment or abiraterone monotherapy. Similarly, the planned melanoma trial administers no treatment alternating with both vemurafenib and cobimetinib. The conceptual paradigms introduced here (cycles of tumor evolution, evolutionary absorbing region, and evolutionary velocity) provide a path forward to reducing the complexity of selecting between 2ⁿ treatment choices at every clinical decision point.

Another trial, currently accruing patients with castrate-sensitive prostate cancer (NCT03511196), attempts to infer underlying dynamics of tumor subpopulations by comparing two clinical biomarkers over time: PSA and testosterone levels. Patients with rising PSA levels without a respective rise in testosterone may indicate an emergent castrate-resistant population (TP). When this occurs, abiraterone is administered to counter this resistant population while leaving the serum testosterone unchanged to bolster the T⁺ population. Even without precise monitoring of underlying subpopulations, changes may be inferred using evolutionary principles.

Especially when developed in collaboration with clinicians, evolutionary models of cellular competition can provide clarity and power despite their simplicity. As a first step, we have illustrated these three adaptive paradigms in a specific case study (prostate cancer) for n = 2 drugs controlling m = 3 phenotypes (testosterone producers, testosterone-dependent, testosterone-independent cells). These techniques can be extended to any multidrug adaptive therapy setting. The first step is to carefully choose competition parameters (i.e., the payoff matrix) of each cell phenotype under consideration and draw the dynamical phase portraits of possible evolutionary trajectories implies by the set of competition parameters. The power of evolutionary game theory is that the relative fitness of each phenotype is often easy to determine, allowing for straightforward construction and analysis of the model.

Each treatment will eventually fail due to resistance. The resistant state is the absorbing stable state of tumor phenotypic composition associated with continuous treatment. Example trajectories for a given treatment can be plotted, predicting the timing to reach this absorbing resistant state. When adding a second treatment in sequence, an absorbing bounded region of phenotype space (Fig. 6) can be drawn, which we term the evolutionary “absorbing region” reachable with two or more treatments. This absorbing state space is bounded by the evolutionary cycle that connects the stable state points from each treatment (note: we consider points in the interior of the simplex and the global equilibrium only). The absorbing region may be a small region (i.e., Fig. 6B) or large (i.e., Fig. 6A). Upon combining all treatments, the state space is expanded, also leading to more options for cycles in new regions of the state space. Intuition would indicate that new drugs should be introduced that have “orthogonal” stable states—stable equilibria far from the equilibria of existing drugs. These orthogonal drugs open up the treatment space, allowing for more options in choosing cycles (Fig. 6D). The closeness of two stable states from two distinct treatments will give an idea of the “orthogonality” of two treatments, aiding treatment selection, the sequential ordering of treatments, and the timing of switching between treatments.

It is also important to note that not treating still allows the tumor to evolve. This no treatment case is often associated with slow evolution due to low selection. During treatment holidays used to avoid resistance to a first treatment, it may be wise to choose a faster (high selection) second treatment that gives a similar resultant tumor composition over the slow evolution of no treatment. Adaptive therapy is a promising step toward ecologically inspired personalized medicine. Optimizing multidrug adaptive therapy is not a straightforward task, but mathematical models are a powerful abstraction of clinical intuition, enabling the generation of new treatment schedules and comparisons with standard of care.

As indicated above, continuous monitoring of evolving populations is important when designing treatment schedules to mitigate the emergence of resistance. This monitoring is facilitated by mathematical models of treatment response combined with putative biomarkers of tumor burden. Adaptive therapy has been attempted clinically using proxy measurements for tumor volume such as PSA for prostate, lactate dehydrogenase for melanoma, thyroglobulin for differentiated thyroid cancer, or calcitonin in patients with medullary thyroid cancer.

Parameterizing a mathematical model using a tumor burden biomarker is a cost-effective method of implementing multidrug adaptive therapies in the clinic, as serial measurements of PSA and lactate dehydrogenase are relatively inexpensive. Monitoring burden biomarkers over time will indicate the efficacy of adaptive therapy, but it will not necessarily elucidate the mechanism, and tumor biopsies may be needed to directly measure tumor response to therapy. These blood biomarkers are crude measurements of tumor burden, but there have been many recent advances in reliable methods of monitoring tumor evolution using circulating tumor DNA as an informative, inherently specific, and highly sensitive biomarker of metastatic cancers (see ref. 57). The concepts introduced in this manuscript can also be applied via more sophisticated methods of measuring or inferring tumor subpopulations, as they become available.

As monitoring methods mature, so must methods of prediction and quantification of evolution (26, 58, 59, 60). After-action reporting utilizing mathematical models parameterized using retrospective clinical data provides a clear opportunity to continuously update resistance management strategies. This after-action reporting betters understanding of the mechanisms of treatment failure, as well as identifies improved strategies. These data demonstrate the feasibility of evolutionary cycles (based on patient-specific parameterization), which control the relapse of resistance. Without this cycling, patients eventually relapse due to the emergence of resistance.

J. Zhang is an employee/paid consultant for Bayer, Dendreon, Clovis Oncology, AstraZeneca and has received speakers bureau honoraria from Sanofi and Merck. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J. West, R.A. Gatenby, J.S. Brown, A.R.A. Anderson

Development of methodology: J. West, J.S. Brown, P.K. Newton, A.R.A. Anderson

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Zhang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. West, L. You, P.K. Newton, A.R.A. Anderson

Writing, review, and/or revision of the manuscript: J. West, J. Zhang, R.A. Gatenby, J.S. Brown, P.K. Newton, A.R.A. Anderson

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A.R.A. Anderson

Study supervision: A.R.A. Anderson

The authors gratefully acknowledge funding from the Physical Sciences Oncology Network (PSON) at the NCI, U54CA193489 (to J. West, J.S. Brown, and A.R.A. Anderson), the European Union's Horizon 2020 research and innovation program, under the Marie Sklodowska-Curie grant agreement No. 690817 (L. You), and the Jayne Koskinas Ted Giovanis (JKTG) Foundation for Health and Policy and the Breast Cancer Research Foundation, private foundations committed to critical funding of cancer research (to P.K. Newton).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.