The value of monitoring to control evolving populations

Optimal Control of a Wright–Fisher Population with Perfect Monitoring.

The optimal control function u¯(xt,t) maximizes the probability that a polymorphism is still present after a time T. In the infinite horizon time limit T→∞, where the optimal control becomes stationary u¯(x), we expect it to also maximize the mean first passage time 〈Tx(0,1)〉 out of the interval (0, 1) (for any fixed starting position x). In Eq. 3, the control variable u appears only linearly. This means that only the two extreme control strengths are ever used to steer the system. This particular type of control is called bang-bang. It follows that the control profile u¯(x) will have the form of a step function with critical frequency xc,σ+u¯={σ,x<xcσ−uc,x≥xc.[4]The only remaining parameter is the critical threshold xc(σ,uc). To find an expression for the objective function 〈T(0,1)〉(xc,σ,uc), we can consider the optimally controlled system as the simplest example of evolution under frequency-dependent, piecewise constant selection σ+u¯(x). The mean first passage time can then be found analytically using standard methods for stochastic processes (24) (SI Text and Fig. S1). At the correct threshold and with strong selective forces (σ≫1 and σ−uc≪−1), the gains are substantial and the polymorphism can be maintained for very long times.

Loss of Control due to Imperfect Monitoring.

The main assumption made so far was that perfect information is available about the state of the system in the form of continuous (in time), synchronous (without delay), and exact (without error) measurements of x. These requirements are impossible to achieve in practice, when monitoring is always imperfect. As we will see, when the assumption of perfect information is relaxed, not only is control over the system lost, but the control profile u¯(x) also ceases to be optimal. Rather than turning to the theory of partially observable Markov decision problems (25), we will use numerical analysis to demonstrate the effect of monitoring with finite resolution in time (relaxing the first condition).

Consider the situation where measurements of the frequency x are given only at discrete times {τi}, whereas no information is available during the intervals of length Δ=τi+1−τi. The immediate question is, given a measurement xi, what control should one apply while waiting for the next measurement? The perfect-information control u¯(xi) is correct only initially, and thus only in the limit Δ→0. However, it is intuitively clear that a naive protocol, applying u¯(xi) during the entire interval Δ, cannot be optimal, because it does not anticipate the dynamics of x under this regime (see the decrease in survival time in Fig. 1B). For example, for 0<xi<xc the initial control is u=0, σ+u>0 and the frequency will, on average, increase and eventually cross the threshold xc. If one could observe the population at that point, the control should be switched to u=−uc until x crosses xc from above. The total result of the naive strategy is to amplify fluctuations due to this overshooting.

Playing-To-Win vs. Playing-Not-To-Lose.

Without a continuous flow of observations as input, a preemptive control protocol u∗(τ,xi) during the interval Δ must be precomputed and then faithfully carried out. In the discrete-time (Wright–Fisher evolution) setting, there are NΔ generation updates until the next measurement and therefore 2NΔ different bang-bang protocols to choose from. However, the example above suggests to search for a preemptive control in a much smaller space, namely within those protocols that start with either u∗=0 or u∗=−uc and then switch, at some later time τc(xi,Δ), to a neutral regime with u∗=−σ (if available). The size of this space is only 2 NΔ (we note that for finding an optimal control one would need to search protocols in the full space). The effect of such a control scheme is to move the population to a safe place and then try to keep it there. There are two important observations: first, this informed control outperforms the naive protocol significantly, especially for intermediate values of Δ (Fig. 1 B and C); second, the safe parking position moves away from the boundary toward x=0.5 for larger values of Δ (Fig. S3). This shift from an aggressive control strategy under perfect information (Δ=0, xc close to a boundary) to a more and more conservative one (aiming for x=0.5 and trying to stay there) can be summarized as playing-to-win vs. playing-not-to-lose.

A similar loss of control can be expected for other types of monitoring imperfections and is a general feature of stochastic optimal control. It is important to note that the perfect-information control problem, and its solution u¯, is a typical starting point for the analysis. The naive control protocol above is indeed optimal for Δ→0, and still a very good option for Δ≪1 (in units of N generations). In most cases, as we will see in the adaptive cancer therapy model below, finding u¯ is challenging in itself and can be a good guidance for finding well-performing control protocols even under imperfect conditions.

A Minimal Model of Drug Resistance in Cancer.

The desired features of a minimal model of drug resistance in cancer include: (i) a variable tumor cell population size N, (ii) at least two cell types, drug-sensitive and drug-resistant, (iii) a carrying capacity K that describes a (temporary) state of tumor homeostasis, and (iv) the possibility for mutation and selection between cell types. Control over the tumor can be applied via a drug that changes the evolutionary dynamics by increasing, for example, the death rate of sensitive cells. We will assume here, as others have done in the context of cancer (11), a well-mixed cell population where the birth (or rather duplication) rate of cells is regulated by the carrying capacity. The dynamics of the model we have chosen here is encapsulated in the following birth and death rates for sensitive and resistant cells:Bi(ns,nr)=(1+gi)ni1+gNK+ansK+μ0(ni¯−ni),Di(ns,nr)=ni(1+Fi(u)), i∈{s,r}, s¯=r, r¯=s,[5]with free growth rates gs=g+a, gr=g, growth advantage a≥0, and mutation rate μ0 between cell types. Fi encodes the effect of the drug (u=1) or its absence (u=0) on cell type i. For N≪K, the absolute growth rates are gi−Fi. A drug effect of the formFs(u)=ufs and Fr(u)=(1−u)fr,[6]renders the drug effective if fs>g+a. The value fr=0 corresponds to drug resistance as such, but fr>g implies that resistant cells thrive under the drug and are drug-addicted. Such an effect has been observed in mice with melanoma carrying oncogene BRAF mutation treated with vemurafenib (29). Altogether, sensitive and resistant cells initially grow exponentially until the total population size N=ns+nr≈K. For some of the protocols that we derive, it is important that the tumor is close to or can reach carrying capacity during treatment (e.g. patients who are inoperable or whose primary tumors are not resectable). At that stage, competition for resources, space, etc. becomes fierce. If sensitive cells have a differential growth advantage a>0 (they might not have to maintain an expensive resistance mechanism), resistant cells will eventually be removed from the tumor or reduced to a small fraction (of size μ0/a). In reality, this scenario might not materialize, as the next mutation could propel the tumor into a new phase of exponential growth.

For the stochastic version of this process we can assume independent and individual birth and death events with the above probabilities per unit time. In analogy to the Wright–Fisher binomial update rule, here we can use a Poisson-like update.ni→n′i=ni+Δni, with Δni=Δni+−Δni−[7]Δni+∼Pois(Bi), Δni−∼Pois(Di), 〈Δni〉=Bi−Di.The total increment Δni follows a Skellam distribution. When we let K→∞ while fixing the combinations γ≡Kg, α≡Ka, μ≡Kμ0, and ϕi=Kfi and setting t=τ/K, the system is described by a Fokker–Planck evolution equation (30) for the distribution P(xs,xr,t) with xi≡ni/K (SI Text). This scaling exercise is mainly important because it allows to relate systems with small K (hundreds to thousands, as required for numerical analysis) to systems with large K (≳108, as present in real cancers).

Optimal Cancer Therapy with Perfect Monitoring.

With the minimal model of drug resistance in cancer introduced above, we can start the program of stochastic optimal control to compute adaptive therapy protocols. The first task is to define the goal of such a program: what is the quantity one aims to optimize? One very important objective is to maximize the (expected) time until the cancer proceeds to the next, possibly lethal stage. This could mean the emergence of a new cell type with a much higher carrying capacity, e.g. with metastatic potential. We will denote this critical event simply with a “driver” event or “metastasis.” The rate of metastasis emergence is a combination of tumor size and the rate ν0 (per cell and generation) for the necessary features to appear via mutation.

Earlier, the optimal control for the Wright–Fisher evolution example turned out to be a piecewise constant function of allele frequency. Here, we need to find a control profile u¯(ns,nr). With perfect information, we would know ns and nr at all times and would base the control decision adaptively on these measurements. In analogy to Eq. 1, the control objective can be expressed asu¯0:T(ns0,nr0)=argmaxu0:Texp(−ν0∑t=0T(nst+nrt)),[8]where the right-hand side is the probability that metastasis has not yet happened after T generations (for ν0K≪1). Here, the control objective is a nonlinear function of the entire trajectory (ns,nr)0:T. As such, it is the simplest manifestation of a so-called risk-sensitive control problem (see e.g. ref. 31). The above formulation assumes a finite (receding) horizon time T and also that control itself is cost-free. In cancer therapy, especially chemotherapy, this is certainly not the case: the side effects of treatment incur a considerable cost in terms of life quality and medical care. The difficulty, however, lies in quantifying these control costs in a manner that would make them comparable to the potential costs considered here. A simple extension of the above control objective to include such a control cost is demonstrated in SI Text and Fig. S4. The recurrence equation for the cost-to-go J(ns,nr,t) for the control objective above is given by (31)J(ns,nr,t)=e−ν0(ns+nr)maxu∈{0,1}〈J(t+1;u)〉,[9]〈J(t′;u)〉≡∑n′s,n′rJ(n′s,n′r,t′) W(n′s,n′r|ns,nr;u),[10]with boundary condition J(T)=1. The (microscopic) transition matrix W is the product of the two Skellam distributions resulting from Eq. 7 (including boundary conditions). With this equation, we can solve the dynamic programming task numerically for moderate values of K (see Fig. S5). For the numerical analysis, we have to introduce an upper bound N˜≫K+K for the population size. The resulting control profiles for a number of different parameter regimes is shown in Fig. 2.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Control of a tumor cell population. (A–D) The optimal control profile under perfect information about ns and nr for different parameters of the cancer model. In the white areas, u¯=0 (no drug), whereas in the gray areas, u¯=1 (with drug). The arrows indicate the deterministic flow. All profiles were calculated via Eq. 9 with T=K/ν generations and K=500 with an absorbing boundary at N=750. The sample trajectories were simulated with K=104 and controlled according to these profiles. The coloring of the trajectories shows the temporal evolution from blue to red. (A) When selection against resistance is stronger than driver emergence, α≫ν, the optimal protocol is to wait until resistant cells are cleared from the system before the drug is applied. (B) For higher driver emergence rates, the drug is applied earlier, which can lead to cycles. (C) For drug-sensitive (ϕs≫γ) and drug-addicted cells (ϕr≫γ) with high mutation (μ≫1), the control in the symmetric case (ϕs=ϕr) is a simple majority rule and very effective. (D) For smaller mutation (μ=1), the optimal strategy first homogenizes the tumor before trying to remove it.

In the case of ϕr=0, resistant cells are unaffected by the drug. If maintenance of the resistance mechanism is costly (α>0), the only way that they can be removed from the population is when selection can act against them. This only happens at N∼K and with u=0 (no drug). If this can take place before the next driver typically appears (if α=Ka≫KNν0∼K2ν0≡ν), then the optimal control protocol is to postpone treatment until the resistant cells are sufficiently cleared from the system so that genetic drift has a chance of removing the remaining part (Fig. 2A). However, this parameter regime of very high selection against resistance and/or very low rate of driver mutation, and therefore this therapy option, is not realistic for cancer. For higher values of ν/α, the optimal strategy is to apply the drug earlier (Fig. 2B). This procedure can lead to cycles of tumor size reduction followed by regrowth, with the overall effect of extending the time until metastasis.

If ϕr>γ (and ϕs>γ+α), resistant cells are actually drug-addicted and thrive only in its presence. Such a situation would be easy to control with perfect information about ns and nr. For example, if mutation between cell types is rapid (μ≫1), a majority rule is optimal (Fig. 2C). For a lower mutation rate, the optimal protocol first tries to amplify one cell type before switching to an environment that is deadly for most cells (Fig. 2D).

The effectiveness of different therapy protocols is compared in Fig. 3 with 1,000 stochastic forward simulations (with K=104) for the parameter setting of Fig. 2D. Whereas no therapy (u=0) and all-out therapy (u=1) both ultimately end with the occurrence of metastasis, adaptive therapy can bring the tumor size down to zero in the majority of cases. In metronome therapy, the drug is applied (withheld) for fixed time intervals τon (τoff). With numerically optimized values for the time intervals, metronome therapy is quite competitive compared with adaptive therapy.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Comparison of cancer therapies. For the parameter setting of Fig. 2D, different therapies are compared via 1,000 forward simulations with K=104. Shown is the fraction of runs that have not yet developed a metastasis mutation as a function of elapsed time. All-out maximum dosage therapy (u≡1) is only slightly better than no therapy (u≡0) in avoiding metastasis. Much better is metronome therapy with τon=τoff=0.1 K generations, with almost 70% success rate (values numerically optimized). Adaptive therapy removes the tumor in close to 90% of runs.

All these control strategies require perfect information, not only in the sense of the earlier Wright–Fisher example (continuous, synchronous, and exact), but also in terms of the inner tumor composition N=ns+nr, which presupposes that sensitive and resistant cells can be distinguished.

Loss of Therapy Efficacy due to Low-Resolution Monitoring.

There are very few cases where the genetic basis for a drug-resistance mechanism is known and can be specifically monitored (29, 32). In most cases the regrowth of the tumor under the drug is observed without understanding the exact biological processes responsible for the resistance. Here we aim to find rational control strategies when only the total tumor cell population size can be monitored. The adaptive therapy protocol that was applied by Gatenby et al. in ref. 16 (coupling the drug concentration to the tumor size) is one example of such a strategy.

Consider the situation where only the total population size N=ns+nr can be (perfectly) monitored, whereas the dynamical laws in Eqs. 5–7 and all parameter values are known. Under these circumstances, the perfect-information optimal control profiles from the last section cannot be used directly. However, there is still valuable information available. The response Nτ→Nτ+Δτ of the tumor size to a control choice over a time interval Δτ can give an indication of the inner tumor composition. As we have seen earlier, the length of time interval Δτ should be shorter than all other intrinsic time scales to enable control. The procedure we follow is to first derive an effective propagator Weff(Nτ+Δτ|Nτ,Nτ−Δτ,uτ−Δτ;uτ) (SI Text) and then repeat the cost-to-go calculation of Eqs. 9 and 10. This propagator takes into account not only the current size Nτ, but also the last measurement Nτ−Δτ and the last control decision uτ−Δτ. It follows from the microscopic W used in Eq. 10 by integrating over the internal degrees of freedom ns and nr at the three time points. Accordingly, the control profile is now a function of (Nτ−Δτ,uτ−Δτ,Nτ). For the parameter values leading to the majority rule in Fig. 2C, the new control profile is shown in Fig. S6. The drug regimen (u=0 or u=1) is maintained as long as the tumor size decreases sufficiently. At the first sign of possible reversal, the regimen is switched.