Optimizing adaptive cancer therapy: dynamic programming and evolutionary game theory

TL;DR

This paper develops a systematic method using dynamic programming and evolutionary game theory to optimize adaptive cancer therapy, aiming to reduce drug usage and improve recovery rates compared to standard treatments.

Contribution

It introduces a novel framework combining dynamic programming and evolutionary game theory to optimize adaptive cancer treatment policies systematically.

Findings

Optimal strategies significantly reduce drug usage.

Adaptive policies increase likelihood of tumor recovery.

Trade-offs between drug amount and recovery time are characterized.

Abstract

Recent clinical trials have shown that the adaptive drug therapy can be more efficient than a standard MTD-based policy in treatment of cancer patients. The adaptive therapy paradigm is not based on a preset schedule; instead, the doses are administered based on the current state of tumor. But the adaptive treatment policies examined so far have been largely ad hoc. In this paper we propose a method for systematically optimizing the rules of adaptive policies based on an Evolutionary Game Theory model of cancer dynamics. Given a set of treatment objectives, we use the framework of dynamic programming to find the optimal treatment strategies. In particular, we optimize the total drug usage and time to recovery by solving a Hamilton-Jacobi-Bellman equation based on a mathematical model of tumor evolution. We compare adaptive/optimal treatment strategy with MTD-based treatment policy. We show that optimal treatment strategies can dramatically decrease the total amount of drugs prescribed as well as increase the fraction of initial tumour states from which the recovery is possible. We also examine the optimization trade-offs between the total administered drugs and recovery time. The adaptive therapy combined with optimal control theory is a promising concept in the cancer treatment and should be integrated into clinical trial design.

Figures

Equations

