Abstract

Keywords

Cancer, Radiotherapy, Anti-Angiogenesis, Sparse Controls, Optimal Control, Pontryagin’s Maximum Principle, SQH Method

Share and Cite:

1. 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 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., [1] [2] [3] . However, while these references discuss models for immunotherapy, the combination of radiation therapy and anti-angiogenic treatment as in [4] 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. [5] and the other one by Ergun et al. [6] . 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 [3] and [4] . In this paper, we investigate a more general cost functional that corresponds to minimizing the tumor volume along the entire time horizon and including L¹- and L²-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

$\stackrel{\dot{}}{p}=p\left(a-\xi \ln \left(p\right)\right),a>\xi >0.$

While the proliferation rate a of the cells is constant, the death rate $\xi \ln \left(p\right)$ 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

$q=\exp \left(\frac{a}{\xi }\right).$

Using this normalized carrying capacity, we obtain

$\stackrel{\dot{}}{p}=\xi p\left(\frac{a}{\xi }-\ln \left(p\right)\right)=\xi p\left[\ln \left(\exp \left(\frac{a}{\xi }\right)\right)-\ln \left(p\right)\right]=-\xi p\ln \left(\frac{p}{q}\right).$ (1)

For $p<q$ the tumor grows ( $\stackrel{\dot{}}{p}>0$ ) until $p=q$ . For $p>q$ the tumor shrinks ( $\stackrel{\dot{}}{p}<0$ ) 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

$\stackrel{\dot{}}{q}=S\left(p,q\right)-I\left(p,q\right).$

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

$\stackrel{\dot{}}{q}=bp-d{p}^{\frac{2}{3}}q,$ (2)

with the birth rate $b>0$ and the death rate $d>0$ . 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

$\stackrel{\dot{}}{q}=b{q}^{\frac{2}{3}}-d{q}^{\frac{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

$\stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right),$

where $\kappa \in \left[\mathrm{0,1}\right]$ . 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

$\stackrel{\dot{}}{p}=-\xi p\ln \left(\frac{p}{q}\right),\stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right).$ (4)

In Figure 1, we show the evolution of this system for initial values with $p\left(0\right)<q\left(0\right)$ and for different $\kappa$ . We see that the dynamics obtained with $\kappa =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\le \kappa \le 0.5$ and, choosing $\kappa$ close to one a slower increase of q can be observed. From this result, it appears that $\kappa =0.5$ and $\kappa =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] .

2.1. 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

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 $\left[\mathrm{0,}{u}_{max}\right]$ and represents the dose of the anti-angiogenic medicine. With the anti-angiogenic elimination parameter $\gamma >0$, we can augment the equation for q in (4) as follows

$\stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right)-\gamma qu.$

Hence, our model for an anti-angiogenetic mono-therapy is given by

$\stackrel{\dot{}}{p}=-\xi p\ln \left(\frac{p}{q}\right),\stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right)-\gamma qu.$ (5)

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.

2.2. 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.

To model this treatment, we introduce the control w, which represents the dose of radiation and takes its values in $\left[\mathrm{0,}{w}_{max}\right]$ . Following a model from Wein et al. [7] , the damage that is done to the tumor by radiation is modelled as follows

$-p\left(t\right)\left(\alpha +\beta {\displaystyle {\int }_0^t}w\left(s\right){\text{e}}^{-\rho \left(t-s\right)}\text{d}s\right)w\left(t\right)\mathrm{,}$

with the radiosensitive parameters $\alpha ,\beta >0$ depending on the treated tissue and the tissue repair rate $\rho >0$ . To simplify the expression above, we introduce the function

$r\left(t\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_0^t}w\left(s\right){\text{e}}^{-\rho \left(t-s\right)}\text{d}s\mathrm{.}$

This is the solution to a linear ODE given by

$\stackrel{\dot{}}{r}=-\rho r+w,\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}r\left(0\right)=0.$

Hence, the term that quantifies the damage done to the tumor can be written as follows

$-\left(\alpha +\beta r\right)pw\mathrm{.}$

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

$-\left(\eta +\delta r\right)qw\mathrm{.}$

Notice that the radiosensitive parameters $\eta ,\delta >0$ have different values, because malignant and healthy tissues have different characteristics.

Summarizing, our controlled model of cancer’s development and treatment is given by

$\begin{array}{l}\stackrel{\dot{}}{p}=-\xi p\ln \left(\frac{p}{q}\right)-\left(\alpha +\beta r\right)pw,\\ \stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right)-\gamma qu-\left(\eta +\delta r\right)qw,\\ \stackrel{\dot{}}{r}=-\rho r+w.\end{array}$ (6)

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] .

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. $p\left(T\right)$ . 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²- and L¹-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

