OPTIMAL CONTROL FOR A MATHEMATICAL MODEL FOR CHEMOTHERAPY WITH PHARMACOMETRICS

. An optimal control problem for an abstract mathematical model for cancer chemotherapy is considered. The dynamics is for a single drug and includes pharmacodynamic (PD) and pharmacokinetic (PK) models. The aim is to point out qualitative changes in the structures of optimal controls that occur as these pharmacometric models are varied. This concerns (i) changes in the PD-model for the eﬀectiveness of the drug ( e


Introduction
Optimal control for mathematical models of cancer chemotherapy has a long history dating back to the 1970s and 1980s and the fundamental work by Eisen [5], Swierniak and Kimmel [10,39] and Swan [36,37].This work continued throughout the 1990s (e.g., see the monograph by Martin and Teo [26] or [40]), and revitalized in the 2000s when novel approaches to cancer therapy such as angiogenic treatments (e.g., [18]) or immunotherapies (e.g., [28]) became available.This research continues strongly into the present time as the numerous references in [32] attest to.
While the literature on this topic is large, comparatively little attention has been given to the importance which pharmacometric models have in the overall process.Such models describe the relations between a drug's dose rate u, its concentration c in the bloodstream (or other central and peripheral compartments), and the actual effects e which the drug has.Specifically, pharmacokinetic models (PK) describe the relations between the drug dose rate and the drug's concentration ("what the body does to the drug") while pharmacodynamic models (PD) describe the effects that the drug has on the disease ("what the drug does to the body").Models for PK are generally described by low-dimensional time-invariant linear systems with real eigenvalues [8]; PD models can be highly nonlinear functional relations of the form ϕ(c) often described by generalised Hill functions [25,30].The question of interest underlying this research is whether, and if so, to what extent, the mathematical models used in the modeling determine the structure of optimal controls.This is not only a mathematically non-trivial topic, but it is also of relevant practical interest [17,19,32].
In this paper, we consider a mathematical model for cancer chemotherapy with a single agent as an optimal control problem.We take the point of view that the total amount of agents which is deemed tolerable for treatment has been determined a priori based on a medical assessment of its side effects.This decision may have been made based on general guidelines or based on side effects the drug has which are not included in the model.The aim then is to use the agent in a "best possible" manner.In optimal control theory this is modelled by minimizing an objective function which tries to quantify terms linked to the overall progress of the disease.Here we simply interpret this as to minimize the tumor volume with the given amount of the agent.In the modeling, we include a 1-compartment pharmacokinetic equation and general terms of the form e = ϕ(c)x (with x denoting a compartment of interest, e.g., tumor cells, immune system, healthy cells, etc.) model pharmacodynamic effects.If the optimal control problem is affine in the control -a common feature of such models -then optimal solutions are typically given by finite concatenations of bang and singular controls.Bang controls describe the common practical treatment approach of piecewise constant therapy protocols that alternate between periods of full-dose treatment and rest periods.Singular controls, on the other hand, correspond to specific time-varying dosing regimens.While less conventional, such controls, which administer the agent at lower than maximum tolerable doses, often are the best choices for the mathematical models.For example, they are an integral portion and completely characterize optimal solutions for anti-angiogenic treatments, both as monotherapy [18] and in combinations with chemo- [27] and radiotherapies [13,21].
In this paper, we derive formulas that characterize singular controls and their optimality within an abstract mathematical model with a generic model for PD.These results are then specified to the PD models commonly used in practice.These include Skipper's linear log-kill model, a Michaelis-Menten type E max -model -probably the most commonly used equation in the pharmacological industry -and sigmoidal models described by generalised Hill-type functions.The results are illustrated with a mathematical model for anti-angiogenic treatment with E max -type pharmacokinetic model.We also compare the solutions with the solutions for the corresponding model when PK is not taken into account, i.e., when drug dosage and concentration are identified.For the latter problem formulation, optimal controls are continuous (agreeing with an interpretations of controls as drug concentrations) and have portions when the controls lie in the interior of the control set.These portions are closely related to the singular controls and their trajectories in the full model analyzed here.For the example under consideration, the dynamics of the concentrations -along an interior control for one model and along a singular control for the other model -are identical.Our observations about the significance of the models for PK and PD will be summarized in a discussion section at the end of the paper.

A Brief discussion of pharmacometrics of drugs
Pharmacometrics is the branch of pharmacology which quantifies how a substance administered to a living organism behaves in the organism.It comprises pharmacokinetics (PK) and pharmacodynamics (PD).Mathematically, this is the study of the relationship between a drug's dose rate, its concentration and its effects on the body (see Fig. 1).Various models have been developed to represent these processes and we briefly discuss the most common ones.

Pharmacokinetics (PK)
Pharmacokinetic equations model the time evolution of a drug's concentration in the blood plasma.The standard model that describes the concentration c = c(t) of the drug in the bloodstream in response to a continuous dose rate u = u(t) is one of exponential growth and decay (see Fig. 2) given by ċ = −ρc + u, c(0) = 0, (2.1)  with ρ the clearance rate of the drug.If u max denotes the maximum dose rate, the maximum achievable concentration c max is given by When administration of the drug is stopped (u ≡ 0), the concentration dissipates at an exponential rate ρ.We note that the maximum concentration is proportional to the maximum dose rate and if this one is multiplied tenfold, then so is the concentration in this model.While such a relation may be reasonable over a certain range of dose rates, it is not if the dose rates become very high.This so-called 1-compartment model is the simplest way to model the relations between the dose rate and the concentration of the drug.In more detailed models, the concentrations in the bloodstream are separated from the concentrations in a peripheral compartment (such as tissue or an organ of interest) and the interactions between those compartments are taken into account (see Fig. 3).This generally is modelled by a 2-dimensional linear system with positive eigenvalues.Higher dimensional compartmental models are used rarely, but arise, for example, if the peripheral compartment is divided further into a "shallow" and "deep" compartment.For example, a 3-compartment model is often used for insulin pumps [34].
In this paper, we only consider a 1-compartment model, but use a slightly more general bilinear equation From a mathematical perspective, this equation causes no additional difficulties in our framework and it reduces to the standard linear model if we choose η = 0.For this reason we use equation (2.3) throughout the paper.The second parameter η allows to regulate the increase in the concentration.When no drugs are administered, the clearance rate is ρ.In this framework the maximum concentration exhibits a Michaelis-Menten type dependence on the maximum dose rate which is given by (2.4) It thus saturates at the value 1 η which gives some biological meaning to this quantity.Furthermore, for the model to be meaningful we need to impose the restriction ρ + ηu max > 0. We also remark in closing this discussion that the Michaelis-Menten shape arises in enzyme kinematics from exactly such a bilinear relationship (e.g., see [4,15]).

Pharmacodynamics (PD)
Pharmacodynamic models describe the effect e which a drug has at concentration c.They are usually given by functional relations of the form e = s(c)x with x denoting the cells in a compartment on which the drug is acting (e.g., tumor cells, immune system, healthy cells, etc.) and the coefficient s(c) representing the concentration dependent killing rate of the drug.There exist several approaches to model the function s = s(c) with a linear relation s = σc based on the linear log-kill hypothesis [35] probably the most common model.Here σ is a positive constant σ.This model is applicable over a certain range of concentrations, but generally is not a valid model for low or high drug concentrations.This is related to the fact that some drugs have little to no effect if their concentrations are too low and, on the other hand, their effectiveness saturates at high concentrations.It is thus reasonable to model the effect of a single chemotherapeutic agent more generally with a function s that depends on the concentration and takes values in some bounded interval [0, E max ].Commonly used forms are the E max model that more accurately describes the intensity of the effect of "fast" acting drugs for high concentrations (see Fig. 4) and sigmoidal functions that also capture the behavior at lower concentrations (see Fig. 5).Mathematically, these models are given by the equations (2.5) Here EC 50 denotes the concentration at which half of the maximum effect E max is realized.Both EC 50 and E max are commonly used parameters in pharmacology.Thus the E max model is described by a Michaelis-Menten (MM) type equation while generalised Hill-type functions are used for the second term.We remark that although this is commonly the case, n need not be an integer.
In practice, drugs are typically administered in the upper ranges of the dose level-concentration curves and here the E max -model is adequate.For this reason, in this paper we focus on this model more, but still carry  out the calculations with the general form s = s(c).Using such a model introduces nonlinearities which not only from the underlying medical aspect, but also from the point of view of the mathematical modeling lead to significant and interesting questions.While an optimal control problem with a linear log-kill model is affine in the controls which leads to bang and singular controls as candidates for optimality, convexity/concavity properties are introduced if a Michaelis-Menten term is used and this changes the structure of optimal controls.

Optimal control for a general model for chemotherapy with pharmacometrics
We formulate a mathematical model for chemotherapy with a single pharmaceutical agent and general dynamics.Pharmacodynamics are modelled by a generic, strictly increasing function s which later will be specified to the models discussed above in Section 2.2.Pharmacokinetics is modelled through a bilinear system which includes the standard linear 1-compartment model as special case η = 0, but makes the equation somewhat more flexible.The optimal control problem that will be considered is control-affine which leads to concatenations of bang and singular controls as the prime candidates for optimality.Based on the necessary conditions for optimality of the Pontryagin maximum principle [29], we compute the singular controls and, for various choices of the PD model s, discuss their optimality.

Problem formulation
We consider a general dynamics of the form ẋ = f (x) + s(c)g(x), x(0) = x 0 . (3.1) The variable c ∈ R + is a non-negative scalar which denotes the concentration of the chemotherapeutic agent and the state x ∈ D ⊂ R n + is a column vector of non-negative numbers which represents cells in various compartments that are included in the modeling.These may be cancer cells, but possibly also other cells related to the vasculature or tumor-immune system interactions.The dynamics is described by smooth vector fields f and g defined on D. The drift vector field f describes the evolution of the (mathematically) uncontrolled system when no drugs are given.(In medicine, this is commonly referred to as the "control".)The control vector field g models the effects which the drug at concentration c has on the system under the pharmacodynamic model s = s(c).We always assume that s : [0, ∞) → [0, ∞), c → s(c), is a monotonically strictly increasing function which satisfies s(0) = 0.
We denote the dose rate of the drug by u and this is the control of the system.We include a 1-compartment model for the PK of the drug given by the bilinear equation with ρ the clearance rate of the drug.Formally, controls are taken as Lebesgue measurable functions with values in a compact interval, u : [0, T ] → [0, u max ].This assumption guarantees the existence of optimal controls and, indeed, for such a model optimal controls need not be piecewise continuous [11].
The variables x and c together describe the state of the system and states are always functions of time t.Generally, however, in order not to burden the notation, we drop the argument t.As a generic rule, Greek letters denote parameters/constants. Lemma 3.1.Given a Lebesgue measurable function with values in [0, u max ], the concentration c is well defined over [0, T ] and for all times satisfies In particular, it holds that ηc(t) < 1 for all times.
Proof.It follows from u ≥ 0 and c(0) = 0 that c(t) ≥ 0 for all t ∈ [0, T ].The equilibrium point of the dynamics for the constant control u(t) ≡ u max is given by c max = umax ρ+ηumax and thus the concentration always lies in the range [0, c max ).
Whereas a solution to equation (3.2) always exists on the full interval [0, T ] on which the control is defined, this need not be the case for the dynamics (3.1).It rather depends on the vector fields f and g whether this will be the case.Since the concentration c is bounded and since s is continuous, the function s(c) is bounded and it follows from well-known results in ordinary differential equations (e.g., see [9]) that the solution x will exist on all of [0, T ] if the vector fields f and g are linearly bounded, i.e., if there exist non-negative constants a and b such that For all the standard models for tumor development (linear, logistic or Gompertzian growth) this will be satisfied.We call a Lebesgue measurable function u : [0, T ] → [0, u max ] an admissible control if the corresponding trajectory (x, c) exists on the full interval [0, T ].Side effects of treatment need to be included in the modeling.Here we take the point of view that an acceptable amount of drugs has been determined a priori (for example, based on the toxicity of the agent) and thus the total amount of agents to be given is limited: The question thus becomes how to use this amount in the "best possible" way.Mathematically, the total amount of drugs to be given is incorporated into the dynamics as an isoperimetric constraint.This is done by adding the trivial dynamics which simply tracks the amount of agents administered to the system.Using † to denote the transpose, we define a new and extended state variable z as the column vector z = (x, c, y) † ⊂ R n+2 + and define the augmented drift and control vector fields F and G as The dynamics can then be expressed in the concise form The aim is to control the system in some sense in the "best possible" way.Now, what exactly is meant by this is up to interpretation.It is generally accepted that some compromise needs to be made between minimizing the tumor cells (i.e., kill as much of the cancer as possible) and side effects.In the approach taken here side effects have already been limited through the total amount A of agents to be administered and therefore we aim to maximize the tumor kill or, equivalently, aim to minimize the tumor cells.Thus we take the objective in the form with α and β n-dimensional row vectors of non-negative weights.These weights will be taken positive for components of x that correspond to tumor cells, but they may be taken as zero for other compartments.We thus consider the following optimal control problem: [CC] Minimize the functional J over all admissible controls u : [0, T ] → [0, u max ] subject to the dynamics described by the equations and the terminal constraint y(T ) ≤ A.
Alternatively, instead of a priori specifying the amount of agents to be used, one could build in a compromise between minimizing tumor cells and side effects into the objective and take it in the form with γ and δ additional non-negative weights.With an objective of this type side effects are modelled indirectly and one can put the emphasis either on the concentration (γ = 0, δ > 0) or the dose rates (γ > 0, δ = 0).Mathematically, the two formulations are almost identical and in this paper, for sake of clarity, we only consider one approach.The general framework based on minimizing the criterion (3.9) is discussed in [22].
We do want to emphasize, however, that in both approaches an L 1 -integral of either the dose rate or the concentration is used to model the total dose (amount of drugs used).This corresponds to standard practice in the pharmaceutical industry of using AUC (area under the curve) and allows us to relate these terms to the side effects of treatment.Other models, such as when quadratic expressions in the control are used, lead to mathematically simpler frameworks, but lack a meaningful interpretation.

Necessary conditions for optimality
Necessary conditions for optimality for the optimal control problem [CC] are given by the Pontryagin maximum principle [29].(For some more recent references on the topic, see [1,2,31].)These conditions are phrased in terms of the so-called (control) Hamiltonian H of the problem.In this function the constraints of the problem are adjoined to the objective with multipliers.According to an interpretation of multipliers as linear functionals, we consequently write multipliers as row vectors and indicate this by an asterisk as superscript at the respective spaces.Thus, introducing multipliers λ ∈ (R n ) * , µ ∈ R * and ν ∈ R * for each of the equations in the dynamics and a multiplier λ 0 ∈ R * for the objective functional, the Hamiltonian takes the following form for problem [CC]: Combining all multipliers into one covector Λ = (λ, µ, ν) ∈ (R n+2 ) * and using the state z and the augmented vector fields F and G, the Hamiltonian simply is the inner product of the multipliers with the objective and the dynamics: If u * is an optimal control for problem [CC] with corresponding trajectory x * = (x * , c * , y * ), then it follows from the Pontryagin Maximum Principle [29] that there exists an absolutely continuous co-vector Λ : [0, T ] → (R n+2 ) * and a constant λ 0 ≥ 0 such that the following conditions are satisfied: (i) [non-triviality condition] the multipliers λ 0 and λ(t) do not vanish simultaneously: here D denotes the Jacobian matrix of the partial derivatives with respect to the state x; (iii) [transversality conditions] λ(T ) = λ 0 α, µ(T ) = 0 and ν(T ) ≥ 0, (iv) [minimization condition] the Hamiltonian function H is minimized point-wise over the control set a.e. by the optimal control along the covector λ and the optimal trajectory x * .That is, (vi) the complementary slackness condition ν(y * (T ) − A) = 0 holds.In particular, by equation (3.14) the multiplier ν is a non-negative constant.
In vector form the adjoint equation and transversality conditions read Proof.If λ 0 = 0, then it follows from the adjoint equations and transversality conditions that λ and µ vanish identically over the interval [0, T ].The non-triviality condition (i) on the multipliers therefore implies that ν > 0 and thus the minimum condition on the Hamiltonian function H gives us that u * (t) ≡ 0. This, however, is not a feasible solution for a well-posed model and we thus exclude it from our considerations.Henceforth we normalize the multiplier λ 0 to λ 0 = 1 and drop it in the notation.Since the Hamiltonian function H is linear in the control, we obtain from the minimization condition that optimal controls satisfy the following condition: (3.17) The function Φ = Φ(t) is the term multiplying the control in the definition of the Hamiltonian H and is called the switching function of the optimal control problem.The name comes from the fact that optimal controls switch between the extreme values 0 and u max of the control set if Φ(τ ) = 0 and Φ(τ ) = 0.The corresponding controls are called bang-bang controls.If the switching function vanishes identically over some non-empty open interval J, then the minimization property by itself does not give any information about the control.In this case, however, also all derivatives of Φ(t) vanish identically on J and this typically determines the corresponding controls.Such controls are called singular.In the biomedical scenario, bang-bang controls correspond to protocols with switchings between maximal dose, often called MTD (maximum tolerated dose), and rest periods.Singular controls, on the other hand, represent time-varying administrations of a drug using lower doses.The question whether optimal controls are bang-bang or singular thus has an immediate interpretation and practical relevance for the structure of optimal treatment protocols.The following result, the generalized Legendre-Clebsch Condition (e.g., see [1,31]) gives a higher order necessary condition for optimality of singular controls. (3.18)

Singular controls for problem [CC]
We compute singular controls for the problem [CC] and evaluate the Legendre-Clebsch condition.For this we need an efficient way to compute the derivatives of the switching function.We recall that the Lie bracket of two vector fields A and B on R n is given by where D denotes the matrix of partial derivatives with respect to the coordinate z.The following formula is verified by a direct calculation (e.g., see [31]).Lemma 3.5.Let z(•) be a solution to the dynamics (3.7) for the control u and let Λ be a solution to the corresponding adjoint equation.Then, for any continuously differentiable vector field L, the derivative of the function is given by For problem [CC] we have the following formulas for the Lie brackets of F and G up to order 2: The switching function for problem [CC] is given by with z = z * evaluated along the optimal controlled trajectory.For simplicity of notation, we drop the time variable t and the asterisks denoting the optimal trajectory in these formulas.The first and second derivatives of the switching function are thus given by The coefficient at the control u in the second derivative is given by and a singular control is of order 1 if this term is not equal to zero.The singular control is then computed by solving equation (3.23) for u * (t).These expressions simplify along a singular control.If an optimal control u * is singular on an open interval I, then it follows for all t ∈ I from the fact that the switching function and its derivative vanish identically on I that Substituting these expressions into the above formulas, we obtain If ν = 0, then it follows from equation (3.25) that we have µ(t) ≡ 0 for t ∈ I. Equation (3.26), respectively the adjoint equation (3.13), give that λ(t), g(x * (t)) ≡ 0 and thus ∂ Φ ∂u (t) ≡ 0, i.e., a singular control is of higher order.Furthermore, since the Hamiltonian vanishes identically, this also implies that Using the above relations it follows that and additional relations follow by differentiating this equation.This procedure generates an over-determined system of equations which generically has no solutions [31] and is difficult to satisfy in particular cases.Nevertheless, this could be possible, but it requires knowledge of f and g to preclude this.
In the following, we make the generic assumption that ν is positive.In this case, both µ and λ(t), g(x * (t)) are negative along an optimal singular control.We thus get the following result from equation (3.27): Proposition 3.6.An extremal singular control u * on an open interval I is of intrinsic order 1 if and only if The strengthened Legendre-Clebsch condition for minimality is satisfied if and only if We have some immediate corollaries that apply to the models for PD given in Section 2.2.Recall that we assume that the function s is strictly increasing and that the multiplier ν is positive.
Corollary 3.7.(a) If η ≥ 0 and s is a strictly concave function, then singular controls are of order one and the strengthened Legendre-Clebsch condition for minimizing the objective J is satisfied.(b) If η > 0 and s is concave, then singular controls are of order one and the strengthened Legendre-Clebsch condition for minimizing the objective J is satisfied.(c) If η = 0, i.e., for the standard linear model of PK, the strengthened Legendre-Clebsch condition for minimizing the objective J is satisfied where s is strictly concave and it is violated where s is strictly convex.
Proof.In case (a) we have that s (c) s (c) < 0 while this quantity is non-positive in case (b).In either case, 2η − (1 − ηc) s (c) s (c) > 0. In case (c), we have that s (c) s (c) > 0 if s is convex and this quotient is negative if s is concave which gives the result.Proof.Direct calculations verify that s is a strictly increasing function which has a unique inflexion point at changing from strictly convex below c 0 to strictly concave above this value.The statement therefore follows from condition (c) in the previous corollary.
We still give the formula for a singular control of order 1.In this case, solving equation (3.23) for the control gives that Using the previously calculated Lie brackets and equation (3.27), this formula takes the form Recall that the dynamics for the concentration is given by and thus the first term in equation (3.30) is the zero-dynamics for the concentration.Hence the evolution of the concentration along a singular control is given by the second term.In the case η = 0, the traditional linear PK model, this dynamics takes the form

.31)
Together with the equations for the state x and the multiplier λ, and away from the hyper-surface in the cotangent bundle defined by λ, g(x) = 0, this differential equation forms a closed system from which the other multipliers and the control have been eliminated.This system can therefore be integrated to obtain a formula for the optimal concentration, but this will need to be done numerically.

