Modeling and Numerical Solution of a Cancer Therapy Optimal Control Problem

A mathematical optimal-control tumor therapy framework consisting of radioand anti-angiogenesis control strategies that are included in a tumor growth model is investigated. The governing system, resulting from the combination of two well established models, represents the differential constraint of a non-smooth optimal control problem that aims at reducing the volume of the tumor while keeping the radioand anti-angiogenesis chemical dosage to a minimum. Existence of optimal solutions is proved and necessary conditions are formulated in terms of the Pontryagin maximum principle. Based on this principle, a so-called sequential quadratic Hamiltonian (SQH) method is discussed and benchmarked with an “interior point optimizer—a mathematical programming language” (IPOPT-AMPL) algorithm. Results of numerical experiments are presented that successfully validate the SQH solution scheme. Further, it is shown how to choose the optimisation weights in order to obtain treatment functions that successfully reduce the tumor volume to zero.


Introduction
Cancer has a growing impact on our society, because it is among the main causes of illness and death worldwide.On account of this, there exist many treatment options as surgery, chemotherapy, radiation therapy, hormonal therapy, immunotherapy and anti-angiogenic treatment.For all these therapies, it is important models consider the dynamics between the tumor volume p and the carrying capacity q.One of the most commonly used models for tumor growth is based on the Gompertzian growth law as follows ( ) ( ) For p q < the tumor grows ( 0 p >  ) until p q = .For p q > the tumor shrinks ( 0 p <  ) again until p q = is reached.
Next, we consider a time-varying carrying capacity q.The basic idea is a combination of stimulatory (S) and inhibitory (I) effects as follows ( ) ( ) , , .q S p q I p q = −  A modelling issue is the choice of S and I, and for this reason we consider the model proposed by Hahnfeldt et al. [5] as follows 2 3 , q bp dp q = −  (2) with the birth rate 0 b > and the death rate 0 d > .This is a well-recognized mathematical model for time-varying carrying capacity.However, it couples the tumor volume variable to the carrying capacity.
On the other hand, a model of time-varying carrying capacity that does not involve the tumor volume explicitly is due to Ergun et al. [6].This model is computationally convenient since p does not appear in the equation.We have 2 4 3 3 .
Based on validation with real data [5] [6], both models appear promising in the quest of an accurate description of tumor growth.For this reason, we consider a combination of the two models ( 2) and (3) as follows ( ) Together with the equation for the tumor growth (1), we obtain the following differential system that models the evolution of the tumor volume and of the carrying capacity of the vasculature.We have ( ) In Figure 1, we show the evolution of this system for initial values with ( ) ( ) and for different κ .We see that the dynamics obtained with 0 κ = , that is with (2), dominates the evolution of our model given by (4).The model (4) exhibits a steep increase at the initial time for the variable q for 0 0.5 and, choosing κ close to one a slower increase of q can be ob- served.From this result, it appears that 0.5 κ = and 1 κ = are representative of two different dynamical behaviour and we shall consider these choices in the numerical experiments.
In the next two sections, we introduce two control mechanisms in (4) that represent the treatment of cancer by anti-angiogenesis and radiotherapy, respectively [3].

Anti-Angiogenesis
The angiogenesis is a process where a growing tumor develops its own blood vessels, which provide the tumor with oxygen and nutrients.The anti-angiogenesis To model this treatment, we introduce a control u that takes its values in [ ] max 0,u and represents the dose of the anti-angiogenic medicine.With the anti-angiogenic elimination parameter 0 γ > , we can augment the equation for q in (4) as follows ( )  Hence, our model for an anti-angiogenetic mono-therapy is given by ( ) The anti-angiogenic treatment influences the carrying capacity of the vascularity q, but as q appears in the equation for p, it also influences the tumor volume p.

Radiotherapy
Radiotherapy is a treatment that uses ionizing radiation to kill cancer cells.For this purpose and to minimize damage on healthy tissues the tumor should be well localized.bq dq bp dp q qu r qw r r w This model is completely specified by giving the values of the parameters appearing in it.These values are specified in Table 1; see [6] and [8].