$\begin{array}{l}\min J\left(\left(p,q,r\right),\left(u,w\right)\right)\\ :=\mathrm{<mi>\; </mi>}\frac{\sigma }{2}{\displaystyle {\int }_0^T}{p}^2\left(t\right)\text{d}t+\frac{\vartheta }{2}{p}^2\left(T\right)+\frac{{\nu }_u}{2}{\Vert u\Vert }_{L^2\left(0,T\right)}^2+{\mu }_u{\Vert u\Vert }_{L^1\left(0,T\right)}+\frac{{\nu }_w}{2}{\Vert w\Vert }_{L^2\left(0,T\right)}^2+{\mu }_w{\Vert w\Vert }_{L^1\left(0,T\right)}\end{array}$ (7)

subject to

$\begin{array}{l}\stackrel{\dot{}}{p}=-\xi p\ln \left(\frac{p}{q}\right)-\left(\alpha +\beta r\right)pw,\\ \stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right)-\gamma qu-\left(\eta +\delta r\right)qw,\\ \stackrel{\dot{}}{r}=-\rho r+w,\end{array}$ (OCAR3)

with the initial conditions

$p\left(0\right)=p_0,\mathrm{<mi>\; </mi>\; <mi>\; </mi>}q\left(0\right)=q_0,\mathrm{<mi>\; </mi>\; <mi>\; </mi>}r\left(0\right)=0,$

and the Lebesgue measurable functions $u(\cdot )$ and $w(\cdot )$ take their values in $\left[\mathrm{0,}{u}_{max}\right]$ and $\left[\mathrm{0,}{w}_{max}\right]$, respectively.

The parameters $\sigma \mathrm{,}\vartheta \mathrm{,}{\nu }_u\mathrm{,}{\mu }_u\mathrm{,}{\nu }_w\mathrm{,}{\mu }_w\ge 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.

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 $\kappa \in \left[\mathrm{0,1}\right]$ and $u=w=0$ has the equilibrium

points $\left({\bar{p}}_1,{\bar{q}}_1,{\bar{r}}_1\right)=\left(0,0,0\right)$ and $\left({\bar{p}}_2,{\bar{q}}_2,{\bar{r}}_2\right)=\left(\bar{p},\bar{q},0\right)$ with $\bar{p}=\bar{q}={\left(\frac{b}{d}\right)}^{\frac{3}{2}}$ .

Proof. Consider

$-\xi p\ln \left(\frac{p}{q}\right)=\mathrm{0\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\Rightarrow \mathrm{<mi>\; </mi>\; <mi>\; </mi>}p=\mathrm{0\; <mi>\; </mi>}\vee \mathrm{<mi>\; </mi>}p=q$

$\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right)=0$