Anti-angiogenic treatment with Michaelis-Menten type pharmacodynamics
We apply these results to a well-studied mathematical model for tumor development under angiogenic signaling formulated by Hahnfeldt, Panigrahy, Folkman and Hlatky [7].In that paper, linear models are used both for PK and PD.
Dropping the PK equation, this model was considered by some of us as an optimal control problem with the aim of minimizing the tumor volume for an a priori given amount of agents.A complete solution in terms of a regular synthesis of optimal controlled trajectories has been given [18,31].Optimal solutions in that case typically are concatenations of a short full dose segment followed by a long piece where the control is singular and end with a short period when no drugs are given.On this interval after effects of the previous administration still reduce the tumor volume.The most important structural element in this solution is the singular arc.When PK is taken into account, the order of this arc increases from 1 to 2 and this leads to optimal solutions which have chattering arcs (Fuller phenomenon) [20] and are not practicable.Nevertheless, knowing the structure of the optimal solution allows us to give excellent and simple sub-optimal approximations [11,12].
In [24] we have considered this model without PK, but with a saturating E max -model for PD.This leads to significant mathematical changes.The nonlinearities introduced into the dependence of the model on the control now make optimal controls continuous in agreement with an interpretation of the controls as concentrations.
Here, based on the formulas developed above, we consider the changes that the addition of a linear model for PK causes.Once more, now optimal controls are concatenations of bang and singular pieces with the singular controls of order 1.This generates similar structures as in the original version.

