Mathematical Details on a Cancer Resistance Model

One of the most important factors limiting the success of chemotherapy in cancer treatment is the phenomenon of drug resistance. We have recently introduced a framework for quantifying the effects of induced and non-induced resistance to cancer chemotherapy (Greene et al., 2018a, 2019). In this work, we expound on the details relating to an optimal control problem outlined in Greene et al. (2018a). The control structure is precisely characterized as a concatenation of bang-bang and path-constrained arcs via the Pontryagin Maximum Principle and differential Lie algebraic techniques. A structural identifiability analysis is also presented, demonstrating that patient-specific parameters may be measured and thus utilized in the design of optimal therapies prior to the commencement of therapy. For completeness, a detailed analysis of existence results is also included.

1. Introduction

The ability of cancer chemotherapies to successfully eradicate cancer populations is limited by the presence of drug resistance. Cells may become resistant through a variety of cellular and micro-environmental mechanisms (Gottesman, 2002). These mechanisms are exceedingly complex and diverse, and remain to be completely understood. Equally complex is the manner in which cancer cells obtain the resistant phenotype. Classically resistance was understood to be conferred by random genetic mutations; more recently, the role of epigenetic phenotype switching was discovered as another mediator of resistance (Pisco et al., 2013). Importantly, both of these phenomena were seen as drug-independent, so that the generation of resistance is functionally separate from the selection mechanism (e.g., the drug). However, experimental studies from the past ten years suggest that drug resistance in cancer may actually be induced by the application of chemotherapy (Sharma et al., 2010; Pisco et al., 2013; Goldman et al., 2015; Doherty et al., 2016; Shaffer et al., 2017).

In view of the multitude of ways by which a cancer cell may become chemoresistant, we have previously introduced a mathematical framework to differentiate drug-independent from drug-dependent resistance (Greene et al., 2019). In that work, we demonstrated that induced resistance may play a crucial role in therapy outcome, and also discussed methods by which a treatment's induction potential may be identified via biological assays. An extension of the work was outlined in the conference paper (Greene et al., 2018a), where a formal optimal control problem was introduced and an initial mathematical analysis was performed. The aim of this work is to formalize the parameter identifiability properties of our theoretical model, to establish the existence of the optimal control introduced in Greene et al. (2018a), and to precisely classify the optimal control structure utilizing the Pontryagin Maximum Principle and differential-geometric techniques. A numerical investigation of both the control structure and corresponding objective is also undertaken as a function of patient-specific parameters, and clinical conclusions are emphasized.

The work is organized as follows. In section 2, we briefly review the mathematical model together with the underlying assumptions. Section 3 restates the optimal control problem, and the Maximum Principle is analyzed in section 4. A precise theoretical characterization of the optimal control structure is summarized in section 5. In section 6, we compare theoretical results with numerical computations, and investigate the variation in control structure and objective as a function of parameters. Conclusions are presented in section 8. We also include additional properties, including details on structural identifiability and existence of optimal controls, in Section 7.

2. Mathematical Modeling of Induced Drug Resistance

We briefly review the model presented in Greene et al. (2019) and analyzed further in Greene et al. (2018a). In that work, we have constructed a simple dynamical model which describes the evolution of drug resistance through both drug-independent (e.g., random point mutations, gene amplification, stochastic state switching) and drug-dependent (e.g., mutagenicity, epigenetic modifications) mechanisms. Drug-induced resistance, although experimentally observed, remains poorly understood. It is our hope that a mathematical analysis will provide mechanistic insight and produce a more complete understanding of this process by which cancer cells inhibit treatment efficacy.

A network diagram of the model under consideration is provided in Figure 1. Specifically, we assume that the tumor being studied is composed of two types of cells: sensitive (x₁) and resistant (x₂). For simplicity, the drug is taken as completely ineffective against the resistant population, while the log-kill hypothesis (Traina and Norton, 2011) is assumed for the sensitive cells. Complete resistance is of course unrealistic, but can serve as a reasonable approximation, especially when toxicity constraints may limit the total amount of drug that may be administered. Furthermore, this assumption permits a natural metric on treatment efficacy that may not exist otherwise (see section 3). The effect of treatment is considered as a control agent u(t), which we assume is a locally bounded Lebesgue measurable function taking values in $\text{R}$₊. Here u(t) is directly related to the applied drug dosage D(t), and in the present work we assume that we have explicit control over u(t). Later, during the formulation of the optimal control problem (section 3), we will make precise specifications on the control set U. Even though an arbitrary dosage schedule is unrealistic as a treatment strategy, our objective in this work is to understand the fundamental mathematical questions associated with drug-induced resistance, so we believe the simplification is justified. Furthermore, our results in section 5 suggest that the applied optimal treatment should take a relatively simple form, which may be approximated with sufficient accuracy in a clinical setting. Sensitive and resistant cells are assumed to compete for resources in the tumor microenvironment; this is modeled via a joint carrying capacity, which we have scaled to one. Furthermore, cells are allowed to transition between the two phenotypes in both a drug-independent and drug-dependent manner. All random transitions to the resistant phenotype are modeled utilizing a common term, ϵx₁, which accounts for both genetic mutations and epigenetic events occurring independently of the application of treatment. Drug-induced deaths are assumed of the form du(t)x₁, where d is the drug cytotoxicity parameter relating to the log-kill hypothesis. Drug-induced transitions are assumed to be of the form αu(t)x₁, which implies that the per-capita drug-induced transition rate is directly proportional to the dosage [as we assume full control on u(t), i.e. pharmacokinetics are ignored]. Of course, other functional relationships may exist, but since the problem is not well-studied, we consider it reasonable to begin our analysis in this simple framework. The above assumptions then yield the following system of ordinary differential equations (ODEs):

$\begin{array}{l}\frac{d{x}_1}{dt}=(1-(x_1+x_2))x_1-(ϵ+\alpha u(t))x_1-du(t)x_1\\ \frac{d{x}_2}{dt}=p_r(1-(x_1+x_2))x_2+(ϵ+\alpha u(t))x_1.& (1)\end{array}$

All parameters are taken as non-negative, and 0 ≤ p_r < 1. The restriction on p_r emerges due to (1) already being non-dimensionalized, as p_r represents the relative growth rate of the resistant population with respect to that of the sensitive cells. The condition p_r < 1 thus assumes that the resistant cells divide more slowly than their sensitive counterparts, which is observed experimentally (Shackney et al., 1978; Lee, 1993; Brimacombe et al., 2009). As mentioned previously, many simplifying assumptions are made in system (1). Specifically, both types of resistance (random genetic and epigenetic) are modeled as dynamically equivalent; both possess the same division rate p_r and spontaneous (i.e., drug-independent) transition rate ϵ. Thus, the resistant compartment x₂ denotes the total resistant subpopulation.

FIGURE 1

Figure 1. Visualization of interactions considered in system (1).

The region

$\begin{array}{ll}\Omega =\{(x_1,x_2)|0\le x_1+x_2\le 1,x_1,x_2\ge 0\}& (2)\end{array}$

in the first quadrant is forward invariant for any locally bounded Lebesgue measurable treatment function u(t) taking values in $\text{R}$₊. Furthermore, if ϵ > 0, the population of (1) becomes asymptotically resistant:

$\begin{array}{ll}\left(\begin{array}{l}{x}_1(t)\hfill \\ x_2(t)\hfill \end{array}\right)\stackrel{t\to \infty }{\to }\left(\begin{array}{l}0\hfill \\ 1\hfill \end{array}\right).& (3)\end{array}$

For a proof, see Theorem 2 in SI A in Greene et al. (2019). Thus in our model, the phenomenon of drug resistance is inevitable. However, we may still implement control strategies which, for example, may increase patient survival time. Such aspects will inform the objective introduced in section 3. For more details on the formulation and dynamics of system (1), we refer the reader to Greene et al. (2019).

3. Optimal Control Formulation

As discussed in section 2, all treatment strategies u(t) result in an entirely resistant tumor: $\bar{x}:=({\bar{x}}_1,{\bar{x}}_2)=(0,1)$ is globally asymptotically stable for all initial conditions in region Ω. Thus, any chemotherapeutic protocol will eventually fail, and a new drug must be introduced (not modeled in this work, but the subject of future study). Therefore, selecting an objective which minimizes tumor volume (x₁ + x₂) or resistant fraction [x₂/(x₁ + x₂)] at a fixed time horizon would be specious for our modeling framework. However, one can still combine therapeutic efficacy and clonal competition to influence transient dynamics and possibly prolong patient life, as has been shown clinically utilizing real-time patient data (Gatenby et al., 2009).

Toxicity as well as pharmacokinetic constraints limit the amount of drug to be applied at any given instant. Thus, we assume that there exists some number M > 0 such that u(t) ≤ M for all t ≥ 0. Any Lebesgue measurable treatment regime u(t) is considered, so that the control set is U = [0, M] and the set of admissible controls ${U}$ is

$\begin{array}{ll}{U}=\{u:[0,\infty )\to \left[0,M\right]|u\text{ is Lebesgue measurable}\}.& (4)\end{array}$

Recall that all cellular populations have been normalized to remain in [0, 1]. We assume that there is a critical tumor volume V_c, at which treatment, by definition, has failed. Our interpretation is that a tumor volume larger than V_c interferes with normal biological function, while x₁ + x₂ ≤ V_c indicates a clinically acceptable state. Different diseases will have different V_c values. For technical reasons needed in section 5 we assume that V_c < 1 − ϵ. This is a mild assumption, since genetic mutation rates ϵ are generally small (Loeb et al., 1974), and our interest is on the impact of induced resistance. Thus

$\begin{array}{ll}{V}_c\in (0,1-ϵ).& (5)\end{array}$

Define t_c as the time at which the tumor increases above size V_c for the first time. To be precise,

$\begin{array}{ll}{t}_c(u):=max\{T|x_1(t)+x_2(t)\le V_c\text{ for all }t\in \left[0,T\right]\}.& (6)\end{array}$

Since all treatments approach the state (0, 1), t_c(u) is well-defined for each treatment u(t). For simplicity, denote t_c = t_c(u) in the remainder of the work. The time t_c is then a measure of treatment efficacy, and our goal is then to find those controls u_* which maximize t_c. Writing in standard form as a minimization problem, we have the following objective:

$\begin{array}{ll}\underset{u\in {U}}{\min }\{J(u)\}=\underset{u\in {U}}{\min }\left\{-{\displaystyle {\int }_0^{t_c}1dt}\right\}.& (7)\end{array}$

We are thus seeking a control $u_{*}(t)\in {U}$ which maximizes t_c, i.e. solves the time-optimal minimization problem (7) restricted to the dynamic state equations given by the system described in (1) and the condition x₁(t) + x₂(t) ≤ V_c for all 0 ≤ t ≤ t_c. Note that the above is formulated (using the negative sign) as a minimization problem to be consistent with previous literature and results related to the Pontryagin Maximum Principle (PMP) (Ledzewicz and Schättler, 2012). Note that maximization is still utilized in section 7.2 and section 4.1, and we believe that the objective will be clear from context. To be consistent with notation utilized later, we denote the system (1) as

$\begin{array}{ll}\dot{x}=f(x)+u(t)g(x),& (8)\end{array}$

where

$\begin{array}{ll}f(x)=\left(\begin{array}{c}(1-(x_1+x_2))x_1-ϵx_1\\ p_r(1-(x_1+x_2))x_2+ϵx_1\end{array}\right),& (9)\end{array}$

$\begin{array}{ll}g(x)=\left(\begin{array}{c}-(\alpha +d)\\ \alpha \end{array}\right)x_1& (10)\end{array}$

and x(t) = (x₁(t), x₂(t)). By continuity of solutions, the time t_c must satisfy the terminal condition (t_c, x(t_c)) ∈ N, where N is the line x₁ + x₂ = V_c in Ω, i.e.,

$\begin{array}{ll}N={\psi }^{-1}(0)\cap \Omega ,& (11)\end{array}$

where

$\begin{array}{ll}\psi (x_1,x_2):=x_1+x_2-V_c.& (12)\end{array}$

With this notation, the path-constraint

$\begin{array}{ll}\psi (x_1(t),x_2(t))\le 0& (13)\end{array}$

must also hold for 0 ≤ t ≤ t_c. Equation (13) ensures that the tumor remains below critical volume V_c for the duration of treatment. Equivalently, the dynamics are restricted to lie in the set Ω_c ⊆ Ω, where

$\begin{array}{ll}{\Omega }_c:=\{(x_1,x_2)|0\le x_1+x_2\le V_c,x_1,x_2\ge 0\},& (14)\end{array}$

for all times t such that t ∈ [0, t_c]. The initial state

$\begin{array}{ll}{x}_0=(x_1^0,x_2^0)& (15)\end{array}$

is also assumed to lie in Ω_c. Except for section 7.1 where we restrict to the case $x_2^0=0$, the remainder of the work allows for arbitrary $x_2^0\in [0,V_c)$.

4. Maximum Principle

We dedicate the present section to characterize the optimal control utilizing the Pontryagin Maximum Principle (PMP). The subsequent analysis is strongly influenced by the Lie-derivative techniques introduced by Sussmann (1982, 1987a,b,c). For an excellent source on both the general theory and applications to cancer biology, see the textbooks by Ledzewicz and Schättler (2012) and Schättler and Ledzewicz (2015).

Before starting our analysis of the behavior and response of system (1) to applied treatment strategies u(t) utilizing geometric methods, we would like to mention that we have not found a reference for existence of optimal controls for a problem such as this, due perhaps to the non-standard character of it (maximization of time, path constraints). For this reason, we have added a self-contained proof of regarding existence in section 7.2.

4.1. Elimination of Path Constraints