$-\rho r=\mathrm{0\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\Rightarrow \mathrm{<mi>\; </mi>\; <mi>\; </mi>}r=0$

From the second equation with $q=p$, we obtain

$p^{\frac{2}{3}}\left(b-d{p}^{\frac{2}{3}}\right)=\mathrm{0\; <mi>\; </mi>}\vee \mathrm{<mi>\; </mi>}p\left(b-d{p}^{\frac{2}{3}}\right)=\mathrm{0\; <mi>\; </mi>\; <mi>\; </mi>}\Rightarrow \mathrm{<mi>\; </mi>\; <mi>\; </mi>}p=q=\mathrm{0\; <mi>\; </mi>}\vee \mathrm{<mi>\; </mi>}\left(b-d{p}^{\frac{2}{3}}\right)=0.$

Altogether, we have $p=q=\mathrm{0\; <mi>\; </mi>}\vee \mathrm{<mi>\; </mi>}p=q={\left(\frac{b}{d}\right)}^{\frac{3}{2}}$ and $r=0$ . □

To show that the equilibrium point $\left(\bar{p}\mathrm{,}\bar{q}\mathrm{,0}\right)$ is locally asymptotically stable, we set $z:=\ln p$, $y:=\ln q$ and obtain

$\begin{array}{l}\stackrel{\dot{}}{z}=-\xi \left(z-y\right)=:f_1\left(z,y\right),\\ \stackrel{\dot{}}{y}=b\exp \left(z-y\right)-d\exp \left(\frac{2z}{3}\right)=:f_2\left(z,y\right).\end{array}$ (8)

For this system, we denote the equilibrium point with $\left(\bar{z}\mathrm{,}\bar{y}\right)$ where $\bar{z}=\ln \bar{p}$ and $\bar{y}=\ln \bar{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

$A_f\left(z,y\right)=\left(\begin{array}{cc}\frac{\partial f_1}{\partial z}& \frac{\partial f_1}{\partial y}\\ \frac{\partial f_2}{\partial z}& \frac{\partial f_2}{\partial y}\end{array}\right)\left(z,y\right)=\left(\begin{array}{cc}-\xi & \xi \\ b{\text{e}}^{z-y}-\frac{2}{3}d{\text{e}}^{\frac{2z}{3}}& -b{\text{e}}^{z-y}\end{array}\right).$

At the equilibrium point $\left(\bar{z}\mathrm{,}\bar{y}\right)$, we obtain

$A_f\left(\bar{z}\mathrm{,}\bar{y}\right)=\left(\begin{array}{cc}-\xi & \xi \\ b-\frac{2}{3}d{\text{e}}^{\frac{2\bar{z}}{3}}& -b{\text{e}}^{z-y}\end{array}\right)\mathrm{.}$

We have that the trace is $\text{tr}{A}_f\left(\bar{z}\mathrm{,}\bar{y}\right)=-\xi -b<0$ and the determinant is

$det{A}_f\left(\bar{z},\bar{y}\right)=\frac{2}{3}d{\text{e}}^{\frac{2\bar{z}}{3}}>0$ . Since the trace is the sum of the eigenvalues and the

determinant is the product of the eigenvalues of $A_f$, we conclude ${\lambda }_1+{\lambda }_2<0$ and ${\lambda }_1{\lambda }_2>0$ . That implies ${\lambda }_1,{\lambda }_2<0$, hence the equilibrium point $\left(\bar{z}\mathrm{,}\bar{y}\right)$ is locally asymptotically stable.

Next, we show that the solution to (6) with $u=w=0$ 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 $\left(\bar{p}\mathrm{,}\bar{q}\right)$ cannot be reached according to the stability discussion above. On this account, we restrict our examination to the following domain

$\mathcal{D}=\left\{\left(p,q\right)\in {ℝ}^2|0<p<\bar{p},\mathrm{<mi>\; </mi>0}<q<\bar{q}\right\}=\left\{\left(p,q\right)\in {ℝ}^2|0<p,q<{\left(\frac{b}{a}\right)}^{\frac{3}{2}}\right\}.$ (9)

Next, we consider the uncontrolled case with $\kappa =1$ and neglect the equation for r since in this case r is not coupled to the $\left(p\mathrm{,}q\right)$ -system. We have

$\stackrel{\dot{}}{p}=-\xi p\ln (pq)$

$\stackrel{\dot{}}{q}=b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}=b{q}^{\frac{2}{3}}\left(1-\frac{d}{b}{\bar{q}}^{\frac{2}{3}}\right)=b{q}^{\frac{2}{3}}\left(1-{\left(\frac{q}{\bar{q}}\right)}^{\frac{2}{3}}\right).$

The equation for q does not depend on p, hence q grows (independently of p) until $q=\bar{q}$ (equilibrium point) is reached. As we start with initial values $\left(p_0\mathrm{,}{q}_0\right)\in \mathcal{D}$ with $p_0<q_0$, 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=\bar{p}=\bar{q}=q$ (equilibrium point) is reached. Hence in this case a solution that starts in $\mathcal{D}$ will never leave $\mathcal{D}$ .

Now, we show the same solution properties for the controlled case with $\kappa =0$ . Consider the following problem

$\begin{array}{l}\stackrel{\dot{}}{p}=-\xi p\ln \left(\frac{p}{q}\right)-\left(\alpha +\beta r\right)pw,\\ \stackrel{\dot{}}{q}=bp-d{p}^{\frac{2}{3}}q-\gamma qu-\left(\eta +\delta r\right)qw,\\ \stackrel{\dot{}}{r}=-\rho r+w.\end{array}$ (10)

Theorem 1 The set $\mathcal{D}$ is positively invariant for the system (10).

Proof. We show that if $\left(p_0\mathrm{,}{q}_0\right)\in \mathcal{D}$ and $u\mathrm{,}w$ are arbitrary admissible controls defined on some interval $\left[\mathrm{0,}T\right]$, then the solution $\left(p(\cdot )\mathrm{,}q(\cdot )\right)$ with the initial condition $\left(p\left(0\right),q\left(0\right)\right)=\left(p_0,q_0\right)$ exists for all $t\in \left[\mathrm{0,}T\right]$ and $\left(p\left(t\right)\mathrm{,}q\left(t\right)\right)\in \mathcal{D}$ .

Consider $\stackrel{\dot{}}{r}=-\rho r+w,\mathrm{<mi>\; </mi>}r\left(0\right)=0$ . Because of $w\in \left[\mathrm{0,}{w}_{max}\right]$, the unique solution exists and is given by

$r\left(t\right)={\displaystyle {\int }_0^t}w\left(s\right){\text{e}}^{-\rho \left(t-s\right)}\text{d}s\ge 0.$

On the boundary set $\left\{\left(p,q\right)\in \mathcal{D}:p=\bar{p},\mathrm{<mi>\; </mi>\; 0}<q<\bar{q}\right\}$, we have

$\stackrel{\dot{}}{p}=-\underset{>0}{\underset{︸}{\xi \bar{p}}}\ln \left(\underset{>1}{\underset{︸}{\frac{\bar{p}}{q}}}\right)-\left(\alpha +\beta r\right)\bar{p}w<0,$

thus the vector field points into $\mathcal{D}$ .

Now, we look at the boundary set $\left\{\left(p,q\right)\in \mathcal{D}:0<p<\bar{p},\mathrm{<mi>\; </mi>}q=\bar{q}\right\}$ . Note that the nullclines for $\stackrel{\dot{}}{q}=0$ are given by

$q=E_{\nu \omega }\left(p,r\right):=\frac{bp}{d{p}^{\frac{2}{3}}+\gamma \nu +\left(\eta +\delta r\right)\omega },$

for the constant controls $u=\nu$ and $w=\omega$ . We rewrite the q-equation as follows

$\begin{array}{c}\stackrel{\dot{}}{q}=bp-d{p}^{\frac{2}{3}}q-\gamma uq-\left(\eta +\delta r\right)qw\\ =\left(E_{\nu \omega }\left(p,r\right)-q\right)\underset{>0}{\underset{︸}{\left(dp-d{p}^{\frac{2}{3}}q-\gamma \nu q-\left(\eta +\delta r\right)q\omega \right)}}.\end{array}$

We obtain that $\stackrel{\dot{}}{q}>0$ for $q<E_{\nu \omega }\left(p,r\right)$ and $\stackrel{\dot{}}{q}<0$ for $q>E_{\nu \omega }\left(p,r\right)$ . $E_{\nu \omega }$ is strictly increasing for fixed r with $E_{\nu \omega }\left(0,r\right)=0$ and

$E_{\nu \omega }\left(\bar{p},r\right)=\frac{b}{b+\gamma \nu +\left(\eta +\delta r\right)\omega }\bar{p}\le \bar{p}.$

Therefore, for all p with $0<p<\bar{p}$, all constant controls $\nu \in \left[0,u_{\max }\right]$, $\omega \in \left[0,w_{\max }\right]$ and for every $r\ge 0$, we have

$E_{\nu \omega }\left(p,r\right)<E_{\nu \omega }\left(\bar{p},r\right)\le \bar{p}=\bar{q}.$

Hence $\stackrel{\dot{}}{q}<0$ for trajectories starting in points $\left(p\mathrm{,}q\right)$ with $0<p<\bar{p}$ and $q=\bar{q}$ .

On the coordinates axes for $p=0$ and $q=0$ the dynamics has singularities,

in the sense that $-\xi pln\left(\frac{p}{q}\right)$ is not defined. Therefore, we consider the lines

$p=xq$ for $x=0$ . Now, it is sufficient to show that $\stackrel{\dot{}}{x}$ is positive for small $x>0$ as follows

$\begin{array}{c}\stackrel{\dot{}}{x}=\frac{\text{d}}{\text{d}t}\left(\frac{p}{q}\right)=-\xi x\ln \left(x\right)-x\left(bx-d{p}^{\frac{2}{3}}-\gamma u-\left(\eta +\delta r\right)w\right)\\ >x\left(-\xi \ln \left(x\right)-bx\right)>0\end{array}$

for x small enough. On the other hand, $\stackrel{\dot{}}{x}$ is negative for large x

$\begin{array}{c}\stackrel{\dot{}}{x}<x\left(-bx+d{\bar{p}}^{\frac{2}{3}}+\gamma u_{\max }+\left(\eta +\delta \underset{r\in \left[0,T\right]}{\max }r\right)w_{\max }\right)\\ =bx\left(1+\frac{\gamma }{b}{u}_{\max }+\frac{\eta +\delta \underset{r\in \left[0,T\right]}{\max }r}{b}{w}_{\max }-x\right)<0\end{array}$

for

$x>1+\frac{\gamma }{b}{u}_{\max }+\frac{\eta +\delta \underset{r\in \left[0,T\right]}{\max }r}{b}{w}_{\max }.$

Notice that $ma{x}_{r\in \left[\mathrm{0,}T\right]}r$ exists, because as we discuss below, our differential equation has an absolutely continuous solution. Hence, the region $\mathcal{D}$ is positively invariant for system (10). □

Now, let us again look at the system (6) with $u,w=0$ . We take a closer look at the q-equation given by

$\stackrel{\dot{}}{q}=\kappa \left(b{q}^{\frac{2}{3}}-d{q}^{\frac{4}{3}}\right)+\left(1-\kappa \right)\left(bp-d{p}^{\frac{2}{3}}q\right).$

For initial values in $\mathcal{D}$ the solutions to this equation with $\kappa =0$ or $\kappa =1$ are positive and tend to zero as the equilibrium point $\left(\bar{p}\mathrm{,}\bar{q}\mathrm{,0}\right)$ is reached (see Lemma 1). Since $\stackrel{\dot{}}{q}$ is the linear combination for $0<\kappa <1$, it is also positive and tends to zero as the equilibrium is reached. The p- and r-equations do not depend on $\kappa$ and therefore we can deduce that for initial values $\left(p_0\mathrm{,}{q}_0\mathrm{,}{r}_0\right)\in \mathcal{D}\times \left\{0\right\}$ the solution cannot be negative for an arbitrary $\kappa \in \left[\mathrm{0,1}\right]$ .

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).

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 $x=\left(p,q,r\right)$ to (6) are in