A brief review of the model
We briefly review the underlying model.In order to avoid undue duplications, however, we refer to the existing literature for both a deeper discussion of the modeling and some of the mathematical details.All of these can all be found, for example, in [31].
In this model by Hahnfeldt et al. [7], the spatial aspects of the underlying consumption-diffusion process that stimulate and inhibit angiogenesis are incorporated into a non-spatial 2-compartment model with the primary tumor volume, p, and the carrying capacity of the vasculature, q, as its principal variables.The dynamics thus consists of two ODEs that describe the evolution of the tumor volume and its carrying capacity and are given by In this formulation, the variable u represents the concentration of the anti-angiogenic agent and its effect is modelled by the linear log-kill term E max uq with E max representing the efficacy of the agent.The carrying capacity q and tumor volume p 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.The term bp represents the stimulation of the vasculature by the tumor and is taken proportional to the tumor volume.The term dp 2 3 q models the inhibition which is represented as an interaction term between the tumor surface and the vasculature.The coefficients b and d are mnemonically labelled for birth and death, respectively.A listing of the parameters and the numerical values used in calculations is given in Table 1.The parameter values for the dynamics are taken from the paper [7].
Here we instead use a Michaelis Menten term to model PD in the form E max c EC50+c q.Mathematically, it is possible to normalize EC 50 to EC 50 = 1 by scaling the concentration c in units of EC 50 , but this is not commonly done in pharmacology and hence this parameter is retained.Adding a standard linear model for PK, we consider the following optimal control problem: [A] Given a fixed amount A of anti-angiogenic agents, among all Lebesgue measurable functions u : [0, T ] → [0, u max ] defined over an interval [0, T ] with the terminal time T free, find the control u (representing the dosage of the drug) which minimizes the tumor size J(u) = p(T ) at the terminal time subject to the dynamics ṗ = −ξp ln p q , p(0) = p 0 , (4.1) q − E max c EC 50 + c q, q(0) = q 0 , (4.2) We note that there are no data available for the parameters occurring in the E max formulation for PD and these values have been chosen to illustrate some of the mathematical features.
Administering anti-angiogenic drugs 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.If the a priori given amount A of inhibitors is too small, it is not possible to reduce the tumor volume for some initial conditions (p 0 , q 0 ).In this case, the mathematically optimal solution for problem [A] is given by T = 0, i.e., all that is possible in such a case is to slow down tumor growth, but the tumor cannot be reduced.This generates mathematically degenerate situations which are easily analyzed, but which we exclude here for simplicity of presentation.We call the data (p 0 , q 0 , A) well-posed if an objective value better than p 0 is achievable and in this case the final time T along an optimal control is positive.
We furthermore make the simplifying assumption that This condition ensures that the effectiveness of the drug on the diagonal (p = q) is high enough to lower the vasculature.Under this assumption, it can be shown [18,31] that optimal controls for well-posed initial data have the following three intuitive properties: (i) optimal controls end with an interval (τ, T ] where the control is u * (t) ≡ 0; (ii) the endpoint of the optimal controlled trajectory lies on the diagonal: p * (T ) = q * (T ); and (iii) all inhibitors have been used up:

