# Mathematical Details on a Cancer Resistance Model

^{1}Department of Mathematics, Clarkson University, Potsdam, NY, United States^{2}Department of Mathematics and Center for Quantitative Biology, Rutgers University, Piscataway, NJ, United States^{3}Department of Electrical and Computer Engineering, Department of Bioengineering, Northeastern University, Boston, MA, United States^{4}Laboratory of Systems Pharmacology, Program in Therapeutic Science, Harvard Medical School, Boston, MA, United States

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*_{1}) and resistant (*x*_{2}). 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*_{1}, 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*_{1}, 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*_{1}, 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):

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*_{2} denotes the total resistant subpopulation.

The region

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:

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: $\stackrel{\u0304}{x}:=({\stackrel{\u0304}{x}}_{1},{\stackrel{\u0304}{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*_{1} + *x*_{2}) or resistant fraction [*x*_{2}/(*x*_{1} + *x*_{2})] 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

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*_{1} + *x*_{2} ≤ *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

Define *t*_{c} as the time at which the tumor increases above size *V*_{c} for the first time. To be precise,

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:

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*_{1}(*t*) + *x*_{2}(*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

where

and *x*(*t*) = (*x*_{1}(*t*), *x*_{2}(*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*_{1} + *x*_{2} = *V*_{c} in Ω, i.e.,

where

With this notation, the path-constraint

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

for all times *t* such that *t* ∈ [0, *t*_{c}]. The initial state

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*_{1} *be the first time such that x*_{*}(*t*) ∈ *N*. *Fix* δ > 0 *such that t*_{1} − δ > 0, *and*

*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*^{0}, *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 $\stackrel{\u0304}{z}$ such that $\stackrel{\u0304}{z}(\tau )=\xi $, τ > *t*_{1} − δ, concatenate $\stackrel{\u0304}{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*_{1} > *t*_{c}, contradicting the global optimality of *x*_{*}.

Recall that *t*_{1} 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*_{1}, *x*_{2}) such that *x*_{1} + *x*_{2} < *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*_{1} + *x*_{2} < *V*_{c}. Furthermore, the constraint (13) has generic order one. In other words,

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*:

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*_{1} 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**. Region in Ω_{c} where *L*_{Y}*V*(*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*

*and the point* ${x}^{*}=({x}_{1}^{*},{x}_{2}^{*})\in N$ *with coordinates*

*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*_{1} + *x*_{2} *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*_{Y}*V*(*x*) denotes the Lie derivative of *V*(*x*) with respect to *Y*. Thus, for *x* ∈ *N*,

Assuming $M>\frac{(1-{V}_{c})(1-{p}_{r})}{d}$, the sign of *L*_{Y}*V*(*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

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*_{1}. Thus, it is feasible for *x* ∈ *N* if ${x}_{1}\in ({x}_{1}^{*},{V}_{c})$.□

**Proposition 4**. *For x* = (*x*_{1}, *x*_{2}) ∈ Ω_{c} *with*

*the Lie derivative L_{Y}V(x) is positive*.

*Proof*. As in Proposition 2, we compute *L*_{Y}*V*(*x*) directly:

where the first inequality utilizes *V* ≤ *V*_{c}, and the second relies on (21)□

**Proposition 5**. *For*

*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

Note that we are assuming here that the maximal dosage *M* satisfies $M>\frac{1-\u03f5}{\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

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

The Maximum Principle then yields the following theorem:

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

1. (λ_{0}, λ(*t*)) ≠ 0 *for all t* ∈ [0, *t*_{c}].

2. λ(*t*) = (λ_{1}(*t*), λ_{2}(*t*)) *satisfies the second-order differential equation*

3. *u*_{*}(*t*) *minimizes H pointwise over the control set U*:

*Thus, the control u*_{*}(*t*) *must satisfy*

*where*

4. *The Hamiltonian H is identically zero along the extremal lift* (*x*_{*}(*t*), *u*_{*}(*t*), λ(*t*)):

*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*(λ_{0}, λ(*t*), *x*_{*}(*t*), *u*_{*}(*t*)) = −λ_{0}. By item (4) in Theorem 6, *H* ≡ 0, which implies that λ_{0} = 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} = 0), as well as characterize the switching function dynamics via the Lie bracket.

**Proposition 8**. *For all x*_{1} ∈ Ω, *x*_{1} > 0, *the vector fields f*(*x*) *and g*(*x*) *are linearly independent*.

*Proof*. Define *A*(*x*) = *A*(*x*_{1}, *x*_{2}) to be the matrix

The determinant of *A* can calculated as

where

As *x*_{1}(*t*) + *x*_{2}(*t*) ≤ 1 for all *t* ≥ 0, κ(*x*(*t*)) ≥ 0, and we see that det*A*(*x*) = 0 in Ω if and only if *x*_{1} = 0, completing the proof.□

The line *x*_{1} = 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*_{1}(*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} = 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} = 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} = 0 and Φ(*t*) = 〈λ(*t*), *g*(*x*_{*}(*t*))〉, Equation (22) reduces to

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 λ_{0} = 1. Indeed, λ(*t*) may be rescaled by λ_{0} > 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

We recall the Lie bracket as the first-order differential operator between two vector fields *X*_{1} and *X*_{2}:

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

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

where

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*_{1} + *x*_{2} < *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*]:

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*))〉:

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

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*_{1}(*t*), *x*_{2}(*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 Φ^{(n)}(*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*_{1}, *x*_{2}) plane where singular arcs are attainable, as demonstrated in the following proposition.

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

*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*_{1}, *x*_{2}) phase plane, we may calculate $\dot{\Phi}(t)$ from equation (43):

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*_{1} + *bx*_{2} − *c* = 0. Thus, define regions in the first quadrant as follows:

Recall that Ω_{c} is simply the region in Ω prior to treatment failure, i.e., 0 ≤ *V* ≤ *V*_{c}, *x*_{1}, *x*_{2} ≥ 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

**Figure 3**. Domain in (*x*_{1}, *x*_{2}) plane. **(A)** Region where γ changes sign. We see that inside the triangular region *x*_{1} + *x*_{2} ≤ 1 of the first quadrant, γ changes sign only along the line *ax*_{1} + *bx*_{2} − *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*_{1} + *bx*_{2} = *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*,

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*

*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*_{1}.

*Thus*, ${L}$ *contains a segment* $\stackrel{\u0304}{L}\subset {L}$ *which is a singular arc. Note that* $\stackrel{\u0304}{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 $\stackrel{\u0304}{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*_{1} + *bx*_{2} = *c*, taking the time derivative of this equation yields *aẋ*_{1} + *bẋ*_{2} = 0, and substituting the expressions (1) we compute *u*_{s} as

where *a, b, c* are given by (37–39) and *x*_{2} and *x*_{1} satisfy *ax*_{1} + *bx*_{2} = *c*. As discussed previously, *x*_{1}(*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*_{1} (i.e. as a *feedback law*) for ${x}_{1}\in (\stackrel{\u0304}{s},\frac{c}{a})$, for some $\stackrel{\u0304}{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}=\stackrel{\u0304}{s}$ where vector field *Y* is parallel to ${L}$. Thus, at this $\stackrel{\u0304}{s}$, the trajectory must leave the singular arc, and enter the region ${\Omega}_{c}^{-}$. As *ẋ*_{2} ≥ 0, trajectories must follow ${L}$ in the direction of decreasing *x*_{1} (see Figure 4). We summarize these results in the following theorem.

**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*_{2} = 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 $\stackrel{\u0304}{\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*_{1}, *x*_{2}) *plane as a segment of the line* ${L}$. *Along this singular arc, the control is given by Equation* (*52*), *where ax*_{1} + *bx*_{2} = *c*. *Therefore, in this case the necessary minimum conditions on u*_{*} *from* (*25*) *can be updated as follows*:

*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*_{2} is increasing, no segment $\stackrel{\u0304}{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 $\stackrel{\u0304}{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* $\stackrel{\u0304}{L}$.

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

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 $\stackrel{\u0304}{\tau}>\tau $ such that $\Phi (\stackrel{\u0304}{\tau})=0,\dot{\Phi}(\stackrel{\u0304}{\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 $\stackrel{\u0304}{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 $\stackrel{\u0304}{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 $\stackrel{\u0304}{L}$ by the definition of ${L}$.

**THEOREM 15**. *The line segment* $\stackrel{\u0304}{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 $\stackrel{\u0304}{L}$ at *x*, i.e. we find the unique solution to

Note that we can invert (50):

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

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

As *L*_{X}γ < 0 and *L*_{Y}γ > 0 on $\stackrel{\u0304}{\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

for $q\in \stackrel{\u0304}{L}$. Let ϕ ∈ ($\text{R}$^{2})* such that

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 $\stackrel{\u0304}{L}$ on some interval containing zero, and we can compute $\dot{\Phi}$ as [using (34)]

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 $\stackrel{\u0304}{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 $\stackrel{\u0304}{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 2*k* − 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

where ad_{Z} is the adjoint endomorphism for a fixed vector field *Z*:

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

along the arc. Note that $\frac{\partial H}{\partial u}=\Phi $, so that the above is simply the *u* coefficient of the 2*k*-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 \stackrel{\u0304}{L}$, 〈λ(*t*), [*g*, [*f, g*]](*x*(*t*))〉 > 0. *Thus, the Legendre-Clebsch condition is violated, and the singular arc* $\stackrel{\u0304}{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:

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

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

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

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

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

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

Theorem 16 then implies that the singular arc is suboptimal, i.e. that $\stackrel{\u0304}{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

In fact, we compute it explicitly:

Then, along any controlled trajectory (*x, u*) defined on [*t*_{0}, *t*_{1}], the integral of ω computes the time *t*_{1} − *t*_{0}:

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 \stackrel{\u0304}{L}$ and the corresponding unique *XY* trajectory connecting these points in ${\Omega}_{c}^{-}$ (note that uniqueness is guaranteed as long as *q*_{1} and *q*_{2} are sufficiently close).

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

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

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

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**. *XY* (solid) and *XYXY* (dashed) trajectories taking *q*_{1} to *q*_{2} 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*

*where* (*YX*)^{n}: = (*YX*)^{n−1}*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*_{p}*Yu*_{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 τ_{1} and τ_{2} 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* ∈ [τ_{1}, τ_{2}]. Thus, we can consider the alternate feasible trajectory which remains on *N* between the points (*S*(τ_{1}), *R*(τ_{1})) and (*S*(τ_{2}), *R*(τ_{2})); 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),

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

i.e., a longer amount of time is achieved by remaining on the boundary *N*. Hence the arc *u*_{p}*Yu*_{p} is sub-optimal if it is feasible, as desired.

**Figure 7**. Comparison of *u*_{p}*Yu*_{p} arc and an arc that remains on *N* (hence *u* ≡ *u*_{p}) between the points [*S*(τ_{1}), *R*(τ_{1})] and [*S*(τ_{2}), *R*(τ_{2})], 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 τ_{2} − τ_{1}, implying that *u*_{p}*Yu*_{p} is sub-optimal.□

Note that in Theorem 17, the switchings must occur across the singular arc $\stackrel{\u0304}{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*:

*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*_{Y}*V*(*x*) must satisfy

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

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**. Parameter values and initial conditions used throughout section 6, unless stated otherwise.

Theorem 17 characterizes the qualitative form of the optimal control:

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*_{p}*Y*, 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**. 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*_{1}) and resistant (*x*_{2}) temporal dynamics. **(B)** Control structure of form *YXu*_{p}*Y*. **(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**. Phase plane corresponding to Figure 8. Trajectory which optimal control is of the form *YXu*_{p}*Y* 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*_{1} > 0.1059 and will increase for *x*_{1} < 0.1059.

**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*_{1}), resistant (*x*_{2}), and volume (*x*_{1} + *x*_{2}) 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**. Optimal time *t*_{c} for each of the computed controls appearing in 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**. 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*_{1}, *t*_{2}] with 0 < *t*_{1} < *t*_{2} < *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*):

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:

Here $0<{x}_{1}^{0}<1$, so that (*x*_{1}(*t*), *x*_{2}(*t*)) ∈ Ω for all *t* ≥ 0. We note that *x*_{2}(0) = 0 is not restrictive, and similar results may derived under the more general assumption 0 ≤ *x*_{2}(0) < 1. The condition *x*_{2}(0) = 0 is utilized both for computational simplicity and since *x*_{2}(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:

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

and *x*(*t*) = (*x*_{1}(*t*), *x*_{2}(*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:

Note that *x*_{1}(*t*), *x*_{2}(*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,

*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

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*_{X}*H* of a *C*^{ω} function *H* by a *C*^{ω} (i.e. real-analytic) vector field *X*:

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*_{Y}*L*_{X}*H* = *L*_{Y}(*L*_{X}*H*), 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*),

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*_{0} appear due to evaluation at *t* = 0. The observation space is then defined as the span of the elements *Y*(*x*_{0}, *U*):

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

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

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*_{0} = (*S*_{0}, 0) we can recursively determine the parameters *p*:

Since *F*_{1} = *F*_{2}, 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

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

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

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

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

For any compact subset *E* of Ω such that ${x}_{0}\in E,\stackrel{\u0304}{x}\notin E$, we associate to each control (and hence to its corresponding trajectory) a time *t*_{E}(*u*) such that

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

**THEOREM 19**. *Define*

*With the above construction, T*_{*} *is finite*.

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

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

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

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

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*^{1}([0, ∞))

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*^{1}([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

for any *t*_{k−1} ∈ [0, ∞), where *x* solves the IVP

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\text{}\mathrm{\text{and}}\text{}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

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:

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}]:

Since this result is independent of *k*, this implies that

The corresponding trajectory *x*(*t, u*) thus never enters *K*, contradicting the the global stability of $\stackrel{\u0304}{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

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

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

**THEOREM 20**. *An absolutely continuous function x*:[*a, b*] ↦ $\text{R}$^{2} *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

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}$^{2} such that

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

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

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

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

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}$^{2}, and ϕ:$\text{R}$ × $\text{R}$^{2} ↦ $\text{R}$ a cost function. We would like to maximize ϕ(*T, x*(*T*)) over admissible controls with initial and terminal constraints:

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}$^{2}, *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

Since *S* ⊂ [0, *T*] × $\text{R}$^{n}, 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

for almost all *t* ∈ [0, *T*_{*}]. Equations (120) imply that

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

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}$^{2}, and that

is continuous. Lastly, we have seen that all solutions remain in the closure $\stackrel{\u0304}{\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**. 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**. 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**. 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**. 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**. 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:

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^{−7} 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**. 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^{−7} for figure η = 0, the remaining figures have a maximum relative error of 5.5727 × 10^{−7} 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*_{1}, *t*_{2}] with 0 ≤ *t*_{1} < *t*_{2} < *t*_{c}. It remains to understand the point of entry [*x*_{1}(*t*_{1}), *x*_{2}(*t*_{1})] and exit [*x*_{1}(*t*_{2}), *x*_{2}(*t*_{2})] (Figure 20A). What is the significance of the times *t*_{1} and *t*_{2} 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*_{1}(*t*_{*}), *x*_{2}(*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*_{1} = 20.03. Is this a coincidence? We would like to understand the dependence of the entry time *t*_{1} 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**. Visualization of model that includes a reverse phenotype transition from resistant to sensitive. *x*_{1} 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. (A)** Example of an arc with feedback control with entry point [*x*_{1}(*t*_{1}), *x*_{2}(*t*_{1})] an exit point [*x*_{1}(*t*_{2}), *x*_{2}(*t*_{2})] the exit point **(B)** Example of an arc that does not slides but reaches the boundary *V* = *V*_{c} at the contact point (*x*_{1}(*t*_{*}), *x*_{2}(*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.

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

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

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

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

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.

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

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

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.

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

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

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

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

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

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

Keywords: drug resistance, chemotherapy, phenotype, optimal control theory, singular controls

Citation: Greene JM, Sanchez-Tapia C and Sontag ED (2020) Mathematical Details on a Cancer Resistance Model. *Front. Bioeng. Biotechnol.* 8:501. doi: 10.3389/fbioe.2020.00501

Received: 30 January 2020; Accepted: 29 April 2020;

Published: 17 June 2020.

Edited by:

Julio R. Banga, Consejo Superior de Investigaciones Científicas (CSIC), SpainReviewed by:

Sebastian Sager, Otto von Guericke University Magdeburg, GermanyJames William Thomas Yates, AstraZeneca, United Kingdom

Copyright © 2020 Greene, Sanchez-Tapia and Sontag. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Eduardo D. Sontag, sontag@sontaglab.org