$X=\left\{x\in L^2\left(0,T\right):0\le p\le \bar{p},\mathrm{<mi>\; </mi>\; 0}\le q\le \bar{q},\mathrm{<mi>\; </mi>\; 0}\le r\le \frac{w_{\max }}{\rho }\right\}$ for all controls $u\in U_{ad}\mathrm{<mo>\; :\; </mo>}=\left\{u\in L^2\left(\mathrm{0,}T\right)\mathrm{<mo>\; :\; </mo>}u\left(t\right)\in \left[\mathrm{0,}{u}_{max}\right]\right\}$ and $w\in W_{ad}:=\left\{w\in L^2\left(0,T\right):w\left(t\right)\in \left[0,w_{\max }\right]\right\}$ . The set $X\times U_{ad}\times W_{ad}$ is weakly sequentially compact (see [11] Theorem 2.11). As our objective $J\mathrm{<mo>\; :\; </mo>}X\times U_{ad}\times W_{ad}\to ℝ$, $\left(x\mathrm{,}u\mathrm{,}w\right)\mapsto J\left(x\mathrm{,}u\mathrm{,}w\right)$ 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 $M\subseteq {ℝ}^n$ .

2) A control set $U\subseteq {ℝ}^m$ .

3) The controlled dynamical system

$\stackrel{\dot{}}{x}=f\left(t,x,u\right),$ (11)