Singular controls for problem [A]
In the framework of the general formulation in Section 3, for the optimal control problem [A] we have that x = (p, q) † , α = 1, β = 0, η = 0, s(c) = Emaxc EC50+c , and the vector fields f , g and their Lie bracket [f, g] are given by The Hamiltonian for problem [A] takes the form and the adjoint equations read The last multiplier ν is a non-negative constant.Suppose an optimal control u * is singular over an open interval I.The switching function Φ is given by Φ = µ + ν and thus µ(t) = const for t ∈ I.By equation (4.11), we then have that If ν = 0, then these relations imply that µ and λ 2 also vanish identically.Since H ≡ 0 as well, this implies that λ 1 also vanishes on I contradicting the terminal condition λ 1 (T ) = 1.Thus the multiplier ν is positive.
For problem [A], equation (3.30) for the singular control reduces to The multipliers can be eliminated from this equation and it is possible to determine the singular control as a feedback function of the state variables (p, q, c) † alone.For, using equation (4.12) and the fact that the Hamiltonian vanishes identically, it follows that (4.14) Figure 6.A cross section through the graph of the singular control u sing (c; p, q) for c = 2.
This equation allows us to eliminate the term ξ λ1 λ2 from equation (4.13).The singular control therefore can be computed as a feedback control in the form Naturally, the general formulas for the singular control derived from Φ ≡ 0 ignore the aspect of whether the resulting control will be admissible, i.e., satisfies the prescribed bounds on the control.This always needs to be analyzed for the specific model under consideration.Here the singular control is a feedback function and thus, if we fix the point x = (p, q) † , the resulting parameter dependent function u sing = u sing (c; p, q) can be analyzed.Direct computations verify that the function c → u sing (c; p, q) is strictly increasing above the diagonal (p > q) and that it is strictly convex above the diagonal and strictly concave below the diagonal: Furthermore, the equations u sing (c; p, q) = 0 and u sing (c; p, q) = u max are quadratic in c and can therefore be solved explicitly, albeit with messy expressions.These properties allow us to determine the region where the singular control is admissible.Figure 6 shows an example of a cross section through the graph of the singular control u sing = u sing (c; p, q).The singular control has a singularity on the diagonal (as it is obvious from Eq. (4.15)) with positive values above the diagonal -the main region of interest -which increase with smaller tumor volumes.The latter feature is reminiscent of the same property for the model with a log-kill model for PD.