The Optimal Control-Treatment Problem
Usually, in the context of optimal control of cancer development models, the objective of the control is to minimize the volume of the tumor at final time, i.e.
( ) p T .See Schättler and Ledzewicz [3] for a detailed discussion of this setting.Now, while we keep this objective, we introduce additional terms that include a reduction of the tumor volume p along the time evolution, and L 2 -and L 1 -norms of the controls u and w.With respect to the side effects of anti-angiogenetic medicine and radiotherapy, it is reasonable to have penalty terms for the corresponding controls.
The parameters , , , , , 0 σ ϑ ν µ ν µ ≥ can be chosen differently to obtain different settings.We refer to this optimal control problem as the (OCAR3) control problem.Notice that in [9] a similar model is considered where only the tumor volume at the final time enters in the controls' objective.

Analysis of the Anti-Angiogenesis and Radio-Therapy Model
In this section, we analyse our cancer development and treatment differential system in (OCAR3).We have the following lemma.
Lemma 1 The model (6) with has the equilibrium points ( ) ( ) , , 0, 0, 0 p q r = and ( ) ( ) , , , ,0 From the second equation with q p = , we obtain Altogether, we have To show that the equilibrium point ( ) , , 0 p q is locally asymptotically stable, we set : ln z p = , : ln y q = and obtain For this system, we denote the equilibrium point with ( ) , z y where ln z p = and ln y q = .We do not consider the third equation for r, because r is not relevant for our analysis since it neither appears in the p or the q equations, nor does p or q appear in the equation for r.Notice that the r-equation is asymptotically stable.
Using Taylor expansion, we obtain the following Jacobi matrix ( ) ( ) At the equilibrium point ( ) , z y , we obtain ( ) We have that the trace is ( ) Since the trace is the sum of the eigenvalues and the determinant is the product of the eigenvalues of f A , we conclude 1 2 0 λ λ + < and 1 2 0 , 0 λ λ < , hence the equilibrium point ( ) , z y is locally asymptotically stable.
Next, we show that the solution to (6) with if the initial values are non-negative.From a medical point of view this is reasonable, because it makes no sense to have negative volumes and capacities.Also values for p and q that are larger than the equilibrium point ( ) , p q cannot be reached according to the stability discussion above.On this account, we restrict our examination to the following domain Next, we consider the uncontrolled case with 1 κ = and neglect the equation for r since in this case r is not coupled to the ( ) , p q -system.We have


The equation for q does not depend on p, hence q grows (independently of p) until q q = (equilibrium point) is reached.As we start with initial values ( ) 0 0 , p q ∈  with 0 0 p q < , the tumor volume initially grows until p q = , but since q is growing (independently of p), we will again have p q < and the tumor continues to grow until p p q q = = = (equilibrium point) is reached.
Hence in this case a solution that starts in  will never leave  .Now, we show the same solution properties for the controlled case with 0 κ = .Consider the following problem ( ) ( ) ln , , .p p p r pw q q bp dp q qu r qw r r w Theorem 1 The set  is positively invariant for the system (10).
Proof.We show that if ( ) 0 0 , p q ∈  and , u w are arbitrary admissible con- trols defined on some interval [ ] 0,T , then the solution , the unique solution exists and is given by  thus the vector field points into  .Now, we look at the boundary set ( ) { } , : 0 , p q p p q q ∈ < < =  . Note that the nullclines for 0 q =  are given by , : , bp q E p r dp r νω γν η δ ω = = + + + for the constant controls u ν = and w ω = .We rewrite the q-equation as follows , .
q bp dp q uq r qw E p r q dp dp q q r q νω γ η δ E νω is strictly increasing for fixed r with

E p r E p r p q
νω νω < ≤ = Hence 0 q <  for trajectories starting in points ( ) , p q with 0 p p < < and q q = .
On the coordinates axes for 0 p = and 0 q = the dynamics has singularities, in the sense that ln p p q ξ   −     is not defined.Therefore, we consider the lines p xq = for 0 x = .Now, it is sufficient to show that x  is positive for small ∈ exists, because as we discuss below, our differential equation has an absolutely continuous solution.Hence, the region  is posi- tively invariant for system (10).□ Now, let us again look at the system (6) with , 0 u w = .We take a closer look at the q-equation given by ( )


For initial values in  the solutions to this equation with 0 κ = or 1 κ = are positive and tend to zero as the equilibrium point ( ) , , 0 p q is reached (see Lemma 1).Since q  is the linear combination for 0 1 κ < < , it is also positive and tends to zero as the equilibrium is reached.The pand r-equations do not depend on κ and therefore we can deduce that for initial values  the solution cannot be negative for an arbitrary κ ∈ .The model is well defined and by applying Caratheodory's Theorem (see for example [10], Theorem 54 Proposition C.3.6) one can show existence and uniqueness of a solution to (6).

