OPTIMAL CONTROL FOR A MATHEMATICAL MODEL FOR ANTI-ANGIOGENIC TREATMENT WITH MICHAELIS-MENTEN PHARMACODYNAMICS

. We analyze the structure of optimal protocols for a mathematical model of tumor anti-angiogenic treatment. The control represents the concentration of the agent and we consider the problem to administer an a priori given total amount of agents in order to achieve a minimum tumor vol- ume/maximum tumor reduction. In earlier work, this problem was studied with a log-kill type pharmacodynamic model for drug eﬀects which does not account for saturation of the drug concentration. Here we study the eﬀect of incorporating a Michaelis-Menten (MM) or E max -type pharmacodynamic model, the most commonly used model in the ﬁeld of pharmacometrics. We compare the formulations of both problems and the resulting solutions. The reformulated problem with E max pharmacodynamics is no longer linear in the control. This results in qualitative changes in the structure of optimal controls which, in line with an interpretation as concentrations, now are continuous while discontinuities exist if the log-kill model is used which is more in line with an interpretation of the control as dose rates. In spite of these qualitative diﬀerences, similarities in the structures of solutions can be observed. Both aspects are discussed theoretically and illustrated numerically.

1. Introduction. We are interested in the question how an a priori specified amount of drugs (which has been determined based on potential side effects) can be administered to achieve a 'best possible' outcome in drug treatments. Clearly, the expression 'best possible' is open to interpretation and in this paper, for sake of specificity, we consider the problem to maximize the tumor reduction achievable with the given amount of drugs. In the modeling of such problems, an essential role is played by the relations between the drug's dosage, u, its concentration in the bloodstream, c, and the actual effects the drug has. These relations, summarily called pharmacometrics, are studied under the terms pharmacokinetics (PK) and pharmacodynamics (PD) in pharmacology. Pharmacokinetic models (PK) describe the relations between the drug dosage and the drug's concentration in the blood stream ("what the body does to the drug") while pharmacodynamic (PD) models describe the effects that the drug has on the disease ("what the drug does to the body"). The underlying question of interest in our research is whether, and if so, to what extent, the mathematical equations used in the modeling determine the structure of optimal controls. In this paper, as a case study, we investigate the effects that the modeling of PD-a topic still far from understood for many drugs-has on optimal controls for a model for anti-angiogenic therapy.
As a general rule, the PK of a drug is well-understood. Low-dimensional, linear differential equations with a small number of compartments (typically ranging from 1 to 3) are sufficient to provide an adequate modeling [4]. When longer time horizons are considered, drug dosage and drug concentrations are commonly identified as PK acts on a much faster time-scale and typically the steady-state value of c is proportional to u. In such a case, the PK of drugs need not be taken into account. The PD of drugs, on the other hand, is a much more intricate and less understood subject. Pharmacodynamic models are generally described by functional relations of the form ϕ(c)x where x represents a compartment of cells on which the drugs are acting (e.g., tumor cells, vasculature, immune system, healthy cells, etc.) and the coefficient ϕ(c) models the effect of the drug at concentration c. Over a limited range of concentrations (not too low and particularly, not too high), linear models based on the log-kill hypothesis [16] of the form ϕ(c) = Gc with G a positive constant are adequate. However, drug effects saturate with increasing concentrations and also may require some minimum concentration levels to be effective at all. In the pharmaceutical industry, Michaelis-Menten type relations of the form ϕ(c) = E max c U C50+c with E max denoting the maximum effect the drug can have and U C 50 the concentration at which half of this effect is realized, called the E max model in pharmacology, are commonly used to model these saturating effects at higher concentrations while sigmoidal functions also model the ineffectiveness of drugs at low concentrations [9,13]. Clearly, these equations can be linearized over certain ranges, but, for example, the resulting slopes (i.e., the coefficients G to be used) will be very different if one linearizes around a point on the upper end of this curve or say the point U C 50 . This may significantly alter the structure of solutions and it seems worthwhile to compare the solutions to such structurally different PD models (log-kill, Michaelis-Menten, sigmoidal models) [6,15].
Here we focus on a comparison between the linear log-kill and the concave E max models. Drugs are typically administered in the upper range of the concentrationeffect curve and an E max model based on Michaelis-Menten terms is commonly used in approaches pursued by the pharmaceutical industry (e.g., see [5,10]). Thus it seemed important to focus on this model in our study. From the point of view of an optimal control formulation, the mathematical structures of these two problem formulations differ significantly: if a linear log-kill model is used, then the resulting Hamiltonian function H for the optimal control problem becomes linear in the control leading to bang-bang and singular controls as well as all kind of possible concatenations as candidates for optimality. Such controls typically have discontinuities and relate to an interpretation of the controls as dose rates. An E max -model, on the other hand, brings in an element of convexity that takes away the linear structure in the control and for such a problem formulation optimal controls are continuous in line with an interpretation of the controls as concentrations. Mathematically, the two problem formulations (with linear, respectively, E max models for PD) have fundamentally different properties as optimal control problems: the first one leads to what mathematically are considered degenerate situations since the matrix ∂ 2 H ∂u 2 of the second derivatives of the Hamiltonian H vanishes identically while in the second one, depending on the sign of the multiplier, this quantity will be positive or negative along optimal solutions. This dichotomy clearly raises the possibility of significant structural changes between the solutions of these optimal control problems as they indeed exist. On the other hand, if the concentrations along the optimal solutions do not vary excessively, one would expect the solutions to these two problem formulations to be close to each other. The question to what extent this indeed is the case is a non-trivial one mathematically with potential practical implications on the underlying models and it is the focus of our research.
In this paper, we analyze a case study for this question. The particular model we have chosen is the problem to administer a given amount of anti-angiogenic agents in order to minimize the tumor volume. The dynamics is based on the well-established model by Hahnfeldt et al. [3] and we have chosen this example simply since a complete solution for the model with a linear pharmacokinetic model is known from previous work of some of the authors [7,15]. Mathematically, while there are many similarities in the analysis, simply because of the similarities in the dynamics, there are also significant differences that are brought in by the nonlinear model in the controls. In fact, we do not yet have an equally complete solution for the problem with an E max pharmacodynamic model, but there are strong correlations for the cases for which we can establish the theoretical structure of solutions for both models and numerical solutions reveal strong similarities. Our emphasis here lies less on giving a complete solution for the E max model, but on comparing some of its results with the solutions for a linear log-kill model.
The paper is organized as follows: We briefly summarize the underlying model in section 2 also introducing well-posed initial data. Necessary conditions for optimality are then given in section 3 where we also show that optimal controls are continuous. Similarly as for the log-kill model, optimal controls alternate between intervals where the control is constant given by a maximum or no dose and intervals where the control is smooth and lies in the interior of the control set. Junctions between these segments are theoretically analyzed in section 4 and, for some biologically realistic initial conditions and parameters, the optimal concatenation structure is established in section 5. Numerical examples of solutions are given in section 6.