Discussion
Clearly, any change in the modeling, be it changes in the models for PK and/or PD, or simply some changes in parameter values, will cause quantitative changes.Conceptually it is more significant if qualitative changes arise as this may alter the overall approach to the problem.We have shown in this and some of our earlier papers [23,24] that, in the absence of a model for PK, this is the case if a linear log-kill model for PD is replaced with an E max -model.Incorporating a linear PK model into the dynamics which uses the E max -model for PD then restores the original structure of optimal controls as concatenations of bang and singular arcs.the yellow trajectory corresponds to giving the total amount in one full dose session while the total dose is divided into two equal parts with a rest period for the blue trajectory.The points when the control switches from administration of the drug to a rest period are marked with a red cross on the corresponding trajectories.Administration of the drug then starts again when the diagonal p = q is reached.In order to motivate the search for optimal controls, we include some comments to the point that even as simple a question as whether it is advantageous to administer a given dose "all at once", so-to-speak, or if it would be more beneficial to split it into smaller sub-doses is not easily decided.Figures 7 and 8, using model [A] as an example, illustrate effects which the clearance rate of the agent (even in a linear model for PK) has on optimal solutions.There seems to be no clear a priori rule that would indicate whether it is better to give agents as a single dose up-front or to divide the total dose into several equal smaller doses.
In our figures, in order to better visualize tumor reductions, we plot p as the vertical variable and q as the horizontal variable.All trajectories start from the initial condition (p 0 ; q 0 ) = (5000; 5000) on the diagonal, give the same total amount of agents, A = 300, and use the parameter values from Table 1 for the dynamics.In these computations, we only compare protocols which give the same total dose A in 1, 2, 3 or 4 equal portions.If the total dose is given all at once, we still let the tumor volume decrease until the minimum tumor volume is realized when the trajectory crosses the diagonal p = q.If the total dose is broken up into equal smaller sub-doses, then at the end of the administration of each sub-dose, a rest period (u = 0) is inserted.The duration of these rest periods, however, varies and the time is chosen until the system reaches the diagonal.During this interval the tumor volume still decreases while it would increase again if the rest period would be longer as the system moves into the region below the diagonal.As such controls are clearly not optimal, we do not consider such cases.In Figure 7 we show the trajectories when either all drugs are given in one session (yellow trajectory) or the total amount is divided into two equal portions (blue trajectory).In the figure on the left the clearance rate ρ is chosen as ρ = 1 while it is chosen as ρ = 0.3 in the panel on the right.While the trajectory which administers all drugs in one segment does better for the larger clearance rate ρ = 1, it is the trajectory which divides the dose into two pieces which does better for ρ = 0.3.The reason lies in the fact that the concentrations decays much slower in the second case and thus a greater benefit (in the sense of a larger tumor reduction) is taken from the rest period.As the concentration remains higher for a longer time, this leads to better tumor reductions.Yet, if the number of portions into which the total dose is divided is increased further, no additional benefit is gained.Figure 8 shows trajectories which correspond to divisions into 3 (green) and 4 (red) segments and these trajectories perform worse than the blue one with 2 pieces, the red one quite a bit worse.Thus, increasing the number of segments further does not lead to better outcomes.The point we want to make with this simple comparison is that the question of how an optimal dose rate looks like, simply cannot be decided by just looking at a few ad hoc choices.
We note that the initial burst in the vasculature seen in the trajectories at their starting point also is an effect of the PK as the concentration first needs to build up to some level before it becomes effective.At the later times when the trajectories reach the diagonal again after the rest periods, concentrations are still high enough so that this phenomenon is not seen again.
We close with a discussion of what motivated the calculations given in this paper.There exist good heuristic reasons for a modeling approach that simplifies the system by ignoring selective (so-to-speak higher order) effects and then incorporate them, as needed, later on.For medical problem formulations this approach is clearly taken in the approach to PK and PD in the literature and we simply are interested to what extent this is justified.In the particular example considered above, model [A] for anti-angiogenic therapy, there indeed exists a close relation between the singular control for the full model (with E max PD and linear PK) and the problem when the PK model (4.3) is dropped from the system.In the latter system, the concentration of the agent is identified with the dose rate and this combined quantity becomes the control.This problem has been analyzed in [24] and we briefly compare these structures without delving into the computations.In order to distinguish between the concentrations in the model, however, we write v for the control (i.e., concentration dosage) in that model and we add a tilde to the corresponding variables.It is shown in [24] that optimal controls v * are continuous and, as long as they lie in the interior of the control set, these controls are the solutions of the equation Formally, this equation is identical with equation (4.12) of this paper if we normalize the extra multiplier as ν = 1 ρ .While the optimal control v * takes values in the interior of the control set, it is differentiable and it follows from implicit differentiation of equation (5.1) that We give two numerically computed examples of extremal controls.Figure 9 shows an example from [24] where the optimal control only has a small interior piece at the end which connects the full dose control with the control u = 0.It follows from necessary conditions for optimality for the problem considered in that paper that optimal controls, as for the problem considered here, must end with a rest period and thus such a connecting piece must exist as optimal controls are continuous.For the model when PK has been taken into account this no longer is required and in this case the corresponding control (not shown) is simply bang-bang with a switch from maximum dose to u = 0 close to the end.This agrees with an interpretation of the controls as dose rates while they represent concentrations in the previous model.Figure 10 shows an example of an extremal where a singular segment is present in the control and it is demonstrated with a simple comparison that this control does better than the bang-bang control.The initial conditions are (p 0 ; q 0 ) = (12 000; 15 000) and here we consider only a small total amount of drugs, A = 46.75.While the differences in the terminal state of the tumor volume are small, they are still significant with the small amount of drugs.Generally, however, as the formula for the singular control has a singularity on the diagonal, optimal trajectories which contain singular arcs are only found when the system reaches a region above the diagonal in the state space which is at a good distance from the diagonal.This is reminiscent of the identical behavior for the optimal controls for the original model from [7] with a linear log-kill PD when the PK model is dropped [18,31].While there exist some clear similarities between all four model (with linear or E max PD and with or without PK), these are not evident from the onset and it required some detailed analysis of necessary conditions for optimality to uncover these for this particular problem formulation.It is a commonly observed procedure in models for drug treatment not to include a model for PK initially and take the point of view that this can be analyzed separately and be added later on.While this is a good heuristic, it nevertheless leads to some awkward ambiguities.Clearly, in such a model the controls represent the concentrations of the agent.Yet, in many models optimal controls will be discontinuous resembling dose rates.This discrepancy is unavoidable and is caused by identifying dose rates and the concentration in the control.Adding a PK model resolves this issue, but if a linear log-kill model is used for PD this will cause difficulties if the optimal controls contain lower-dose singular controls (see, for example, [11]).If instead an E max -model is used for PD, such interpretational discrepancies disappear as the controls now are continuous in the model without PK and the addition of the PK-model clearly separates dose rates and concentrations.In a certain sense, the model formulation considered in this paper, i.e., the formulation with a E max model for PD and linear PK, is the clearest from an interpretation point of view: concentrations are piecewise continuously differentiable functions generated by possibly piecewise continuous controls.The chattering phenomenon which is possible for models with a linear log-kill PD disappears.Overall, from the point of view of interpreting the results and relating them to the underlying problem, this modeling approach presented here therefore seems to be the most satisfactory one.