$$\{\begin{array}{l}{\dot{x}}_D=x_D({\psi }_D-<\psi >)\\ {\dot{x}}_G=x_G({\psi }_G-<\psi >)\\ {\dot{x}}_V=x_V({\psi }_V-<\psi >)\end{array}$$

$$\{\begin{array}{l}q=\frac{x_V}{x_V+x_D}\\ p=x_G\end{array}\text{or }\text{leavevmode}\text{leavevmode}\{\begin{array}{l}{x}_D=(1-q)(1-p)\\ x_G=p\\ x_V=(1-p)q\end{array}$$

$$\begin{cases}\dot{q}=q(1-q)\Big{(}\frac{b_{v}}{n+1}\sum\limits_{k=0}^{n}p^{k}-c\Big{)}\\ \dot{p}=p(1-p)\Big{(}\frac{b_{a}}{n+1}-(b_{v}-c)q\Big{)}\end{cases}$$

$$\frac{b_a}{n+1}<b_v-c<cn.$$

$$\{\begin{array}{ll}p(t)<r_b,& \text{if therapy succeeds;}\\ p(t)>1-f_b,& \text{if therapy fails.}\end{array}$$

$$\Delta=\Big{\{}(q,p)\in[0,1]\times[0,1]:p<r_{b}\text{ or }p>1-f_{b}\Big{\}}.$$

$$\begin{cases}\dot{q}(t)=q(t)\Big{(}1-q(t)\Big{)}\Big{(}\frac{b_{v}}{n+1}\sum\limits_{k=0}^{n}p^{k}(t)-c\Big{)},\\ \dot{p}(t)=p(t)\Big{(}1-p(t)\Big{)}\Big{(}\frac{b_{a}}{n+1}-(b_{v}-c)q(t)-\bm{d(t)}\Big{)};\\ q(0)=q_{0},\leavevmode\nobreak\ p(0)=p_{0}.\end{cases}$$

$$T(q_{0},p_{0},d(\cdot))=\min\Big{\{}t\in{\mathbb{R}}_{+}\leavevmode\nobreak\ |\leavevmode\nobreak\ \bm{(}q(t),p(t))\in\Delta,\leavevmode\nobreak\ q(0)=q_{0},\leavevmode\nobreak\ p(0)=p_{0}\Big{\}}.$$

$$g(q,p)=\{\begin{array}{ll}+\infty ,& \text{if }p>1-f_b,\\ 0,& \text{otherwise.}\end{array}$$

$$J(q_{0},p_{0},d(\cdot))=\int\limits_{0}^{T}(d(s)+\sigma)ds+g\Big{(}q\big{(}T\big{)},p\big{(}T\big{)}\Big{)}.$$

$$u(q_0,p_0)=\underset{d(\cdot )}{\inf }J(q_0,p_0,d(\cdot ))$$

$$\min\limits_{d\in[0,d_{max}]}\left\{\nabla u(q,p)\cdot\left(\begin{array}[]{c}\dot{q}(q,p,d)\\ \dot{p}(q,p,d)\\ \end{array}\right)+d+\sigma\right\}=0,\qquad(q,p)\in([0,1]\times[0,1])/\Delta.$$

$$\{\begin{array}{ll}u(q,p)=0,& \text{ if }p<r_b;\\ u(q,p)=+\infty ,& \text{ if }p>1-f_b.\end{array}$$

$$J(q_{0},p_{0},d(\cdot))\,=\,\int\limits_{0}^{T}K\left(d(s)\right)\,ds\,+\,g\Big{(}q\big{(}T\big{)},p\big{(}T\big{)}\Big{)},$$

$$u(q_{0},p_{0})\,=\,\int\limits_{0}^{\tau}K\left(d^{*}(t)\right)\,ds\,+\,u\Big{(}q(\tau),p(\tau)\Big{)}$$

$$\nabla u(q_{0},p_{0})\cdot\left(\begin{array}[]{c}\dot{q}(q_{0},p_{0},d^{*}_{0})\\ \dot{p}(q_{0},p_{0},d^{*}_{0})\\ \end{array}\right)+d^{*}_{0}+\sigma=0.$$

$$\min\limits_{d\in[0,d_{max}]}\left\{\leavevmode\nobreak\ \Big{[}1-u_{p}p(1-p)\Big{]}d\;+\;u_{q}q(1-q)\Big{(}\frac{b_{v}}{n+1}\sum\limits_{k=0}^{n}p^{k}-c\Big{)}+u_{p}p(1-p)\Big{(}\frac{b_{a}}{n+1}-q(b_{v}-c)\Big{)}+\sigma\leavevmode\nobreak\ \right\}\;=\;0.$$

$$d^{*}=\begin{cases}d_{max},&\quad\text{if }\Big{(}1-u_{p}p(1-p)\Big{)}<0;\\ 0,&\quad\text{otherwise}.\end{cases}$$

$$\begin{cases}\dot{y}(t)=f\Big{(}y(t),d(t)\Big{)}\\ y(0)=x,\end{cases}$$

$${\tilde{x}}_d=x+\tau f(x,d).$$

$$u(x)=\underset{d}{\min }\left\{\tau K(d)+u\left({\tilde{x}}_d\right)\right\}+o(\tau ).$$

$${\tilde{x}}_d=\frac{||x_{2d}-{\tilde{x}}_d||}{||x_{1d}-x_{2d}||}{x}_{1d}+\frac{||x_{1d}-{\tilde{x}}_d||}{||x_{1d}-x_{2d}||}{x}_{2d}\text{ and }{\tau }_d=\frac{||x-{\tilde{x}}_d||}{||f(x,d)||}.$$

$$U(x)=\underset{d\in \{0,d_{max}\}}{\min }\left\{\frac{||x-{\tilde{x}}_d||}{||f(x,d)||}K(d)+\frac{||x_{2d}-{\tilde{x}}_d||}{||x_{1d}-x_{2d}||}U(x_{1d})+\frac{||x_{1d}-{\tilde{x}}_d||}{||x_{1d}-x_{2d}||}U(x_{2d})\right\},$$

Full text

METHODS: Given a set of treatment objectives, we use the framework of dynamic programming to find the optimal treatment strategies. In particular, we optimize the total drug usage and time to recovery by solving a Hamilton-Jacobi-Bellman equation based on a mathematical model of tumor evolution.

RESULTS: We compare adaptive/optimal treatment strategy with MTD-based treatment policy. We show that optimal treatment strategies can dramatically decrease the total amount of drugs prescribed as well as increase the fraction of initial tumour states from which the recovery is possible. We also examine the optimization trade-offs between the total administered drugs and recovery time.

CONCLUSIONS: The adaptive therapy combined with optimal control theory is a promising concept in the cancer treatment and should be integrated into clinical trial design.

Keywords: adaptive therapy; optimal treatment policy; evolutionary game theory; Hamilton-Jacobi-Bellman equation; heterogeneity

1 Introduction

Intratumoral heterogeneity is being increasingly recognized as a cause of metastasis, progression and resistance to therapy 1. While genetic instability, a hallmark of malignancy 2, can result in this heterogeneity, it is being increasingly understood that eco-evolutionary factors, like selection and clonal interference, can also drive and maintain it3, 4.

While sequencing technologies have enabled increasingly in-depth quantitative understanding of the genetic heterogeneity, relatively little experimental work has sought to directly quantify the eco-evolutionary interactions involved. As more studies come to light showing the efficacy of treatments based on eco-evolutionary trial designs, this lack of quantification is coming into focus.

In line with standard, cell-autonomous growth-based theories, conventional chemotherapy is given to patients at the maximum tolerated doses (MTD): the highest doses that most patients can safely tolerate. Although the MTD-based chemotherapy offers advantages in survival compared to no therapy, cures remain elusive, and side effects can be severe. In addition to the toxicity, it is known that relapse is nearly inevitable due to the emergence of therapeutic resistance: a process driven by Darwinian evolutionary dynamics in which the MTD-based chemotherapy kills off the chemotherapy-sensitive cells and chemo-refractory cells eventually dominate in the tumor. While it is unknown whether these resistant cells are present before therapy, or acquire resistance mutations during therapy, it is the process of variation and selection under standard therapy that drives the inevitable failure in the patient.

Metronomic chemotherapy has been proposed as a possible alternative to the MTD strategy5, 6. Metronomic chemotherapy is given in an on/off fashion at frequent time intervals according to a set periodic schedule. The idea behind this method is to give less overall therapy, thereby increasing tolerability, and allowing time for therapy sensitive cells to regrow, allowing for resensitization of the tumor. Frustratingly, the results of clinical trials of metronomic chemotherapy have been “variable” 7 and many of them have not demonstrated significant efficacy 8, 9, 10 as compared to standard therapy. Recent studies argue that metronomic chemotherapy highly depends on timing, and the right scheduling can improve the results of its usage 11.

Based on the hypothesis that disease dynamics depend on the evolution of tumor heterogeneity as modulated by competition between subtypes, the idea of using adaptive therapy (AT) has been proposed 12. AT is much like metronomic therapy, with an important difference. AT administers doses of therapy according to the current state of tumor growth and its anticipated evolutionary changes (trajectory). These can be estimated using direct (e.g., taking biopsies) or indirect (e.g., antigen testing, mathematical modeling) methods. Therefore, unlike the MTD-based or metronomic protocols, AT does not have preset schedule and it adjusts the doses and timing before a tumor becomes chemotherapy-resistant, prolonging time to this event. Recently, the adaptive strategies have shown promise in pre-clinical trials of breast cancer 13 and a phase 2 clinical trial in metastatic castrate-resistant prostate cancer 14.

These two recent successes 13, 14 in adaptive therapy have been based on mathematical modeling of tumor evolution under therapy using a dynamical systems approach based on Evolutionary Game Theory (EGT) 15, 16. This formalism explicitly considers interactions between sub-populations and models their fitness in frequency dependent terms. EGT has been used to theoretically consider many scenarios in cancer before, including therapy scheduling and timing in prostate cancer 17, 18, 14, 19; the use of tumor microenvironment targeting therapy in glioblastoma 20; the trade-off between healthy tissue and cancer in multiple myeloma 21, 22; and drug resistance in general 23, 24, 25. These theoretical studies, combined with the recent empiric realizations, suggest significant opportunities to improve therapy by using this evolutionarily enlightened approach. Nevertheless, therapeutic decisions in general practice are currently not based on this knowledge and continue to use the MTD paradigm.

Assuming that an oncologist has perfect information about the current state of a tumor, and a faithful mathematical model that can predict its trajectory (these, of course, are very strong assumptions), it is not clear how he/she should adjust the schedule and doses. Based on a stage of the disease and patient’s needs, the therapy can have different final goals: maximization of patient’s duration of life, ensuring the best possible quality of the rest of life, decreasing probability of new metastases appearing, decreasing time/cost of the treatment, etc.. Unfortunately, an oncologist can usually only focus on one or two of these goals, having some reasonable constraints on the secondary parameters. Thus, an important step toward optimizing AT is to define an objective of the therapy and “translate” it into mathematical language. The next step is to quantify how good each particular strategy is with respect to that chosen objective. Optimizing this objective is a mathematical goal which can be addressed by the tools of optimal control theory.

Optimal control theory, a branch of mathematics typically applied to problems in engineering, can be applied to a wide class of problems arising in oncology 26. The first application of optimal control theory in cancer was done by Swan and Vincent 27 who found the optimal treatment strategy for multiple myeloma with the objective to minimize the total amount of drugs used applying the Pontryagin Minimum Principle (PMP) 28. Since then, others have used the PMP to different optimal cancer treatment problems: a chemotherapy optimization under evolving drug resistance 29, 30, 31, optimal scheduling of a vessel disruptive agent 32, MAPK inhibitors 33 input in cancer treatment, minimization amount of drugs prescribed in tumor Cimmune model 34, finding compromise between drug toxicity and tumor repression for the myeloma bone disease 35, and many others.

While these approaches have offered benefits in their ability to formally optimize problems written as dynamical systems, the PMP method has several limitations. First, PMP yields only a necessary condition for an optimum, and any locally optimal trajectory of the control system satisfies PMP. Local optimality means that the trajectory is optimal when comparing it with its small perturbations, but there may well be a different trajectory that is even better (globally optimal, compared with all possible trajectories). Secondly, PMP provides a time-dependent (open-loop) control: given an initial state, the method provides an optimal treatment strategy as a function of time – therefore a treating oncologist has to follow it regardless of the changing state of the tumor. However, if the underlying model has been perturbed or includes some noise (like a tumor acquiring mutations, say), the control cannot adapt to these unexpected changes.

A different approach to analysis of an optimal control problem is a feedback (closed-loop) control point of view. Using the Hamilton-Jacobi-Bellman (HJB) equations, one can obtain controls that depend on the current state of the dynamical system (current distribution of sub-populations of cancer cells) rather than only the current time36, 37. In this case, the treating oncologist’s decisions can be adjusted if something unexpected has happened with the trajectory. Moreover, the HJB equations guarantees that the resulting treatment feedback strategies are globally optimal. Despite these advantages of the HJB over PMP method, there are few works 38, 39 which use the feedback control paradigm to find an optimal treatment strategy.

Here, we apply the HJB approach to solve for optimal treatment strategies for a model of lung cancer proposed by Kaznatcheev and colleagues 40. In that paper the authors introduce an evolutionary game (system of replicator-type equations) that models the dynamics of three sub-populations of tumor cells. The article highlights the importance of a good scheduling in the polyclonal regime, when the game has cyclic dynamics. The article has an example of two different scheduling strategies with the same set of initial parameters that lead the system to opposite outcomes: putative recovery, versus putative death of a patient. While several qualitatively different treatment schedules are presented, optimal therapy is not discussed. Given the growing interest in EGT in clinical applications14 and recent work connecting these models using direct in vitro parameterization 41, we believe the optimization of therapies based on such models will become increasingly important and the HJB-based approach will be used far more often in the future.

2 Methods

In this article we focus on a model of cancer evolution that has been proposed in Kaznatcheev et al.40 and summarized in Box 1. This model considers interactions between three different sub-populations of cancer cells playing a modified version of the public goods social dilemma: glycoltyic cells (GLY), vacscular overproducers (VOP) and cells called defectors (DEF) which use both strategies to “cheat” on the others. GLY cells are anaerobic and produce lactic acid. Both VOP and DEF cells are aerobic. In addition, VOP cells spend extra energy to produce VEGF (a protein that improves the vasculature, benefiting both VOP and DEF). Based on the replicator model from EGT15, 16 summarized in equations (1) and (3), the evolution of the tumor is described by tracking the changing proportions of GLY, VOP and DEF cells in the full population. The patient is viewed as recovered when the GLY proportion falls below some low threshold $r_{b}$. (Below this recovery barrier, the validity of replicator-based model is harder to justify and we assume that the GLY cells are essentially extinct.) Conversely, we assume that GLY cells suppress other tumor cells and a patient dies if the total proportion of aerobic cells (VOP and DEF sub-populations combined) falls below some low threshold (a failure barrier) $f_{b}$.

For a range of parameter values (4), this model predicts a heterogeneous regime†††Other tumor regimes (fully angiogenic and glycolyctic) also exist outside of this parameter range. They are less interesting from the point of view of treatment strategies, but we still consider them for the sake of completeness in Section 6.3 of Supplementary Materials. in cancer evolution with coexistent and oscillating proportions of GLY, VOP, and DEF. Without any treatment, these sub-populations follow cyclic dynamics and a patient never recovers; see Figure 1.

Following Section 4.1. in Kaznatcheev et al. 40 we consider a cell-type-targeting therapy that preferentially penalizes the fitness of GLY cells; see formula (7) in Box 2. During the treatment, a doctor defines the timing of therapy and its intensity. This time-dependent intensity $d(t)$ can vary between [math] (no therapy) and $d_{max}>0$ (the maximum tolerated dose, MTD). Two extreme case ($d(t)=0$ versus $d(t)=d_{max}$ for all times $t$) are illustrated in Figures 1 and 1 respectively. In the latter case, GLY cells become extinct and the patient recovers quickly, however it is natural to ask whether this treatment strategy is optimal in some sense (i.e. could the recovery be much delayed if the patient received therapy less often or at a lower intensity?). In the following section we will show that the MTD-based treatment can lead to an avoidably high cumulative amount of therapy (see Figure 2(c) or might even fail to achieve a recovery in situations where AT-based treatment would otherwise have succeeded (see Figure 6), much like the early results from Zhang et al. in metastatic prostate cancer 14.

One natural objective function to minimize is the total amount of therapy administered over the course of treatment (which in this case could be a surrogate for both toxicity and cost). This can be quantified as $D=\int\limits_{0}^{T}d(t)dt,$ where the total time of treatment $T$ is dependent on the initial cancer subpopulation fractions and on our chosen therapy policy $d(\cdot)$. However, this objective is problematic for two reasons. First, the minimum of $D$ is clearly attained without any therapy (taking $d(t)\equiv 0$ implies $D=0$ and $T=+\infty$), even though the dynamics become cyclic and the recovery is never achieved. Second, if we constrain our minimization to only those $d(\cdot)$ that lead to recovery, an optimal treatment policy does not exist. Instead, there is a sequence of treatment policies that lead to successively smaller $d$ values but with an unbounded increase in corresponding treatment times $T$. The idea of such policies is simple: travel along the therapy-free trajectories of Figure 1 for most of the time, but use short bursts of therapy only when the drugs are most effective. To approach the optimally small $d$, one would need to use shorter and shorter bursts, resulting in policies that are hard to implement in practice and would require unrealistically long treatment times $T$, but would yield a situation like a chronic disease, where while the tumor is never cured, it is always controlled.

In order to get a meaningful optimal policy we will penalize the treatment time by a time penalty $\sigma>0$. The total time spent on the treatment, including time when a patient skips the doses, is an important factor by itself. Much longer time spent on a treatment causes worse quality of life and additional costs for a patient. Therefore, our objective is to minimize the sum of a therapy cost $D$ and a treatment time cost $\sigma T$, while guaranteeing the treatment results in recovery. For every choice of $\sigma>0,$ the resulting treatment policies are thus Pareto-optimal with respect to $D$ and $T$.

In this paper we consider an objective function that is equal to $\int\limits_{0}^{T}d(t)\,dt\,+\,\sigma T$, if a patient recovers, and is equal to $+\infty$ otherwise. The value function, $u$, is defined as the minimum of this objective function (over the set of treatment policies), and any policy $d(\cdot)$ that realizes this minimum is called optimal.

Due to the structure of this optimization problem, one can show that virtually all optimal treatment policies are bang-bang: at any given time $t$, they either administer therapy at MTD-rate ($d(t)=d_{max}$) or administer no drugs at all ($d(t)=0$); see Section 6.1 in Supplementary Materials. For such policies, the objective function becomes a weighted sum of the total therapy time $\widetilde{T}$ (when $d(t)\equiv d_{max}$) and the total treatment time $T$, with $d_{max}$ and $\sigma$ as the corresponding weights. Moreover, this allows for a simple visual representation of any such policy: splitting the full state space into two parts (MTD dynamics vs. no-therapy dynamics), shown in yellow and blue respectively in all figures throughout this paper, and simplifies application of therapy into something familiar to clinicians: on or off.

3 Results

3.1 Quantifying the benefits of optimal treatment strategies

In this section we apply the optimal control theory to an example considered in Kaznatcheev et al.40. The details of the optimal control problem formulation are summarized in Box 3. We take the same model parameters as in Figure 2 of Kaznatcheev et al.40: $b_{a}=2.5$, $b_{v}=2$, $c=1$, $n=4,$ strength of MTD $d_{max}=3$, and initial state $(q_{0},p_{0})=(0.6,0.9)$, which by formula (2) is equivalent to $(x_{D},x_{G},x_{V})=(0.04,0.9,0.06)$. We also use $\sigma=0.01$ to incorporate a time-penalty absent in the original model. We take $r_{b}=f_{b}=10^{-1.5}$ as recovery and failure barriers‡‡‡This change in parameter values is meant to decrease the computational cost of our numerical approach (see section 6.2 in Supplementary Materials.) The original $r_{b}=f_{b}=10^{-4}$ from Kaznatcheev et al.40 would require computations on a finer mesh.

In Figure 2 we compare the treatment cost (10) and treatment time (8) of trajectories corresponding to four different treatment strategies.

The first two strategies we consider are similar to those modeled in Kaznatcheev et al.40: 2(a) is an example of a bad policy that may cause a failure by stopping the therapy prematurely, while 2(b) is a good policy based on ad-hoc adjustment of the start time for the therapy. We also illustrate a “MTD-based”’ policy 2(c), which is analogous to the standard of care using MTD as long as possible. Even though both 2(b) and 2(c) lead to recovery, neither of these is optimal (with the MTD-based approach resulting in excessive amount of drugs, captured by the higher cost). The policy minimizing our objective function can be found by solving the HJB equation; we illustrate it in 2(d).

The corresponding therapy on/off regions and the resulting vector field are shown in Figure 3. The zoomed version shows that trajectories can be prevented from crossing the failure barrier by using the MTDs just before crossing. In fact, a chattering control (with intermittent and sufficiently frequent use of MTDs) would be sufficient to guarantee this.

The level sets of $u$ in Figure 3 show that value functions need not be smooth. Since the gradient of $u$ is used to determine the optimal course of action (therapy on or off), there can actually be more than one optimal policy for initial states on a shockline (where that gradient is undefined). We show an example of such trajectories for an initial point $(x_{D},x_{G},x_{V})\approx(0.417,0.311,0.272)$ (solid green and red lines in Figure 3). Both trajectories yield the same cost of $2.764$. Moreover, non-smoothness of the value function often poses a challenge for method’s based on Pontryagin Maximum Principle (PMP) 28 even if the initial state is not on a shockline. For example, perturbing the initial state to a nearby $(x_{D},x_{G},x_{V})=(0.35,0.3,0.35)$, denoted by a cross in Figure 3, one sees two locally optimal trajectories and PMP might yield either of these depending on the initial guess. The green one is however inferior to the globally optimal red trajectory, which is always recovered by solving the HJB equation.

3.2 Optimization trade-offs: total administered drugs versus time to recovery

The optimal therapy-on regions are clearly dependent on specific values of all model parameters. Here we explore their dependence on $\sigma$ and $d_{max}$. Recall that, for every policy leading to recovery, the overall cost of treatment is a sum of the “therapy cost” (i.e., the total amount of drugs administered, $D=\int\limits_{0}^{T}d(t)dt$) and the treatment-time cost $\sigma T$. Since the optimal control is bang-bang, this can be re-written as a weighted sum of the time-till-recovery $T$, and the total drug therapy time $\tilde{T}\leq T$. That is, for controls based on repeated therapy-off/MTD-level-therapy switches, we can re-write the overall cost as $J=d_{max}\tilde{T}+\sigma T,$ with the ratio between the weights $(\sigma/d_{max})$ representing the “relative importance” of $T$ and $\tilde{T}$ for the optimization. But the functional role of these weights is quite different: while $\sigma$ can be chosen to reflect our preferences, the MTD-rate $d_{max}$ is dictated by the medical reality, which will be patient and drug specific. By varying $d_{max}$ while keeping $(\sigma/d_{max})$ constant, we can study the role played by the MTD-level in determining optimal policies under a fixed relative preference between the objectives. In Figure 4 we conduct this experiment for two different initial states and the same set of game parameters ($b_{a}=2.5$, $b_{v}=2$, $c=1$, $n=4$), demonstrating that the value of $d_{max}$ strongly influences the optimal policies and the shape of the “therapy-on” yellow regions.

On the other hand, for any fixed/biological $d_{max}$ value, we can vary $\sigma$ to study how the trade-off between $\tilde{T}$ and $T$ affects the optimization. This experiment is conducted for the same two initial states in Figure 5. As we can see, smaller $\sigma$ entails larger total time. Intuitively, this happens since it becomes “cheaper” to pause the therapy until we reach a “better” state to administer the drugs. Larger $\sigma$ leads to a shorter time-to-recovery $T$, but also an increase in the total amount of administered drugs $D=d_{max}\tilde{T}$ and a larger therapy-on region (shown in yellow) in the state space.

3.3 “Incurable” states and periodic trajectories under MTD treatment.

One might think that, despite being sub-optimal, an aggressive MTD-based strategy is at least always fully reliable and the resulting trajectories are guaranteed to reach the recovery zone from every initial configuration, as shown in Figure 2(c). Indeed, if $\frac{b_{a}}{n+1}\leq d_{max}$, the MDT-based policy ($d(t)\equiv d_{max}$) guarantees that $\dot{p}$ is always negative; see equation (7). But with $\frac{b_{a}}{n+1}>d_{max}$ the recovery might not be attained with the constant use of MTDs (even if some other treatment policies are successful!).

Consider, for example, the following set of parameters: $b_{a}=4$, $b_{v}=2$, $c=1$, $n=4$; $r_{b}=f_{b}=10^{-1.5}$; $d_{max}=0.3$, $\sigma=0.03$ and an initial state $(x_{D},x_{G},x_{V})=(0.02,0.8,0.18).$ Under these parameters the MTD-based therapy has a periodic trajectory§§§This is easy to prove by redefining the parameter $b_{a}:=b_{a}-(n+1)d_{max}>0$ and reducing the MTD-based case to the periodic behavior of the original uncontrolled system (3) in the parameter regime (4).; see Figure 6(b). Since the treatment time is infinite, the cost (10) of such a policy is $+\infty$. (In reality this would lead to the emergence of drug resistance and eventual failure, but this biological situation is not modeled in Kaznatcheev et al.40.)

We can see that neither of two extreme strategies (“no-drugs-at-all” in Figure 6(a) and the MTD-based “drugs-all-the-time” in Figure 6(b)) can bring the trajectory to the recovery zone. However, their adaptive combination can still achieve the objective. We show a trajectory corresponding to the optimal policy in Figure 6(c). With a larger failure zone (e.g., $r_{b}=f_{b}=10^{-1}$), a previously successful MTD-based treatment might even result in death (Figure 6(d)), while the adaptive strategy still leads to recovery (Figure 6(e)).

For a fixed treatment policy, we define its corresponding “incurable” area to be a set of states starting from which it is impossible to cross the recovery barrier. For example in Figure 6, the incurable area of the MTD-based policy includes the state $(x_{D},x_{G},x_{V})=(0.02,0.8,0.18)$ (when $r_{b}=f_{b}=10^{-1.5}$ or $r_{b}=f_{b}=10^{-1}$). However, this state is not in the (dramatically smaller) incurable area of the adaptive/optimal policy; see Figure 7.

Starting from any incurable configuration, one could similarly pose a different control problem of maximizing the time until crossing the failure barrier. While we do not address it here, we note that the HJB approach would be quite suitable to find optimal treatment policies for this problem as well.

Of course, the “incurable areas” are highly dependent on model parameters. In section 6.4 of Supplementary Materials, we show that they can grow due to an increase in the MTD rate $d_{max}$ or a decrease in vascularization benefits $b_{v}$.

4 Discussion

By now it is widely accepted that cancer is an evolutionary process, and that variation and selection drive the emergence of drug resistance. While this new knowledge is driving cancer research forward, it has largely not yet affected clinical practice, with the majority of clinical protocols relying on MTD-based approaches, which invariably fail in the setting of most metastatic disease. This is changing with the advent of adaptive therapy – therapeutic strategies specifically designed with a changing regimen prescribed: one that adapts to an evolving tumor 42.

To optimally design AT protocols, the underling dynamics of the tumor growth and treatment response must be known. While the methods of learning these dynamics are still in their infancy, the last decade has seen a flurry of activity using evolutionary game theory and other evolutionary models for a range of prescribed dynamics. These works have shown qualitative changes in tumor behavior in a range of treatment scenarios, and importantly, demonstrated that the order 43, sequence 40 and timing 17 of therapy can drastically change the outcomes. As we come closer to the reality of evolutionarily designed therapeutic trials in the mainstream, it is important then that we develop methods not just to make our outcomes better, but also to formally optimize them.

However, before using the mathematical tools for optimization, one also needs to choose a specific quantifiable criterion for comparing the outcomes. Once that criterion is selected and the underlying mathematical model is sufficiently accurate, the best treatment strategy can be found by the techniques of optimal control theory. In this paper, we show how this can be done for one particular heterogeneous cancer model previously described in Kaznatcheev et al.40. Given a set of model parameters ($b_{v},\leavevmode\nobreak\ b_{a},\leavevmode\nobreak\ c,\leavevmode\nobreak\ n$) and treatment/recovery parameters ($d_{max},r_{b},f_{b}$), we can find an optimal treatment policy for any initial distribution of cells ($q_{0},\leavevmode\nobreak\ p_{0}$). We show that the optimal treatment policy can have multiple regimes: always on $d^{*}(\cdot)\equiv d_{max}$, always off $d^{*}(\cdot)\equiv 0$, and involving several contiguous treatment periods. For the latter, the challenge is to accurately approximate the on/off “switching curves” in the state space. We show that the definition of optimal treatment policy is heavily dependent on a parameter $\sigma$, a “cost” term, describing the relative importance of minimizing the total amount of drugs administered versus the total time to recovery. We further show that, for some parameter regimes, there are “incurable regions” in the state space – the starting configurations that will not lead to a recovery regardless of the chosen treatment strategy – suggesting an alternative therapy (or goal) should be considered. Moreover, for some other starting configurations, the “always on” strategy might not lead to a recovery even if the recovery is actually achievable with some on/off hybrid strategies.

Just as any other model, the approach in Kaznatcheev et al.40 is based on simplifying assumptions (e.g., only the subpopulation fractions are important, and no novel types can arise), which limit its practical applicability. But our message is broader, and we use this specific model primarily to illustrate the general optimization approach. With adaptive therapy trials now under way showing early promise 14, we propose that such methods will become increasingly integrated into trial design discussions, with increasingly detailed cancer evolution equations or even in data-driven/equation-free framework.

Of the two main approaches of optimal control theory, the Pontryagin Maximum Principle (PMP) has been much more widely used in cancer treatment research up till now. In contrast, our approach here is based on dynamic programming and the numerical methods for Hamilton-Jacobi-Bellman (HJB) equations. The higher computational cost of these methods is balanced by several important practical considerations. First, they yield a policy in feedback-form and are thus more robust to modeling/measurement errors. Second, they always return the globally optimal treatment strategies and avoid some of the pitfalls well-know for the PMP-based methods (e.g., see Figure 3). With the advent of efficient numerical methods, we posit that the HJB equations will be soon playing a larger role in treatment optimization.

There are several obvious directions for extending our approach. First, the ability to optimize outcomes for a range of criteria will open new avenues to quantify physician/patient discussions concerning trade-offs and personalized therapy that have previously been only qualitative. In the current paper, computing optimal policies for different values of $\sigma$ can be viewed as a small step in this direction. But there are many other possible optimization criteria of practical interest. Methods for approximating all Pareto-optimal policies are available44 but are usually more computationally challenging. For probabilistic cancer evolution models, one can also choose between optimizing different characteristics of the same random quantity. (E.g., minimize the average time-to-recovery versus maximizing the probability of recovery in the next year). Finally, one can also use the choice of criterion to promote robustness by systematically treating possible measurement/modeling errors as perturbations chosen by some adversarial player. Such “games-against-nature”, as described recently by Stankova and colleagues 42, can be similarly treated by solving Hamilton-Jacobi-Isaacs equations.

Another important limitation of the above techniques is our assumed full knowledge of the system state. E.g., we assumed that the exact subpopulation fractions are known at every point in time and can be used to decide whether to administer the drugs. In practice, one can periodically obtain an approximation of these quantities (e.g., based on a repeat biopsy), but most of the time the decisions must be made based on some less invasive measurements (e.g., based on PSA-levels in castrate-resistant prostate cancer model14 or in cell free circulating DNA as recently shown by Khan and colleagues 45). A rigorous treatment of such optimization challenges will require the framework of partially-observable controlled processes; see e.g., Davis and Varaiya46.

With cancer evolution playing a larger and larger role in our thinking about therapy, and adaptive/evolutionary therapy coming to the fore, methods for optimizing these approaches will grow in importance. We suggest that the HJB method be used in these scenarios, and that clinicians and modelers begin discussions of further method development along these lines.

Acknowledgement

MG acknowledges the support from The Institute for Data and Decision Analytics at The Chinese University of Hong Kong, Shenzhen during his visit when this research was completed.

JGS would like to thank the NIH Case Comprehensive Cancer Center support grant P30CA043703 and the Calabresi Clinical Oncology Research Program, National Cancer Institute of Award Number K12CA076917.

AV would like to thank the Simons Foundation for its fellowship support and the National Science Foundation (award DMS-1738010) for supporting development of numerical methods for Hamilton Jacobi equations. A part of this work was performed during a sabbatical leave at Princeton/ORFE, and AV is grateful to ORFE Department for its hospitality.

5 Conflict of Interest

The authors declare no conflict of interest.

6 Supplementary Materials

6.1 Deriving a Hamilton-Jacobi-Bellman equation

We start by explaining the logic of dynamic programming that yields the HJB PDE (12) (Box 3, Section 3.1), and the “bang-bang” property of optimal treatment policies.

Recall that the evolving composition of cancer sub-populations can be fully defined by $(q(t),p(t))$; see formulas (2) and (7). The process is tracked until we cross either a recovery or failure barrier; i.e., until the trajectory leaves \Omega=\Big{(}[0,1]\times[0,1]\Big{)}\backslash\Delta, with the terminal set $\Delta$ defined in formula (6). For an arbitrary initial state $(q_{0},p_{0})\in\Omega,$ the goal is to choose our treatment policy to minimize the integral of an instantaneous cost $K(d(t))=d(t)+\sigma$ up to the terminal time $T=T(q_{0},p_{0},d(\cdot)).$ I.e., the total cost of starting at $(q_{0},p_{0})$ and using a policy $d(\cdot)$ is

[TABLE]

where $g$ is the terminal cost specified on $\Delta$ in formula (9). The value function $u(q_{0},p_{0})$ is the result of minimizing $J$ over all available treatment policies, and we say that the policy $d^{*}(\cdot)$ is optimal if $u(q_{0},p_{0})=J(q_{0},p_{0},d^{*}(\cdot))$.

Bellman’s Optimality Principle 47 is the key idea of dynamic programming. It states that, if we move along any optimal trajectory, a remaining (yet to be traversed) part of that trajectory is in itself optimal from our current configuration/state. In terms of the above model,

[TABLE]

should hold for every sufficiently small $\tau>0.$ Assuming that the value function $u(q,p)$ and $d^{*}(t)$ are smooth, one can use Taylor series and take the limit $\tau\rightarrow 0$ to obtain

[TABLE]

Here $d^{*}_{0}=d^{*}(0)$ is the optimal initial rate of therapy starting from $(q_{0},p_{0})$ and $(\dot{q},\dot{p})$ are specified by the right hand side of the ODEs in (7). Since (15) does not involve $d^{*}(t)$ for any $t>0$, it is now natural to switch to a feedback control perspective based on a state-dependent (rather than explicitly time-dependent) optimal control $d^{*}_{0}(q,p)$. Since the latter is a priori unknown, a Hamilton-Jacobi PDE (12) is obtained by minimizing over all available control values $d\in[0,d_{max}]$ and demanding that (15) should hold at every $(q,p)\in\Omega.$ Additional boundary conditions $u=g$ are specified on $\Delta$ by (13).

The above derivation is merely formal since the value function $u$ is typically non-smooth. Indeed, (12) rarely has classical solutions, and if one considers Lipschitz-continuous weak solutions (by demanding that the PDE should hold wherever $\nabla u$ is defined), one immediately loses the uniqueness. Additional test conditions introduced by Crandall and Lions48 are employed to pick out a viscosity solution – the unique weak solution coinciding with the value function of the original control problem 36. Convergence to this viscosity solution is also a requirement for all numerical methods for HJB equations used in control-theoretic applications.

Using the dynamics (7) specific to our model, the HJB equation (12) can be re-written as follows:

[TABLE]

The linear $d$-dependence of the minimized expression allows us to find the minimizer in closed form:

[TABLE]

Therefore, an optimal treatment policy takes only extreme values – either [math] or $d_{max}$. This is usually called the bang-bang property.

Using (17) in practice would require knowing $u_{p}$ at every point $(q,p)$. In principle $u_{p}$, can be computed along an optimal trajectory (backwards, from the recovery barrier to our initial state $(q_{0},p_{0})$ without solving the PDE on the entire $\Omega$. This is in a sense the main idea of Pontryagin Maximum Principle (PMP) 28. This method will work (and will be much more computationally attractive) as long as $u$ remains smooth along the optimal trajectory. Unfortunately, the PMP has no way of identifying whether a backward-traced trajectory passes through a shockline (where $\nabla u$ is undefined). In practice, this would result in obtaining locally (rather than globally) optimal treatment policies; see the example in Figure 3c. We thus focus on solving the full HJB equation, yielding a variational formula for $d^{*}(q,p).$

6.2 Numerical methods for the Hamilton-Jacobi-Bellman equation

We obtain an approximate solution to HJB equations on a regular triangulated mesh over the $(x_{D},x_{G},x_{V})$ space; see Figure 8. Since at any moment of time $x_{V}(t)\equiv 1-x_{G}(t)-x_{D}(t)$, it is enough to consider an ODE system for two sub-populations $x_{D}(t)$ and $x_{G}(t)$. To simplify the notation, we will write

[TABLE]

where $y(t)=(x_{D}(t),x_{G}(t))$, $x=(x_{D}(0),x_{G}(0))$ and $f(\cdot)$ denotes the right-hand side of (7) in $(x_{D},x_{G})$ coordinates using the transformation (2).

Our approximation scheme is based on a first-order accurate semi-Lagrangian discretization 49. Starting at a meshpoint $x$ and using control $d$, we assume that the rate of change is constant for a small time $\tau$, yielding a new state

[TABLE]

Assuming that the running cost is also constant over that small time interval, one can rewrite Bellman’s optimality principle as

[TABLE]

Since $\tilde{x}_{d}$ is usually not a meshpoint, $u\left(\tilde{x}_{d}\right)$ is approximated by interpolation using the neighboring meshpoint values. (This is the key idea of all semi-Lagrangian techniques.) While there are many ways to choose $\tau,$ we select it for each $d$ value individually to guarantee that $\tilde{x}_{d}$ lies on a mesh line and only two neighboring values are needed for the interpolation; a similar approach has been used in 50, 51. More specifically, suppose that a vector $f(x,d)$ anchored at $x$ lies within a triangle $xx_{1d}x_{2d}$; see Figure 8(a). Then $\tilde{x}_{d}$ lies on a segment $x_{1d}x_{2d}$ with

[TABLE]

Recalling that only extreme rates $d$ can be optimal due to the bang-bang property, we obtain a coupled system of discretized equations:

[TABLE]

which must hold for each meshpoint $x\in\Omega$. The boundary conditions are handled by setting $U=0$ when $x_{G}<r_{b}$ and $U=+\infty$ when $x_{G}>1-f_{b}$.

In our implementation, the above coupled system of discretized equations is handled by Gauss-Seidel iterations, with an additional speed up through alternating meshpoint orderings (in a “Fast Sweeping” fashion)52, 53, 54. Another alternative would be to decouple the system dynamically – by selecting larger $\tau_{d}$ adaptively so that “already known” mesh values would be sufficient for updating the still-tentatively-known $U$ values. The latter “Ordered Upwind” approach has been primarily used in problems with geometric dynamics 51, 55, 56 and offers advantages when optimal trajectories frequently change directions. In the future, it would be interesting to extend it (as well as its two-scale hybrids with sweeping 57) to therapy optimization problems, particularly for the case of small $\sigma$.

From implementation purposes, it is easier to conduct computations in a Cartesian coordinate system (Figure 8(a)), which is equivalent to a regular triangular mesh (Figure 8) by a linear transformation.

To ensure the accuracy of the value function in Figure 3, we have used $n=9000$ meshpoints along one side of the GLY-VOP-DEF triangle, yielding $N=37\leavevmode\nobreak\ 988\leavevmode\nobreak\ 686$ meshpoints in $\Omega$ with $r_{b}=f_{b}=10^{-1.5}.$ The algorithm terminates when the difference between value functions in sequential iterations falls below $10^{-5}$, which required a total of 62 iterations (sweeps) for this example.

In the future, we hope to reduce the computational cost by using a higher-order accurate semi-Lagrangian discretization 49 on a coarser mesh, employing “Ordered Upwind” techniques to reduce or eliminate the coupling in the discretized system.

6.3 Fully angiogenic and glycolyctic tumors

In section 3 we have focused on optimizing treatment policies for polyclonal tumors. Under “therapy-off” policy any trajectory of a polyclonal tumour has periodic dynamics and model parameters satisfy (4). Two other types of tumours are also possible under the model40: fully angiogenic and fully glycolyctic.

A tumor has a fully angiogenic regime if \max\Big{(}\frac{b_{a}}{n+1},cn\Big{)}<b_{v}-c. If the model (3) satisfies this condition all cells tend to switch to VOP type and the trajectory converges to the recovery zone40. In some sense, the fully angiogenic regime is less interesting case for our analysis because even without any therapy a patient will recover. However, the optimal control analysis still might be useful when, for example, time-penalty is high (a patient wants to recover as soon as possible) and some amount of drugs can be applied to accelerate the recovery, see Figure 9.

A tumor has a fully glycolyctic regime if $\frac{b_{a}}{n+1}>b_{v}-c$, all cells tend to be GLY type cells and trajectories converge to the failure zone from any initial state. Even if a treatment policy gives some short-term results, the trajectory will turn towards the failure zone once the therapy is stopped. Nevertheless, crossing the recovery barrier means full recovery under assumptions of the model40. We consider an example of optimal policy for fully glycolyctic tumour in Figure 9.

6.4 Incurable areas changing under parameter variation

In Figure 10 we examine the changes in the optimal/minimal incurable area due to variations in the MTD rate and model parameters.