is given by the function $f\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}T\right]\times M\times U\to {ℝ}^n$, $\left(t\mathrm{,}x\mathrm{,}u\right)\mapsto f\left(t\mathrm{,}x\mathrm{,}u\right)$ . We assume

that the partial derivative $\frac{\partial f}{\partial x}\left(t\mathrm{,}x\mathrm{,}u\right)$ is continuous as a function of all variables.

4) The class $\mathcal{U}$ of admissible controls is taken to be a set of piecewise continuous functions u defined on a compact interval $\left[\mathrm{0,}T\right]\subseteq ℝ$ with values in the control set U.

Definition 1 The pair $\left(x(\cdot )\mathrm{,}u(\cdot )\right)$ is called admissible if it is a solution to the differential Equation (11) and if $u\left(t\right)\in U$ for all $t\in \left[\mathrm{0,}T\right]$ .

The objective of the control $u\in \mathcal{U}$ is to minimize the following functional

$J\left(x\mathrm{,}u\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_0^T}L\left(s\mathrm{,}x\left(s\right)\mathrm{,}u\left(s\right)\right)\text{d}s+g\left(x\left(T\right)\right)\mathrm{.}$

Where $g\mathrm{<mo>\; :\; </mo>}{ℝ}^n\to ℝ$ is continuously differentiable and $L\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}T\right]\times M\times U\to ℝ$ is continuous in $\left(t\mathrm{,}x\mathrm{,}u\right)$, differentiable in x for fixed $\left(t\mathrm{,}u\right)\in ℝ\times U$, and the derivative with respect to x is continuous.