Conclusion
We have analyzed a general version of an optimal control problem for cancer chemotherapy with a single drug when various pharmacodynamic (PD) and pharmacokinetic (PK) models are included in the modeling.If PK is not included in the modeling, the controls represent drug concentrations and optimal controls change continuously between full dose protocols and no dose rest periods along smooth connecting pieces that take values in the interior of the control set.Once PK is taken into account, the role of the control changes from concentrations to dose rates and optimal controls are typically discontinuous and they may contain singular arcs which represent intermediate, time-varying dose rates.The qualitative properties of optimal controls not only depend on whether or not PK is taken into account, but also on specifics of the functional relation which is

Figure 1 .
Figure 1.A schematic representation of the interactions of pharmacometric models.

Figure 4 .
Figure 4. E max model for the PD of a drug.

Figure 5 .
Figure 5. Sigmoidal model for the PD of a drug.

Definition 3 . 3 .
Let Γ = (Λ, z * , u * ) be an extremal lift for the optimal control problem [CC] which is singular on an open interval I ⊂ [0, T ], i.e., the switching function Φ vanishes identically on I.The corresponding control u * is said to be singular on I and the response of the system z * is called a singular arc.The singular control is said to be of intrinsic order k over I if the first 2k − 1 derivatives of the switching function do not depend on the control u while the coefficient at the control u in the derivative of order 2k of the switching function does not vanish.