Optimal Solutions and the Pontryagin Maximum Principle
In this section, we discuss existence of optimal controls and their characterization by optimality conditions.
The existence of an optimal solution to the (OCAR3) control problem can be shown as follows.We know that all solutions ( ) , , x p q r = to (6) are in ( ) is weakly sequentially compact (see [11] Theorem 2.11).As our objective : , , , , x u w J x u w  is convex and continuous (see [12] Theorem 2.14) and thus weakly lower semicontinuous (see [11] Theorem 2.12).Therefore we apply the direct method in calculus of variation and obtain the existence of an optimal solution to our problem (OCAR3).
Next, we discuss the characterization of optimal solutions in the framework of the Pontryagin maximum principle (PMP).For this purpose and for ease of exposition, we consider a general controlled differential model and control setting with the following properties.We have 1) An open and connected state space 3) The controlled dynamical system ( ) , , , is given by the function , , , , t x u f t x u  .We assume that the partial derivative ( ) is continuous as a function of all variables.
4) The class  of admissible controls is taken to be a set of piecewise continuous functions u defined on a compact interval [ ] ⋅ is called admissible if it is a solution to the differential Equation (11) and if ( ) The objective of the control u ∈  is to minimize the following functional

T J x u L s x s u s s g x T
, and the derivative with respect to x is continuous.
Our optimal control problem is now given as follows

T J x u L s x s u s s g x T x f t x t u t x x u t U t T
Notice that the optimal control problem (OCAR3) that was introduced in Section 1 is of this form.
Definition 2 The Hamiltonian function :  for the optimal control problem (12) is defined as follows In our case we have 0 0 λ > and we can assume 0 1 λ = without loss of gene- rality.This situation is called normal extremal lift; see [13]. Theorem ) with the terminal conditions ( ) ( )

Numerical Optimization
In this section, we deal with the numerical implementation of our control framework that belongs to the class of optimize-before-discretize methods.However, in the case of a discretize-before-optimize approach, one could consider the method proposed in [15].In the first part of this section, we introduce our optimization scheme.In the second part, we refer to the operating mode of the IPOPT solver that we use together with the programming language AMPL for benchmarking our SQH scheme.
The main idea of the SQH scheme is the straightforward pointwise minimization of the Hamiltonian function in a way that has been first proposed in [16] [17].Notice that, in this approach, the pointwise update of the control is performed on all grid points and only thereafter the corresponding state is computed.However, this pointwise update of the control may result in large changes of the value of the state variable that makes the proposed approach less robust.
This problem is also discussed in [16] where this issue is left open.On the other hand, we have the quadratic regularisation method proposed in [18], where the Hamiltonian function is augmented with a weighted quadratic term that penalizes deviations from the previous control value to keep the updates of the control sufficiently small.Further, in [18] every pointwise update of the control is followed by a global update of the state variable, which makes this approach very time consuming.In the SQH scheme, we combine the advantages of the two schemes performing a pointwise update of the augmented Hamiltonian on all grid points and recalculating the state variable after the control update.The augmented Hamiltonian is given by ˆ, , , , , , , ,   ˆ:  , , , , , , , We remark that by increasing  a sufficient descent of the cost functional in each iteration can be achieved, see ([19], Lemma 4.1).
While we refer to [19] for a detailed analysis of convergence of the SQH scheme, in the following algorithm we present the implementation of this method.For this purpose, we denote     To determine these therapies, an optimal control problem was formulated considering a cost functional including the tumor volume and L 1 -and L 2 -penalty terms for the controls.After the proof of existence of minimizers, the necessary optimality conditions that characterize these minimizers were deduced in the framework of the Pontryagin maximum principle.Based on this PMP framework, the SQH method was used for numerical solution.This algorithm was used to solve the optimal cancer therapy problem with different experimental settings.Furthermore, optimal solutions obtained by the SQH algorithm were compared with the optimal solution obtained by the IPOPT solver together with the programming language AMPL.This comparison showed that the SQH method is faster by a factor 10 than IPOPT.In a final series of experiments it was shown that it is actually possible to choose the optimisation parameters in such a way to reduce the volume of the tumor and the related carrying capacity to zero.

3 )
Calculate x  from (10) corresponding to v  and calculate τ κ < : STOP and return k v .Else go to 2).

Figure 7 .
Figure 7. Third experimental setting: optimal solution with 2 u µ = and While the proliferation rate a of the cells is constant, the death rate