Our optimal control problem is now given as follows

$\begin{array}{l}min\mathrm{<mi>\; </mi>}J\left(x\mathrm{,}u\right)={\displaystyle {\int }_0^T}L\left(s,x\left(s\right),u\left(s\right)\right)\text{d}s+g\left(x\left(T\right)\right)\\ \text{s}\text{.t}\text{.}\stackrel{\dot{}}{x}=f\left(t,x\left(t\right),u\left(t\right)\right),\mathrm{<mi>\; </mi>}x\left(0\right)=x_0,\\ \mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}u\left(t\right)\in U\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}\forall t\in \left[0,T\right].\end{array}$ (12)

Notice that the optimal control problem (OCAR3) that was introduced in Section 1 is of this form.

Definition 2 The Hamiltonian function $H\mathrm{<mo>\; :\; </mo>}ℝ\times {ℝ}^n\times {ℝ}^n\times {ℝ}^m\to ℝ$ for the optimal control problem (12) is defined as follows

$H\left(t\mathrm{,}x\mathrm{,}\lambda \mathrm{,}u\right)={\lambda }_{0}L\left(t\mathrm{,}x\mathrm{,}u\right)+{\lambda }^{\text{T}}f\left(t\mathrm{,}x\mathrm{,}u\right)\mathrm{,}$

with ${\lambda }_0\in ℝ$ and $\lambda \in {ℝ}^n$ .

In our case we have ${\lambda }_0>0$ and we can assume ${\lambda }_0=1$ without loss of generality. This situation is called normal extremal lift; see [13] .

Theorem 2 (Pontryagin’s maximum principle) Let $\left(x^{\mathrm{<mo>\; *\; </mo>}}(\cdot )\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}(\cdot )\right)$ be admissible for the optimal control problem (12). If $\left(x^{\mathrm{<mo>\; *\; </mo>}}(\cdot )\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}(\cdot )\right)$ is optimal, then in every point $t\in \left[\mathrm{0,}T\right]$ in which $u^{\mathrm{<mo>\; *\; </mo>}}(\cdot )$ is continuous, we have

$H\left(t\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{\lambda }^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)=\underset{u\in U}{min}H\left(t\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{\lambda }^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}u\right)\mathrm{,}$ (13)

where ${\lambda }^{\mathrm{<mo>\; *\; </mo>}}(\cdot )$ is the solution to the adjoint equation

${\stackrel{\dot{}}{\lambda }}^{\mathrm{<mo>\; *\; </mo>}}=-H_x\left(t\mathrm{,}{x}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{\lambda }^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\mathrm{,}{u}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\right)\mathrm{.}$