We begin our analysis by separating interior controls from those determined by the path-constraint (13) (equivalently, x ∈ N). The following theorem implies that outside of the one-dimensional manifold N, the optimal pair (x_*, u_*) solves the same local optimization problem without the path and terminal constraints. More precisely, the necessary conditions of the PMP (see section 4.2) at states not on N are exactly the conditions of the corresponding maximization problem with no path or terminal constraints.

THEOREM 1. Suppose that x_* is an optimal trajectory. Let t₁ be the first time such that x_*(t) ∈ N. Fix δ > 0 such that t₁ − δ > 0, and

$\begin{array}{ll}\xi =x_{*}(t_1-\delta ).& (16)\end{array}$

Define z(t): = x_*(t)|_{t ∈ [0,t1−δ]}. Then the trajectory z is a local solution of the corresponding time maximization problem t_c with boundary conditions x(0) = x⁰, x(t_c) = ξ, and no additional path constraints. Hence at all times t, the path z (together with the corresponding control and adjoint) must satisfy the corresponding unconstrained Pontryagin Maximum Principle.

Proof. We first claim that z satisfies the path-constrained maximization problem with boundary conditions $x(0)=x^0,x(t_c)=\xi$. This is a standard dynamic programming argument: if there exists a trajectory $\bar{z}$ such that $\bar{z}(\tau )=\xi$, τ > t₁ − δ, concatenate $\bar{z}(t){|}_{t\in \left[0,\tau \right]}$ with x_*(t)|_{t ∈ [τ,tc]} at t = τ to obtain a feasible trajectory satisfying all constraints. This trajectory then has total time τ + δ + t_c − t₁ > t_c, contradicting the global optimality of x_*.

Recall that t₁ was the first time such that x_*(t) ∈ N. Since z is compact, we can find a neighborhood of z that lies entirely in {x|x ∉ N}. As the Maximum Principle is a local condition with respect to the state, this completes the proof.□

Theorem 1 then tells us that for states x = (x₁, x₂) such that x₁ + x₂ < V_c, the corresponding unconstrained PMP must be satisfied by any extremal lift of the original problem. (Recall that an extremal lift of an optimal trajectory is obtained by adding the Lagrange multipliers, or adjoint variables, to the control and state; see details in Definition 2.2.4, page 95, Chapter 2 of Ledzewicz and Schättler, 2012). We have now demonstrated that the optimal control consists of concatenations of controls obtained from the unconstrained necessary conditions and controls of the form (18). In the next section, we analyze the Maximum Principle in the region x₁ + x₂ < V_c. Furthermore, the constraint (13) has generic order one. In other words,

$\begin{array}{ll}{L}_g\psi =\nabla \psi ·g\ne 0.& (17)\end{array}$

Therefore, the feedback control (also known as the constrained control) can be found by differentiating the function (12) to insure that trajectories remain on the line N:

$\begin{array}{ll}{u}_p(x_1,x_2)=\frac{1}{d}\frac{(1-(x_1+x_2))(x_1+p_r{x}_2)}{x_1}.& (18)\end{array}$

Its existence however does not imply its feasibility, which is discussed below. Specifically, u_p can be shown to be a decreasing function of x₁ which is feasible on the portion of N satisfying $x_1^{*}\le x_1\le V_c$, where $x_1^{*}$ is given in (20), and infeasible elsewhere. This is proven in Proposition 3, and the geometric structure is depicted in Figure 2. Propositions 4 and 5 provide characterizations of the volume dynamics in certain regions of phase space, and are included here for completeness.

FIGURE 2

Figure 2. Region in Ω_c where L_YV(x) is guaranteed to be positive. That is, applying the maximal dosage M results in an increasing cancer population in the yellow-shaded region of phase-space.

Proposition 2. Suppose that the maximal dosage M satisfies

$\begin{array}{ll}M>\frac{(1-V_c)(1-p_r)}{d}.& (19)\end{array}$

and the point $x^{*}=(x_1^{*},x_2^{*})\in N$ with coordinates

$\begin{array}{l}{x}_1^{*}=\frac{p_r(1-V_c)V_c}{dM-(1-V_c)(1-p_r)},\\ x_2^{*}=V_c{\displaystyle (1-\frac{p_r(1-V_c)}{dM-(1-V_c)(1-p_r)})}.& (20)\end{array}$

Denote by Y(x) = f(x) + Mg(x) the vector field corresponding to the maximal allowed dosage M [here, f and g are the functions defined in (9), (10)]. The Lie derivative, for any x ∈ N, of the volume function V(x) = x₁ + x₂ with respect to Y is

(a) positive if $x_1<x_1^{*}$,

(b) zero at $(x_1^{*},x_2^{*})$, and

(c) negative if $x_1>x_1^{*}$.

Proof. We verify the above claims with a direct calculation. Let L_YV(x) denotes the Lie derivative of V(x) with respect to Y. Thus, for x ∈ N,

$\begin{array}{l}{L}_{Y}V(x)=\nabla V(x)·Y(x)\\ \\ =\left(\begin{array}{l}1\\ 1\end{array}\right)·\left(\begin{array}{c}\left[1-V_c-ϵ-(\alpha +d)M\right]x_1\\ \left[ϵ+\alpha M-p_r(1-V_c)\right]x_1+p_r(1-V_c)V_c\end{array}\right)\\ =\left[1-V_c-ϵ-(\alpha +d)M\right]x_1\\ +\left[ϵ+\alpha M-p_r(1-V_c)\right]x_1\\ +p_r(1-V_c)V_c\\ \\ =\left[(1-V_c)(1-p_r)-dM\right]x_1+p_r(1-V_c)V_c.\end{array}$

Assuming $M>\frac{(1-V_c)(1-p_r)}{d}$, the sign of L_YV(x) is as in the statement of the proposition.□

Proposition 2 implies that if the allowable dosage is large enough (Equation 9), treatment can at least decrease the tumor in certain regions of phase space. If this condition was not met, then the applied drug would generally be ineffective in reducing the tumor volume V, and hence not be utilized in a clinical scenario.

Proposition 3. Let x be a point on the line N. The feedback control u_p is unfeasible if $x_1\in (0,x_1^{*})$, and is feasible if $x_1\in (x_1^{*},V_c)$

Proof. For x ∈ N we compute

$\begin{array}{l}{u}_p(x)=\frac{(1-V_c)(1-p_r)}{d}+\frac{(1-V_c)p_r{V}_c}{d{x}_1}\ge 0.\end{array}$

It is straightforward to check that u_p > M if $x_1<x_1^{*}$. In addition, the feedback control, when restricted to points in N, is a decreasing function with respect to x₁. Thus, it is feasible for x ∈ N if $x_1\in (x_1^{*},V_c)$.□

Proposition 4. For x = (x₁, x₂) ∈ Ω_c with

$\begin{array}{ll}{x}_2>\frac{dM-(1-V_c)}{p_r(1-V_c)}{x}_1,& (21)\end{array}$

the Lie derivative L_YV(x) is positive.

Proof. As in Proposition 2, we compute L_YV(x) directly:

$\begin{array}{l}{L}_{Y}V(x)=(1-(x_1+x_2))(x_1+p_r{x}_2)-dM{x}_1\\ \ge (1-V_c)(x_1+p_r{x}_2)-dM{x}_1\\ =\left[(1-V_c)-dM\right]x_1+p_r(1-V_c)x_2\\ >\left[(1-V_c)-dM\right]x_1+p_r(1-V_c)\frac{dM-(1-V_c)}{p_r(1-V_c)}{x}_1\\ =0,\end{array}$

where the first inequality utilizes V ≤ V_c, and the second relies on (21)□

Proposition 5. For

$\begin{array}{l}M>\frac{1-ϵ}{\alpha +d},\end{array}$

trajectories corresponding to the maximal dosage M have a decreasing sensitive cellular population.

Proof. For u(t) ≡ M, the corresponding sensitive trajectory is given by

$\begin{array}{l}\dot{x_1}=(1-(x_1+x_2))x_1-ϵx_1-(\alpha +d)M{x}_1\\ <(1-(x_1+x_2))x_1-ϵx_1-(1-ϵ)x_1\\ =-(x_1+x_2)x_2\le 0\end{array}$

Note that we are assuming here that the maximal dosage M satisfies $M>\frac{1-ϵ}{\alpha +d}$.□

4.2. Maximum Principle and Necessary Conditions at Interior Points

Necessary conditions for the optimization problem discussed in section 3 without path or terminal constraints are derived from the Pontryagin Maximum Principle (Pontryagin, 1987; Ledzewicz and Schättler, 2012). The corresponding Hamiltonian function H is defined as

$\begin{array}{ll}H({\lambda }_0,\lambda ,x,u)=-{\lambda }_0+\langle \lambda ,f(x)\rangle +u\Phi (x,\lambda ),& (22)\end{array}$

where λ₀ ≥ 0 and λ ∈ $\text{R}$². Here 〈·, ·〉 denotes the standard inner product on $\text{R}$² and, since the dynamics are affine in the control u, Φ(x, λ) is the switching function:

$\begin{array}{ll}\Phi (x,\lambda )=\langle \lambda ,g(x)\rangle .& (23)\end{array}$

The Maximum Principle then yields the following theorem:

THEOREM 6. If the extremal (x_*, u_*) is optimal, there exists λ₀ ≥ 0 and a covector (adjoint) $\lambda :\left[0,t_c\right]\to {({\text{R}}^2)}^{*}$, such that the following hold:

1. (λ₀, λ(t)) ≠ 0 for all t ∈ [0, t_c].

2. λ(t) = (λ₁(t), λ₂(t)) satisfies the second-order differential equation

$\begin{array}{ll}\begin{array}{l}\dot{\lambda }(t)=\left(\begin{array}{cc}2x_1+x_2+ϵ-1& p_r{x}_2-ϵ\\ x_1& p_r(2x_2+x_1-1)\end{array}\right)\lambda (t)\\ +u(t)\left(\begin{array}{cc}\alpha +d& -\alpha \\ 0& 0\end{array}\right)\lambda (t)\end{array}& (24)\end{array}$

3. u_*(t) minimizes H pointwise over the control set U:

$\begin{array}{l}H({\lambda }_0,\lambda ,x_{*}(t),u_{*}(t))={\displaystyle \underset{v\in U}{min}}H({\lambda }_0,\lambda ,x_{*}(t),v).\end{array}$

Thus, the control u_*(t) must satisfy