2.
Optimal control for anti-angiogenic treatment with Michaelis-Menten type pharmacodynamics. We consider a dynamical system for tumor development under angiogenic signaling based on the equations by Hahnfeldt, Panigrahy, Folkman and Hlatky [3]. In this model, the spatial aspects of the underlying consumption-diffusion process that stimulate and inhibit angiogenesis are incorporated into a nonspatial 2-compartment model with the primary tumor volume, p, and the carrying capacity of the vasculature, q, as its principal variables. The dynamics consists of two ODEs that describe the evolution of p and q. This model has been extensively analyzed in the literature and we only briefly describe its premises, but refer the reader to [3] or [15] for a detailed development of the mathematical model. Tumor growth is described using a Gompertzian model of the forṁ with ξ a constant tumor growth parameter. The carrying capacity and tumor volume are balanced for p = q while the tumor volume shrinks for inadequate vascular support (q < p) and increases if this support is plentiful (q > p). Different from traditional approaches, the carrying capacity is not considered constant, but becomes a state variable whose evolution is governed by a balance of stimulatory and inhibitory effects. In the model by Hahnfeldt et al. [3], these terms are taken in the formq = bp − dp with bp representing the stimulation of the vasculature by the tumor (which is taken proportional to the tumor volume) and dp 2 3 q modeling the corresponding inhibition (which is represented by an interaction term between the tumor surface and the vasculature). These functional forms reflect important aspects of angiogenesis, namely that tumor-derived inhibitors act more systemically, whereas tumor-derived stimulators act more locally to the individual secreting tumor site [3, pg. 4771]. Finally, µq is a standard death term. The terms ξ, b, d, µ and G in these equations are constant parameters with µ ≥ 0 and all others are positive.
The variable u represents the concentration of the anti-angiogenic agent and its effect is modeled by the linear log-kill term Guq with G representing the efficacy of the agent. In this paper, we replace this pharmacodynamic model with a Michaelis-Menten term G u U C50+u q which incorporates saturating effects (see Fig. 1). This is the model widely in use in the pharmacological industry since it more realistically models drug effects in the upper range of the dose-concentration curve, u ≥ U C 50 , the typical operating range for drug administration. Mathematically, one could normalize U C 50 by scaling the control u in terms of U C 50 . This, however, is not done in pharmacology in order to keep a scale of real doses and thus we also retain this redundant parameter here. Then the dynamics for the evolution of the vasculature takes the following form:q = bp − dp In steady-state, in this model it is assumed that dosage and concentration of the anti-angiogenic agent are proportional and thus we do not include the standard linear pharmacokinetic model, i.e., we identify drug dosage with drug concentrations. In this sense the integral T 0 u(t)dt, where T is the terminal time of the therapy, represents the total amount of drugs given. We are interested in how to administer this fixed amount of drugs in order to achieve the best possible tumor reductions. Also, we limit the maximum dosage that is administered, 0 ≤ u ≤ u max . Then, given an arbitrary control strategy u : [0, T ] → [0, u max ] (which for mathematical reasons we take as a Lebesgue measurable function since this choice guarantees the existence of a solution), the objective J = J(u) is a functional of this control strategy u(·) that gives us a measure for the corresponding outcome. Here we simply use the final tumor volume as this measure, i.e., J(u) = p(T ). We thus consider the following optimal control problem: Mathematically, it is of advantage to adjoin an additional variable y with the dynamicsẏ = u and initial condition y(0) = 0. This turns the isoperimetric constraint (4) into the terminal constraint y(T ) ≤ A. Introducing the vector z = (p, q, y) † and the vector fields f , g 1 and g 2 below, then the dynamics takes the forṁ Mathematically, there is a significant difference between this model and the one with a log-kill PD in that the dynamics is no longer linear in u. The vector field g 2 is constant and is used to count the total dosage of the drugs, i.e., this leads to a term which is linear in u regardless of what model is used for PD. The latter enters as a factor at g 1 , but the vector fields g 1 and g 2 will be the same regardless of the specific model for PD that is used.
In the model formulation, the initial conditions p 0 and q 0 and the amount A of available agents are fixed, but arbitrary. The aim is to determine the structure of optimal controls not just for one particular value, but to obtain classifications of all possible solutions depending on these values. This naturally leads to a variety of cases. The administration of anti-angiogenic drugs directly leads to a reduction of the carrying capacity q of the vasculature, but only indirectly affects the tumor volume p. Indeed, the tumor volume increases regardless of the control whenever p < q and it decreases whenever p > q. It is therefore possible that, for some initial conditions (p 0 , q 0 ), the amount A of inhibitors is too small to achieve a tumor reduction. In such a case, the mathematically optimal solution for problem [OC] is simply given by T = 0. All that is possible in such a case is to slow down tumor growth, but the tumor cannot be reduced. This generates exceptional cases, which are easily analyzed, but also lead to degenerate necessary conditions for optimality that we wish to avoid here. We therefore make the following definition: We say the initial data (p 0 , q 0 , A) are ill-posed for the optimal control problem [OC] if for no admissible control it is possible to reach a point (p, q) with p < p 0 . In such a case, the optimal solution for the problem [OC] is given by T = 0. The initial data (p 0 , q 0 , A) are said to be well-posed if an objective value better than p 0 is realizable. In this case, the final time T along the optimal control is positive.
We briefly recall some dynamical systems properties of equations (1) and (3) for constant controls. Here we ignore the constraint (4) on the total dosage and, since the controls are constant, we can simply rely on the results for the model with a linear log-kill PD form (e.g., see [11] or [15,Sect. 5 and v max = umax U C50+umax . We recall the following facts: for a constant control v, the system has a unique globally asymptotically stable equilibrium point at provided that b − µ > Gv. For b − µ ≤ Gv, the system no longer has an equilibrium point and all trajectories converge to the origin as t → ∞. In this case, theoretically, eradication of the tumor is possible with an unlimited supply of inhibitors. We henceforth make this somewhat simplifying assumption.
[A]: We assume that Lemma 2.2. Suppose assumption [A] is satisfied and u * is an optimal control for problem [OC] defined over the interval [0, T ] with T > 0 (i.e., for well-posed initial data). Then optimal controls end with a no dose segment over a final interval (τ, T ], the endpoint of the optimal trajectory lies on the diagonal, p * (T ) = q * (T ), and all inhibitors have been used up, The proof of this Lemma is identical with the one of Lemma 5.2.1 in [15]. The key point is that assumption [A] allows to lower the tumor volume if the state of the system is on the diagonal. For, since b−µ G < umax U C50+umax , along the control u ≡ u max we have at any point p = q on the diagonal thatq < 0. Hence trajectories corresponding to a full dose control always cross the diagonal D 0 = {(p, q) ∈ D : p = q} from the region D − = {(p, q) ∈ D : p < q} below the diagonal into the region D + = {(p, q) ∈ D : p > q} above the diagonal. Thus, if for some time τ we have that p(τ ) = q(τ ) and y(τ ) < A, then we can still use the control u ≡ u max until all inhibitors have been exhausted and lower the tumor volume. Thus optimal trajectories necessarily end with a no dose segment when the diagonal is reached and all inhibitors have been exhausted.
We also recall that ifp =q = b−µ d 3 2 denotes the equilibrium point for the uncontrolled system (i.e., v = u = 0), then the region is positively invariant for the control system, i.e., if (p 0 , q 0 ) ∈ D and u is an arbitrary admissible control defined on some interval [0, T ] ⊂ [0, ∞), then the solution (p(·), q(·)) to the corresponding dynamics with initial condition (p(0), q(0)) = (p 0 , q 0 ) exists for all t ∈ [0, T ] and lies in D, (p(t), q(t)) ∈ D. Indeed, as is shown in Lemma 5.3.1 in [15], the sub-region N = (p, q) ∈ D : bp > µ + dp 2 3 q is also positively invariant and as initial points outside of this region are medically unrealistic, in this paper we always assume that the initial condition lies in this region: We note that, in order to better visualize tumor reductions in our figures, we graph p as the vertical and q as the horizontal variable.
3. Necessary conditions for optimality for well-posed data. The fundamental first-order necessary conditions for optimality for problem [OC] are given by the Pontryagin maximum principle [12] (for some more recent references on optimal control, see [1,2,14]). With λ = (λ 1 , λ 2 , λ 3 ) : [0, T ] → (R 3 ) * a 3-dimensional covector (which we write as a row vector) the Hamiltonian function H for the optimal control problem is given by If u * is an optimal control for problem [OC] defined over an interval [0, T ] with corresponding trajectory (p * , q * ), then there exist a non-negative constant λ 0 and an absolutely continuous co-vector, λ : [0, T ] → (R 3 ) * , such that (1) not all λ i , i = 0, 1, 2, 3, are zero, (2) λ 1 and λ 2 satisfy the adjoint equationṡ with terminal conditions λ 1 (T ) = λ 0 , λ 2 (T ) = 0 and λ 3 is a non-negative constant, and (3) the optimal control u * minimizes the Hamiltonian function H along the multiplier λ and the optimal trajectory z * almost everywhere on [0, T ] with minimum value 0, i.e., Controlled trajectories (z, u) for which there exist multipliers λ such that these conditions are satisfied are called extremals and the triples (z, u, λ) including the multipliers are called extremal lifts (to the cotangent bundle). Extremals for which λ 0 = 0 are called abnormal while those with a positive multiplier λ 0 are called normal.
Lemma 3.1. If the initial data (p 0 , q 0 , A) are well-posed, then extremals are normal.
Proof. If λ 0 = 0, then it follows from the adjoint equation that λ 1 and λ 2 vanish identically. Hence the non-triviality of the multipliers gives us that λ 3 = 0 and thus u * (t) ≡ 0 by the fact that the Hamiltonian H vanishes identically along the optimal solution. Hence the initial data are ill-posed.
Proof. Since extremals are normal, it follows from the adjoint equation that λ 1 and λ 2 cannot vanish simultaneously. At the terminal time we have that λ 1 (T ) > 0, λ 2 (T ) = 0 andλ 2 (T ) = −ξλ 1 (T ) < 0 and thus both multipliers are positive for t < T and t close enough to T . Suppose there exists a time τ ≥ 0 such that λ 1 and λ 2 are positive on (τ, T ), but one of them vanishes at time τ . If λ 1 (τ ) = 0, then λ 2 (τ ) > 0 and, since the region N is positively invariant under the control system, Contradiction. Thus we must have λ 1 (τ ) > 0 and λ 2 (τ ) = 0. But this leads to the same contradiction:λ Hence λ 1 and λ 2 are positive for t < T . Suppose now that λ 3 = 0. The function u → u U C50+u is strictly increasing for u ≥ 0 and since λ 2 is positive, the Hamiltonian H attains its minimum at u max and the optimal control is constant, u ≡ u max . This is only possible if the entire (p, q)-trajectory except for the terminal point lies below the diagonal. But then the initial conditions are ill-posed.
We henceforth make assumption [A] and consider only well-posed initial data; we also normalize the third multiplier setting λ 3 = 1.
Theorem 3.3. Let u * be an optimal control with corresponding trajectory z * = (p * (t), q * (t), y * (t)) † and multiplier λ and let Then we have that whereũ is the positive solution to the equation In particular, optimal controls are continuous.
Proof. This argument is the same as in the proof of Theorem 3.1 in [8]. We only briefly summarize it giving formulas we shall need later on. Optimal controls minimize the function over the control set [0, u max ]. The first two derivatives of Υ are given by and If Ψ(t) ≤ 0, then Υ is strictly increasing and the minimum over the control set [0, u max ] is taken on at u * (t) = 0. If Ψ(t) > 0, then Υ is strictly convex on (−U C 50 , ∞) with the global minimum at its stationary point. This stationary point is the solution to equation (15). If Ψ(t) < 1, then the stationary point is negative and thus Υ is still strictly increasing over the control set [0, u max ] with the minimum attained for u = 0. As long as the stationary point lies inside the interval , the resulting control is admissible while the minimum is attained for u = u max once the stationary point exceeds u max as Υ is now strictly decreasing over the control set. This gives us equation (14). It follows from the implicit function theorem that the stationary point is smooth as a function of t. While Ψ(t) is positive, we can represent the control in the form u * (t) = max{0, min{ũ(t), u max }}, and thus it follows that u * is continuous. For Ψ(t) ≤ 0, the optimal control is given by u * ≡ 0 which also is the optimal control for 0 ≤ Ψ(t) ≤ 1. Hence optimal controls remain continuous as Ψ becomes nonnegative.
4. Analysis of optimal junctions. For well-posed data, it follows from the transversality conditions that Ψ(T ) = 0 and thus optimal controls end with an interval (τ, T ] where u * (t) ≡ 0. However, the precise sequence of segments when the control lies on the boundary or in the interior needs to be determined and this is the topic of this section. We call any time τ where the form of the optimal control changes from values inside the control interval to a boundary value a junction time and the corresponding point in the state-space a junction. While the optimal control takes values in the interior of the control set it is smooth. Using the dynamics and adjoint equations, and suppressing the variable t, it follows from implicit differentiation of equation (15) that Hence the derivative of the control at a junction has the same sign as λ 2 b − λ 1 ξ. Along the interior control, we have by the stationary point condition that and thus It therefore follows from the fact that H ≡ 0 that Setting x = p q , we obtain that Hence the derivative of the control at a junction has the same sign as the term in braces in equation (24). For a ≥ 0 let L a = (p, q) : p > 0, q > 0, x = p q , bx(ln x − 1) + µ + dp 2 3 + a = 0 .
For a = 0 (i.e., for the uncontrolled system), this curve is the same as the singular base curve for the optimal control problem with a log-linear kill term which has been analyzed in [7]. A detailed derivation of the geometric properties of the curve L 0 is given in Proposition 5.3.3 in [15]. Essentially, L 0 is a closed loop in the first quadrant which starts and ends at the origin and reaches its largest p value on the diagonal. For increasing positive values of a, the loops L a become smaller and for a 1 < a 2 , the loop L a2 lies inside the loop L a1 . We are interested in the value m = G umax U C50+umax 2 and we denote the corresponding loop by L m . We illustrate the geometric shape of these sets in Fig. 2, but refer the reader to [15] for the details of a general mathematical derivation of these properties.
The curves L 0 and L m determine the regions where an optimal control can change from a value in the interior of the control set to a boundary value. If this happens, then the derivative at the junction time must be non-negative if the control changes from u = 0 to the interior control or from the interior control to u = u max and it must be non-positive if the control changes from u = u max to the interior control or from the interior control to u = 0. Hence the first type of junctions cannot happen where this derivative is negative while the second type are impossible if this derivative is positive.
The curve L 0 and the diagonal D 0 = {(p, q) : p = q} form boundary curves for the regions where optimal controls can change from the interior control to u = 0 or vice versa. Let R 1 denote the region outside of the loop L 0 and above the diagonal (D + = {(p, q) : p > q}), R 2 the region inside of the loop L 0 and above the diagonal,  Proof. Under assumption [A], we have that b > µ (based on the biological meaning of the parameters, this condition is always satisfied) and thus the expression defining L 0 is negative on the diagonal near p = 0 inside the loop. Since this expression changes sign only along the curve L 0 , it is negative inside and positive outside of this loop. Therefore it follows from equations (20) and (24) thatu is negative in the regions R 2 and R 4 and positive in the regions R 1 and R 3 . Thus junctions from 0 to the interior control cannot lie in R 2 ∪ R 4 and junctions from the interior control to 0 cannot lie in R 1 ∪ R 3 . Similarly, junctions between the full dose control and the interior control can be analyzed. Let S 1 denote the region outside of the loop L m and above the diagonal, S 2 the region inside of the loop L m and above the diagonal, S 3 the region inside of the loop L m and below the diagonal and S 4 the region outside of the loop L m and below the diagonal. On the diagonal and near the origin the expression defining L m has the same sign as and depending on whether this quantity is positive or negative, different switching scenarios arise. In this paper, we focus on the case when this expression is negative and make the following assumption: [B]: We assume that

Combining assumptions [A] and [B], we therefore have the following condition [AB]
: Proposition 4.2. If condition [AB] is satisfied, then, along an optimal controlled trajectory, there are no junctions from u = u max to the interior controlũ in S 1 ∪ S 3 and there are no junctions from the interior controlũ to u = u max in S 2 ∪ S 4 .
Proof. The argument is the same as for the control u ≡ 0. The expression defining L m is negative inside of the loop and positive outside of it. Hence the derivativeu at the upper control limit u max is positive in the regions S 1 and S 3 and negative in the regions S 2 and S 4 .

5.
Analysis of optimal concatenations. In order to analyze possible concatenations of optimal segments, we need some observations about the flow of trajectories along the constant controls u ≡ 0, u ≡ u max and along the controlũ. We only consider the region D. Along the control u ≡ 0 all trajectories cross the diagonal transversally from the region D + (above the diagonal) into the region D − (below the diagonal). It follows from assumption [A] that trajectories of the system cross the diagonal transversally in the opposite direction, i.e., from D − into D + , along u ≡ u max . This also holds for the interior control.
Lemma 5.1. If the optimal control takes a value in the interior of the control set while the system is on the diagonal, then the corresponding trajectory crosses the diagonal transversally from D − into D + .
Proof. It is convenient to use the projective coordinate x = p q . We need to show that the derivative at x = 1 is positive, i.e.,ẋ |x=1 > 0, along the controlũ. It follows from the dynamics thaṫ Henceẋ Along an extremal control we have that H ≡ 0 and thus, at a point on the diagonal, If the optimal control is given by the stationary controlũ, then we also have from the stationary point condition that and thus on the diagonal verifying the result.
The region S 2 is not positively invariant under the control u = u max . The dynamics has a strong differential-algebraic character with q the fast and p the slow variable. As a result, for a constant control u, trajectories quickly move towards the null-cline N u of the system and then follow this null-cline. The null-cline is given byq = 0, i.e., bx = µ + dp For u = 0, this null-cline lies in D − , but it shifts into D + for sufficiently large constant u. Figure 3 shows a representative example for this scenario for the same parameter values that we used in our numerical computations shown in section 6 (given in Table 1 This point lies in D + and the null-cline enters the loop L m transversally at (p m , q m ). This is easily verified by a direct computation in the coordinates (x, p). If the initial condition p 0 lies above the maximum p-value on the loop L m , i.e., if then trajectories corresponding to the full dose control will cross over from D − into D + at points above L m and will transversally enter this loop at a point very close to (p m , q m ) as shown in Fig. 3.
There exists a δ > 0 such that the region P = S 2 ∩{(p, q) : p ≤ p m +δ} is positively invariant under trajectories corresponding to both the full dose control u max and the controlũ.
Proof. We already have seen that trajectories which start at a point on the diagonal and correspond to the full dose control or to the controlũ enter the region S 2 . Furthermore, above the diagonal, and regardless of the control being used, the tumor volume always decreases. Thus trajectories cannot exit the region P through the horizontal line p = p m + δ. Exit from the region P is therefore only possible through points on the loop L m below p m . Since the null cline N m enters the loop L m transversally, it follows that near-by trajectories corresponding to the constant control u max also enter this loop transversally. Hence, for δ small enough, such trajectories starting on the loop L m with values p ≤ p m + δ enter S 2 . If the interior control is used, thenq will have a larger value than for the control u max . Sinceṗ does not depend on the control, it follows that at such a point the dynamics for the controlũ points to the right of the dynamics for the full dose control. Hence the system also enters S 2 .
Using these results we can limit the switching sequences of optimal controls for initial conditions that represent a typical medical scenario.
Proposition 5.1. Suppose condition [AB] is satisfied and let (p 0 , q 0 , A) be wellposed initial data with (p 0 , q 0 ) ∈ D + and p 0 ≥ p max . Then, while the system lies below the diagonal, optimal controls are at most concatenations of one full dose control segment followed by a segment for the controlũ.
Proof. We only outline the reasoning, but skip the technical details. The optimal control cannot start with u = 0. For, if this were the case and u = 0 on [0, t 1 ], then both p(t 1 ) > p 0 and q(t 1 ) > q 0 . At time t 1 the system has the same amount of inhibitors available, but is in a worse position. It is intuitively clear, but requires verification, that such a control cannot be optimal. Essentially, using the same control shifted to the interval [0, T − t 1 ] without the initial no dose segment will generate a better value. This leaves only full dose segments and segments along the interior controlũ as possible starting controls. Using the same reasoning as above, optimal controls cannot switch to u = 0 while the state lies below the diagonal. Since the initial condition satisfies p 0 ≥ p max , the trajectory will lie in the region S 4 (outside of the loop L m and below the diagonal) where, by Proposition 4.2, no switchings from the interior control to full dose are optimal. Thus, if the control switches to the controlũ, it cannot switch back to full dose while the trajectory lies below the diagonal.
. We note that there need not be a segment when the control is full dose. It is allowed, and quite possible and even likely in many scenarios, that the optimal control is given by the controlũ throughout the time interval when the trajectory lies below the diagonal. An example is given in section 6.
Proposition 5.2. If the optimal controlled trajectory (x * , u * ) enters the region S 2 ∩ {(p, q) : p ≤ p m + δ} at some time τ , then the remaining switching sequence over the interval [τ, T ] is at most a concatenation of a full dose segment (u * ≡ u max ) followed by a segment where the control is given by u * (t) =ũ(t) followed by a final segment when no agents are given (u * ≡ 0): In this sequence, the full dose segment may be absent.
Proof. The region S 2 ∩ {(p, q) : p ≤ p m + δ} is positively invariant under both the full dose u max and the controlũ and no switchings from the interior to full dose controls are optimal while the trajectory lies inside this region. By Proposition 4.1, neither are switchings from u = 0 to the interior control. Hence the remaining switching sequence only allows for the specified concatenation sequence: optimal controls use up all anti-angiogenic agents, possibly with a junction from full dose to the interior control. The point where all inhibitors have been used up lies above the diagonal and thus further tumor reductions still occur along a no dose segment until the trajectory reaches the diagonal. At that point, the minimum possible tumor volume is attained. This defines the terminal time and optimal value for the problem.
Theorem 5.3. Suppose condition [AB] is satisfied and let (p 0 , q 0 , A) be well-posed initial data with (p 0 , q 0 ) ∈ D + and p 0 ≥ p max . Then optimal controls are at most concatenations of the following form: Proof. Since the initial conditions are well-posed, optimal controlled trajectories (x * , u * ) cross the diagonal at some time τ ∈ [0, T ]. It follows from Proposition 5.1 that the control u * is at most a concatenation of a full dose segment followed by a segment along the controlũ over the interval [0, τ ]. Let p d = p * (τ ) and first suppose the optimal trajectory crosses the diagonal along the constant control u max . Since p d > p max , by Proposition 4.2 the optimal control can only switch to the interior control once the corresponding trajectory x * has entered the loop L m . Since this trajectory closely follows the nullcline N m , the trajectory enters the region S 2 ∩ {(p, q) : p ≤ p m + δ} and thus, by Proposition 5.2, the remaining control sequence contains a switch to the control u * (t) =ũ(t) followed by a final no dose segment after all agents have been used up. Thus, in this case the overall concatenation sequence is of the form full dose → interior control → no dose .
If the optimal controlled trajectory crosses the diagonal while the controlũ is used, then, in principle, the full concatenation sequence (31) is possible. For, in this case, in principle, the optimal control could switch back to the full dose control u max in the region S 1 . (The control can only switch to u = 0 for the terminal segment.) The control can then switch back to the interior control only once the trajectory lies in the positive invariant region S 2 ∩ {(p, q) : p ≤ p m + δ} and the remaining sequence is as above. This leads to the full concatenation sequence specified. Corollary 1. Suppose condition [AB] is satisfied and let (p 0 , q 0 , A) be well-posed initial data with (p 0 , q 0 ) ∈ D + and p 0 ≥ p max . If the optimal controlled trajectory crosses the diagonal with a full dose control, then optimal controls are at most concatenations of the following form: full dose → interior control → no dose.
(32) 6. Numerical examples. We give two numerical examples of extremal controlled trajectories that illustrate the two main behaviors of such pairs: one where the diagonal is crossed using the full dose control u max , the other when the diagonal is crossed along the controlũ with values in the interior of the control set. These examples are merely intended to illustrate these mathematical structures. The parameter values used in our calculations are given in Table 1. Because of a symmetry in the equations, without loss of generality the parameter U C 50 can be normalized to U C 50 = 1 and the remaining control parameters u max and A are simply given in terms of multiples of this base unit. The values used here represent just some generic situations. The quotient u U C50+u is dimensionless and, since we have not included an equation for the PK of the agent in the model, we do not list the dimensions of u as either concentrations or dose rates. With an additional linear model for PK, these relations are given by simple conversion factors. The parameters which define the uncontrolled dynamics are based on the biologically validated model in [3]. If the amount of inhibitors is large enough, then the full dose trajectory can reach the loop L m and it follows from our theoretical results that there exists an extremal controlled trajectory whose control follows a full dose → interior control → no dose switching sequence. This extremal control can be computed through a direct procedure. For a reasonably large ε 0 > 0 and 0 ≤ ε ≤ ε 0 define the 1parameter family F ε = {(z ε , u ε )} of controlled trajectories described below and let V (ε) = p ε (T ε ) denote the value of the trajectory corresponding to ε, i.e., the tumor volume at the terminal time T ε for the trajectory z ε .
1. For ε = 0, the control u 0 is constant given by the full dose, u 0 (t) ≡ u max , until all inhibitors are exhausted at time T = A umax . Then the control is concatenated with a no dose segment, u = 0, until the trajectory crosses the diagonal at time T 0 > T . At this point, the tumor volume is recorded as the value for ε = 0, i.e., V (0) = p 0 (T 0 ). 2. For ε > 0 and small, the full dose control is used over the interval [0, T −ε] and at that time the control is changed to the control u ε (t) =ũ(t). At the junction point, the multiplier λ 2 is computed from the stationary point condition (22) while setting the interior control to u max and then λ 1 is computed from the condition that the Hamiltonian (9) vanishes. Simultaneously integrating the dynamics and adjoint equations forward defines the controlũ. This control is used until either all remaining inhibitors have been used up or until the value of the interior control reaches 0.
(a) If all inhibitors are used up before the controlũ reaches 0, then still concatenate the control with a no dose segment until the diagonal is reached and record the tumor volume on the diagonal as the value V (ε) = p ε (T ε ). Note that such controls are not optimal as they are discontinuous. (b) There exists a unique value ε opt when the controlũ reaches the lower limit 0 exactly at the time when all inhibitors are used up. At this point, again concatenate the control with a no dose segment until the diagonal is reached and record the tumor volume on the diagonal as the value V (ε opt ) = p εopt (T εopt ). This control is the only extremal control in the family and we verify the necessary condition for optimality that the terminal value for the multiplier λ 1 is positive. (c) If the interior control reaches 0 before all inhibitors are exhausted, then use up the remaining inhibitors through another full dose segment. Then once again concatenate the control with a no dose segment until the diagonal is reached and record the tumor volume at the end-point as V (ε) = p ε (T ε ). As in the first subcase, such controls are discontinuous and thus are not optimal.
In this family, only the controlled trajectory (z εopt , u εopt ) is an extremal as all the other controls are discontinuous. The minimum value for V is attained at ε opt and this control indeed is locally optimal. (It is almost clear from the construction that this controlled trajectory can be embedded into a field of broken extremals of the same type and this guarantees local optimality [14,Chapter 6].) Figure 4 shows the graph of the value function for this 1-parameter family F ε for the initial data (p 0 , q 0 , A) = (4600, 5800, 12). The sharp minimum is typical for the shape of the value V . Figure 5 shows the optimal control and its corresponding trajectory. The optimal control is characterized by a sharp drop of its value from the maximum dose to no dose and clearly emulates a bang-bang control with one switch from full to no dose. The junction times are t 1 = 7.780 and t 2 = 8.0750 with the terminal time given by T = 8.0754 so that the no dose segment is barely visible on the graph. This example of an optimal solution correlates with optimal solutions which arise in the model with a log-kill pharmacodynamic model of the form Guq when the tumor volumes are small, the singular control in that model has saturated and the singular control is no longer admissible [7]. In that model, however, it is the more common scenario that optimal controls are concatenations of a full dose segment followed by a segment when the control is singular (and takes values in the interior of the control set) and then end with a no dose segment as all inhibitors are exhausted. While the example given above in some sense follows this structure as well, the short interval when the control lies in the interior of the control set much more resembles a bang-bang control. Analogues of the situation when the control is singular over long intervals arise for the model considered in this paper when optimal controls cross the diagonal along the controlũ. A numerical illustration of this situation is given in Figure 6 for the initial conditions (p 0 , q 0 ) = (12000, 15000) and A = 2. Comparing the graph of the control with a typical solution to the original model, it again is evident that the interior control in this case very much mimics even the shape of the singular control for the original model. Actually, computing the  The graphs are shown with U C 50 = 50, U max = 1.5U C 50 = 75 and A = 12U C 50 = 600. The optimal controlled trajectory crosses the diagonal along a full dose segment and the optimal concatenation sequence is full dose →ũ → no dose. The initial full dose segment is shown in black, the intermediate interior control segment in blue and the final no dose segment in black; junctions are marked with an asterisk. The loop L m is shown as the dashed red curve.
controlũ here is significantly harder as the value ofũ on the diagonal is completely determined by the stationary point condition and the fact that the Hamiltonian vanishes. This determines the multiplier λ 2 on the diagonal, but λ 1 is left free.
Here shooting techniques involving the transversality conditions need to be used.
In this case, we do not yet have a firm theoretical basis for the structure of optimal solutions and the result given in Theorem 5.3 only provides an upper bound on the concatenation structure of controls. In numerical computations, we mostly see a lower concatenation sequence like it is shown in Fig. 6. The graphs are shown with U C 50 = 50, U max = 1.5U C 50 = 75 and A = 2U C 50 = 100. The optimal controlled trajectory crosses the diagonal using intermediate doses along the controlũ and the optimal concatenation sequence is of the shortened formũ → no dose. The optimal controlled trajectory is shown in red and as a comparison the corresponding full dose trajectory is shown in blue. Here there is a significant difference in the final tumor volumes.

7.
Conclusion. In this paper, we initiated a case study that compares optimal solutions for the problem to minimize the tumor volume with an a priori given amount of inhibitors for a mathematical model of anti-angiogenic treatment when two qualitatively very different pharmacodynamic models for the effect of the drug are considered. For the system with a linear PD (log-kill assumption), a complete solution for this optimal control problem has been given in [7,15]. However, a linear PD is valid only in a specific limited range and does not account for a saturation of the concentration, an important and relevant phenomenon. In this paper, we started the analysis of solutions for the same dynamics if an E max pharmacodynamic model is used. Mathematically, the two problem formulations lead to very different optimal control problems: the log-kill model for PD gives a Hamiltonian function which is linear in the control generating discontinuous optimal controls while the E max -model generates convexity properties that give continuous optimal controls. In spite of these differences, for both formulations we observe a similar pattern of typical optimal controls starting with an initial full dose period followed by lowered intermediate doses and ending with a no dose period where the tumor volume still decays because of after effects. As for the mathematical analysis, the reformulated model under study here exhibits some simplifications since optimal controls are continuous, but on the other hand also significant difficulties arise.
If the optimal control takes values in the interior of the control set as the system crosses the diagonal, then the value for the control at the crossing point on the diagonal is determined. This severely limits standard numerical techniques to compute such extremals as this phenomenon induces a removable singularity on the diagonal