For a proof see ( [14] , 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 $\left(\left(p^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{q}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{r}^{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{,}\left(u^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{w}^{\mathrm{<mo>\; *\; </mo>}}\right)\right)$ be an optimal solution to (OCAR3). Then in every point $t\in \left[\mathrm{0,}T\right]$ in which $\left(u^{\mathrm{<mo>\; *\; </mo>}}(\cdot )\mathrm{,}{w}^{\mathrm{<mo>\; *\; </mo>}}(\cdot )\right)$ is continuous, we have that the Hamiltonian

$\begin{array}{l}H\left(t,\left(p,q,r\right),\lambda ,\left(u,w\right)\right)\\ =\frac{\sigma }{2}{p}^2+\frac{{\nu }_u}{2}{u}^2+\frac{{\nu }_w}{2}{w}^2+{\mu }_u\left|u\right|+{\mu }_w\left|w\right|-{\lambda }_1\left(\xi p\ln \left(\frac{p}{q}\right)+\left(\alpha +\beta r\right)pw\right)\\ +{\lambda }_2\left(\kappa M_1\left(q\right)+\left(1-\kappa \right)M_2\left(p,q\right)\right)-{\lambda }_2\left(\gamma qu+\left(\eta +\delta r\right)qw\right)+{\lambda }_3\left(-\rho r+w\right)\end{array}$

at $\left(t\mathrm{,}{p}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{q}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{r}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\lambda }^{\mathrm{<mo>\; *\; </mo>}}\right)$ is minimized by $\left(u^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{w}^{\mathrm{<mo>\; *\; </mo>}}\right)$ in $\left[\mathrm{0,}{u}_{\text{max}}\right]\times \left[\mathrm{0,}{w}_{\text{max}}\right]$, where ${\lambda }^{\mathrm{<mo>\; *\; </mo>}}(\cdot )$ is the solution to the adjoint equation

$\begin{array}{l}{\stackrel{\dot{}}{\lambda }}_1^{*}=-\sigma p^{*}+{\lambda }_1^{*}\left(\xi \ln \left(\frac{p^{*}}{q^{*}}\right)+\xi +\left(\alpha +\beta r^{*}\right)w^{*}\right)-{\lambda }_2^{*}\left(1-\kappa \right)\left(b-\frac{2}{3}d{\left(p^{*}\right)}^{-\frac{1}{3}}{q}^{*}\right)\\ {\stackrel{\dot{}}{\lambda }}_2^{*}=-{\lambda }_1^{*}\xi \frac{p^{*}}{q^{*}}-{\lambda }_2^{*}[\kappa \left(\frac{2}{3}b{\left(q^{*}\right)}^{-\frac{1}{3}}-\frac{4}{3}d{\left(q^{*}\right)}^{\frac{1}{3}}\right)\\ \mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}+\left(1-\kappa \right)\left(-d{\left(p^{*}\right)}^{\frac{2}{3}}\right)-\gamma u^{*}-\left(\eta +\delta r^{*}\right)w^{*}]\\ {\stackrel{\dot{}}{\lambda }}_3^{*}={\lambda }_1^{*}\beta p^{*}{w}^{*}+{\lambda }_2^{*}\delta q^{*}{w}^{*}+{\lambda }_3^{*}\rho \end{array}$ (14)

with the terminal conditions

${\lambda }_1^{*}\left(T\right)=\vartheta p^{*}(T)$

${\lambda }_2^{*}\left(T\right)=0$

${\lambda }_3^{*}\left(T\right)=0.$

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 [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

$\begin{array}{l}{K}_{ϵ}\left(t,\left(p,q,r\right),\lambda ,\left(u,w\right),\left(\hat{u},\hat{w}\right)\right)\\ :=H\left(t,\left(p,q,r\right),\lambda ,\left(u,w\right)\right)+ϵ\left({\left(u\left(t\right)-\hat{u}\left(t\right)\right)}^2+{\left(w\left(t\right)-\hat{w}\left(t\right)\right)}^2\right),\end{array}$

where $ϵ>0$ and $K_{ϵ}\mathrm{<mo>\; :\; </mo>}{ℝ}_0^{+}\times {ℝ}^3\times {ℝ}^3\times U\times U\to ℝ$ with $U=\left[\mathrm{0,}{u}_{max}\right]\times \left[\mathrm{0,}{w}_{max}\right]$ . 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 $x:=\left(p,q,r\right)$ and $v=\left(u,w\right)$ .

Algorithm 1 (SQH method)

1) Choose $ϵ>0$, $\hat{\kappa }>0$, $\hat{\sigma }>1$, $\hat{\zeta }\in \left(\mathrm{0,1}\right)$, $\hat{\eta }\in \left(\mathrm{0,}\infty \right)$, $v^0$, compute $x^0$ from (10) corresponding to $v^0$ and ${\lambda }^0$ from (14) corresponding to $x^0$ and $v^0$, set $k\leftarrow 0$ .

2) Minimise $K_{ϵ}$ pointwise

$\tilde{v}=\underset{v\in U}{\arg \min }{K}_{ϵ}\left(t,x^k,{\lambda }^k,v,v^k\right).$

3) Calculate $\tilde{x}$ from (10) corresponding to $\tilde{v}$ and calculate

$\hat{\tau }:={\Vert \tilde{u}-u^k\Vert }_{L^2\left(0,T\right)}^2+{\Vert \tilde{w}-w^k\Vert }_{L^2\left(0,T\right)}^2$ .

4) If

$J\left(\tilde{x}\mathrm{,}\tilde{v}\right)-J\left(x^k\mathrm{,}{v}^k\right)\mathrm{<mo>\; >\; </mo>}-\hat{\eta }\hat{\tau }$ : Choose $ϵ\leftarrow \hat{\sigma }ϵ$

Else if

$J\left(\tilde{x}\mathrm{,}\tilde{v}\right)-J\left(x^k\mathrm{,}{v}^k\right)\le -\hat{\eta }\hat{\tau }$ : Choose $ϵ\leftarrow \hat{\zeta }ϵ$, $x^{k+1}\leftarrow \tilde{x}$, $v^{k+1}\leftarrow \tilde{v}$, calculate ${\lambda }^{k+1}$ from (14) corresponding to $x^{k+1}$ and $v^{k+1}$ and $k\leftarrow k+1$ .

5) If $\hat{\tau }<\hat{\kappa }$ : STOP and return $v^k$ .

Else go to 2).