$\begin{array}{ll}{u}_{*}(t)=\{\begin{array}{ll}0\hfill & \Phi (t)>0,\hfill \\ M\hfill & \Phi (t)<0.\hfill \end{array}& (25)\end{array}$

where

$\begin{array}{ll}\Phi (t):=\Phi (x_{*}(t),\lambda (t)).& (26)\end{array}$

4. The Hamiltonian H is identically zero along the extremal lift (x_*(t), u_*(t), λ(t)):

$\begin{array}{ll}H({\lambda }_0,\lambda (t),x_{*}(t),u_{*}(t))\equiv 0.& (27)\end{array}$

Proof. Most statements of Theorem 6 follow directly from the Maximum Principle, so proofs are omitted. In particular, items (1), (2), and the first part of (3) are immediate consequences (Ledzewicz and Schättler, 2012). Equation (25) follows directly since we minimize the function H, which is affine in u (see Equation 22). The Hamiltonian vanishes along (x_*(t), u_*(t), λ(t)) since it is independent of an explicit time t dependence and the final time t_c is free, the latter being a consequence of the transversality condition.□

Proposition 7. For all t ∈ [0, t_c], the adjoint λ(t) corresponding to the extremal lift (x_*(t), u_*(t), λ(t)) is nonzero.

Proof. This is a general result relating to free end time problems. We include a proof here for completeness. Suppose that there exists a time t ∈ [0, t_c] such that λ(t) = 0. By (22), the corresponding value of the Hamiltonian is H(λ₀, λ(t), x_*(t), u_*(t)) = −λ₀. By item (4) in Theorem 6, H ≡ 0, which implies that λ₀ = 0. This contradicts item (1) in Theorem 6. Hence, λ(t) ≠ 0 on [0, t_c].□

4.3. Geometric Properties and Existence of Singular Arcs

We now undertake a geometric analysis of the optimal control problem utilizing the affine structure of system (8) for interior states (i.e., controls which satisfy Theorem 6). We call such controls interior extremals, and all extremals in this section are assumed to be interior. The following results depend on the independence of the vector fields f and g, which we use to both classify the control structure for abnormal extremal lifts (extremal lifts with λ₀ = 0), as well as characterize the switching function dynamics via the Lie bracket.

Proposition 8. For all x₁ ∈ Ω, x₁ > 0, the vector fields f(x) and g(x) are linearly independent.

Proof. Define A(x) = A(x₁, x₂) to be the matrix

$\begin{array}{ll}\begin{array}{l}A(x)=(f(x)g(x))\\ =\left(\begin{array}{cc}(1-(x_1+x_2)-ϵ)x_1& -(\alpha +d)x_1\\ p_r(1-(x_1+x_2))x_2+ϵx_1& \alpha x_1\end{array}\right).\end{array}& (28)\end{array}$

The determinant of A can calculated as

$\begin{array}{ll}detA(x)=\alpha x_1^2\kappa (x)+p_r(\alpha +d)x_2{x}_1\kappa (x)+ϵd{x}_1^2& (29)\end{array}$

where

$\begin{array}{ll}\kappa (x):=1-(x_1+x_2).& (30)\end{array}$

As x₁(t) + x₂(t) ≤ 1 for all t ≥ 0, κ(x(t)) ≥ 0, and we see that detA(x) = 0 in Ω if and only if x₁ = 0, completing the proof.□

The line x₁ = 0 is invariant in Ω, and furthermore the dynamics in the set are independent of the control u(t). Conversely, $x_1^0>0$ implies that x₁(t) > 0 for all t ≥ 0. We concern our analysis only in this latter case, and so without loss of generality, f(x) and g(x) are linearly independent in the region of interest Ω_c.

We begin by showing that abnormal extremal lifts are easily characterized. We recall that an extremal lift is abnormal if λ₀ = 0, i.e., if the Hamiltonian is independent of the objective.

THEOREM 9. Abnormal extremal lifts at interior points, i.e., extremal lifts corresponding to λ₀ = 0, are constant and given by the maximal (M) or minimal (0) dosage.

Proof. Assume that u_* switches values at some time t. From (25), we must have that Φ(t) = 0. Since λ₀ = 0 and Φ(t) = 〈λ(t), g(x_*(t))〉, Equation (22) reduces to

$\begin{array}{ll}H(t)=\langle \lambda (t),f(x_{*}(t))\rangle =0.& (31)\end{array}$

Thus, λ(t) is orthogonal to both f(x_*(t)) and g(x_*(t)). Since f and g are linearly independent (Proposition 8), this implies that λ(t) = 0. But this contradicts Proposition 7. Hence, no such time t exists, and u_*(t) is constant. The constant sign of Φ thus corresponds to u = 0 or u = M (see Equation 25).□

The control structure for abnormal extremal lifts is then completely understood via Theorem 9. To analyze the corresponding behavior for normal extremal lifts, without loss of generality we assume that λ₀ = 1. Indeed, λ(t) may be rescaled by λ₀ > 0 to yield an equivalent version of Theorem 6. We thus assume that the Hamiltonian H(t) evaluated along (λ(t), x_*(t), u_*(t)) is of the form

$\begin{array}{ll}H(t)=-1+\langle \lambda (t),f(x_{*}(t))\rangle +u_{*}(t)\Phi (t)\equiv 0.& (32)\end{array}$

We recall the Lie bracket as the first-order differential operator between two vector fields X₁ and X₂:

$\begin{array}{ll}\left[X_1,X_2\right](z)=D{X}_2(z)X_1(z)-D{X}_1(z)X_2(z),& (33)\end{array}$

where, for example, DX₂(z) denotes the Jacobian of X₂ evaluated at z. As f and g are linearly independent in Ω, there exist γ, β ∈ C^∞(Ω) such that

$\begin{array}{ll}\left[f,g\right](x)=\gamma (x)f(x)+\beta (x)g(x),& (34)\end{array}$

for all x ∈ Ω. Explicitly, we compute γ and β:

$\begin{array}{ll}\gamma (x)=-\frac{(\alpha +d)x_1^2}{detA(x)}(a{x}_1+b{x}_2-c),& (35)\end{array}$

$\begin{array}{l}\beta (x)=\frac{x_1^2}{\det A(x)}(\alpha (1-p_r)\kappa (x)(\kappa (x)-ϵ)+ϵd(x_1+p_r{x}_2\\ +\kappa (x)-ϵ)),& (36)\end{array}$

where

$\begin{array}{ll}a=\alpha \left((1-p_r)+\frac{d}{\alpha +d}\right),& (37)\end{array}$

$\begin{array}{ll}b=\alpha (1-p_r)+d{p}_r,& (38)\end{array}$

$\begin{array}{ll}c=\alpha (1-p_r)+ϵd.& (39)\end{array}$

Clearly, for parameter values of interest (recall 0 < p_r < 1), a, b, c > 0. The assumption (5) guarantees that β(x) > 0 on 0 < x₁ + x₂ < V_c.

From (25), the sign of the switching function Φ determines the value of the control u_*. As λ and x_* are solutions of differential equations, Φ is differentiable. The dynamics of Φ can be understood in terms of the Lie bracket [f, g]:

$\begin{array}{ll}\dot{\Phi }(t)=\frac{d}{dt}\langle \lambda (t),g(x_{*}(t))\rangle & (40)\end{array}$

$\begin{array}{ll}=\gamma (x_{*}(t))\langle \lambda (t),f(x_{*}(t))\rangle +\beta (x_{*}(t))\Phi (t).& (41)\end{array}$

The last lines of the above follow from (34) as well as the linearity of the inner product. We are then able to derive an ODE system for x_* and Φ. Equation (32) allows us to solve for 〈λ(t), f(x_*(t))〉:

$\begin{array}{ll}\langle \lambda (t),f(x_{*}(t))\rangle =1-u_{*}(t)\Phi (t).& (42)\end{array}$

Substituting the above into (41) then yields the following ODE for Φ(t), which we view as coupled to system (8) via (25):

$\begin{array}{ll}\dot{\Phi }(t)=\gamma (x_{*}(t))+\left(\beta (x_{*}(t))-u_{*}(t)\gamma (x_{*}(t))\right)\Phi (t).& (43)\end{array}$

The structure of the optimal control at interior points may now be characterized as a combination of bang-bang and singular arcs. We recall that the control (or, more precisely, the extremal lift) u_* is singular on an open interval I ⊂ [0, t_c] if the switching function Φ(t) and all its derivatives are identically zero on I. On such intervals, Equation (25) does not determine the value of u_*, and a more thorough analysis of the zero set of Φ(t) is necessary. Indeed, for a problem such as ours, aside from controls determined by the path constraint ψ(x₁(t), x₂(t)) ≤ 0, singular arcs are the only candidates for optimal controls that may take values outside of the set {0, M}. Conversely, times t where Φ(t) = 0 but Φ⁽ⁿ⁾(t) ≠ 0 for some n ≥ 1 denote candidate bang-bang junctions, where the control may switch between the vertices 0 and M of the control set U. Note that the parity of the smallest such n determines whether a switch actually occurs: n odd implies a switch, while for n even u_* remains constant. Equation (43) allows us to completely characterize the regions in the (x₁, x₂) plane where singular arcs are attainable, as demonstrated in the following proposition.

Proposition 10. Singular arcs are only possible in regions of the (x₁, x₂) plane where γ(x) = 0. Furthermore, as x₁(t) > 0 for all t ≥ 0, the region {x ∈ $\text{R}$² | γ (x) = 0} ∩ Ω is the line

$\begin{array}{ll}a{x}_1+b{x}_2-c=0,& (44)\end{array}$

where a, b, c are defined in (37–39).

Proof. As discussed prior to the statement of Proposition 10, a singular arc must occur on a region where both Φ(t) and $\dot{\Phi }(t)$ are identically zero (as well as all higher-order derivatives). Denoting by x_*(t) the corresponding trajectory in the (x₁, x₂) phase plane, we may calculate $\dot{\Phi }(t)$ from equation (43):

$\begin{array}{ll}\dot{\Phi }(t)=\gamma (x_{*}(t)).& (45)\end{array}$

Note we have substituted the assumption Φ(t) = 0. Clearly we must also have that γ(x_*(t)) = 0, thus implying that $x_{*}(t)\in {\gamma }^{-1}(0)$, as desired. The last statement of the proposition follows immediately from Equation (35).□

Proposition 10 implies that singular solutions can only occur along the line ax₁ + bx₂ − c = 0. Thus, define regions in the first quadrant as follows:

$\begin{array}{ll}{\Omega }_c^{+}:=\{x\in \Omega |\gamma (x)>0\},& (46)\end{array}$

$\begin{array}{ll}{\Omega }_c^{-}:=\{x\in \Omega |\gamma (x)<0\},& (47)\end{array}$

$\begin{array}{ll}{L}=\{x\in \Omega |\gamma (x)=0\}.& (48)\end{array}$

Recall that Ω_c is simply the region in Ω prior to treatment failure, i.e., 0 ≤ V ≤ V_c, x₁, x₂ ≥ 0. From (35), Ω_c is partitioned as in Figure 3B. From (35) and (37–39), ${L}$ is a line with negative slope −b/a. Furthermore, necessary and sufficient conditions for ${L}$ to lie interior to Ω_c are $\frac{c}{a},\frac{c}{b}\le V_c$. From (37)–(39), this occurs if and only if

$\begin{array}{l}ϵ\le \min \{\frac{\alpha }{\alpha +d}-\frac{1-V_c}{d}\left(\alpha (1-p_r)+\frac{\alpha d}{\alpha +d}\right),\\ p_r-\frac{1-V_c}{d}\left(\alpha (1-p_r)+d{p}_r\right)\}.& (49)\end{array}$

FIGURE 3

Figure 3. Domain in (x₁, x₂) plane. (A) Region where γ changes sign. We see that inside the triangular region x₁ + x₂ ≤ 1 of the first quadrant, γ changes sign only along the line ax₁ + bx₂ − c = 0. For this line to be interior to Ω_c as depicted, we must be in the parameter regime indicated in (49). X and Y vector fields corresponding to vertices of control set U. For singular controls to lie in U, X and Y must point to opposite sides along ${L}$. (B) Same as in (A), but with α = 0.

As we have assumed that ϵ is small, and that V_c ≈ 1, this inequality is not restrictive, and we assume it is satisfied for the remainder of the work. We note an important exception below: when α = 0 the inequality is never satisfied with ϵ > 0; for such parameter values, line ${L}$ is horizontal (Figure 3B). We note that this does not change the qualitative results presented below. Of course, other configurations of the line ax₁ + bx₂ = c and hence precise optimal syntheses may exist, but we believe the situation illustrated in Figure 3A is sufficiently generic for present purposes.

With the existence of singular arcs restricted to the line γ = 0 by Proposition 10, we now investigate the feasibility of such solutions. Recall that the treatment u(t) must lie in the control set U = [0, M], for some M > 0 corresponding to the maximally tolerated applied dosage. Defining the vector field X(x) and Y(x) as the vector fields corresponding to the vertices of U,

$\begin{array}{ll}\begin{array}{l}X(x):=f(x),\\ Y(x):=f(x)+Mg(x),\end{array}& (50)\end{array}$

a singular control takes values in U at $x\in {L}$ if and only if X(x) and Y(x) point in different directions along ${L}$. More precisely, the corresponding Lie derivatives L_Xγ(x) and L_Yγ(x) must have opposite signs (see Figure 3A). The following proposition determines parameter values where this occurs.

Proposition 11. Suppose that α > 0, so that drug has the potential to induce resistance. Also, let the maximally tolerated dosage M satisfy

$\begin{array}{l}M>\frac{\alpha +d}{\alpha (\alpha +d)+\alpha d}(d\left(\frac{\alpha }{\alpha +d}-ϵ\right)+ϵd(p_r-\alpha )\\ -2\alpha d(1-p_r)).& (51)\end{array}$

Then the following hold along ${L}$:

1. L_Xγ < 0,

2. L_Yγ < 0 as $(x_1,x_2)\to \left(0,\frac{c}{b}\right)$ in Ω,

3. L_Yγ > 0 at $(x_1,x_2)=\left(\frac{c}{a},0\right)$, and

4. L_Yγ is monotonically decreasing as a function of x₁.

Thus, ${L}$ contains a segment $\bar{L}\subset {L}$ which is a singular arc. Note that $\bar{L}$ is precisely the region in ${L}$ where L_Yγ is positive.

Proof. The proof is purely computational.□

Note that if inequality (51) is not satisfied, then singular arcs are not in the domain Ω_c.

The geometry of Proposition 11 is illustrated in Figure 4. Thus, assuming α > 0 and M as in (51), singular arcs exist along the segment $\bar{L}\subset {L}$. Furthermore, the corresponding control has a unique solution u_s, which may be computed explicitly. Indeed, as the solution must remain on the line ${L}$, or equivalently, ax₁ + bx₂ = c, taking the time derivative of this equation yields aẋ₁ + bẋ₂ = 0, and substituting the expressions (1) we compute u_s as

$\begin{array}{ll}{u}_s(t)=\frac{\left(1-(x_1(t)+x_2(t))\right)\left(a{x}_1(t)+p_{r}b{x}_2(t)\right)+ϵ(b-a)x_1(t)}{2\alpha (1-p_r)d{x}_2(t)},& (52)\end{array}$

where a, b, c are given by (37–39) and x₂ and x₁ satisfy ax₁ + bx₂ = c. As discussed previously, x₁(t) > 0 for $x_1^0>0$, so this formula is well-defined. Proposition 11 implies that it is possible to simplify Equation (52) as a function of x₁ (i.e. as a feedback law) for $x_1\in (\bar{s},\frac{c}{a})$, for some $\bar{s}>0$, but since its value will not be needed, we do not provide its explicit form. Note that the maximal dose M is achieved precisely at $x_1=\bar{s}$ where vector field Y is parallel to ${L}$. Thus, at this $\bar{s}$, the trajectory must leave the singular arc, and enter the region ${\Omega }_c^{-}$. As ẋ₂ ≥ 0, trajectories must follow ${L}$ in the direction of decreasing x₁ (see Figure 4). We summarize these results in the following theorem.

FIGURE 4

Figure 4. Geometry of vector fields X and Y with α > 0 and M satisfying (51). As in Proposition 11, this can be understood via the corresponding Lie derivatives of γ. Note that near x₂ = 0, X, and Y point to opposite sides of L, while at $(x_1,x_2)=\left(0,\frac{c}{b}\right)$, both X and Y point away from γ > 0. The line $\overline{\mathit{\text{L}}}$ is the unique singular arc in Ω_c.

THEOREM 12. If α > 0, and M satisfies (51), a singular arc exists in the (x₁, x₂) plane as a segment of the line ${L}$. Along this singular arc, the control is given by Equation (52), where ax₁ + bx₂ = c. Therefore, in this case the necessary minimum conditions on u_* from (25) can be updated as follows:

$\begin{array}{ll}{u}_{*}(t)=\{\begin{array}{ll}0\hfill & \Phi (t)>0,\hfill \\ M\hfill & \Phi (t)<0,\hfill \\ u_s(t),\hfill & \Phi (t)\equiv 0\text{ for }t\in I,\hfill \end{array}& (53)\end{array}$

where I is an open interval. Recall again that this is the optimal control at points interior to Ω_c.

Proof. See the discussion immediately preceding Theorem 12.□

In the case α = 0, the line ${L}$ is horizontal, and as x₂ is increasing, no segment $\bar{L}\subseteq {L}$ is admissible in phase space. Thus, the interior controls in this case are bang-bang; for a visualization (see Figure 3B).

THEOREM 13. If α = 0, there are no singular arcs for the optimal time problem presented in section 3. Thus, the interior control structure is bang-bang.

Outside of the singular arc $\bar{L}$, the control structure is completely determined by (25) and (43). The precise result, utilized later for the optimal synthesis presented in section 5, is stated in the following theorem. We first introduce a convenient (and standard) notation. Let finite words on X and Y denote the concatenation of controls corresponding to vector fields X (u ≡ 0) and Y (u ≡ M), respectively. The order of application is read left-to-right, and an arc appearing in a word may not actually be applied (e.g. XY denotes an X arc followed by a Y arc or a Y arc alone).

THEOREM 14. Consider an extremal lift Γ = ((x, u), λ). Trajectories x remaining entirely in ${\Omega }_c^{+}$ or ${\Omega }_c^{-}$ can have at most one switch point. Furthermore, if $x\in {\Omega }_c^{+}$, then the corresponding control is of the form YX. Similarly, $x\in {\Omega }_c^{-}$ implies that u = XY. Hence multiple switch points must occur across the singular arc $\bar{L}$.

Proof. If τ is a switching time, so that Φ(τ) = 0, Equation (43) allows us to calculate $\dot{\Phi }(\tau )$ as

$\begin{array}{ll}\dot{\Phi }(\tau )=\gamma (x(\tau )).& (54)\end{array}$

Thus, in ${\Omega }_c^{+}$ where γ > 0, $\dot{\Phi }(\tau )>0$, and hence Φ must increase through τ. The expression for the control (25) then implies that a transition from a Y-arc to an X-arc occurs at τ (i.e., a YX arc). Furthermore, another switching time cannot occur unless x leaves ${\Omega }_c^{+}$, since otherwise there would exist a $\bar{\tau }>\tau$ such that $\Phi (\bar{\tau })=0,\dot{\Phi }(\bar{\tau })<0$ which is impossible in ${\Omega }_c^{+}$. Similarly, only XY-arcs are possible in ${\Omega }_c^{-}$.□

The structure implied by Theorem 14 is illustrated in Figure 4. Note that inside the sets ${\Omega }_c^{+},{\Omega }_c^{-}$, and $\bar{L}$, extremal lifts are precisely characterized. Furthermore, the results of section 4.1 (and particularly Equation 18) yield the characterization on the boundary N. What remains is then to determine the synthesis of these controls to the entire domain Ω_c, as well as to determine the local optimality of the singular arc $\bar{L}$. The latter is addressed in the following section.

4.4. Optimality of Singular Arcs

We begin by proving that the singular arc is extremal, i.e. that it satisfies the necessary conditions presented in section 4.2 (note that it is interior by assumption). This is intuitively clear from Figure 4, since X and Y point to opposite sides along $\bar{L}$ by the definition of ${L}$.

THEOREM 15. The line segment $\bar{L}\subset {L}$ is a singular arc.

Proof. We find an expression for u = u(x) such that the vector f(x) + u(x)g(x) is tangent to $\bar{L}$ at x, i.e. we find the unique solution to

$\begin{array}{ll}{L}_{f+ug}(\gamma )=0& (55)\end{array}$

Note that we can invert (50):

$\begin{array}{ll}\begin{array}{l}f(x)=X(x)\\ g(x)=\frac{1}{M}(Y(x)-X(x))\end{array}& (56)\end{array}$

so that $f+ug=(1-\frac{u}{M})X+\frac{u}{M}Y$. Thus,

$\begin{array}{l}{L}_{f+ug}(\gamma )=\left(1-\frac{u}{M}\right)L_X\gamma +\frac{u}{M}{L}_Y\gamma \end{array}$

Setting the above equal to zero, and solving for u = u(x) yields

$\begin{array}{ll}u(x)=M\frac{L_X\gamma (x)}{L_X\gamma (x)-L_Y\gamma (x)}& (57)\end{array}$

As L_Xγ < 0 and L_Yγ > 0 on $\overline{\mathit{\text{L}}}$ by Proposition 11, we see that 0 < u(x) < M. We must also verify that the associated controlled trajectory (57) is extremal by constructing a corresponding lift. Suppose that x(t) solves

$\begin{array}{l}\dot{x}=f(x)+u(x)g(x),\\ x(0)=q,\end{array}$

for $q\in \bar{L}$. Let ϕ ∈ ($\text{R}$²)* such that

$\begin{array}{l}\langle \varphi ,g(q)\rangle =0,\langle \varphi ,f(q)\rangle =1.\end{array}$

Let λ(t) solve the corresponding adjoint Equation (24) with initial condition λ(0) = ϕ. Then the extremal lift Γ = ((x, u), λ) is singular if Φ(t) = 〈λ(t), g(x(t))〉 ≡ 0. By construction of u(x), the trajectory remains on $\bar{L}$ on some interval containing zero, and we can compute $\dot{\Phi }$ as [using (34)]

$\begin{array}{l}\dot{\Phi }(t)=\langle \lambda (t),[f,g](x(t))\rangle \\ =\gamma (x(t))\langle \lambda (t),f(x(t)\rangle +\beta (x(t))\langle \lambda (t),g(x(t))\rangle \\ =\beta (x(t))\Phi (t),\end{array}$

Note that we have used (43) and the fact that γ = 0 by our choice of u. Since Φ(0) = 0 by hypothesis, this implies that Φ(t) ≡ 0, as desired.□

The above then verifies that $\bar{L}$ is a singular arc. Note that an explicit expression for u = u(x) was given in (52), which can be shown to be equivalent to (57).

Having shown that the singular arc $\bar{L}$ is extremal, we now investigate whether it is locally optimal for our time-optimization problem. The singular arc is of intrinsic order k if the first 2k − 1 derivatives of the switching function are independent of u and vanish identically on an interval I, while the 2kth derivative has a linear factor of u. We can compute [this is standard for control-affine systems (8)] that

$\begin{array}{ll}{\Phi }^{2k}(t)=\langle \lambda (t),{\text{ad}}_f^{2k}(g)(x(t))\rangle +u(t)\langle \lambda (t),\left[g,{\text{ad}}_f^{2k-1}(g)\right](x(t))\rangle ,& (58)\end{array}$

where ad_Z is the adjoint endomorphism for a fixed vector field Z:

$\begin{array}{ll}{\text{ad}}_Z(V)=\left[Z,V\right],& (59)\end{array}$

and powers of this operator are defined as composition. Fix an extremal lift Γ = ((x, u), λ) of a singular arc of order k. The Generalized Legendre-Clebsch condition (also known as the Kelley condition) (Ledzewicz and Schättler, 2012) states that a necessary condition for Γ to satisfy a minimization problem with corresponding Hamiltonian H is that

$\begin{array}{ll}{(-1)}^k\frac{\partial }{\partial u}\frac{d^{2k}}{d{t}^{2k}}\frac{\partial H}{\partial u}({\lambda }_0,\lambda (t),x(t),u(t))\ge 0& (60)\end{array}$

along the arc. Note that $\frac{\partial H}{\partial u}=\Phi$, so that the above is simply the u coefficient of the 2k-th time derivative of the switching function (multiplied by (−1)^k). The order of the arc, as well as the Legendre-Clebsch condition, are addressed in Theorem 16.

THEOREM 16. The singular control is of order one. Furthermore, for all times t such that $x(t)\in \bar{L}$, 〈λ(t), [g, [f, g]](x(t))〉 > 0. Thus, the Legendre-Clebsch condition is violated, and the singular arc $\bar{L}$ is not optimal.

Proof. Along singular arcs we must have $\Phi (t),\dot{\Phi }(t),\ddot{\Phi }(t)\equiv 0$, and we can compute these derivatives using iterated Lie brackets as follows:

$\begin{array}{ll}\begin{array}{l}\Phi (t)=\langle \lambda (t),g(x(t))\rangle ,\\ \dot{\Phi }(t)=\langle \lambda (t),\left[f,g\right](x(t))\rangle ,\\ \ddot{\Phi }(t)=\langle \lambda (t),\left[f+ug,\left[f,g\right]\right](x(t))\rangle .\end{array}& (61)\end{array}$

The final of the above in (61) can be simplified as

$\begin{array}{ll}\ddot{\Phi }(t)=\langle \lambda (t),\left[f,\left[f,g\right]\right](x(t))\rangle +u(t)\langle \lambda (t),\left[g,\left[f,g\right]\right](x(t))\rangle \equiv 0,& (62)\end{array}$

which is precisely (58) for k = 1. Order one is then equivalent to being able to solve this equation for u(t). Thus, 〈λ(t), [g, [f, g]](x(t))〉 > 0 will imply that the arc is singular of order one. We directly compute 〈λ(t), [g, [f, g]](x(t))〉 = 〈λ(t), [g, ad_f(g)](x(t))〉. Using Equation (34) and recalling properties of the singular arc [γ = 0 and the remaining relations in (61), as well as basic “product rule” properties of the Lie bracket], we can show that

$\begin{array}{ll}\left[g,\left[f,g\right]\right]=(L_g\gamma )f-\gamma \left[f,g\right]+(L_g\beta )g.& (63)\end{array}$

Recall that for an extremal lift along the arc $\bar{{L}}$,

$\begin{array}{ll}\begin{array}{l}\langle \lambda (t),g(x(t))\rangle \equiv 0,\\ \langle \lambda (t),\left[f,g\right](x(t))\rangle \equiv 0\\ \langle \lambda (t),f(x(t))\rangle \equiv 1.\end{array}& (64)\end{array}$

The first two of the above follow from $\Phi ,\dot{\Phi }\equiv 0$, and the third is a consequence of H ≡ 0 [see (22)]. Equations (63) and (64) together imply that

$\begin{array}{ll}\begin{array}{l}\langle \lambda (t),\left[g,\left[f,g\right]\right](x(t))\rangle =L_g\gamma \langle \lambda (t),f(x(t))\rangle -\gamma \langle \lambda (t),\left[f,g\right](x(t))\rangle \\ +L_g\beta \langle \lambda (t),g(x(t))\rangle \\ =L_g\gamma (x(t))\\ =\frac{1}{M}(L_Y\gamma (x(t))-L_X\gamma (x(t))).\end{array}& (65)\end{array}$

The last equality follows from the representation in (56). As L_Yγ > 0 and L_Xγ < 0 along $\overline{\mathit{\text{L}}}$ (Proposition 11), 〈λ(t), [g, [f, g]](x(t))〉 > 0, as desired. Furthermore,

$\begin{array}{ll}-\langle \lambda (t),\left[g,\left[f,g\right]\right](x(t))\rangle <0,\text{ or equivalently}& (66)\end{array}$

$\begin{array}{ll}{(-1)}^1\frac{\partial }{\partial u}\frac{d^2}{d{t}^2}\frac{\partial H}{\partial u}<0,& (67)\end{array}$

showing that (60) is violated (substituting k = 1). Thus, $\bar{L}$ is not optimal.□

Theorem 16 then implies that the singular arc is suboptimal, i.e. that $\bar{L}$ is “fast” with respect to the dynamics. In fact, we can compare times along trajectories using the “clock form,” a one-form on Ω. As one-forms correspond to linear functionals on the tangent space, and f and g are linearly independent, there exists a unique ω ∈ (TΩ)^* such that

$\begin{array}{ll}{\omega }_x(f(x))\equiv 1,{\omega }_x(g(x))\equiv 0.& (68)\end{array}$

In fact, we compute it explicitly:

$\begin{array}{ll}{\omega }_x=\frac{g_2(x)d{x}^1-g_1(x)d{x}^2}{det(f(x),g(x))}.& (69)\end{array}$

Then, along any controlled trajectory (x, u) defined on [t₀, t₁], the integral of ω computes the time t₁ − t₀:

$\begin{array}{l}{\displaystyle {\int }_x\omega }={\displaystyle {\int }_{t_0}^{t_1}{\omega }_{x(t)}}(\dot{x}(t))\text{\hspace{0.05em}d}t\\ ={\displaystyle {\int }_{t_0}^{t_1}{\omega }_{x(t)}}(f(x(t))+u(t)g(x(t))))\text{\hspace{0.05em}d}t\\ ={\displaystyle {\int }_{t_0}^{t_1}{\omega }_{x(t)}}(f(x(t))\text{\hspace{0.05em}d}t+{\displaystyle {\int }_{t_0}^{t_1}u}(t){\omega }_{x(t)}(g(x(t))))\text{\hspace{0.05em}d}t\\ ={\displaystyle {\int }_{t_0}^{t_1}}\text{d}t\\ =t_1-t_0.& (70)\end{array}$

We can then use ω and Stokes' Theorem to compare bang-bang trajectories with those on the singular arc. See Figure 5 below for a visualization of a singular trajectory connecting $q_1,q_2\in \bar{L}$ and the corresponding unique XY trajectory connecting these points in ${\Omega }_c^{-}$ (note that uniqueness is guaranteed as long as q₁ and q₂ are sufficiently close).

FIGURE 5

Figure 5. Both XY and singular trajectories taking q₁ to q₂.

Let t_S denote the time spent along the singular arc, t_X the time spent along the X arc, and t_Y the time spent along the Y arc. Denote by Δ the closed curve traversing the X and Y arcs positively and the singular arc negatively, with R as its interior. As X and Y are positively oriented (they have the same orientation as f and g), Stokes' Theorem yields

$\begin{array}{ll}{t}_X+t_Y-t_S={\displaystyle {\int }_{\Delta }}\omega ={\displaystyle {\int }_R}d\omega & (71)\end{array}$

Taking the exterior derivative yields the two-form dω see Chapter 2 of (Ledzewicz and Schättler, 2012):

$\begin{array}{ll}d\omega =-\frac{\gamma }{det(f,g)}.& (72)\end{array}$

As the determinant is everywhere positive (see the proof of Proposition 8), and R lies entirely in γ < 0, the integral on the right-hand side of (71) is positive, so that we have

$\begin{array}{ll}{t}_S<t_X+t_Y& (73)\end{array}$

Thus, time taken along the singular arc is shorter than that along the XY trajectory, implying that the singular arc is locally suboptimal for our problem (recall that we want to maximize time). Since local optimality is necessary for global optimality, trajectories should never remain on the singular arc for a measurable set of time points. This reaffirms the results of Theorem 16. A completely analogous statement holds for YX trajectories in the region γ > 0. We can also demonstrate, utilizing the same techniques, that increasing the number of switchings at the singular arc speeds up the trajectory (see Figure 6). This again reinforces Theorem 16, and implies that trajectories should avoid the singular arc to maximize the time spent in Ω_c.

FIGURE 6

Figure 6. XY (solid) and XYXY (dashed) trajectories taking q₁ to q₂ in the region γ > 0. The time difference between the two trajectories can again be related to the surface integral in the region R, where γ < 0. The XY trajectory can then be seen to be slower in comparison.

5. Characterization of Optimal Control

The results of sections 4.1, 4.2, 4.3, and 4.4 may now be combined to synthesize the optimal control introduced in section 3.

THEOREM 17. For any α ≥ 0, the optimal control to maximize the time to reach a critical time is a concatenation of bang-bang and path-constraint controls. In fact, the general control structure takes the form

$\begin{array}{ll}{(YX)}^n{u}_{p}Y& (74)\end{array}$

where (YX)ⁿ: = (YX)ⁿ⁻¹YX for n ∈ $\text{N}$, and the order should be interpreted as left to right. Here u_p is defined in (18).

Proof. Formula (74) is simply a combination of the results presented previously. Note that singular arcs are never (locally) optimal, and hence do not appear in the equation. We also observe that X arcs are not admissible once the boundary N has been obtained, as an X arc always increases V. A Y arc may bring the trajectory back into int(Ω_c), but a YX trajectory is no longer admissible, as the switching structure in ${\Omega }_c^{-}$ is XY (Theorem 14).

The only aspect that remains is to show that once N is reached, the only possible trajectories are either u_p given by (18) or Y, with at most one switching occurring between the two. That is, a local arc of the form u_pYu_p is either sub-optimal or non-feasible (equivalently, outside of the control set U). Suppose that such an arc is feasible, i.e., that for all such points in phase space, 0 ≤ u_p ≤ M [recall that u_p is defined via feedback in (18)]. Denote by τ₁ and τ₂ the times at which the switch onto and off of Y occurs, respectively. Since u_p decreases with S, feasibility implies that u_p(t) ≤ M for all t ∈ [τ₁, τ₂]. Thus, we can consider the alternate feasible trajectory which remains on N between the points (S(τ₁), R(τ₁)) and (S(τ₂), R(τ₂)); see Figure 7 for an illustration. Call τ the time for such a trajectory. Then, using the clock-form ω and the positively-oriented curve Δ which follows N first and Y (in the reverse direction) second, we obtain similarly to (71),

$\begin{array}{ll}\tau -({\tau }_2-{\tau }_1)=-{\displaystyle {\int }_R}\frac{\gamma }{det(f,g)},& (75)\end{array}$

where R: = int(Δ). Recalling that γ < 0 in R (see Figure 4), the previous equation implies that

$\begin{array}{ll}\tau >{\tau }_1-{\tau }_2,& (76)\end{array}$

i.e., a longer amount of time is achieved by remaining on the boundary N. Hence the arc u_pYu_p is sub-optimal if it is feasible, as desired.

FIGURE 7

Figure 7. Comparison of u_pYu_p arc and an arc that remains on N (hence u ≡ u_p) between the points [S(τ₁), R(τ₁)] and [S(τ₂), R(τ₂)], assuming that u_p remains feasible (that is, u_p ∈ [0, M]). Note that γ < 0 in the area of interest, and that a switching of a Y to an X arc is prohibited via the Maximum Principle. Thus, the only possibility is the curve illustrated, which leaves the boundary N for a Y arc before u_p becomes infeasible.

The previous argument has one subtle aspect, as we used results from the Maximum Principle on the boundary set N, where technically it does not apply. However, the above still remains true, since we may approximate the boundary line V = V_c with a curve interior to Ω_c which remains feasible. By continuity, the time along such a curve can be made arbitrarily close to τ, and hence is still greater than τ₂ − τ₁, implying that u_pYu_p is sub-optimal.□

Note that in Theorem 17, the switchings must occur across the singular arc $\bar{L}$, if it exists (recall that it is not admissible if α = 0). The control u_p is determined along the boundary of Ω_c, and provides the synthesis between interior and boundary controls.

We finally include a technical result, which eliminates the optimality of the constrained (boundary) control u_p in certain cases.

Proposition 18. Assume that the maximal dose M is as in Proposition 2:

$\begin{array}{ll}M>\frac{(1-V_c)(1-p_r)}{d}& (77)\end{array}$

If the optimal control becomes maximal in ${\Omega }_c^{-}$ (i.e., u = M in this region), then the control cannot take the boundary value u_p (Equation 18) on an interval. Equivalently, an optimal control cannot end in the form Yu_p.

Proof. Note that if u_* = Y and reaches N at the point x, then the Lie derivative L_YV(x) must satisfy

$\begin{array}{ll}{L}_{Y}V(x)\ge 0& (78)\end{array}$

as V must be increasing along the Y vector field, since it reaches N. But by Proposition 2, this implies that

$\begin{array}{l}{x}_1\le x_1^{*}\end{array}$

Proposition 3 then implies that u_p is unfeasible in this region, completing the proof.□

6. Numerical Results

In this section, we provide numerical examples of the analytical results obtained in previous sections. All figures in this section were obtained using the GPOPS-II MATLAB software (Patterson and Rao, 2014). Parameters and initial values are given in Table 1 shown below, unless stated otherwise.

TABLE 1

Table 1. Parameter values and initial conditions used throughout section 6, unless stated otherwise.

Theorem 17 characterizes the qualitative form of the optimal control:

$\begin{array}{ll}{u}_{*}={(YX)}^n{u}_{p}Y,& (79)\end{array}$

where n is the number of interior switches, u_p the sliding control (18), and X and Y denote the lower and upper corner controls u = 0 and u = M, respectively. We begin by computing sample controls (see Figures 8, 10). Note that the optimal control in Figure 8B takes the form YXu_pY, while that of Figure 10B is an upper corner control Y. The phase plane dynamics corresponding to Figure 8 are also provided in Figure 9. In both cases the cytotoxic parameter was fixed at d = 0.05, while the induced rate of resistance α varies between α = 0.005 in Figure 8 and α = 0.1 in Figure 10. Note that for the smaller value of α (Figure 8), a longer period of treatment success is observed, as the time to treatment failure is approximately 70 time units; compare this with t_c = 24.2 in Figure 10. This result is intuitive, as the treatment less likely to induce resistance is able to be more effective when optimally applied.

FIGURE 8

Figure 8. Numerical solution of the optimal control problem with d = 0.05, α = 0.005, and the remainder of parameters as in Table 1. (A) Sensitive (x₁) and resistant (x₂) temporal dynamics. (B) Control structure of form YXu_pY. (C) Volume dynamics. Note that the trajectory remains on the line V = V_c for most times, with corresponding control u = u_p.

FIGURE 9

Figure 9. Phase plane corresponding to Figure 8. Trajectory which optimal control is of the form YXu_pY with parameter values as in Table 1 except for α = 0.005 and d = 0.05. The yellow dot in the figure represents the $(x_1^{*},x_2^{*})$ point at which Y(x) is tangent to the sliding surface. Here, $(x_1^{*},x_2^{*})=(0.1059,0.7941)$. As proven in Proposition 2, for points on the line N, the tumor volume will decrease along the Y(x) direction if x₁ > 0.1059 and will increase for x₁ < 0.1059.

FIGURE 10

Figure 10. Numerical solution of optimal control problem with d = 0.05, α = 0.1, and the remainder of parameters as in Table 1. (A) Sensitive (x₁), resistant (x₂), and volume (x₁ + x₂) temporal dynamics. (B) Control structure of form Y, i.e., an entirely upper corner control. (C) Phase plane dynamics, plotted with relevant vector fields.

The generality of the previous statement is investigated in Table 2 and Figures 11, 12. The computed optimal times t_c suggest that when the cytotoxicity of the drug (d) is small, higher induction rates (α) actually increase treatment efficacy. For example, for d = 0.001 treatment response increases as α increases (Figure 12A). This could be explained from the fact that sensitive cells have a higher growth rate than resistant cells (assumption p_r < 1). Thus, when the chemotherapeutic drug has a low effectiveness (small d) a larger α value actually helps to reduce the sensitive population size, and therefore extends the time t_c at which the tumor volume exceeds its critical value V_c.

TABLE 2

Table 2. Optimal time t_c for each of the computed controls appearing in Figure 11.

FIGURE 11

Figure 11. Optimal control structures for different α and d values. The blue curve is the computed optimal control, while the red curve is the feedback control along on the boundary of N, which may or may not be optimal or even feasible.

FIGURE 12

Figure 12. Variation in t_c as a function of α. (A) Treatment success time t_c for d = 0.001 with varying α values. (B) Functional dependence of t_c on α for different d parameters. Note that for small d, t_c increases as a function of α, but that this trend is reversed if d is further increased.

The situation is reversed when we consider larger values of d because in this case it would take more time for the tumor to grow to its critical volume V_c if the drug effectiveness is large enough; see for example the row d = 0.5 in Table 2, and the corresponding purple curve in Figure 12B. Figure 12B provides the critical time as a function of α for multiple cytotoxicities d; note the qualitative change in t_c as d increases.

Examining Figure 11 and Table 2 also suggests that as d increases, the feedback control u_p becomes optimal on an interval [t₁, t₂] with 0 < t₁ < t₂ < t_c. More numerical results are provided in section 7.3.

7. Additional Results

7.1. Structural Identifiability

For completeness, we discuss the identifiability of system (1). As our focus in this work has been on control structures based on the presence of drug-induced resistance, we rely on the ability to determine whether, and to what degree, the specific chemotherapeutic treatment is generating resistance.

Ideally, we envision a clinical scenario in which cancer cells from a patient are cultured in an ex vivo assay (for example, see Silva et al., 2017) prior to treatment. Parameter values are then calculated from treatment response dynamics in the assay, and an optimal therapy regime is implemented based on the theoretical work described below. Thus, identifying patient-specific model parameters, specially the induced-resistance rate α, is a necessary step in determining the control structures to apply. In this section, we address this issue, and prove that all parameters are structurally identifiable, as well as demonstrate a specific set of controls that may be utilized to determine α. A self-contained discussion is presented; for more details on theoretical aspects, see Sontag (2017) and the references therein. Other recent works related to identifiability in the biological sciences (as well as practical identifiability) can be found in Eisenberg and Jain (2017) and Villaverde et al. (2016).

We first formulate our dynamical system, and specify the input and output variables. The treatment u(t) is the sole input. Furthermore, we assume that the only clinically observable output is the net tumor volume V(t):

$\begin{array}{ll}V(t):=x_1(t)+x_2(t).& (80)\end{array}$

That is, we do not assume real-time measurements of the individual sensitive and resistant sub-populations. We note that in some instances, such measurements may be possible; however for a general chemotherapy, the precise resistance mechanism may be unknown a priori, and hence no biomarker with the ability to differentiate cell types may be available.

Treatment is initiated at time t = 0, at which we assume an entirely sensitive population:

$\begin{array}{ll}{x}_1(0)=x_1^0,x_2(0)=0.& (81)\end{array}$

Here $0<x_1^0<1$, so that (x₁(t), x₂(t)) ∈ Ω for all t ≥ 0. We note that x₂(0) = 0 is not restrictive, and similar results may derived under the more general assumption 0 ≤ x₂(0) < 1. The condition x₂(0) = 0 is utilized both for computational simplicity and since x₂(0) is generally small (assuming a non-zero detection time, and small drug-independent resistance parameter ϵ; see Greene et al., 2019 for a discussion).

As formulated in section 7.2.1, the above then allows us to formulate our system (1) in input/output form, where the input u(t) appears affinely:

$\begin{array}{ll}\begin{array}{l}\dot{x}(t)=f(x(t))+u(t)g(x(t)),\\ x(0)=x_0,\end{array}& (82)\end{array}$

where (as defined on Equations (9) and (10)) f and g are

$\begin{array}{ll}f(x)=\left(\begin{array}{c}(1-(x_1+x_2))x_1-ϵx_1\\ p_r(1-(x_1+x_2))x_2+ϵx_1\end{array}\right),& (83)\end{array}$

$\begin{array}{ll}g(x)=\left(\begin{array}{c}-(\alpha +d)\\ \alpha \end{array}\right)x_1,& (84)\end{array}$

and x(t) = (x₁(t), x₂(t)). As is standard in control theory, the output is denoted by the variable y, which in this work corresponds to the total tumor volume:

$\begin{array}{ll}\begin{array}{l}y(t,p):=h(x(t),u(t),p)\\ =x_1(t)+x_2(t).\end{array}& (85)\end{array}$

Note that x₁(t), x₂(t) depend on both the input u(t) and parameters p. A system in form (82) is said to be uniquely structurally identifiable if the map (u(t), p) ↦ (u(t), y(t, p)) is injective almost everywhere (Meshkat and Seth, 2014; Eisenberg and Jain, 2017), where p is the vector of parameters to be identified. In this work,

$\begin{array}{ll}p=(x_1^0,d,\alpha ,ϵ,p_r).& (86)\end{array}$

Local identifiability and non-identifiability correspond to the map being finite-to-one and infinite-to-one, respectively. Our objective is then to demonstrate unique structural identifiability for model system (82) [or equivalently (1)], and hence recover all parameter values p from only measurements of the tumor volume y. We also note that the notion of identifiability is closely related to that of observability; for details Anguelova (2004); Sontag (1979) are good references.

To analyze identifiability, we utilize results appearing in, for example (Hermann and Krener, 1977; Wang and Sontag, 1989; Sontag and Wang, 1991), and hence frame the issue from a differential-geometric perspective. Our hypothesis is that perfect (hence noise-free) input-output data is available and in particular, for differentiable controls, that we can compute y and its derivatives. We thus, for example, make measurements of

$\begin{array}{ll}\begin{array}{l}y(0)=h(x(0)),\\ \dot{y}(0)={\displaystyle \frac{d}{dt}{|}_{t=0}h(x(t))}\end{array}& (87)\end{array}$

for appropriately chosen inputs, and relate their values to the unknown parameter values p. If there exist inputs u(t) such that the above system of equations may be solved for p, the system is identifiable. The right-hand sides of (87) may be computed in terms of the Lie derivatives of the vector fields f and g in system (82). We recall the definition of Lie differentiation L_XH of a C^ω function H by a C^ω (i.e. real-analytic) vector field X:

$\begin{array}{ll}{L}_{X}H(x):=\nabla H(x)·X(x).& (88)\end{array}$

Here the domain of both X and H is the first-quadrant triangular region Ω, seen as a subset of the plane, and the vector fields and output function are C^ω on an open set containing Ω (in fact, they are given by polynomials, so they extend as analytic functions to the entire plane). Iterated Lie derivatives are well-defined, and should be interpreted as function composition, so that for example L_YL_XH = L_Y(L_XH), and $L_X^{2}H=L_X(L_{X}H)$.

More formally, let us introduce the observable quantities corresponding to the zero-time derivatives of the output y = h(x),

$\begin{array}{ll}Y(x_0,U)={\displaystyle \frac{d^k}{d{t}^k}{|}_{t=0}h(x(t)),}& (89)\end{array}$

where U ∈ R^k is the value of the control u(t) (without loss of generality, a polynomial of degree k − 1) and its derivatives evaluated at t = 0: U = (u(0), u′(0), ..., u^(k−1)(0)). Here k ≥ 0, indicating that the k^th-order derivative Y may expressed as a polynomial in the components of U (Sontag and Wang, 1991). The initial conditions x₀ appear due to evaluation at t = 0. The observation space is then defined as the span of the elements Y(x₀, U):

$\begin{array}{ll}{F}_1:={\text{span}}_{\text{R}}\{Y(x_0,U)|U\in {\text{R}}^k,k\ge 0\}.& (90)\end{array}$

Conversely, we also define span of iterated Lie derivatives with respect to the output h and vector fields f(x) and g(x):

$\begin{array}{ll}{F}_2:={\text{span}}_{\text{R}}\left\{L_{i_1}\dots L_{i_k}h(x_0)|(i_1,\dots i_k)\in {\{g,f\}}^k,k\ge 0\right\}.& (91)\end{array}$

Wang and Sontag (1989) proved that F₁ = F₂, so that the set of “elementary observables” may be considered as the set of all iterated Lie derivatives F₂. Hence, identifiability may be formulated in terms of the reconstruction of parameters p from elements in F₂. Parameters p are then identifiable if the map

$\begin{array}{ll}p\mapsto \left(L_{i_1}\dots L_{i_k}h(x_0)|(i_1,\dots i_k)\in {\{g,f\}}^k,k\ge 0\right)& (92)\end{array}$

is one-to-one. For the remainder of this section, we investigate the mapping defined in (92).

Computing the Lie derivatives and recalling that x₀ = (S₀, 0) we can recursively determine the parameters p:

$\begin{array}{ll}\begin{array}{l}{x}_1^0=h(x_0),\\ d=-\frac{L_{g}h(x_0)}{x_1^0},\\ \alpha =\frac{L_g^{2}h(x_0)}{d{x}_1^0}-d,\\ ϵ=\frac{L_f{L}_{g}h(x_0)}{d{x}_1^0}+1-x_1^0,\\ p_r=\frac{x_1^0}{1-x_1^0}+\frac{L_g{L}_{f}h(x_0)}{\alpha x_1^0(1-x_1^0)}-\left(1+\frac{d}{\alpha }\right)\left(1-\frac{x_1^0}{1-x_1^0}\right).\end{array}& (93)\end{array}$

Since F₁ = F₂, all of the above Lie derivatives are observable via appropriate treatment protocols. For an explicit set of controls and corresponding relations to measurable quantities [elements of the form (89)], see Greene et al. (2019). Thus, we conclude that all parameters in system (1) are identifiable, which allows us to investigate optimal therapies dependent upon a priori knowledge of the drug-induced resistance rate α.

7.2. Existence Results

For the problem presented in section 3, we are going to verify that the supremum of times t_c(u) for $u\in {U}$ [with t_c(u) as defined in Equation (6)] is obtained by some $u_{*}\in {U}$, i.e., that an optimal control exists. This involves two distinct steps: (1) proving that the supremum is finite, and (2) that it is obtained by at least one admissible control. The following two subsections verify these claims.

7.2.1. Finiteness of the Supremum

We prove that

$\begin{array}{ll}{\displaystyle \underset{u\in {U}}{sup}}{t}_c(u)<\infty & (94)\end{array}$

for the control system introduced in section 3. The result depends crucially on (3), and the fact that the globally asymptotically stable state (0, 1) is disjoint from the dynamic constraint x ∈ Ω_c (see Equation (13)). That is, V_c < 1 is necessary for the following subsequent result to hold, and generally an optimal control will not exist if V_c = 1 or if the path constraint (13) is removed.

Our control system has the form

$\begin{array}{ll}\dot{x}=f(x)+u(t)g(x),& (95)\end{array}$

where x ∈ Ω, $u\in {U}$, and the vector fields f, g:Ω → $\text{R}$² are continuously differentiable. Note that the above vector field is affine (and thus continuous) in the control u. Fix the initial condition

$\begin{array}{ll}x(0)=x_0,& (96)\end{array}$

with x₀ ∈ Ω. Recall that all solutions of (95) and (96) approach the fixed point $\bar{x}:=(0,1)\in \Omega$. That is, for all $u\in {U}$,

$\begin{array}{ll}{x}_u(t)\stackrel{t\to \infty }{\to }\overline{x}.& (97)\end{array}$

Note that we explicitly denote the dependence of the trajectory on the control u, and the above point $\bar{x}$ is independent of the control u.

For any compact subset E of Ω such that $x_0\in E,\bar{x}\notin E$, we associate to each control (and hence to its corresponding trajectory) a time t_E(u) such that

$\begin{array}{ll}{t}_E(u)=max\{T|x_u(t)\in E\text{ for all }t\le T\}.& (98)\end{array}$

The above is well-defined (as a maximum) for each control u, since by assumption x₀ ∈ E and each trajectory asymptotically approaches $\bar{x}\notin E$, x_u is continuous, and E is compact.

THEOREM 19. Define

$\begin{array}{ll}{T}_{*}={\displaystyle \underset{u\in {U}}{sup}}{t}_E(u).& (99)\end{array}$

With the above construction, T_* is finite.

Proof. Consider the sets K, V ⊂ $\text{R}$², with V being an open neighborhood of the steady state $\bar{x}=(0,1)$ and K a compact set in $\text{R}$² such that

$\begin{array}{c}(0,1)\in V⊊K\text{ and }K\cap \{(x_1,x_2)\in {\text{R}}^2:x_1,x_2\ge 0\\ \text{and }0\le x_1+x_2\le V_c\}=\varnothing .\end{array}$

By contradiction, suppose that T_* is not finite, so we can find a sequence of controls ${\left\{v_k\right\}}_{k=1}^{\infty }$ in ${U}$ satisfying

$\begin{array}{ll}{d}_{\infty }\left(x(t,v_k),K\right)\ge ϵ\text{ for all }t\le t_k,\text{with }{t}_k\to \infty .& (100)\end{array}$

where d_∞ denotes the supremum metric and, for each k ∈ $\text{N}$, x(t, v_k) is the solution of the IVP:

$\begin{array}{ll}\begin{array}{l}\dot{x}=f(x)+v_k(t)g(x),\\ x(0)=x_0,\end{array}& (101)\end{array}$

Our aim is to find a control $u\in {U}$ such that x(t, u), solution of system (101), does not enter K for any t > 0. Recall that by the Banach-Alaoglu theorem, the ball

$\begin{array}{ll}B(L^{\infty }([0,\infty )))=\{u\in L^{\infty }([0,\infty )):\Vert u{\Vert }_{\infty }\le M\}& (102)\end{array}$

is a compact set on the weak^* topology and metrizable. Thus, the sequence ${\left\{v_k\right\}}_{k=1}^{\infty }$ must have a weak^*−convergent subsequence ${\left\{u_j\right\}}_{j=1}^{\infty }$ which converges to a control u ∈ L^∞([0, ∞)). In other words, for every ψ ∈ L¹([0, ∞))

$\begin{array}{ll}{\displaystyle \underset{j\to \infty }{lim}}{\displaystyle {\int }_{[0,\infty )}}\psi u_{j}d\mu ={\displaystyle {\int }_{[0,\infty )}}\psi ud\mu ,& (103)\end{array}$

where μ is the usual Lebesgue measure. This means that the sequence ${\left\{u_j\right\}}_{j=1}^{\infty }$ converges to u with respect to the weak^* topology on L^∞([0, ∞)) as the dual of L¹([0, ∞)).

We next prove that ${lim}_{j\to \infty }||x(t,u)-x(t,u_j){||}_{\infty }=0$ for all t ∈ [t_k−1, t_k] and all k ∈ $\text{N}$. In order to do so define

$x_{k-1}=x_0+{\displaystyle {\int }_0^{t_{k-1}}}\left[f(x(s))+u(s)g(x(s))\right]ds$

for any t_k−1 ∈ [0, ∞), where x solves the IVP

$\begin{array}{ll}\begin{array}{l}\dot{x}=f(x)+v_k(t)g(x),\\ x(t_{k-1})=x_{k-1}.\end{array}& (104)\end{array}$

Notice that the fact that the equilibrium (0, 1) is globally asymptotically stable on $\left\{(x_1,x_2)\in {\text{R}}^2:x_1,x_2\ge 0\mathrm{\text{and}}0<x_1+x_2\le V_c\right\}$ implies that x_k−1 is well-defined for any k ∈ $\text{N}$.

The integral form of (104) is given by

$\begin{array}{ll}F(t,x,v_k)=x_{k-1}+{\displaystyle {\int }_{t_{k-1}}^t}\left[f(x)+v_k(s)g(x)\right]ds.& (105)\end{array}$

With the help of the t_k's from (100) and assuming (without loss of generality) that t_k increases as k goes to infinity, we write the set [0, ∞) as the countable union of finite closed intervals:

$[0,\infty )={\displaystyle \bigcup _{k\in \text{N}}}\left[t_{k-1},t_k\right]\text{ where }{t}_0=0.$

Let w_{j, k} and w denote the functions u_j and u restricted to the interval [t_k−1, t_k], respectively. Thus, the sequence ${\left\{w_{j,k}\right\}}_{j=1}^{\infty }$ converges weakly* to w on [t_k−1, t_k]:

$\begin{array}{l}{\displaystyle \underset{j\to \infty }{lim}}||x(t,w)-x(t,w_{j,k}){||}_{\infty }\\ ={\displaystyle \underset{j\to \infty }{lim}}||F(t,x,w)-F(t,x,w_{j,k}){||}_{\infty }& (106)\end{array}$

$\begin{array}{ll}=\underset{j\to \infty }{\lim }{\Vert {\displaystyle {\int }_{t_{k-1}}^{t}w}(s)g(x)ds-{\displaystyle {\int }_{t_{k-1}}^t{w}_{j,k}}(s)g(x)ds\Vert }_{\infty }& (107)\end{array}$

$\begin{array}{ll}\begin{array}{l}=\underset{j\to \infty }{\lim }{\Vert {\displaystyle {\int }_{t_{k-1}}^t[w_{j,k}(s)-w(s)]g(x)ds}\Vert }_{\infty }\hfill \end{array}& (108)\end{array}$

$\begin{array}{ll}=0\text{ for all }t\in \left[t_{k-1},t_k\right]\text{.}& (109)\end{array}$

Since this result is independent of k, this implies that

$\begin{array}{l}{d}_{\infty }(x(t,u),K)=\underset{j\to \infty }{\lim }{d}_{\infty }(x(t,u_j),K)\ge ϵ\text{ for all}\\ t\in [t_{k-1},t_k],\text{independently of k}\in \text{N}.& (110)\end{array}$

The corresponding trajectory x(t, u) thus never enters K, contradicting the the global stability of $\bar{x}$. Hence, T_* must be finite, as desired.□

For the system and control problem defined in sections 2 and 3, the above theorem implies that ${sup}_{u\in {U}}{t}_c(u)$ is finite by taking E = Ω_c.

7.2.2. Supremum as a Maximum

Here we provide a general proof for the existence of optimal controls for systems of the form (95), assuming the set of maximal times is bounded above, which we have proven for our system in section 7.2.1. For convenience, we make the proof as self-contained as possible (one well-known result of Filippov will be cited), and state the results in generality which we later apply to the model of induced resistance. Arguments are adapted primarily from the textbook of Bressan and Piccoli (2007).

Consider again general control systems as in section 7.2.1. Solutions (or trajectories) of (95) will be defined as absolutely continuous functions for which a control $u\in {U}$ exists such that (x(t), u(t)) satisfy (95) a.e., in their (common) domain [a, b].

It is easier and classical to formulate existence with respect to differential inclusions. That is, define the multi-function

$\begin{array}{ll}F(x)=\{f(x)+\omega g(x)|\omega \in U\}.& (111)\end{array}$

Thus, the control system (95) is clearly related to the inclusion

$\begin{array}{ll}\dot{x}\in F(x).& (112)\end{array}$

The following theorem (see Filippov, 1967 for a proof) makes this relationship precise.

THEOREM 20. An absolutely continuous function x:[a, b] ↦ $\text{R}$² is a solution of (95) if and only if it satisfies (112) almost everywhere.

We first prove a lemma demonstrating that the set of trajectories is closed with respect to the sup-norm ||·||_∞ if all the sets of velocities F(x) are convex.

LEMMA 21. Let x_k be a sequence of solutions of (95) converging to x uniformly on [0, T]. If the graph of (t, x(t)) is entirely contained in Ω, and all the sets F(x) are convex, then x is also a solution of (95).

Proof. By the assumptions on f, g, the sets F(x) are uniformly bounded as (t, x) range in a compact domain, so that x_k are uniformly Lipschitz, and hence x is Lipschitz as the uniform limit. Thus x is differentiable a.e., and by Theorem 20, it is enough to show that

$\begin{array}{ll}\dot{x}(t)\in F(x(t))& (113)\end{array}$

for all t such that the derivative exists.

Assume not, i.e., that the derivative exists at some τ, but ẋ(τ)∉F(x(τ)). Since F(x(τ)) is compact and convex, and ẋ(τ) is closed, the Hyperplane Separation Theorem implies that there exists a hyperplane separating F(x(τ)) and ẋ(τ). That is, there exists an ϵ > 0 and a (without loss of generality) unit-vector p ∈ $\text{R}$² such that

$\begin{array}{ll}\langle p,y\rangle \le \langle p,\dot{x}(\tau )\rangle -3ϵ,& (114)\end{array}$

for all y ∈ F(x(τ)). By continuity, there exists δ > 0 such that for |x′ − x(τ)| ≤ δ

$\begin{array}{ll}\langle p,y\rangle \le \langle p,\dot{x}(\tau )\rangle -2ϵ,& (115)\end{array}$

for all y ∈ F(x′). Since x is differentiable at τ, we can choose τ′ > τ such that

$\begin{array}{ll}\begin{array}{l}{\displaystyle |\frac{x({\tau }^{\prime })-x(\tau )}{{\tau }^{\prime }-\tau }-\dot{x}(\tau )|<ϵ,}\\ |x(t)-x(\tau )|<\delta ,\end{array}& (116)\end{array}$

for all t ∈ [τ, τ′]. Equation (116) and uniform convergence then implies that, as p is a unit vector,

$\begin{array}{ll}\langle p,\frac{x_k({\tau }^{\prime })-x_k(\tau )}{{\tau }^{\prime }-\tau }\rangle \stackrel{k\to \infty }{\to }\langle p,\frac{x({\tau }^{\prime })-x(\tau )}{{\tau }^{\prime }-\tau }\rangle \ge \langle p,\dot{x}(\tau )\rangle -ϵ.& (117)\end{array}$

On the other hand, since ẋ(t) ∈ F(x′) for t ∈ [τ, τ′], Equation (115) implies that for k sufficiently large,

$\begin{array}{ll}\langle p,\frac{x_k({\tau }^{\prime })-x_k(\tau )}{{\tau }^{\prime }-\tau }\rangle =\frac{1}{{\tau }^{\prime }-\tau }{\displaystyle {\int }_{\tau }^{{\tau }^{\prime }}}\langle p,\dot{x}(t)\rangle \text{d}t\le \langle p,\dot{x}(\tau )\rangle -2ϵ.& (118)\end{array}$

Clearly, (117) and (118) contradict one another, so that (113) must be true, as desired.□

We now restate the optimal control problem associated to (95). Let S denotes the set of admissible terminal conditions, S ⊂ $\text{R}$ × $\text{R}$², and ϕ:$\text{R}$ × $\text{R}$² ↦ $\text{R}$ a cost function. We would like to maximize ϕ(T, x(T)) over admissible controls with initial and terminal constraints:

$\begin{array}{ll}\begin{array}{c}{\displaystyle \underset{u\in {U},T\ge 0}{max}}\varphi (T,x(T,u)),\\ x(0)=x_0,(T,x(T))\in S.\end{array}& (119)\end{array}$

We now state sufficient conditions for such an optimal control to exist.

THEOREM 22. Consider the control system (95) and corresponding optimal control problem (119). Assume the following:

1. The objective ϕ is continuous.

2. The sets of velocities F(x) are convex.

3. The trajectories x remain uniformly bounded.

4. The target set S is closed.

5. A trajectory satisfying the constraints in (119) exists.

6. S is contained in some strip [0, T] × $\text{R}$², i.e. the set of final times (for free-endpoint problems) can be uniformly bounded.

If the above items are all satisfied, an optimal control exists.

Proof. By assumption, there is at least one admissible trajectory reaching the target set S. Thus, we can construct a sequence of controls u_k:[0, T_k] ↦ U whose corresponding trajectories x_k satisfy

$\begin{array}{l}{x}_k(0)=x_0,\\ (T_k,x_k(T_k))\in S,\\ \varphi (T_k,x(T_k))\stackrel{k\to \infty }{\to }\underset{u\in {U},\overline{T}\ge 0}{\sup }\varphi (\overline{T},x(\overline{T},u)).& (120)\end{array}$

Since S ⊂ [0, T] × $\text{R}$ⁿ, we know that T_k ≤ T for all k. Each function x_k can then be extended to the entire interval [0, T] by setting x_k(t) = x_k(T_k) for t ∈ [T_k, T].

The sequence x_k is uniformly Lipschitz continuous, as f is uniformly bounded on bounded sets. This then implies equicontinuity of ${\left\{x_k\right\}}_{k=1}^{\infty }$. By the Arzela-Ascoli Theorem, there exists a subsequence x_nk such that T_nk → T_*, T_* ≤ T, and x_nk → x_* uniformly on [0, T_*].

Lemma 21 implies that x_* is admissible, so that there exists a control u_*:[0, T_*] ↦ U such that

$\begin{array}{ll}{\dot{x}}_{*}(t)=f(t,x_{*}(t),u_{*}(t))& (121)\end{array}$

for almost all t ∈ [0, T_*]. Equations (120) imply that

$\begin{array}{ll}\begin{array}{l}{x}_{*}(0)=x_0\\ (T_{*},x_{*}(T_{*}))={\displaystyle \underset{n_k\to \infty }{lim}}\varphi (T_{n_k},x_{n_k}(T_{n_k}))\in S.\end{array}& (122)\end{array}$

Note that the second of (122) relies on S being closed. Continuity of ϕ and (120) implies that

$\begin{array}{ll}\varphi (T_{*},x_{*}(T_{*}))={\displaystyle \underset{n_k\to \infty }{lim}}\varphi (T_{n_k},x_{n_k}(T_{n_k}))={\displaystyle \underset{u\in {U},T_{*}\ge 0}{sup}}\varphi (T_{*},x(T_{*},u)).& (123)\end{array}$

Thus, u_* is optimal, as desired.□

For the model of drug-induced resistance, the control set U is the compact set U = [0, M], and for such control-affine systems, convexity of F(x) is implied by the convexity of U. Existence of a trajectory satisfying the constraints is clear; for example, take u(t) ≡ 0. Our objective is to maximize the time to not escape the set N. Note that N is a closed subset of $\text{R}$², and that

$\begin{array}{ll}\varphi (\bar{T},x(\bar{T},u))=\bar{T}.& (124)\end{array}$

is continuous. Lastly, we have seen that all solutions remain in the closure $\bar{\Omega }$, so that |x(t)| ≤ 1 for all $u\in {U}$ and hence solutions are uniformly bounded. Existence is then reduced to Item 6 in the previous theorem. Since the supremum of time t was shown to be finite, Theorem 22 together with Theorem 19 imply that the optimal control for the problem presented in section 3 exists.

7.3. Further Numerical Experiments

In this subsection, we present further numerical experiments (see section 6). Specifically, we study how the values of the relative resistant growth rate and critical volume influence the control structure. We also consider a regularized objective, which suggests that our numerical methods are converging to (at least local) solutions of the optimal control problem.

We first investigate the control structure and treatment outcome as a function of d for a fixed α; these results are presented in Figures 13, 14. Here α = 0.005 is fixed and d is varied on the interval [0.001, 0.1]. Figure 13 presents three of these controls; although none of the controls is of the form YXY, the figure suggests that there may exist a d^* ∈ (0.02062, 0.0207959) where the solution trajectory may intersect the boundary line N only at one point and subsequently switches into a Y arc, thus providing the existence of a YXY control. Figure 14 suggests that increasing d for a fixed α increases the overall effectiveness of the treatment for all values of α, and that decreasing the induction rate α allows for longer tumor control. However, for small values of d, increasing α may provide a better treatment outcome (see, for example, the intersection of the yellow and purple curves in Figure 14).

FIGURE 13

Figure 13. Computed optimal controls for α = 0.005 and (A) d = 0.0206, (B) d = 0.020624489795918, and (C) d = 0.207959. Note that the control in (A) takes the form Y, while that in (B,C) is of the form YXu_p.

FIGURE 14

Figure 14. Variation in t_c as a function of d. (A) t_c response for varying d values. Note that treatment efficacy generally increases with increasing d. (B) α = 0.1.

We also investigated how the shape of the optimal control changes for different values of the resistant growth fraction (p_r) and/or the critical tumor volume (V_c). We run several simulations for V_c ∈ {0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9} and p_r ∈ {0.2, 0.3, 0.5, 0.7, 0.85, 0.9, 0.95, 0.98, 0.99}. We found that when the reproduction rate of resistant cells is close to the reproduction rate of sensitive cells (p_r near 1), the best strategy is to not give any drug at the beginning of treatment. This is perhaps to prolong the appearance of fast-growing resistance cells which cannot be eliminated with treatment. A representative set of the controls for these simulations are shown in Figure 15.

FIGURE 15

Figure 15. Optimal control structures for different V_c and p_r values. The blue curve is the computed optimal control, while the red curve is the feedback control along on the boundary of N, which may or may not be optimal or even feasible.

We further simulated the following parameter sets: V_c ∈ {0.3, 0.6}, d ∈ {0.01, 0.05, 0.1, 0.5, 0.75, 1} and p_r ∈ {0.2, 0.3, 0.5, 0.7, 0.85, 0.9, 0.95, 0.98, 0.99}. Figure 16 shows some of the controls for these simulations for the case when V_c = 0.6, while Figure 17 shows some of the controls for the case V_c = 0.3. In both figures, we observe that independently of the value of resistant growth rate p_r, if the chemotherapeutic drug has a low effectiveness (d small) then the best strategy is to give the maximum possible drug dosage during treatment. However, when d increases past d = 0.1, the control structure changes qualitatively. When V_c = 0.6 and the resistant reproduction rate is close to the reproduction rate of sensitive cells, the best strategy is to start with no drug treatment while for case V_c = 0.3 (independently of the value of p_r) the best strategy is to give the maximum drug dosage from the start.

FIGURE 16

Figure 16. Optimal control structures for V_c = 0.6, and different d and p_r values. The blue curve is the computed optimal control, while the red curve is the feedback control along on the boundary of N, which may or may not be optimal or even feasible.

FIGURE 17

Figure 17. Optimal control structures for V_c = 0.3, and different d and p_r values. The blue curve is the computed optimal control, while the red curve is the feedback control along on the boundary of N, which may or may not be optimal or even feasible.

Before ending this section, we would like to mention that to verify the performance of the numerical software, we approached the original problem by a sequence of regularized problems, which is done by adding a quadratic term to the Lagrangian. More precisely, we considered the perturbed performance index:

$\begin{array}{ll}{J}_{\eta }[u]=-{\displaystyle {\int }_0^{t_c}\left[1-\frac{(1-\eta )}{2}{u}^2(t)\right]}dt\text{ for }\eta \in [\text{0},\text{1}].& (125)\end{array}$

Notice that Equation (125) represents a family of performance indexes parameterized by η. The original performance index corresponds to η = 1. Furthermore, for η ≠ 1 the optimal control problem is regular and solvers such as GPOPS-II (used here) or SNOPT should provide accurate solutions. Thus, to test the accuracy to the case η = 1, we investigated the corresponding control structure in the limit η → 1. An example of different controls, for η values 0, 0.5, 0.7, 0.9, 0.95, 0.999, 0.99999, and 1, are shown on Figure 18. For each case we obtained different relative errors: the largest relative error of 4.0338 × 10⁻⁷ occurs for η = 0, with the remaining values of η having smaller relative errors. From the values η = 0.95, η = 0.999 and η = 0.99999 in Figure 18 we can see that as η → 1 the computed control approaches the solution to the original problem (case η = 1).

FIGURE 18

Figure 18. Different perturbed controls for α = 0.005 and d = 0.05. Here, from (A–H), the value of η is 0, 0.5, 0.7, 0.9, 0.95, 0.999, 0.99999, and 1, respectively. The maximum relative error is of 4.0338 × 10⁻⁷ for figure η = 0, the remaining figures have a maximum relative error of 5.5727 × 10⁻⁷ or smaller.

8. Conclusions

In this work, we have provided a rigorous analysis of the optimal control problem introduced in Greene et al. (2018a). That is, we have formally applied optimal control theory techniques to understand treatment strategies related to a model of induced drug resistance in cancer chemotherapy introduced in Greene et al. (2019). Although the model is relatively simple, it has recently been found to be highly successful in matching experimental data (Gevertz et al., 2019; Johnson et al., 2020), which we believe justifies the careful analysis presented here. An optimal control problem is then presented which maximizes a specific treatments therapy window. A formal analysis of the optimal control structure is performed utilizing the Pontryagin Maximum Principle and differential-geometric techniques. Optimal treatment strategies are realized as a combination of bang-bang and path-constrained arcs, and singular controls are proved to be sub-optimal. Numerical results are presented which verify our theoretical results, and demonstrate interesting and non-intuitive treatment strategies. We have also shown that a drug's level of resistance induction is identifiable, thus allowing for the possibility of designing therapies based on individual patient-drug interactions (see section 7.1).

Under the assumption that sensitive cells have a higher growth rate than resistant cells, our results (section 6) indicate that when using a chemotherapeutic drug with low cytotoxicity, the time at which the tumor volume exceeds its critical value t_c would be larger when the transition rate of the drug is high (see for example Table 2, on cases d = 0.001 and d = 0.01, as α has larger values the end time t_c becomes larger). The situation is reversed when we consider larger values of drug effectiveness because in this case it would take more time for the tumor to grow to its critical volume whenever the drug effectiveness is large enough. Also, our simulations indicate that it is optimal to apply the maximal dosage M subsequent to sliding along the boundary V = V_c (e.g., Figure 9), prior to treatment failure.

Clearly, further analysis is required in order to understand this phenomenon, and its implications for clinical scenarios. Although our model considers only an idealized scenario where resistance is unavoidable, we see that induced resistance dramatically alters therapy outcome, which underscores the importance of understanding its role in both cancer dynamics and designing chemotherapy regimes.

Other questions remain open for future work:

♢ Several studies indicate that drug-tolerance is a phenotypic property that appears transiently under the presence of the drug (Goldman et al., 2015). A next step to this research is to incorporate a reverse transition rate (from resistant to sensitive cells) that represents this phenotype-switching (see Figure 19).

♢ For controls where the trajectory remains on the boundary V = V_c (u_p), the feedback control is optimal during a time interval [t₁, t₂] with 0 ≤ t₁ < t₂ < t_c. It remains to understand the point of entry [x₁(t₁), x₂(t₁)] and exit [x₁(t₂), x₂(t₂)] (Figure 20A). What is the significance of the times t₁ and t₂ with respect to parameter values?

♢ Do there exist conditions, once the trajectory reaches V_c, under which the optimal trajectory no longer slides? Is it possible that at the time t_* the point [x₁(t_*), x₂(t_*)] is a contact point (Figure 20B)? Some numerical results suggest that such a contact point may exist and give rise to a YXY control structure (Figure 13).

♢ We have shown that an optimal control can switch at most once in each of the regions ${\Omega }_c^{+}$ and ${\Omega }_c^{-}$. Numerically we did not observe any bang-bang controls of the form YXY, although its existence was strongly suggested. The existence of a bang-bang junction in ${\Omega }_c^{-}$ is therefore of interest.

♢ For all examples plotted in Figure 11 with d ≥ 0.1, the entry time occurs approximately at the same value t₁ = 20.03. Is this a coincidence? We would like to understand the dependence of the entry time t₁ and on parameters α, d, p_r, M, and/or ϵ.

♢ We would like to extend models to include multiple, possibly non-cross resistant, cytotoxic agents. Indeed, clinical practice generally includes multiple agents applied concurrently and sequentially, and we plan on investigating strategies when different types of drugs may be applied. For example, what control strategies arise when a targeted therapy exists which targets the resistant sub-population? What order should the agents be applied, and for how long? Are intermediate doses now optimal? Mathematically, all of these questions may be studied, and the results may be clinically relevant.

FIGURE 19

Figure 19. Visualization of model that includes a reverse phenotype transition from resistant to sensitive. x₁ denotes the sensitive cancerous cell population, y_i the drug-induced resistant cancerous cell population, and y_s the non-drug-induced resistant cell population.

FIGURE 20

Figure 20. (A) Example of an arc with feedback control with entry point [x₁(t₁), x₂(t₁)] an exit point [x₁(t₂), x₂(t₂)] the exit point (B) Example of an arc that does not slides but reaches the boundary V = V_c at the contact point (x₁(t_*), x₂(t_*)).

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This research was supported in part by NSF grants 1716623 and 1849588.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dr. Anil Rao for technical suggestions regarding the optimization formulation and the use of GPOPS-II. This manuscript has been released as a pre-print at bioRxiv (Greene et al., 2018b).

References

Anguelova, M. (2004). Nonlinear Observability and Identifiability: General Theory and a Case Study of a Kinetic Model for S. cerevisiae. Chalmers University of Technology.

Google Scholar

Bressan, A., and Piccoli, B. (2007). Introduction to mathematical control theory. AIMS Ser. Appl. Math. Philadelphia

Google Scholar

Brimacombe, K. R., Hall, M. D., Auld, D. S., Inglese, J., Austin, C. P., Gottesman, M. M., et al. (2009). A dual-fluorescence high-throughput cell line system for probing multidrug resistance. Assay Drug Dev. Technol. 7, 233–249. doi: 10.1089/adt.2008.165

PubMed Abstract | CrossRef Full Text | Google Scholar

Doherty, M., Smigiel, J., Junk, D., and Jackson, M. (2016). Cancer stem cell plasticity drives therapeutic resistance. Cancers 8:8. doi: 10.3390/cancers8010008

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenberg, M. C., and Jain, H. V. (2017). A confidence building exercise in data and identifiability: Modeling cancer chemotherapy as a case study. J. Theor. Biol. 431, 63–78.

PubMed Abstract | Google Scholar

Filippov, A. F. (1967). Classical solutions of differential equations with multi-valued right-hand side. SIAM J. Control. 5, 609–621.

Google Scholar

Gatenby, R. A., Silva, A. S., Gillies, R. J., and Frieden, B. R. (2009). Adaptive therapy. Cancer Res. 69, 4894–4903. doi: 10.1158/0008-5472.CAN-08-3658

CrossRef Full Text | Google Scholar

Gevertz, J. L., Greene, J. M., and Sontag, E. D. (2019). Validation of a mathematical model of cancer incorporating spontaneous and induced evolution to drug resistance. bioRxiv. doi: 10.1101/2019.12.27.889444

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, A., Majumder, B., Dhawan, A., Ravi, S., Goldman, D., Kohandel, M., et al. (2015). Temporally sequenced anticancer drugs overcome adaptive resistance by targeting a vulnerable chemotherapy-induced phenotypic transition. Nat. Commun. 6:6139. doi: 10.1038/ncomms7139

PubMed Abstract | CrossRef Full Text | Google Scholar

Gottesman, M. (2002). Mechanisms of cancer drug resistance. Annu. Rev. Med. 53, 615–627. doi: 10.1146/annurev.med.53.082901.103929

PubMed Abstract | CrossRef Full Text | Google Scholar

Greene, J., Sanchez-Tapia, C., and Sontag, E. (2018b). Mathematical details on a cancer resistance model. bioRxiv [preprint]. doi: 10.1101/475533

CrossRef Full Text | Google Scholar

Greene, J., Sanchez-Tapia, C., and Sontag, E. D. (2018a). Control structures of drug resistance in cancer chemotherapy. Proc. IEEE Conf. Decis. Control. doi: 10.1109/CDC.2018.8618653

CrossRef Full Text | Google Scholar

Greene, J. M., Gevertz, J. L., and Sontag, E. D. (2019). Mathematical approach to differentiate spontaneous and induced evolution to drug resistance during cancer treatment. JCO Clin. Cancer Inform. 3, 1–20. doi: 10.1200/CCI.18.00087

PubMed Abstract | CrossRef Full Text | Google Scholar

Hermann, R., and Krener, A. (1977). Nonlinear controllability and observability. IEEE Trans. Automatic Control 22, 728–740. doi: 10.1109/TAC.1977.1101601

CrossRef Full Text | Google Scholar

Johnson, K. E., Howard, G. R., Morgan, D., Brenner, E., Gardner, A. L., Durrett, R. E., et al. (2020). Integrating multimodal data sets into a mathematical framework to describe and predict therapeutic resistance in cancer. bioRxiv [preprint]. doi: 10.1101/2020.02.11.943738

CrossRef Full Text | Google Scholar

Ledzewicz, U., and Schättler, H. (2012). Geometric Optimal Control. Theory, Methods and Examples, 1st Edn. New York, New York: Springer. doi: 10.1007/978-1-4614-3834-2

CrossRef Full Text | Google Scholar

Lee, W.-P. (1993). The role of reduced growth rate in the development of drug resistance of hob1 lymphoma cells to vincristine. Cancer Lett. 73, 105–111. doi: 10.1016/0304-3835(93)90251-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Loeb, L. A., Springgate, C. F., and Battula, N. (1974). Errors in DNA replication as a basis of malignant changes. Cancer Res. 34, 2311–2321.

PubMed Abstract | Google Scholar

Meshkat, N., and Seth, S. (2014). Identifiable reparametrizations of linear compartment models. J. Symbolic Comput. 63, 46–67.

Google Scholar

Patterson, M. A., and Rao, A. V. (2014). GPOPS-II: A matlab software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming. ACM Trans. Math. Softw. 41:1. doi: 10.1145/2558904

CrossRef Full Text | Google Scholar

Pisco, A. O., Brock, A., Zhou, J., Moor, A., Mojtahedi, M., Jackson, D., et al. (2013). Non-darwinian dynamics in therapy-induced cancer drug resistance. Nat. Commun. 4:2467. doi: 10.1038/ncomms3467

PubMed Abstract | CrossRef Full Text | Google Scholar

Pontryagin, L. S. (1987). Mathematical Theory of Optimal Processes. New York, NY; London, UK; Paris, Montreux, Tokyo: Gordon and Breach Science Publishers.

Google Scholar

Schättler, H., and Ledzewicz, U. (2015). Optimal Control for Mathematical Models of Cancer Therapies. New York, NY: Springer. doi: 10.1007/978-1-4939-2972-6

CrossRef Full Text | Google Scholar

Shackney, S. E., McCormack, G. W., and Cuchural, G. J. (1978). Growth rate patterns of solid tumors and their relation to responsiveness to therapy: an analytical review. Ann. Intern. Med. 89, 107-121. doi: 10.7326/0003-4819-89-1-107

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaffer, S. M., Dunagin, M. C., Torborg, S. R., Torre, E. A., Emert, B., Krepler, C., et al. (2017). Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature 546:431. doi: 10.1038/nature22794

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, S., Lee, D., Li, B., and Quinlan, M. E. A. (2010). A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell 141, 69–80. doi: 10.1016/j.cell.2010.02.027

CrossRef Full Text | Google Scholar

Silva, A., Silva, M. C., Sudalagunta, P., Distler, A., Jacobson, T., Collins, A., et al. (2017). An ex vivo platform for the prediction of clinical response in multiple myeloma. Cancer Res. 77, 3336–3351.

PubMed Abstract | Google Scholar

Sontag, E. D. (1979). On the observability of polynomial systems, I: Finite-time problems. SIAM J. Control Optimization. 17, 139–151. doi: 10.1137/0317011

CrossRef Full Text | Google Scholar

Sontag, E. D. (2017). Dynamic compensation, parameter identifiability, and equivariances. PLoS Comput. Biol. 13:e1005447.

PubMed Abstract | Google Scholar

Sontag, E. D., and Wang, Y. (1991). “I/O equations for nonlinear systems and observation spaces,” in Decision and Control, 1991., Proceedings of the 30th IEEE Conference on (IEEE), 720–725.

Google Scholar

Sussmann, H. (1982). “Time-optimal control in the plane,” in Feedback Control of Linear and Nonlinear Systems, eds D. Hinrichsen and A. Isidori (Berlin, Heidelberg: Springer), 244–260. doi: 10.1007/BFb0006833

CrossRef Full Text | Google Scholar

Sussmann, H. (1987a). Regular synthesis for time-optimal control of single-input real analytic systems in the plane. SIAM J. Control Optim. 25, 1145–1162. doi: 10.1137/0325062

CrossRef Full Text | Google Scholar

Sussmann, H. (1987b). The structure of time-optimal trajectories for single-input systems in the plane: The C^∞ nonsingular case. SIAM J. Control Optim. 25, 433–465. doi: 10.1137/0325025

CrossRef Full Text | Google Scholar

Sussmann, H. (1987c). The structure of time-optimal trajectories for single-input systems in the plane: the general real analytic case. SIAM J. Control Optim. 25, 868–904. doi: 10.1137/0325048

CrossRef Full Text | Google Scholar

Traina, T. A., and Norton, L. (2011). “Log-kill hypothesis,” in Encyclopedia of Cancer, ed M. Schwab (Berlin, Heidelberg: Springer), 2074–2075. doi: 10.1007/978-3-642-16483-5_3409

CrossRef Full Text | Google Scholar

Villaverde, A. F., Barreiro, A., and Papachristodoulou, A. (2016). Structural identifiability of dynamic systems biology models. PLoS Comput. Biol. 12:e1005153. doi: 10.1371/journal.pcbi.1005153

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., and Sontag, E. D. (1989). On two definitions of observation spaces. Syst. Control Lett. 13, 279–289.

Google Scholar