Theorem 3 . 4 .
Let u * be an optimal control for problem [CC] which is singular of intrinsic order k over an open interval I ⊂ [0, T ] with corresponding extremal lift Γ = (Λ, z * , u * ).Then it holds that

25 ) and ( 1 −
ηc) s (c) λ, g(x) ≡ ρµ.(3.26) Case (a) of this corollary directly applies to the E max -model, s(c) = Emaxc EC50+c , while case (b) applies to the standard linear log-kill function s(c) = sc.The situation is more complicated for a Hill-type sigmoidal function, s(c) = Emaxc n EC n 50 +c n , as s has an inflexion point.Corollary 3.8.Let η = 0 and consider the function s(c) = Emaxc n EC n 50 +c n .Then the strengthened Legendre-Clebsch condition for minimizing J is satisfied where

Figure 7 .
Figure 7. Trajectories for piece-wise constant controls using ρ = 1 (left) and ρ = 0.3 (right):the yellow trajectory corresponds to giving the total amount in one full dose session while the total dose is divided into two equal parts with a rest period for the blue trajectory.The points when the control switches from administration of the drug to a rest period are marked with a red cross on the corresponding trajectories.Administration of the drug then starts again when the diagonal p = q is reached.

. 2 )Figure 9 .
Figure 9.An optimal control (left) and corresponding trajectory (right) for the optimal control problem when PK is not taken into account.The initial condition is (p 0 ; q 0 ) = (4600; 5800) and A = 600.The initial full dose segment is shown in black, the intermediate interior control segment is shown in blue and the final no dose segment is shown in black.Junctions (i.e., the points on the trajectory which correspond to changes in the control between interior and boundary values) are marked with an asterisk.The dashed red curve L m is part of a loop that bounds a region where junctions between full dose and interior controls are allowed.

Figure 10 .
Figure 10.Comparison of a controlled trajectory which has a singular arc (red curve) with a trajectory which follows a full dose-no dose structure (blue curve) for model [A].The two junction points from full dose to the singular control and from the singular control to no dose are shown with a red asterisk on the red curve and the junction point from full dose to no dose is marked by a blue asterisk on the blue curve.The numerical differences in the final tumor values for this example are small and the graph on the right gives a blow-up of the trajectories near the terminal point.
16))Triples Γ = (Λ, z * , u * ) consisting of a controlled trajectory (z * , u * ) and a multiplier Λ for which all the conditions (i)-(vi) of the Pontryagin maximum principle stated above are satisfied are called extremals.They are called normal if the multiplier λ 0 corresponding to the objective is positive and abnormal if this multiplier is zero.
Lemma 3.2.If the optimal control u * does not vanish identically, then extremals for problem [CC] are normal, i.e., the multiplier λ 0 is positive.

Table 1 .
Variables and parameter values used in numerical simulations.