We remark that the update of $\tilde{v}=\left(\tilde{u},\tilde{w}\right)$ is given by

$\tilde{u}\left(t\right)=\min \left\{u_{\max },\max \left\{0,\frac{-{\mu }_u+{\lambda }_2\left(t\right)\gamma q\left(t\right)+2ϵu\left(t\right)}{{\nu }_u+2ϵ}\right\}\right\},$

$\tilde{w}\left(t\right)=\min \left\{w_{\max },\max \left\{0,\frac{-{\mu }_w+{\lambda }_1\left(t\right)\left(\alpha +\beta r\left(t\right)p\left(t\right)+{\lambda }_2\left(t\right)\left(\eta +\delta r\left(t\right)\right)q\left(t\right)-{\lambda }_3\left(t\right)+2ϵw\left(t\right)\right)}{{\nu }_w+2ϵ}\right\}\right\}.$

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., [20] . The way we choose to perform this task is by using the modeling language AMPL [21] and the solver IPOPT [22] . 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 $T=10$ (unless otherwise stated) and set $N=1000$ . The two-dimensional control, denoted by $\left(u\mathrm{,}w\right)$, takes values in the set $U=\left[0,u_{\max }\right]\times \left[0,w_{\max }\right]$ with $u_{\max }=15$ and $w_{\max }=1$ . For the discretization of U we choose ${\mathcal{N}}_u=50$, ${\mathcal{N}}_w=10$ . As initial guess for the controls, we take the zero function. Further, we initialize our state variables with $p_0=8000$, $q_0=10000$ and $r_0=0$ . The parameters in Algorithm 1 are chosen as $\hat{\kappa }={10}^{-7}$, $\hat{\sigma }=50$, $\hat{\zeta }=0.15$, $\hat{\eta }={10}^{-7}$ and the initial guess $ϵ=0.1$ .

For the first experiment setting, we choose additionally $\kappa =0.5$ and the L¹- and L²-penalty parameters are given by ${\mu }_u={\mu }_w=0$ and ${\nu }_u={\nu }_w=1$, 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 $p\left(T\right)=22.2847$, that is, a reduction of more than two orders of magnitude of the initial volume. However, as the carrying capacity of the vasculature with $q\left(T\right)=325.253$ is larger than $p\left(T\right)$ 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 $w_{\max }=1$ in the entire interval $\left[\mathrm{0,}T\right]$ . 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 ${\mu }_u$ and ${\mu }_w$ are varying while the other parameters are kept fixed. We choose ${\nu }_u={\nu }_w=0.01$, the

rest of the setting remains as in our first example. We see, that the increase of the penalty parameters ${\mu }_u$ and ${\mu }_w$ results in higher values for p, q and J. This is reasonable since the increase of ${\mu }_u$ and ${\mu }_w$ describes higher costs for the controls.

In the second numerical experiment, we consider $\kappa =1$ and the L¹-penalty parameters ${\mu }_u={\mu }_w=0.01$ . The L²-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 $p\left(T\right)=15.1133$, which is smaller than in the first

setting, but now the carrying capacity of the vasculature is larger with $q\left(T\right)=518.716$, 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 L¹-cost towards this task.

Indeed, control costs of L¹-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 $T=7$ (7 days), and the following control bounds $u_{\max }=10$, $w_{\max }=8$ . The other parameters are chosen as follows: $\sigma =1$, $\vartheta =1$, $\kappa =0.5$ . We focus on the L¹ weights while choosing ${\nu }_u={\nu }_w=0.01$ for the L² costs of the controls. This is our third experimental setting that is organised as follows.

In our first experiment of this series, we set ${\mu }_u=0.1$ and ${\mu }_w=1$ . 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 ${\mu }_u$ with the purpose to force it become “sparse”, that is, to become zero for some interval of time. For this

reason, we keep all weights equal but slightly increase ${\mu }_u=0.3$ . 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 ${\mu }_u$ and ${\mu }_w$, as we do in a third experiment, choosing ${\mu }_u=1$ and ${\mu }_w=6$ . The corresponding results are depicted in Figure 6. Notice that this result is qualitatively similar to the first one with ${\mu }_u=0.1$ and ${\mu }_w=1$ : we do not get a striking better result. Finally, we choose ${\mu }_u=2$ and ${\mu }_w=10$ . In this case, the L¹ 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.

7. Conclusions

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.

To determine these therapies, an optimal control problem was formulated considering a cost functional including the tumor volume and L¹- and L²-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.

Acknowledgements

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.

Conflicts of Interest

The authors declare no conflicts of interest.

References