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 to balance the benefits of each treatment with its negative side effects. Therefore a natural mathematical approach to cancer therapy is to consider a mathematical model of the time evolution of tumor that includes the action of the therapy as a control mechanism with the purpose to minimize the tumor volume, while keeping at a minimum the negative side effects on the healthy cells. In order to find an optimal therapy, we formulate an optimal control problem that requires to minimize the tumor volume in a given time horizon and to maximize the health-related quality of life of the patient.
Concerning previous works on optimal control in drug therapy, we refer to, e.g.,    . However, while these references discuss models for immunotherapy, the combination of radiation therapy and anti-angiogenic treatment as in  and in our work is not considered. Another new aspect of our cancer therapy model is the combination of two tumor growth systems, one proposed by Hahnfeld et al.  and the other one by Ergun et al.  . We discuss our model in Section 2, where we illustrate the inclusion of control mechanisms and provide values of the model’s parameters that result from real data.
The typical optimization objective used in cancer therapy models that are considered in the literature is to minimize the tumor volume at the terminal time; see  and  . In this paper, we investigate a more general cost functional that corresponds to minimizing the tumor volume along the entire time horizon and including L1- and L2-costs of the treatment modelled by the control functions. In Section 2.3, we discuss this new cost functional and formulate our optimal control problem. In Section 3, we provide a detailed analysis of our tumor development and treatment model and discuss equilibrium points of this dynamical system and the positivity of solutions. In Section 4, we prove existence of optimal solutions to our optimal control problem and their characterization by the optimality conditions in the framework of the Pontryagin maximum principle (PMP). In Section 5, we illustrate the sequential quadratic Hamiltonian (SQH) method. Also in this section, we consider the operating mode of the “interior point optimizer―a mathematical programming language” (IPOPT-AMPL) solver that we use as a benchmark for our optimization method, before we focus on our SQH scheme. In Section 6, we present results of numerical experiments that successfully validate the effectiveness of our numerical optimization procedure that computes the same results as the IPOPT-AMPL scheme, while requiring considerably less time. Further, we use our SQH solution scheme to show that it is possible to set the optimisation weights in such a way to obtain treatment functions that successfully reduce the tumor volume to zero. A section of conclusion completes this work.
2. Mathematical Modeling of Cancer Development and Treatment
We investigate a new mathematical model for cancer development and treatment resulting from a combination of two complementary mathematical models. Both 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
While the proliferation rate a of the cells is constant, the death rate increases with a growing tumor volume p. The value q, where the proliferation rate equals the death rate, is called the carrying capacity given by
Using this normalized carrying capacity, we obtain
For the tumor grows ( ) until . For the tumor shrinks ( ) again until 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
A modelling issue is the choice of S and I, and for this reason we consider the model proposed by Hahnfeldt et al.  as follows
with the birth rate and the death rate . 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.  . This model is computationally convenient since p does not appear in the equation. We have
Based on validation with real data   , 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
where . 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 , 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 and, choosing close to one a slower increase of q can be observed. From this result, it appears that and 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  .
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
Figure 1. Solutions to (4) for different κ.
therapy is an indirect treatment since it does not fight the tumor cells directly but influences the tumor’s micro-environment, in particular the vasculature. The lack of oxygen and nutrients will force the tumor to shrink.
To model this treatment, we introduce a control u that takes its values in and represents the dose of the anti-angiogenic medicine. With the anti-angiogenic elimination parameter , 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 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.
To model this treatment, we introduce the control w, which represents the dose of radiation and takes its values in . Following a model from Wein et al.  , the damage that is done to the tumor by radiation is modelled as follows
with the radiosensitive parameters depending on the treated tissue and the tissue repair rate . To simplify the expression above, we introduce the function
This is the solution to a linear ODE given by
Hence, the term that quantifies the damage done to the tumor can be written as follows
Now, we have to take into account that the radiation has also a damaging effect on the healthy tissues. Specifically, the damage on the carrying capacity of the vascularity q is given by
Notice that the radiosensitive parameters have different values, because malignant and healthy tissues have different characteristics.
Summarizing, our controlled model of cancer’s development and treatment is given by
This model is completely specified by giving the values of the parameters appearing in it. These values are specified in Table 1; see  and  .
2.3. 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. . See Schättler and Ledzewicz  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 L2- and L1-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.
We define our optimal control problem with anti-angiogenesis and radiotherapy as follows
Table 1. Model parameters.
with the initial conditions
and the Lebesgue measurable functions and take their values in and , respectively.
The parameters can be chosen differently to obtain different settings. We refer to this optimal control problem as the (OCAR3) control problem. Notice that in  a similar model is considered where only the tumor volume at the final time enters in the controls’ objective.
3. 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 and has the equilibrium
points and with .
From the second equation with , we obtain
Altogether, we have and . □
To show that the equilibrium point is locally asymptotically stable, we set , and obtain
For this system, we denote the equilibrium point with where and . 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 , we obtain
We have that the trace is and the determinant is
. Since the trace is the sum of the eigenvalues and the
determinant is the product of the eigenvalues of , we conclude and . That implies , hence the equilibrium point is locally asymptotically stable.
Next, we show that the solution to (6) with is always non-negative, 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 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 and neglect the equation for r since in this case r is not coupled to the -system. We have
The equation for q does not depend on p, hence q grows (independently of p) until (equilibrium point) is reached. As we start with initial values with , the tumor volume initially grows until , but since q is growing (independently of p), we will again have and the tumor continues to grow until (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 . Consider the following problem
Theorem 1 The set is positively invariant for the system (10).
Proof. We show that if and are arbitrary admissible controls defined on some interval , then the solution with the initial condition exists for all and .
Consider . Because of , the unique solution exists and is given by
On the boundary set , we have
thus the vector field points into .
Now, we look at the boundary set . Note that the nullclines for are given by
for the constant controls and . We rewrite the q-equation as follows
We obtain that for and for . is strictly increasing for fixed r with and
Therefore, for all p with , all constant controls , and for every , we have
Hence for trajectories starting in points with and .
On the coordinates axes for and the dynamics has singularities,
in the sense that is not defined. Therefore, we consider the lines
for . Now, it is sufficient to show that is positive for small as follows
for x small enough. On the other hand, is negative for large x
Notice that exists, because as we discuss below, our differential equation has an absolutely continuous solution. Hence, the region is positively invariant for system (10). □
Now, let us again look at the system (6) with . We take a closer look at the q-equation given by
For initial values in the solutions to this equation with or are positive and tend to zero as the equilibrium point is reached (see Lemma 1). Since is the linear combination for , it is also positive and tends to zero as the equilibrium is reached. The p- and 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  , Theorem 54 Proposition C.3.6) one can show existence and uniqueness of a solution to (6).
4. 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 to (6) are in
for all controls and . The set is weakly sequentially compact (see  Theorem 2.11). As our objective , is convex and continuous (see  Theorem 2.14) and thus weakly lower semicontinuous (see  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 .
2) A control set .
3) The controlled dynamical system
is given by the function , . 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 with values in the control set U.
Definition 1 The pair is called admissible if it is a solution to the differential Equation (11) and if for all .
The objective of the control is to minimize the following functional
Where is continuously differentiable and is continuous in , differentiable in x for fixed , and the derivative with respect to x is continuous.
Our optimal control problem is now given as follows
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
with and .
In our case we have and we can assume without loss of generality. This situation is called normal extremal lift; see  .
Theorem 2 (Pontryagin’s maximum principle) Let be admissible for the optimal control problem (12). If is optimal, then in every point in which is continuous, we have
where is the solution to the adjoint equation
For a proof see (  , 2.4.2).
By applying PMP to our optimal control problem (OCAR3), we obtain the following necessary conditions for an optimal solution.
Proposition 3 Let be an optimal solution to (OCAR3). Then in every point in which is continuous, we have that the Hamiltonian
at is minimized by in , where is the solution to the adjoint equation
with the terminal conditions
5. 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  . 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   . 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  where this issue is left open. On the other hand, we have the quadratic regularisation method proposed in  , 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  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
where and with . We remark that by increasing a sufficient descent of the cost functional in each iteration can be achieved, see (  , Lemma 4.1).
While we refer to  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 and .
Algorithm 1 (SQH method)
1) Choose , , , , , , compute from (10) corresponding to and from (14) corresponding to and , set .
2) Minimise pointwise
3) Calculate from (10) corresponding to and calculate
: Choose , , , calculate from (14) corresponding to and and .
5) If : STOP and return .
Else go to 2).
We remark that the update of is given by
To benchmark our novel method with a well-known solution scheme, we remark that the numerical approximation of an optimal control problem as (OCAR3) can interpreted as a nonlinear optimization problem that can be solved by using a nonlinear programming approach; see, e.g.,  . The way we choose to perform this task is by using the modeling language AMPL  and the solver IPOPT  . We refer to the resulting optimization framework as the IPOPT-AMPL solver.
6. Numerical Experiments
This section is devoted to the investigation of the effectiveness of our numerical optimization procedure. We present results of experiments with the (OCAR3) optimization problem solved by the SQH method given by Algorithm 1. Further, to assess the computational efficiency of our scheme, we compare results of simulation with those obtained with the IPOPT-AMPL solver.
For all numerical experiments of this section, we consider (unless otherwise stated) and set . The two-dimensional control, denoted by , takes values in the set with and . For the discretization of U we choose , . As initial guess for the controls, we take the zero function. Further, we initialize our state variables with , and . The parameters in Algorithm 1 are chosen as , , , and the initial guess .
For the first experiment setting, we choose additionally and the L1- and L2-penalty parameters are given by and , respectively. By using Algorithm 1 with this setting, we obtain the optimal solution showed in Figure 2. An almost identical solution is obtained with the IPOPT-AMPL solver. The tumor volume reduces to , that is, a reduction of more than two orders of magnitude of the initial volume. However, as the carrying capacity of the vasculature with is larger than the tumor will grow again after the treatment.
In Table 2, we further compare the results obtained with the SQH and IPOPT-AMPL schemes showing that the resulting values of the objective and the states at the terminal time T are similar. In both cases, the control w is a constant function taking the maximum value in the entire interval . However, the SQH method only needs a 1/50 of the CPU time needed by the IPOPT-AMPL scheme.
In Table 3, we consider a setting, where the parameters and are varying while the other parameters are kept fixed. We choose , the
Figure 2. Optimal solution obtained with the IPOPT-AMPL and SQH schemes for the first experimental setting.
Table 2. Comparison of the results of the SQH and IPOPT-AMPL methods for the first experimental setting.
Table 3. Results for increasing L1 penalty parameters.
rest of the setting remains as in our first example. We see, that the increase of the penalty parameters and results in higher values for p, q and J. This is reasonable since the increase of and describes higher costs for the controls.
In the second numerical experiment, we consider and the L1-penalty parameters . The L2-penalty parameters and the other parameters are the same as for the first experimental setting. With this setting and using Algorithm 1, we obtain the solution displayed in Figure 3.
Since the value of the control u at the beginning is bigger than in the first experimental setting, the states p and q are initially decreasing faster than before. The tumor volume reduces to , which is smaller than in the first
Figure 3. Optimal solution obtained by the IPOPT-AMPL and SQH schemes for the second experimental setting.
setting, but now the carrying capacity of the vasculature is larger with , so again the tumor will grow after the treatment. As above, IPOPT-AMPL provides a comparable solution and the results presented in Table 4 also confirm this fact. Again the computation time for the SQH method is less than the one needed by IPOPT.
In the experiments above, our focus was the comparison of the SQH method with the IPOPT-AMPL scheme. For this reason, less attention has been put on the ability of our optimal control formulation to deliver effective treatment. In the following experiments, we would like to show that it is indeed possible to find an optimisation setting that results in control functions that are able to reduce the volume and carrying capacity to zero at final time. In particular, we show the importance of the L1-cost towards this task.
Indeed, control costs of L1-type are considered in the literature for their ability to promote sparsity of controls. However, this feature is usually validated with a single control function, whereas in our case two control functions are considered that act on a nonlinear coupled system. In fact, as results of our experiment show, the choice of the weights of the costs of the control is a delicate issue.
For our experiments, we choose a time horizon (7 days), and the following control bounds , . The other parameters are chosen as follows: , , . We focus on the L1 weights while choosing for the L2 costs of the controls. This is our third experimental setting that is organised as follows.
In our first experiment of this series, we set and . The results of the SQH scheme with this setting are depicted in Figure 4. Notice that the volume and the carrying capacity become quickly zero. We also see that the control w (radiotherapy) is acting only in the first period of the treatment, while the anti-angiogenesis dose is u remains non-zero for the whole treatment. Therefore it seems natural to increase the weight with the purpose to force it become “sparse”, that is, to become zero for some interval of time. For this
Table 4. Comparison of results with the SQH and the IPOPT-AMPL schemes for the second experimental setting.
Figure 4. Third experimental setting: optimal solution with and .
reason, we keep all weights equal but slightly increase . The results of this experiment are shown in Figure 5. As expected, we can see a reduction of the control function u towards zero. However, it remains non-zero, and it has the effect that the control w is less sparse. Even more important, with this setting the carrying capacity remains non-zero. This result may suggest that we should increase both and , as we do in a third experiment, choosing and . The corresponding results are depicted in Figure 6. Notice that this result is qualitatively similar to the first one with and : we do not get a striking better result. Finally, we choose and . In this case, the L1 weights are large enough to get sparse control functions, but this is achieved at the cost of a non-zero and increasing carrying capacity; see Figure 7. Therefore our control framework allows to make the appropriate tuning of the optimisation weights in order to obtain promising successful treatments as with the first and third experiments.
A mathematical cancer therapy model was presented and investigated. This model resulted from the combination of two existing models for the simulation of cancer development and included two therapy mechanisms representing radiation and anti-angiogenesis inhibitors.
Figure 5. Third experimental setting: optimal solution with and .
Figure 6. Third experimental setting: optimal solution with and .
Figure 7. Third experimental setting: optimal solution with and .
To determine these therapies, an optimal control problem was formulated considering a cost functional including the tumor volume and L1- and L2-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.
We thank the Referee for helpful comments and pointers to the literature.
This publication was funded by the German Research Foundation (DFG) and the University of Würzburg in the funding programme Open Access Publishing.