ROBUSTNESS OF OPTIMAL CONTROLS FOR A CLASS OF MATHEMATICAL MODELS FOR TUMOR ANTI-ANGIOGENESIS

Abstract. We describe optimal protocols for a class of mathematical models for tumor anti-angiogenesis for the problem of minimizing the tumor volume with an a priori given amount of vessel disruptive agents. The family of models is based on a biologically validated model by Hahnfeldt et al. [9] and includes a modification by Ergun et al. [6], but also provides two new variations that interpolate the dynamics for the vascular support between these existing models. The biological reasoning for the modifications of the models will be presented and we will show that despite quite different modeling assumptions, the qualitative structure of optimal controls is robust. For all the systems in the class of models considered here an optimal singular arc is the defining element and all the syntheses of optimal controlled trajectories are qualitatively equivalent with quantitative differences easily computed.


Introduction.
A growing tumor needs a steady supply of oxygen and nutrients for cell duplication.Initially, during avascular growth, it is provided through the surrounding environment.As the tumor becomes larger these mechanisms become inadequate and tumor cells enter the dormant stage of the cell cycle.As a consequence, vascular endothelial growth factors (VEGF) are released that stimulate the formulation of new blood vessels and capillaries in order to supply the tumor with needed nutrients.This process is called tumor angiogenesis.Tumor antiangiogenesis is a treatment approach for cancer that aims at depriving the tumor of this vasculature.Ideally, without an adequate support network, the tumor shrinks.
Anti-angiogenic treatment was already proposed in the early seventies by J. Folkman [7], but became medically possible only with the discovery of the inhibitory mechanisms of the tumor in the nineties [4,12].It brings in external anti-angiogenic 2 HEINZ SCH ÄTTLER URSZULA LEDZEWICZ AND BENJAMIN CARDWELL agents to disrupt the growth of endothelial cells which form the lining of newly developing blood vessels and capillaries.Rather than targeting the fast duplicating, genetically unstable and continuously mutating tumor cells, this indirect treatment approach targets the genetically stable endothelial cells.As a consequence, no clonal resistance to angiogenic inhibitors has been observed in experimental cancer [1].Since developing drug resistance all too often is the limiting factor in conventional chemotherapy treatments, tumor anti-angiogenesis has been called a new hope for the treatment of tumors [13].Although these high hopes have not been realized in practice, there still is strong interest and active research on tumor antiangiogenesis as a method that normalizes the vasculature [10,11] and thus, when combined with traditional treatments like chemotherapy or radiotherapy, enhances the efficiency of these procedures.
In this paper, we formulate a class of mathematical models for tumor antiangiogenesis as optimal control problems.Specifically, we consider the problem of how to schedule an a priori given amount of anti-angiogenic (e.g., vessel disruptive) agents in order to minimize the tumor volume.(For similar formulations with a modified objective where the tumor volume is minimized over a fixed therapy horizon, see the work of Swierniak et al. [23,24].)This problem will be analyzed for a class of mathematical models that include and are based on a model that was developed and biologically validated by Hahnfeldt, Panigrahy, Folkman and Hlatky [9], a group of researchers then at Harvard Medical School.The principal state variables are the primary tumor volume, p, and the carrying capacity of the vasculature, q.The latter is a measure for the tumor volume sustainable by the vascular network.The dynamics describes the interactions between these variables.The tumor volume p changes according to some growth function dependent on the variable carrying capacity q and the q-dynamics consists of a balance of stimulatory and inhibitory effects.Hahnfeldt et al. carry out an asymptotic analysis of the underlying consumption-diffusion model that leads to the form for these terms proposed in [9] (see section 2).However, there exists some freedom in the modeling and by now the original formulation has undergone several modifications with a main one by Ergun, Camphausen and Wein [6].While significant modeling changes are made in the dynamics for the vascular support in this model, the solutions to the optimal control problem we consider are in fact qualitatively identical.Here we more generally consider a family of models that in a certain sense interpolate the dynamics for the vascular support between these two models.We show that optimal controlled trajectories for all these models are characterized by the fact that there exists a unique optimal "path" in (p, q)-space along which the optimal tumor reductions are realized.Optimal controls steer the system to this path as quickly as possible (using maximum dosages) and then follow the optimal path determined by a so-called singular arc (using specific partial dosages) until all inhibitors run out.We compare these optimal paths for the various models.
2. Mathematical Modeling.The models we consider are minimally parameterised and population based.They incorporate the spatial aspects of the underlying consumption-diffusion processes that stimulate and inhibit angiogenesis into a nonspatial 2-compartment model with the primary tumor volume, p, and its carrying capacity, q, as variables.Intuitively, the latter defines the ideal tumor volume sustainable by the vascular network.It is related to the volume of endothelial cells and for this reason we also call it the endothelial support of the tumor for short.
The growth of the primary tumor volume p can generally be modelled in a form ṗ = ξpF p q (1) where ξ denotes a tumor growth parameter and F is a function of the scalar variable x = p q .The function F is assumed to have basic properties that reflect cell growth.Typically F : (0, ∞) → R, x → F (x), satisfies F (1) = 0 and is taken to be twice continuously differentiable and strictly decreasing.Given the definition of the variables, these are natural conditions to impose: since q is the carrying capacity, for p = q the endothelial support and tumor volume are balanced and thus p should not change whereas the tumor volume should shrink for inadequate endothelial support (p > q) and increase if ample support is available (p < q).It also is reasonable to assume that these processes are more pronounced the smaller the quotient x becomes.Many commonly used growth models, such as Gompertzian and logistic growth functions, fit this description.In this paper, as in the original models in [9] and [6], we employ a Gompertzian growth model with F (x) = − ln x.Clearly, other growth functions like, for example, logistic models, are equally realistic.But for the optimal control problems considered here we already vary the dynamics for the carrying capacity of the vasculature and keeping a general model for tumor growth leads to too broad a picture.For this reason we restrict our analysis to one specific model here.In the original variables we thus have that ṗ = −ξp ln p q . (2) If one defines y = ln(p), then this simply corresponds to an exponential growth model on y with a forcing term defined by the carrying capacity q, ẏ = −ξy +ξ ln(q).There exist classical studies in the medical literature that support this assumption for various tumors like, for example, breast cancer [20,21].For this reason, and although it clearly has shortcomings for small tumor volumes, this equation is widely used for modeling cancer growth.
The dynamics for the carrying capacity q consists of a balance between stimulatory and inhibitory effects.Its basic structure can be written in the form q = S(p, q) − I(p, q) − µq where I and S, respectively, denote endogenous inhibition and stimulation terms and the term µq that has been separated describes the net balance between endothelial cell proliferation and loss to the endothelial cells through natural causes (death etc.).These effects are small when compared with the stimulation and inhibition exerted by the tumor and µ typically is a small constant that is set to 0 in many models.It is included here, but with the understanding that its value is small in absolute value.The components that define the dynamics are the functional forms for the endogenous inhibition and stimulation terms.Hahnfeldt et al. [9] carry out a spatial analysis of the underlying consumption-diffusion model that leads to the following two principal conclusions for the relations between endogenous inhibition and stimulation: 1.The inhibitor will impact endothelial cells in a way that grows like volume of cancer cells to the power 2 3 .The exponent 2  3 arises through the interplay of the surface of the tumor through which the inhibitor needs to be released with the volume of endothelial cells that mainly determine the carrying capacity q.
2. The inhibitor term tends to grow at a rate of q α p β faster than the stimulator term where α + β = 2 3 .Based on the first conclusion in [9] the inhibitor term is taken in the form with d a constant, the "death" rate.But there exists freedom in the choice of α and β and this has become a source for other models considered in the literature.In the original work [9] the natural assumption is made to take the stimulation term proportional to tumor volume, with b a constant, here mnemonically labelled as the "birth" rate.This corresponds to the choice α = 1 and ).But based on the underlying derivation in [9] other choices are possible and, for example, taking α = 0 and β = 2  3 results in the equally simple form S(p, q) = bq when the stimulation is proportional to the carrying capacity.The second choice generates a considerably simpler q-dynamics for the model in which q factors.Since p and q are expected to be closely related in steady-state, one expects similar behavior for these models.In [5] d'Onofrio and Gandolfi fully analyze both dynamics with special attention to the effects of periodic treatment regimes.
A third model that more strongly makes the steady state assumption of relating p with q was formulated by Ergun, Camphausen and Wein in [6].In fact, identifying p and q in the q-dynamics, in this paper the following forms are used for inhibition and stimulation I(q) = dq and S(q) = bq This choice eliminates any p-dependence from the dynamics for the endothelial support.The authors' motivation for this approximation lies in a different balance for the substitution of stimulation and inhibition that slows down the q-dynamics.
In the original model [9] this dynamics is very fast and as a result the steady-state is reached quickly.In the modification in [6] the inhibition term is only taken proportional to the tumor radius (the exponent 2 3 in the first conclusion of the analysis from [9] is replaced with 1  3 ) and thus the premises of this model are no longer consistent with the implications of the analysis in [9].As a possible justification for the choice bq 2 3 for the stimulation term it could be argued that the necrotic core of the tumor does not contribute to the stimulation of the vasculature and thus the exponent 2/3 could be interpreted as scaling down the stimulation effects from the tumor's volume p to its surface area p 2 3 and then, as the argument is made in [6], replacing p with q in the steady-state analysis.The second premise of [9], 3 , is retained in the modification in [6] and, with p and q interchangeable, this term thus becomes dq 4 3 .The main advantage of replacing p with q lies in a significant mathematical simplification.On the other hand, this step eliminates a direct link between tumor volume p and endothelial support q and thus the overall consequences of this modeling change might appear drastic.Nevertheless, as we shall see, its implications on the structure of optimal controls are not.
In this paper we consider a formulation that interpolates between these two approaches, i.e., as in [6] retains the exponent 1  3 to model the impact of inhibitors Model inhibition I(p, q) stimulation S(p, q) Reference (H 1 ) dp [6] Table 1.Models for inhibition and stimulation on the vasculature, but does not replace p with q.This leads to the following inhibition and stimulation terms with which is consistent with the second premise of [9].More generally, however, we shall consider the stimulation term of the form with θ a parameter.The choice θ = 1 thus corresponds to the term chosen in [9] while θ = 2 3 is consistent with the modification made in [6].Note that for θ = 1 we have I S = d b qp − 2 3 and thus α + β = 1 3 violating the second modeling premise of [9].Thus these models can be thought of as interpolating between the models of [9] and [6].We summarize and label the q-dynamics of all four models in Table 1.
3. Anti-Angiogenic Treatment as an Optimal Control Problem.We now add a control u that represents treatment with an anti-angiogenic agent to the model.The variable u denotes the angiogenic dose rate and the control set U is chosen as a compact interval, U = [0, a], with a denoting an a priori set maximum dosage.As a first approximation we do not include pharmacokinetic equations for the anti-angiogenic agents and thus identify the dosage of the agents with their concentration in the plasma.Hence the effects on the carrying capacity are simply modelled by subtracting a term γqu from the dynamics for q with the constant coefficient γ representing the anti-angiogenic killing parameter.
One of the advantages of anti-angiogenic treatment is that there are basically no serious side-effects.But most anti-angiogenic agents are so-called "biological drugs" that need to be grown in a lab and hence are limited and very expensive.The problem of how to administer a fixed amount of inhibitors to achieve the "best possible" effect thus arises naturally.A formulation suggested by Ergun et al. in [6] and then taken up by us in [15,16,17,18] is to maximize the tumor reduction achievable with an a priori given amount of angiogenic inhibitors, Adding this isoperimetric constraint as a third variable, mathematically this becomes the following optimal control problem: [OC]: For a free terminal time T , minimize the value p(T ) subject to the dynamics q = S(p, q) − I(p, q) − µq − γqu, q(0) = q 0 , ( over all Lebesgue measurable functions u : [0, T ] → [0, a] for which the corresponding trajectory satisfies y(T ) ≤ A.
It is clear from the problem formulation that the variables p and q need to be positive.For the models under consideration this indeed is satisfied for any admissible control function (see, e.g., [5,15,17]) and need not be postulated as an additional constraint.It is also clear from the dynamics for p, (1), that p increases for p < q regardless of what control is being used.As a result, for some degenerate initial conditions (p 0 , q 0 ) it is possible that the (mathematically) optimal solution to problem [OC] is given by T = 0.This situation arises when the amount of available inhibitors simply is too small to reach a point that would have a lower p-value than p 0 .In such a case it is not possible to decrease the tumor volume with the available amount of inhibitors and thus the mathematically "optimal" solution for problem [OC] simply is to do nothing and take T = 0.It is still possible to slow down the tumor's growth, for example, by giving the full dose u = a until all inhibitors run out (but this need not be the best way of doing this, [17]).Hence the formulation [OC] includes a number of subcases which we simply want to exclude here for simplicity of presentation.We thus make the following definition.Definition 3.1.We say the initial condition (p 0 , q 0 , A) is ill-posed for the system under consideration if for no admissible control it is possible to reach a point (p, q) with p < p 0 .In this case the optimal solution for the problem [OC] is given by T = 0.The initial condition (p 0 , q 0 , A) is well-posed if the final time T along the optimal control is positive.
Clearly, whether or not a given initial condition (p 0 , q 0 , A) is well-posed depends on the specific system under consideration (growth function, inhibition and stimulation terms, values of the parameters etc.).The notion of an ill-posed initial condition is inherent to the dynamics and simply reflects the fact that inadequate amounts of anti-angiogenic agents cannot reduce the tumor volume.If, rather than limiting these amounts a priori as it is done in (10), this integral would be added to the objective, then formally this issue does not arise.But for certain initial conditions it nonetheless exists in the form of optimal solutions that give unrealistically high amounts of anti-angiogenic agents.For our problem any initial condition that satisfies p 0 ≥ q 0 is well-posed and henceforth we only consider well-posed initial conditions (p 0 , q 0 , A).
The mathematical analysis of the problem starts with the application of the Pontryagin Maximum Principle (e.g., see [22,2,3]) which provides first-order necessary conditions for optimality of a control u.For a row-vector λ = (λ 1 , λ 2 , λ 3 ) ∈ (R 3 ) * , define the Hamiltonian H = H(λ, p, q, u) as If u * is an optimal control defined over the interval [0, T ] with corresponding trajectory (p * , q * , y * ), then there exists an absolutely continuous co-vector, λ : [0, T ] → (R 3 ) * , such that the following conditions hold: (a) λ 3 is constant, and λ 1 and λ 2 satisfy the adjoint equations with transversality conditions (b) for almost every time t ∈ [0, T ] the optimal control u * (t) minimizes the Hamiltonian along (λ(t), p * (t), q * (t)) over the control set [0, a] with minimum value given by 0.
Remark: Controlled trajectories for which there exists a multiplier λ such that these conditions are satisfied are called extremals.We gave a simplified version of the conditions of the Maximum Principle where an extra non-negative multiplier λ 0 associated with the objective to be minimized is positive and has been normalized to λ 0 = 1, so-called normal extremals.In the case when λ 0 = 0 (abnormal extremals), it can easily be seen that any corresponding control must be constant, a trivial case that need not be of any interest to us here.
Since the Hamiltonian H is linear in u, it follows that the minimizing control typically is given by one of the boundary values u = 0 or u = a, the so-called bang controls.Defining the switching function Φ as we have that If Φ(t) = 0, then any value u ∈ [0, a] satisfies the minimum condition.However, if Φ(t) ≡ 0 on an open interval I, then also all derivatives of Φ(t) must vanish on I and this typically allows to compute the control.Controls of this kind are called singular [2].Optimal controls then need to be synthesized from these candidates through an analysis of the switching function.For example, if Φ(τ ) = 0, but Φ(τ ) = 0, then the control switches between u = 0 and u = a depending on the sign of Φ(τ ).In general, however, it is a challenging problem to determine the optimal structure of concatenations between bang and singular controls.
For the models (H 0 ) and (H 1 ), see Table 1, considered in [9] and the modification (E) from [6] we have given complete solutions to the optimal control problem [OC] in the papers [15,17,18].While optimal controls are bang-bang with just one switching for the system (H 0 ), both for the original formulation (H 1 ) and the modification (E) optimal controls typically start with a full dose segment, u = a, that steers the system to an optimal singular arc, a specific curve in (p, q)-space that gives the response to a singular control.Optimal trajectories then follow this singular arc using the time-dependent singular control until all available anti-angiogenic agents have been exhausted.At the end there still is an interval corresponding to the control u = 0 when due to after effects an additional reduction in the tumor volume occurs until the minimum is realized as the trajectory crosses the diagonal, p = q.Although significant modifications have been made in the modeling assumptions for these two models, the syntheses of optimal controlled trajectories are qualitatively identical.
Figs. 1 and 2 illustrate and compare the optimal syntheses for the models (H 1 ) and (E) for a set of parameter values that are taken from [9].While the equations that define the optimal singular arcs naturally differ (see [15,17]), the concatenation sequences for optimal trajectories are identical and in both cases it is the singular arc that determines the overall optimal structures.We shall show below that this singular arc is preserved in all the models (I θ ) and that its geometry is virtually identical with the one for the model (H 1 ).We also would like to point out that the assumption of a Gompertzian tumor growth model in ( 2) is important for these structures since it guarantees the existence of the singular arc.For a classical logistic growth model it can be shown that singular arcs do not exist [15,23,24] and optimal controls are bang-bang.Interestingly, the number of switchings is small and even for this case the structure of optimal controlled trajectories can be thought of as being the same one as it is given here, but with the interpretation that the length of the singular arc has shrunk to zero. 4. Singular Control and Singular Arc for the Models (I θ ).We write the state of the system as z = (p, q, y) T and express the dynamics in the form where In this notation the switching function becomes the inner product of the multiplier λ = (λ 1 , λ 2 , λ 3 ) and the vector field g along the solution z(t), Φ(t) = λ(t), g(z(t)) .More generally, the derivative of a function where h is any continuously differentiable vector field, can easily be computed directly and is given by with [f, h] denoting the Lie bracket of the vector fields f and h given by and Df the matrix of the partial derivatives of f .We therefore have that and for the model (I θ ) direct calculations verify that these brackets are given by and It is a necessary condition for minimality of the singular control, the Legendre-Clebsch condition [2], that If this quantity is negative, i.e., if the strengthened Legendre-Clebsch condition holds, then we can formally solve the equation Φ(t) ≡ 0 (see ( 26)) for the singular control as this model we have that The switching function Φ(t) = λ 3 − λ 2 (t)γq * (t) vanishes along a singular arc.It follows that λ 3 cannot vanish.For, otherwise λ 2 must vanish as well and then by ( 16) also λ 1 vanishes identically along the singular arc.Since λ 1 and λ 2 are solutions to a homogeneous linear differential equation this contradicts λ 1 (T ) = 1.It then can be shown as in [17] that λ 3 actually must be positive.Hence λ 2 is positive along a singular arc.Thus the strengthened Legendre-Clebsch condition is satisfied and singular controls will be locally optimal.
The vector fields g, [f, g] and [g, [f, g]] are linearly independent everywhere and therefore we can express the vector field [f, [f, g]] as a linear combination of this basis, say Hence But λ vanishes against g and [f, g] along the singular arc, 5. Comparing the structure with the optimal syntheses for the models (H 1 ) and (E), the same features are easily recognized.While the actual geometric shape of the singular arc more closely resembles the one for the model (H 1 ), the trajectories corresponding to the control u = 0 much closer follow trajectories similar to those for the model (E).

Conclusion.
Assuming a Gompertzian tumor growth model, for a wide class of mathematical models for tumor anti-angiogenesis, including models that make significant differences in their modeling assumptions, there exists a unique path in (p, q)-space along which optimal tumor reductions are achieved.This path is given by a so-called singular arc whose precise analytical structure can be determined explicitly, but of course depends on the functions entering the dynamics of the models.These solutions are fully robust with respect to parameter variations.In our simulations above we have used parameter values from [9] that were biologically validated, but the underlying theoretical results do not need to make any assumptions about these numerical values and are generally valid.Naturally, the feedback control that defines the optimal dosages along the singular arc is not a feasible treatment strategy.However, for both the models (H 1 ) and (E) there exist excellent, and themselves fully robust suboptimal approximations to the optimal controls given by piecewise constant controls [14,19].In fact, these approximations also are highly insensitive towards the initial value q 0 of the carrying capacity and only show strong dependence on the initial tumor volume p 0 .For example, for model (H 1 ) and parameter values that are specified in [19], Fig. 6 compares the minimum tumor volumes that are achievable with the optimal control (red curve), with constant full (blue) or half (green) dose treatments, and with another constant dose treatment when this dose is given by the averaged value of the optimal control (black curve).In these computations the initial tumor volume was held constant at p 0 = 12, 000 [mm 3 ] and the initial value q 0 for the carrying capacity was varied from q 0 = 2, 500 [mm 3 ] to q 0 = 12, 000 [mm 3 ].For the optimal and averaged optimal dose the changes in the minimum tumor volume as the initial condition for the carrying capacity varies remain within 1% of the best value.Given the fact that q 0 is an idealistic quantity that cannot be measured accurately, this feature of the optimal protocols would seem to be of practical importance.Also note that the averaged optimal control, a simple, easily computed and realizable treatment protocol, already provides excellent approximations to the best possible values [19].Knowing the structure of the mathematically optimal controls for the various models, as some were established in this paper, allows to judge these much simpler and realistic protocols.

Figure 2 .
Figure 2. Optimal synthesis for the model (E)

Figure 6 .
Figure 6.Value functions for various control strategies for problem (H 1 ) as a function of the initial value q 0 of the carrying capacity for fixed initial tumor volume p 0 = 12, 000 [mm 3 ]