1. Introduction

Cancer is one of the most common causes of death in industrialized countries. The complex process by which normal cells are transformed into cancer cells is called carcinogenesis and results from progressive abnormalities in the genetic material of the transformed cells. Malignant tumors (cancers) have a specific capacity to invade and destroy the underlying mesenchyme (local invasion). The tumor cells need nutrients via the bloodstream and produce a range of proteins that stimulate the growth of blood vessels into the tumor, thus allowing continuous growth to occur. For vascularisation to occur, the nearest vessel or capillary needs to become destabilised so that the endothelial cells lining the vessel can loosen from their neighbours, and migrate through the extracellular matrix towards the tumor. Only after a tumor has recruited its own blood supply can it expand in size [2]. Tumors do this via the production of angiogenic factors secreted into local tissues and stroma; this process has been termed the angiogenic switch. The new vessels are not well formed and are easily damaged so that the invading tumor cells may penetrate these and lymphatic vessels. Tumor fragments may be carried in these vessels to local lymph nodes or to distant organs where they may produce secondary tumors. This process of angiogenesis (blood vessel formation from the existing vascular network) is one of the hallmarks of cancer [20].

After observing these phenomena, Judah Folkman suggested the substantial potential of tumor angiogenesis as a therapeutic target [15]. Since in normal healthy adults the process of angiogenesis is very limited, it should, at least in theory, be possible to inhibit tumor angiogenesis without affecting normal tissues. Anti-angiogenic therapies have become one of the most promising approaches in anti-cancer drug development and successful preclinical research data are leading to clinical trials based on different strategies [14]. Approaches currently under evaluation for inhibiting angiogenesis may either be direct (targeting cell surface-bound proteins/receptors) or indirect (targeting growth factor molecules). Because angiogenesis is a complex process with multiple, sequential, and interdependent steps, this complexity creates many potential targets for inhibition. Therefore, an anti-angiogenic effect can be achieved by targeting angiogenic stimulators, angiogenic receptors, extracellular matrix proteins, extracellular matrix proteolysis, control mechanisms of angiogenesis, or the endothelial cells directly. Despite the fact that these approaches put forward innovative ideas for successful cancer treatment, at present there are a number of problems in clinical trials on humans that require very attentive studies and critical interpretations. Compounds that perform quite well in preclinical studies fail to give similar results in patients with cancer. There are more than a few reasons that can explain the presence of such differences between preclinical and clinical outcomes [7]. An important constrain in efficient anti-angiogenic therapy is the accessibility of tumors to anti-angiogenic agents. The genetic instability and high mutation rate of tumor cells is responsible, in part, for the frequent emergence of acquired drug resistance to conventional cytotoxic anticancer therapy (see e.g. [22]). However, vascular endothelial cells, like bone marrow cells, are genetically stable and have a low mutation rate. Unfortunately, contrary to hopes, that anti-angiogenic therapy would be a strategy to bypass drug resistance [21], two types of resistance have been observed. First, evasive resistance which includes revascularization as a result of upregulation of alternative pro-angiogenic signals, and second, intrinsic resistance which includes rapid adaptive responses observed by the absence of any beneficial effect of anti-angiogenic agents [1].

Nowadays anti-angiogenic therapy is considered rather as an essential component of multidrug cancer therapy ([17,43]) especially combined with chemotherapy. Although tumor eradication in such combined therapy may be still the primary goal, the chaotic structure of the angiogenically-created network leads to another target for anti-angiogenic agents, namely using angiogenic inhibitors to normalize the abnormal vasculature (the so-called pruning effect) and thus facilitate drug delivery [32], [11]. Continuous treatment with angiogenic inhibitors ultimately leads to a decrease in tumor blood flow and to a decreased tumor uptake of co-administrated cytotoxic drugs. In periodic therapy the main goal of anti-angiogenic agents is to normalize tumor vasculature. A number of anti-angiogenic clinical trials currently in progress have been designed to compare the effects of a particular cytotoxic agent alone with the effects of the same agent in combination with an angiogenesis inhibitor. These effects of combination therapy, which have also been observed for the combination of radiation therapy and angiogenesis inhibitors [43], could play a significant role in the clinical evaluation and effects of angiogenesis inhibitors. It is also worth mentioning that anti-angiogenic therapy was found to be efficient for slowly growing tumors, which are difficult to target by classical chemotherapy.

The administration of cytotoxic drugs often results in significant side effects which may reflect either the primary anti-proliferative action of the drug, some less well understood but predictable toxicological effect, or may be entirely idiosyncratic. Whereas side-effects of chemotherapy are already relatively well investigated after many years of application, we still do not know much about the side-effects of anti-angiogenic therapy. Anti-angiogenic agents do not require a very high dose to fulfil their function, but obvious possible complications might be related to menstruation, diabetes and wound healing and the long-term effects of therapy require attention.

Pharmacokinetic factors contribute towards mechanisms of resistance; for example, it is important to realize that for many anticancer drugs the administered form of the drug is not necessarily the active form. Generally, pharmacokinetic effects should be taken into account in scheduling anticancer drugs. While cytotoxic drugs mostly have a half-life time of about a few hours, the half life-time of anti-angiogenic agents may vary over a wide range, from 15 minutes (e.g. angiostatin) up to 20 days (bevacizumab), see e.g. [6,44]. Drug resistance may lead to lack of response at the time of treatment or, following an initial response, the tumor regrows. On regrowth, a decision may be made whether to repeat the same regimen or to switch to a second line therapy. This decision is usually based on the initial response to the drug and to the specific drug-free interval [7].

We concentrate here on the class of two-compartmental models proposed by Hahnfeldt et al. [19] with two control variables representing effects of two anticancer modalities and multiple delays introduced in these control variables to take into account pharmacokinetic/pharmacodynamic (PK/PD) effects and additional requirements resulting from clinical recommendation (for example, a delay in use of a cytotoxic agent sufficient for pruning vessels by an anti-angiogenic agent). Our study was inspired by recently reported results of clinical trials with two angiogenic inhibitors characterized by different half-lives combined with chemotoxic agents (see e.g. [45]).

The question which of the different goals of a combined therapy, mentioned before, could be reached in a finite treatment horizon could be answered, at least theoretically, by the analysis of the controllability of the dynamical systems used as models of the processes of tumor growth in the presence of vascularization. Controllability is a qualitative property of dynamical control systems and its meaning, roughly speaking, is the following property: a dynamical system is controllable if it is possible to steer it from an arbitrary initial state to an arbitrary final state using the set of admissible controls. In the existing literature there are many different definitions of controllability strongly depending on the class of dynamical control systems (see e.g. [23] and references therein). In the present paper, we consider constrained local controllability problems for second-order finite-dimensional semilinear stationary dynamical systems described by a set of two ordinary differential state equations with multiple delays in control variables. The local nature of the conditions requires to first drive the system to the neighbourhood of the desired final state. One way to achieve this is to design therapy protocols which are optimal in the sense of minimization of this neighbourhood. Thus an important part of this paper is devoted to the problem of optimization of control for the model discussed and to numerical experiments which demonstrate the efficiency of the methods.

---

2. Two-compartmental models of combined therapy and their properties

Hahnfeldt et al. model [19] is based on experimental data from anti-angiogenic therapy and non therapy trials of Lewis lung tumors in mice. Roughly speaking the main idea of this class of models is to incorporate the spatial aspects of the diffusion of factors that stimulate and inhibit angiogenesis into a non-spatial two-compartmental model for cancer cells and vascular endothelial cells. If $\, p$, denotes the size of the cancer cell population (proportional to a tumor volume) and $\, q\,$ a parameter describing the size of the vascular network, then such growth could be expressed by a Gompertz-type growth equation. A second equation describes vascular network growth, and includes stimulators of angiogenesis, inhibitory factors secreted by tumor cells and natural mortality of the endothelial cells (characterized by constant parameters $b$, $d$ and $\mu$ respectively). In this model $\, \xi\,$ denotes the proliferation ability of the cells. Inhibitory factors are proportional to $p^{2/3}$ (i.e. to a surface dimension), because they concentrate nearby the area of the active surface between the tumor and the vascular network. The effect of therapy in such models can be included in the form of control actions entering the system as multipliers in the bilinear terms. Since anti-angiogenic agents disturb directly only the vascular network, the control variable $u$ is present only in the first equation. The second variable $v$ related to chemotherapy appears in both equations [38]. The coefficients $\varphi, \eta, \gamma$ are non-negative constants (conversion factors) that relate the dosages of anti-angiogenic $u$ and cytotoxic $v$ agents ($\varphi$ is much greater than $\eta$).

$\dot p = - \xi\, p \ln \left( \frac{p}{q} \right) - \varphi\, v\, p,$

(1)

$\dot q = b\, p - d\, q \, p^{2/3} - \mu q - \gamma\, u\, q - \eta\, v\, q .$

(2)

Similar behavior could be obtained if Gompertz-type growth is substituted by a logistic type:

$\dot p = \xi\, p \left( 1-\frac{p}{q}\right) - \varphi\, v\, p .$

(3)

The modification of this model proposed in [9], also satisfies Hahnfeldt's suggestions described above with the only difference which may be called vascularity-based stimulation in contrast to tumor-based stimulation in the original Hahnfeldt model.

$\dot q = b\, q - d\, q\, p^{2/3} - \mu\, q - \gamma\, u\, q - \eta\, v\, q .$

(4)

Combining the models of tumor growth and the associated models of vascular network growth we obtain a set of two-compartmental models, the properties of which have been compared in [40]. All these models when uncontrolled (without therapy) have the same equilibrium point defined by the same values of both variables $p$ and $q$:

$p^* = q^* = ((b - \mu)/d)^{3/2}$

(5)

This equilibrium point is both locally and globally asymptotically stable [9].

For a constant dosage of antitumor drugs in combined therapy this result enables to find such continuous protocols which lead to asymptotic eradication of the vascular network and in turn of the tumor. In this case the values of $p$ and $q$ in equilibrium are not the same [42], but still they are closely related by a linear map. For example, a condition for constant doses of anti-angiogenic ($u_c$) and cytotoxic ($v_c$) drugs ensuring complete asymptotic removal of the tumor for model (1), (4) or (3), (4) is given by

$u_c + v_c \eta/ \gamma = b/ \gamma \Longrightarrow p^*, q^* \rightarrow 0 .$

Yet another simplification proposed in [13] also satisfies the assertions proposed in [19] but the dynamics of vessel carrying support is independent of the size of the tumor:

$\dot q = b\, q^{2/3} - d\, q^{4/3} - \gamma\, u\, q - \eta\, v\, q .$

(6)

This model does not contain the natural mortality factor. Although this term has been present in the previously discussed models, all simulation results presented by the authors are obtained for $\mu=0$. The reason is the relatively small impact of this term for the dynamics of the considered system. Thus to simplify analysis and to have the possibility of comparison of results for all the models mentioned above we will omit this term in our further theoretical consideration. Nevertheless, in numerical simulations presented in section 4 we introduce again a small non-zero $\mu$. The reason is that in section 4 we compare results for different tumor growth models (Gompertz type, logistic type, without delays, with delays) and different objective functionals but with one model of vascular network dynamics given by (2).

This leads to the simpler form of the equilibrium which is also relevant for the model, discussed in [13]:

$p^* = q^* = (b/d)^{3/2} .$

(7)

Constant or periodic therapies which ensure tumor eradication, discussed e.g. in [9], have an important drawback: they should be applied for a long therapy horizon. Realistic control problems related to combined anticancer therapy should be formulated as finite horizon control problems, and in [8] results of simulation for simple protocols of continuous and periodic therapy for finite treatment horizons are presented. Parameters proposed by Hahnfeldt et al. [19] were used in order to implement each model under similar conditions. In the periodic protocol anti-angiogenic treatment has been implemented as the starting therapy to take into account that the vascular network should be normalized before chemotherapy.

One way of checking, at least theoretically, whether there exist protocols enabling reachability of such different final targets to be reached is to find conditions for controllability of the models under discussion, Using logarithmic transformation of state variables we can transform the models presented above into semilinear equations. As mentioned before, for practical reasons, we omit the natural mortality factor represented by parameter $\mu$. Defining:

$\begin{array}{lll} x = \ln (p/p^*), \;\;\;\;& y = \ln (q/q^*), \;\;\;\;& x^*=y^*=0, \\ \tau = \xi t, \;\;\;\;& \vartheta = b / \xi, \;\;\;\;& \bar{\vartheta} = (b d)^{1/2}/\xi, \\ \sigma = \gamma / \xi, \;\;\;\;& \epsilon = \varphi/\xi, \;\;\;\;& \zeta = \eta / \xi, \\ x'= dx / d \tau, \;\;\;\;& y' = d y / d \tau, \;\;\;\;& \end{array}$

(8)

we are led to the following system for model (1), (2):

$x'(\tau) = y(\tau) - x(\tau) - \epsilon v(\tau),$

(9)

$y'(\tau) = \vartheta ( e^{x(\tau) - y(\tau)} - e^{(2/3)x(t)} ) - \sigma u(\tau) - \zeta v(\tau) .$

(10)

If the Gompertz-type growth of the tumor is substituted by the logistic-type one (3) then equation (9) has the form:

$x'(\tau) = 1 - e^{x(\tau)-y(\tau)} - \epsilon v(\tau),$

(11)

Similarly, if the dynamics of vasculature capacity is modeled by (4) as in [9] or by (6) as in [13], then (10) should be substituted by

$y'(\tau) = \vartheta ( 1- e^{(2/3)x(\tau)}) - \sigma u(\tau) - \zeta v(\tau),$

(12)

or

$y'(\tau) = \bar{\vartheta} ( e^{(-1/3)y(\tau)} - e^{(1/3)y(\tau)}) - \sigma u(\tau) - \zeta v(\tau),$

(13)

respectively.

Semilinear stationary finite-dimensional control systems are described by the following ordinary differential state equation, where from now the current time is denoted again by $t$ :

$\underline{x}'(t) = A\underline{x}(t) + F(\underline{x}(t)) + B\underline{u}(t)$

(14)

with initial conditions $\underline{x}(0)$, where the state $\underline{x}(t) \in \mathbb{R}^n$ and the control $\underline{u}(t) \in \mathbb{R}^m$, $A$ is a $n \times n$ dimensional constant matrix, and $B$ is a $n \times m$ dimensional constant matrix. Moreover, let us assume that the nonlinear mapping $F: X \rightarrow X$ is continuously differentiable near the origin such that $F(0)=0$, and let $X$ and $U$ denote state and control spaces, respectively. In practice, admissible controls are always required to satisfy certain additional constraints. We assume that the set of values of controls $U_c \subset U$ is a given closed and convex cone with nonempty interior and vertex at zero.

The associated linear dynamical equation for semilinear dynamical system (14) is defined as:

$\underline{z}'(t) = C \underline{z}(t) + B\underline{u}(t) , \underline{z}(0) = 0,$

(15)

where

$C = A + F\underline{_x}(0)$

is a $n \times n$ dimensional constant matrix. In the case of the models discussed above, the state vector is $\underline{x} = [x, y]^{tr}$, the control vector is $\underline{u} = [u, v]^{tr}$ $\underline{z}$ is the state of the associated linear system, and $F\underline{_x}$ is a Jacobian matrix. Here, and later on $(.)^{tr}$ denotes transposition of a vector or a matrix.

As we have already mentioned, one of the main difficulties in planning such combined therapies is related to the pharmacodynamic/pharmacokinetic (PK/PD) properties of the drugs. In [42] we have proposed to model these effects by including different time delays for different agents. The delays in the models may be introduced also to illustrate the idea of vessel pruning which demands the administration of chemotherapy with a delay with respect to anti-angiogenic agents. To include delays in controls we may modify equations (9) and (10). We consider delays in chemotherapy protocols, which is justified for example if we combine Sunitinib (angiogenic inhibitor) with Cisplatin, or in both types of agents if we combine two different anti-angiogenic agents e.g. Erlotinib (Tarceva) and Bevacizumab (Avastin) with Cisplatin or Paclitaxel.

Thus (14) should be substituted by

$\underline{x}'(t) = A\underline{x}(t) + F(\underline{x}(t)) + b_0\, u(t)+b_1\, v(t-h)$

(16)

or by

$\underline{x}'(t) = A\underline{x}(t) + F(\underline{x}(t)) + b_0\, u(t)+ b_1\, v(t-h)+ b_2\, u(t-h_1))$

(17)

and (15) is replaced by

$\underline{z}'(t) = C \underline{z}(t) + b_0 u(t) + b_1 v(t-h))$

(18)

or by

$\underline{z}'(t) = C \underline{z}(t) + b_0\, u(t) + b_1\, v(t-h)) + b_2\, u(t-h_1)$

(19)

Accordingly, including delay(s) in (10) leads to

$y'(t) = \vartheta\, (e^{x(t)-y(t)}- e^{(2/3)x(t)}) - \sigma \, u(t) -\zeta\, v(t-h))$

(20)

or

$y'(t) = \vartheta (e^{x(t)-y(t)}- e^{(2/3)x(t)} ) - \sigma u(t) -\zeta v(t-h) - \sigma_1 u(t-h_1),$

(21)

respectively. Equation (9) should be transformed to

$x'(t) = y(t) - x(t) - \epsilon\, v(t-h) .$

(22)

It is important to remember that for equations with delays the initial condition should contain not only conditions for $x$ and $y$ but also initial functions for $v$ in the interval $[-h, 0)$ and (possibly) for $u$ in $[-h_1, 0)$.

---

3. Local relative controllability of models with multiple delays

For the semilinear dynamical system (14), it is possible to define many different concepts of controllability. We shall focus our attention on the so-called constrained controllability in the time interval $[0, T]$. In order to do that, first of all let us introduce the notion of the attainable set at time $T>0$ from zero initial conditions, denoted shortly by $K_T(U_c)$ and defined as follows:

$K_T(U_c) = \{\, \underline{x} \in X \, | \, \underline{x} = \underline{x}(T, \underline{u}), \; \underline{u}(t) \in U_c \, \}$

(23)

where $\underline{x}(t, \underline{u})$ for $t > 0$ is the unique solution of the differential state equation (14) with zero initial conditions and a given admissible control. Under the assumptions stated for the nonlinear term $F$, such a solution always exists. Now, using the concept of the attainable set, we recall the well-known definitions of constrained controllability in $[0, T]$ for a semilinear dynamical system.

Definition 3.1. The dynamical system (14) is said to be $U_c$-locally controllable in $[0, T]$ if the attainable set $K_T(U_c)$ contains a neighborhood of zero in the space $X$.

Definition 3.2. The dynamical system (14) is said to be $U_c$-globally controllable in $[0, T]$ if $K_T(U_c) = X$.

The main result is the following sufficient condition for constrained local controllability of the semilinear dynamical system (14) which will be used to study controllability of the models of combined anticancer therapy.

Theorem 3.3. ([24]) Suppose that (ⅰ) $F(0) = 0$, (ⅱ) $U_c \in U$ is a closed and convex cone with vertex at zero, and (ⅲ) the associated linear control system (15) is $U_c$-globally controllable in $[0, T]$. Then the semilinear stationary dynamical control system (14) is $U_c$-locally controllable in $[0, T]$.

To verify the assumption about constrained global controllability of the linear time invariant dynamical system, we may use the following Theorem 3.4.

Theorem 3.4. ([3]) Suppose the set $U_c$ is a cone with vertex at zero and nonempty interior in the space $\mathbb{R}^m$. Then the associated linear dynamical control system (15) is $U_c$-globally controllable in $[0, T]$ if and only if

1. it is controllable without any constraints, i.e.,

$rank[B, CB, C^2B, ..., C^{n-1}B] = n,$

(24)

2. there is no real eigenvector $\underline{w} \in \mathbb{R}^n$of the matrix $C^{tr}$ satisfying the inequality

$\underline{w}^{tr} B\underline{u} \leq 0 \;\; \forall \;\underline{u} \in U_c.$

(25)

Since in the case of our models the control variables are bounded from above to drive the system in the neighbourhood of the origin we require additionally that no eigenvalue of $C$ has a positive real part [3]. Theorem 3.3 could be proved using the generalized open mapping theorem. Condition 2. in Theorem 3.4 could be also interpreted in the following way: For each real eigenvector $\underline{w} \in \mathbb{R}^n$ of the matrix $C^{tr}$ there exist controls $\underline{u} \in U_c$ such that $\underline{w}^{tr}B\underline{u}$ changes its sign. Moreover, for single input systems this condition is equivalent to the requirement that matrix $C$ has only complex eigenvalues (see Corollary from [3]).

In the case of systems with delays in control variables there exist more possible definitions of controllability which may be used. This variety is related to different understanding of the notion of state of a dynamical system in this case. The most frequently used are relative and absolute controllabilities. Attainable sets could be defined similarly as in (23). The main difference is that the set of admissible controls is a cone in the linear space $L^{\infty}([0, T], U_c)$. The definitions of local and global relative controllability are now the same as Definitions 3.1 and 3.2, respectively, , only the admissible controls should be interpreted differently.

Before formulating counterparts of Theorems 3.3 and 3.4 for the relative controllability of semilinear systems with delays we denote $D=[b_0, b_1, b_2]$. Then we have:

Theorem 3.5. ([24]) Suppose that (ⅰ) $F(0) = 0$, (ⅱ) $U_c \in U$ is a closed and convex cone with vertex at zero, (ⅲ) the associated linear control system (19) is $U_c$-globally relatively controllable in $[0, T]$. Then the semilinear stationary dynamical control system (17) is $U_c$-locally relatively controllable in $[0, T]$.

To verify the assumption (ⅲ) about constrained global relative controllability of the linear time invariant dynamical system., we may use the following Theorem 3.6.

Theorem 3.6. ([24]) Suppose the set $U_c$ is a cone with vertex at zero and nonempty interior in the space $\mathbb{R}^m$. Then the associated linear dynamical control system (19) is $U_c$-globally relatively controllable in $[0, T]$ with $T>\max \{h, h_1 \}$ if and only if

1. it is controllable without any constraints, i.e.,

$rank[D, CD, C^2D, ..., C^{n-1}D] = n,$

(26)

2.there is no real eigenvector $\underline{w} \in R^n$ of the matrix $C^{tr}$ satisfying the inequality

$\underline{w}^{tr}D\underline{u} \leq 0 \;\; \forall \; \underline{u} \in U_c.$

(27)

As before, in the case of bounded control coordinates we require that there are no eigenvalue of C with a positive real part.

The proofs of these theorems are based on the generalized open mapping theorems (see [24]). Relative controllability of the system guarantees that $\underline{x}$ reaches a desired point at the given final time $T$ but does not assure that it will stay at this point after this time even if the control will be equal to zero. This latter problem may be solved by absolute controllability which is related to a concept of a complete state of the systems with delays. Nevertheless, this type of controllability will not be discussed in this paper.

The constrained local controllability of the models of combined anticancer therapy presented in the previous section without delays and with a single delay was checked by us in [42]. Now we shortly recall these results and extend them for the case of multiple delays in control. To focus attention we present results for the Hahnfeldt et al. model. The admissible controls are assumed to be positive, hence the set of admissible controls is a positive cone $U_c$ in the space $\mathbb{R}^2$. Taking into account the general form of the semi-linear dynamic system (14) we have for model (9) and (10):

$A = \left[\begin{array}{cc}-1 & 1 \\ 0 & 0 \end{array} \right], \;\; F(x, y) = \left[\begin{array}{c} 0 \\ \vartheta (e^{x-y}-e^{(2/3)x}) \end{array} \right], \;\; B = \left[\begin{array}{cc} 0 &-\epsilon \\-\sigma &-\zeta \end{array} \right], \;\; \underline{x} = \left[\begin{array}{c} x \\ y \end{array} \right].$

Thus we have

$F(0, 0) = \left[\begin{array}{c} 0 \\ 0 \end{array} \right], \;\; F\underline{_x}(0, 0) = \left[\begin{array}{cc} 0& 0 \\ \vartheta/3 &-\vartheta \end{array} \right], \;\; C = A +F\underline{_x}(0, 0) = \left[\begin{array}{cc}-1& 1 \\ \vartheta/3 &-\vartheta \end{array} \right].$

(28)

It is worth to note that the associated linear system will be the same for both Gompertz-type (9) and logistic-type (11) growth equations. We use Theorem 3.4 presented in this section. The characteristic polynomial $P(s)$ for matrix $C^{tr}$ has the form

$P(s) = det(sI-C^{tr}) = det \left[\begin{array}{cc} s+1& \vartheta/3 \\ 1 & s+\vartheta \end{array} \right] = s^2 + s(1+\vartheta) + \frac{2}{3} \vartheta.$

Hence, we have $\Delta(\vartheta) = \vartheta^2 - \frac{2}{3} \vartheta + 1 > 0.$ This means that there are always two real eigenvalues leading to the conclusion that in the case of single input (i.e. monotherapy) the sufficient condition of local constrained controllability is not satisfied. For controllability verification in the case of two control variables (combined therapy) we have

$rank \left[B \, , \, CB \right] = 2 = n.$

The eigenvalues have the form

$s_1 = 0.5 \left(-1-\vartheta - \sqrt{\Delta(\vartheta)} \right)\, < 0, s_2 = 0.5 \left(-1-\vartheta + \sqrt{\Delta(\vartheta)} \right)\, < 0,$

and the corresponding real eigenvectors are

$\underline{w}_1 = \left[\begin{array}{c}-1 \\ (\vartheta + s_1)^{-1} \end{array} \right], \;\; \underline{w}_2 = \left[\begin{array}{c}-1 \\ (\vartheta + s_2)^{-1} \end{array} \right] .$

Thus we have

$\begin{array}{rcl} \underline{w}_1^{tr} B \underline{u} = -(\vartheta+s_1)^{-1} \sigma u + (\epsilon - (\vartheta+s_1)^{-1} \xi )v, \\ \underline{w}_2^{tr} B \underline{u} = -(\vartheta+s_2)^{-1} \sigma u + (\epsilon - (\vartheta+s_2)^{-1} \xi)v, \\ \end{array}$

(29)

We can check that there exists a combination of admissible controls such that the expressions (29) will change their signs. Therefore the sufficient condition of local constrained controllability is satisfied for the combined therapy. The conditions of local controllability do not change if we model cancer population growth by a logistic-type equation instead of the Gompertz-type because of the same linear approximation of both equations.

Now we can study the effect of time delays in control variables on the controllability conditions. If we limit our discussion to relative controllability we can use Theorems 3.5 and 3.6 to study conditions of local relative controllability of the model (20), (22). Since in this case matrix $C$ is the same as in the relative model without delays and $D=B$, the conditions are satisfied analogously as before and the conclusions are the same. The only difference is that the treatment time $T$ should be greater than the time delay $h$, because in the opposite case the system is governed by only one control variable and the controllability conditions are not satisfied. In the case of multiple delays (model (22), (21)), the use of Theorems 3.5 and 3.6 leads to the conclusion that in the case when $T>h_1>h$ or $h_1>T>h$ local relative controllability of the model is guaranteed, but if $T < h < h_1$ the sufficient conditions are not satisfied.

---

4. Optimization of treatment protocols

To our knowledge, [13] was the first paper in which optimal protocols for combined anticancer therapy (anti-angiogenic agents combined with radiotherapy) were discussed. The authors used a simplified model of angiogenesis (see section 2) combined with an LQ model to measure the effects of radiation therapy. Though optimization methods were not applied systematically, the findings in [13] suggest that optimal strategies for free final time combine bang-bang and singular controls. Ledzewicz and Schättler [25] presented a complete solution in the form of an optimal synthesis for the control problem related to anti-angiogenic therapy for this model, and obtained a similar optimal strategy containing singular arcs for the original Hahnfeldt et al. model [26].

Meanwhile, different results are obtained in [41] for the D'Onofrio-Gandolfi model (see section 2) in the case of a fixed time of anti-angiogenic therapy. The most important conclusion is that intermediate drug doses (singular arcs) are not optimal and that the optimal protocol switches between maximal dose and no drug intervals (bang-bang control). Singular arcs are not feasible, since there are no finite intervals of constant solutions to the adjoint equations. Similar properties were found for the Hahnfeldt et al. model with logistic tumor growth [40]. In [28] the broad class of models from this family was analysed and the results from [25,26,40,41] were confirmed as special cases. Suboptimal strategies for the original Hahnfeldt et al. model for minimization of tumor volume with anti-angiogenic therapy using bang-bang optimal controls were described in [27]. Simple suboptimal protocols for models with and without a linear pharmacokinetic equation are presented in [29] and [30]. The big advantage is that these protocols realize tumor volume dynamics close to the optimal ones. Similar research including optimal singular arcs is described in [31]. For piecewise constant dosage protocols, a very good approximation to optimal solutions may be obtained. However, small doses have no significant effect on tumor development, but on the other hand a too high dosage is not efficient enough to justify its enormous cumulative cost. Preliminary results about optimal controls for a mathematical model that combines anti-angiogenic therapy with a chemotherapeutic killing agent were presented in [39] and [38] for the D'Onofrio-Gandolfi model and the fixed treatment horizon. Once more, the optimal strategy had no singular arcs. Moreover, similarly as in [40] the control objective took into account not only the final size of the tumor but also the final size of the vascular network. In [12] a problem of synthesis of optimal controls for the same family of models as in [28] is discussed, and the multicontrol problem for the original Hahnfeldt model with free treatment time was solved numerically. In [31] optimal and suboptimal protocols for this class of problems were compared.

In this section, we present optimization results for the Hahnfeldt et al. model using an objective function that balances the terminal values of the tumor and the vasculature. Control delays can occur in both control variables. In the non-delayed case and for free terminal time we compare numerical results for a Gompertz-type growth function (see [31]) with those obtained for a logistic-type growth function. For fixed terminal time we obtain a strong numerical control chattering when the Gompertz growth function is used. Hence, for fixed terminal time we use only the logistic growth function and present numerical results both in the non-delayed and the delayed case.

The computations are based on a two-stage procedure. First, the optimal control problem is discretized which results in a large-scale non-linear programming (NLP) problem. This NLP problem can be conveniently formulated with the help of the Applied Modeling Programming Language AMPL created by Fourer et al. [16]. AMPL is linked to the Interior-Point optimization solver IPOPT developed by Wächter and Biegler [46]. We use $5000-10000$ grid points and the trapezoidal rule as integration method. Alternatively, the control package NUDOCCCS implemented by Büskens [4] (cf. also [5]) provides a highly efficient method for solving the discretized control problem, because it allows to implement higher order integration methods. In a second step, the switching times are optimized directly using the arc-parametrization method in combination with NUDOCCCS. This method requires that a singular control can be determined in feedback form. The resulting optimization problem is called the Induced Optimization Problem (IOP); cf. [33,34].

---

4.1. Optimal control problem for the Hahnfeldt model with control delays

Now it will be more convenient to return to the model in the form (1), (2), (3). More precisely, we consider the dynamical system with the state variables $p$ (tumor volume) and $q$ (vasculature), the control variables $u$ (anti-angiogenic agent) and $v$ (chemotoxic agent), and control delays $h_1 \geq 0$ in $u$ and $h \geq 0$ in $v$:

$\begin{array}{l} \dot p(t) = f(p(t), q(t)) - \varphi\, p(t)\, v(t-h), \\ \dot q(t) = b\, p(t) - q(t)\, (d\, p(t)^{2/3} + \mu + \gamma_1\, u(t) + \gamma_2\, u(t-h_1) + \eta \, v(t-h) )\, ). \end{array}$

(30)

Initial conditions for $p$ and $q$ are given by

$p(0) =p_0, \;\;\;\; q(0) = q_0 .$

(31)

Due to the delays in the control variables $u$ and $v$, initial control functions have to be prescribed:

$\begin{array}{llllll} u(\tau) = 0 \;\; & \;\; \mbox{ for } \;\;& \;\;-h_1 \leq \tau < 0 \, ;\;\; &\;\; v(\tau) = 0 \;\; & \;\; \mbox{ for } \;\;&\;\; -h \leq \tau < 0 \, . \end{array}$

(32)

The initial conditions for the controls are void in the case $h=h_1=0$ of no delays. We consider two types of growth functions:

Gompertz growth function  $f(p, q) = - \xi \, p\, \ln (p/q)$ (as in (1)),

Logistic growth function  $f(p, q) = \xi \, p\, (1-p/q)$ (as in (3)).

To measure the total amount of control agents used, we introduce two artificial state variables $w$ and $z$ which are defined by the differential equations

$\begin{array}{rl} \dot w(t) = u(t), & \;\;\; w(0) = 0, \\ \dot z(t) = v(t), & \;\;\; z(0) = 0, \end{array}$

(33)

and prescribe the control constraints

$\begin{array}{lll} 0 \leq u(t) \leq u_{\rm max}, &\;\; 0 \leq v(t) \leq v_{\rm max} & \;\;\forall \;\, 0 \leq t \leq T, \\ \int_0^T u(t) dt = w(T) \leq w_{\rm max} \, , & \;\;\int_0^T v(t) dt = z(T) \leq z_{\rm max}\, . \end{array}$

(34)

The objective function is the weighted sum of the terminal values $p(T)$ and $q(T)$,

$J(u) = p(T) + \alpha \, q(T) \;\;( \alpha \geq 0 ).$

(35)

Then the optimal control problem $[OC]$ consists in determining measurable (piecewise continuous) control functions $u, v: [0, T] \rightarrow \mathbb{R}$ that minimize the objective (35) subject to the dynamical equations (30), (33), initial conditions (31), (32) and the control constraints (34).

In [31] only the case $\alpha=0$ in the objective (35) was treated. The numerical results in Sections 4.2 and 4.3 show that the terminal value $q(T)$ can be considerably reduced by choosing appropriate positive weights $\alpha > 0$. Essentially, we shall take the parameters from [31], but double the value of $\varphi$, choose a nonzero value for $\mu$ and adopt new parameters associated with the delayed control variables:

$\begin{array}{lllll} \xi = 0.084, &\;\; b = 5.85, &\;\; d = 0.00873, & \;\;\gamma_1 = 0.15, \\ \gamma_2 = 0.1, &\;\; \varphi = 0.2, &\;\; \eta = 0.05, &\;\;\mu = 0.02 . \end{array}$

(36)

In the non-delayed case with $h = h_1=0$ we set $\gamma_2 = 0$. Note that the control $v$ enters only as a delayed control $v(t-h)$. Thus we can expect that the optimal control will shift the non-delayed control $v(t)$ to the left by $h$ time units (Figure 5).

Let us briefly discuss the necessary optimality conditions of a Maximum Principle as they were recently derived in [18] for optimal control problems with multiple control and state delays. We denote by $u_d$ the delayed control variable with $u_d(t) = u(t-h_1)$ and by $v_d$ the delayed control variable with $v_d(t) = v(t-h)$. Denoting the state vector by $\underline{x}=(p, q, w, z) \in \mathbb{R}^4$ and the adjoint variable by $\lambda=(\lambda_p, \lambda_q, \lambda_w, \lambda_z) \in \mathbb{R}^4$, the Hamiltonian of the delayed control problem is given by

$\begin{array}{rcl} H(\underline{x}, \lambda, u, v, u_d, v_d) &=& \lambda_p ( f(p, q) - \varphi\, p\, v_d ) + \lambda_w u + \lambda_z v \\ && + \lambda_q ( b\, p - q\, (d\, p^{2/3} + \gamma_1\, u + \gamma_2\, u_d + \eta \, v_d ) ). \end{array}$

(37)

Since there is no delay in the state variables, the adjoint equations $\dot \lambda(t) = - H\underline{_x}$ do not contain the advanced time:

$\begin{array}{rcl} \dot{\lambda}_p(t) &=& - \lambda_p(t)( f_p(p(t), q(t)) - \varphi v(t-h) ) - \lambda_q(t)(b- \frac{2}{3} q(t)\, d\, p(t)^{-1/3} ), \\ \dot{\lambda}_q(t) &=& - \lambda_p(t) f_q(p(t), q(t)) \\ && + \lambda_q(t)\, (d\, p(t)^{2/3} + \mu + \gamma_1 u(t) + \gamma_2 u(t-h_1) + \eta v(t-h)), \\ \dot{\lambda}_w(t) &=& 0, \\ \dot{\lambda}_z(t) &=& 0. \end{array}$

(38)

Here, the subscripts $p$ and $q$ in $f$ denote partial derivatives of $f(p, q)$. In view of the objective (35) and constraints (34), the transversality conditions for the adjoint variable are

$\lambda_p(T) = 1, \, \lambda_q(T) = 0, \; \lambda_w(T)(w(T)-w_{\rm max}) = 0, \, \lambda_z(T)(z(T)-z_{\rm max}) = 0.$

(39)

Our computations show that $w(T)=w_{\rm max}$ and $z(T)=z_{\rm max}$ always holds for the chosen data. Then the last two equations in (39) mean that $\lambda_w(T)$ and $\lambda_z(T)$ are undetermined. The optimal control $u(t)$ minimizes the sum of Hamiltonians (cf. [18])

$\begin{array}{l} H(\underline{x}(t), \lambda(t), u, v(t), u(t-h_1), v(t-h)) \\ + \chi_{[0, T-h_1]} (t+h_1)\, H(\underline{x}(t+h_1), \lambda(t+h_1), u(t+h_1), v(t+h_1), u, v(t-h+h_1)) \end{array}$

with respect to $u \in [0, u_{\rm max}]$, where $\chi_{[0, T-h_1]}$ denotes the characteristic function of the interval $[0, T-h_1]$. Likewise, the optimal control $v(t)$ minimizes the sum

$\begin{array}{l} H(\underline{x}(t), \lambda(t), u(t), v, u(t-h_1), v(t-h)) \\ + \chi_{[0, T-h]}(t+h)\, H(\underline{x}(t+h), \lambda(t+h), u(t+h), v(t+h), u(t-h_1+h), v) \end{array}$

with respect to $v \in [0, v_{\rm max}]$. Since both controls appear linearly in the Hamiltonian, we are led to the switching functions

$\begin{array}{l} \phi_u(t) = -\lambda_q(t) \gamma_1 q(t) + \lambda_w(t) - \chi_{[0, T-h_1]} (t+h_1)\lambda_q(t+h_1)\, q(t+h_1) \gamma_2, \\ \phi_v(t) = \lambda_z(t) - \chi_{[0, T-h]}(t+h) ( \lambda_p(t+h) p(t+h) \varphi + \lambda_q(t+h) q(t+h) \eta ), \end{array}$

(40)

which determine the minimizing controls by the control law

$\underline{u}(t) = \left\{ \begin{array}{lcl} 0, & \mbox{if} & \phi\underline{_u}(t) > 0 \\ {\rm singular}, & \mbox{if} & \phi\underline{_u}(t) = 0 \;\; \forall \; t \in I_s \subset [0, T] \\ \underline{u}_{\rm max}, & \mbox{if} & \phi\underline{_u}(t) < 0 \end{array} \right\}, \underline{u} \in \{u, v\}.$

(41)

In the next sections, the sign conditions in this control law will be checked numerically for all computed solutions.

---

4.2. Optimal control solution without delays for free final time $T$.

In the following, we shall compare solutions for the Gompertz-type growth function [31] with solutions for the logistic growth function. The following initial conditions and control bounds for $u$ and $v$ are taken from [31]:

$p(0) = 12000, \; q(0) = 15000, \;\; u_{\rm max} = 75, \; w_{\rm max} = 300, \;\; v_{\rm max} = 2, \; z_{\rm max} = 10.$

---

4.2.1. Optimal control for Gompertz-type growth function $f(p, q) = - \xi \, p\, ln (p/q)$

In this section, we minimze the objective $J(u) = p(T)$ with $\alpha=0$ in the objective (35). Since the terminal time $T$ is free, the singular control $u$ can be obtained in feedback form [12],

$u = u_{\rm sing}(p, q) = \frac{1}{\gamma}\left( \, \xi \ln\left(\frac{p}{q}\right) + b \, \frac{p}{q} + \frac{2}{3} \xi \, \frac{d}{b}\, \frac{q}{p^{1/3}} - (\mu + d\, p^{2/3} ) \right) + \frac{\varphi - \eta}{\gamma}\, v,$

(42)

provided that the control $v(t)$ is piecewise constant on a singular arc of $u$.

The optimal controls have the following structure:

$(u(t), v(t)) = \left\{ \begin{array}{lcl} (u_{\rm max}, 0) & \mbox{for} & 0 \leq t < t_1 \\ (u_{\rm sing}(p(t), q(t)), 0) & \mbox{for} & t_1 \leq t < t_2 \\ (u_{\rm sing}(p(t), q(t)), v_{\rm max}) & \mbox{for} & t_2 \leq t \leq t_3 \\ (0, v_{\rm max}) & \mbox{for} & t_3 < t \leq T \end{array} \right\}.$

(43)

The control $u(t)$ is a bang-singular-bang control with switches at $t_1$ and $t_3$, whereas $v(t)$ is a bang-bang control with only one switch at $t_2$. Using this control structure with preassigned control values, we can directly optimize the switching times $t_1, t_2, t_3$ and the terminal time $T$ to minimize the terminal value $p(T)$ subject to the constraints $w(T) = w_{\rm max} = 300$ and $z(T) = z_{\rm max} = 10$. We solve this Induced Optimization Problem (IOP) by the arc-parametrization method [33] and the control package NUDOCCCS and obtain the following numerical results:

$\begin{array}{rclrclrcl} p(T) &=& 1246.00, \;\; & q(T) & = & 1700.56, \;\;& T &=& 5.5415, \\ t_1 &=& 0.090502, \;\; & t_2 &=& 0.54153, \;\; & t_3 &=& 5.3806. \end{array}$

The numerical results differ considerably from those in [31], since here we have $\varphi = 0.2$ and $\eta =0.05$ in contrast to $\varphi = 0.1$ and $\eta = 0$ in [31]. It is noteworthy that second-order sufficient conditions are satisfied for the IOP, since the projected Hessian $PH$ of the Lagrangian is computed as the positive definite $2\times 2$ matrix

$PH = \left( \begin{array}{cc} 99.40 & -89.07\\ -89.07 & 91.72 \end{array} \right).$

The optimal controls $(u, v)$ and the state variables $(p, q)$ are shown in Figure 1. Note that in view of (43) the control $u(t)$ has a positive jump at the switching time $t_2 = 0.54153$ of the control $v(t)$; cf. Figure 1(a). By applying the methods of synthesis analysis in [37], it can be shown that the control and state trajectories in Figure 1 provide a strict strong minimum. As in [31], the chemotherapeutic agent $v(t)$ should not be administered ab initio but only in a terminal part of the treatment period. In the literature, this effect is commonly referred to as "pruning".

It is of practical interest that the bang-singular-bang control (43) can be approximated by the following simpler control with piecewise constant values,

$(u(t), v(t)) = \left\{ \begin{array}{lcl} (u_{\rm max}, 0) & \mbox{for} & 0 \leq t < t_1 \\ (u_c, 0) & \mbox{for} & t_1 \leq t < t_2 \\ (u_c, v_{\rm max}) & \mbox{for} & t_2 \leq t \leq t_3 \\ (0, v_{\rm max}) & \mbox{for} & t_3 < t < T \end{array} \right\},$

(44)

where the switching times and final time are fixed to

$t_1 = 0.1, \;\; t_2 = 0.5, \;\; t_3 = 5.0, \;\; T = 5.5 .$

The constant value $u(t)=u_c$ in the interval $[t_1, t_3]$ is treated as an optimization variable. The arc-parametrization method [33] and the control package NUDOCCCS furnish the following numerical results:

$p(T) = 1265.49, \;\; q(T) = 3517.13, \;\; u_c = 59.6939.$

Though (44) is a crude approximation of the optimal control (43), the functional value $p(T) = 1265.49$ is only slightly bigger than the optimal value $p(T) = 1246.00$.

---

4.2.2. Optimal control for logistic-type growth function $f(p, q) = \xi \, p\, (1-p/q)$

Since the terminal time $T$ is free, we could obtain a formula for the singular control in feedback form. However, for the logistic growth function we only obtain bang-bang controls as was predicted in [38].

Solving the discretized control problem with $N= 10000$ grid points we find $v(t) \equiv 2$ and the following bang-bang control:

$\label{control-T-free-logistic-vmax=2-zmax=10} u(t) = \left\{ \begin{array}{lcl} 0 & \mbox{for} & 0 \leq t < t_1 \\ u_{\rm max} & \mbox{for} & t_1 \leq t \leq t_2 \\ 0 & \mbox{for} & t_2 < t \leq T \end{array} \right\} .$

(45)

The IOP with respect to the switching times $t_1, t_2$ and the terminal time $T$ yields the solution

$\begin{array}{rclrclrcl} p(T) &=& 1112.45, \;\; & q(T) & = & 3800.32, & && \\ t_1 &=& 0.31499, \;\; & t_2 &=& 4.31499. \;\; & T &=& 5.0. \end{array}$

The optimal control $u$ and the state variables $(p, q)$ are shown in Figure 2. The optimal value $p(T) = 1112.45$ is smaller than the value $p(T) = 1246.00$ for the Gompertz growth function. This is due to the fact that the logistic growth function produces a larger decrease in the tumor volume $p(t)$. However, the optimal terminal value $q(T) = 3800.32$ is much larger than the value $q(T)= 1700.56$ for the Gompertz growth function.

The solution in Figure 2 satisfies theSSC for bang-bang controls in [34], Chapter 7. Namely, SSC hold for the IOP, since the projected Hessian of the Lagrangian is the positive number $15.02$. Moreover, the following strict bang-bang property holds:

$\begin{array}{lllll} \phi_u(t) > 0 & \forall \;\; 0 \leq t < t_1, & \dot\phi_u(t_1) < 0, & \phi_u(t) < 0 & \forall \;\; t_1 < t < t_2, \\ \dot\phi_u(t_2) > 0, & \phi_u(t) > 0 & \forall \;\; t_2 < t \leq T, & \phi_v(t) < 0 & \forall \;\; 0 \leq t \leq T . \end{array}$

(46)

The strict bang-bang property for the control $u$ can be seen in Figure 2 (a). In order to reduce the rather high value $q(T) =3800.32$ we consider the objective $J(u) = p(T) + 0.5 q(T)$. Again, we find $v(t) \equiv 2$, while $u(t)$ is a bang-bang control with only one switch at $t_1$; cf. Figure (3). The numerical results are

$p(T)= 1134.06, \; q(T) = 559.412, \;\; t_1 = 1, \; T = 5.$

It is remarkable that the terminal value of $q$ is reduced from $q(T)=3800.32$ to $q(T) = 559.412$, whereas $p(T) = 1134.06$ is only slightly increased versus $p(T) = 1112.45$.

The numerical results show that theSSC [33,34] for bang-bang controls are satisfied.

---

4.3. Optimal control solution for fixed final time$T = 16$

We have chosen such a control horizon because we want to include control delays and compare with results without delay. Since the maximal delay is greater than 10 and the computations for the model without delays, and free final time lead to control horizons between 5 and 6, such choice is justified.

In this section, we consider only the logistic growth function $f(p, q) = \xi p (1-p/q),$ because the discretization method for the Gompertz growth function produces a strong numerical control chattering, from which it is hard to detect the correct control structure. In view of the rather large delay $h_1=10.6$ in the control $u$, which will be considered later, we have chosen the final time $T=16 > h_1$. The initial conditions

$p(0) = 12000, \;\; q(0) = 15000$

are as in the previous sections, but we choose different control bounds to accommodate to the much larger time horizon:

$u_{\max} = 40, \; w_{\max} = 400, \;\; v_{\max} = 2, \; z_{\max} = 10.$

Here, the ratio $w_{\max} / u_{\max} = 10$ roughly corresponds to the ratio $w_{\max} / u_{\max} = 4$ in section 4.2 which is multiplied by the ratio $\approx 2.5$ of terminal times.

---

4.3.1. Optimal control solution without delays

We are interested in a simple control structure which also produces a rather small terminal value $q(T)$. For that reason, we choose the objective (35) with $\alpha = 0.2$ and thus we minimize

$J(u) = p(T) + 0.2 \, q(T).$

The discretization approach yields the following control structure:

$(u(t), v(t)) = \left\{ \begin{array}{lcl} (0, 0) & \mbox{for} & 0 \leq t < t_1 \\ (u_{\rm max}, 0) & \mbox{for} & t_1 \leq t < t_2 \\ (u_{\rm max}, v_{\rm max}) & \mbox{for} & t_2\leq t \leq T \end{array} \right\}.$

(47)

and the numerical results

$p(T) = 1130.06, \; q(T) = 987.466, \;\; t_1 = 6, \; t_2 = 11, \; T=16.$

The controls and switching functions and state trajectories are displayed in Figure 4.

The remarkable fact about the controls $u(t)$ and $v(t)$ is that both controls are zero on the initial interval $[0, t_k], \, k=1, 2.$ The numerical results show that the controls satisfy the strict bang-bang property

$\begin{array}{lllll} \phi_u(t) > 0 & \forall \;\; 0 \leq t < t_1, & \dot\phi_u(t_1) < 0, & \phi_u(t) < 0 & \forall \;\; t_1 < t \leq T, \\ \phi_v(t) > 0 & \forall \;\; 0 \leq t < t_2, & \dot\phi_v(t_2) < 0, & \phi_v(t) < 0 & \forall \;\; t_2 < t \leq T. \end{array}$

In this problem, the first-order sufficient conditions [34] are satisfied, since the switching times $t_1=6$ and $t_2=12$ are determined by the boundary conditions $w(T) = 400$ and $z(T)=10$. Therefore, the solution displayed in Figure 4 provides a strict strong minimum.

Inspecting the switching function $\phi_v$, one recognizes that it is bending upwards on a terminal interval and is approaching zero. Then by increasing the weight $\alpha = 0.2$ in the objective to values $\alpha \geq 0.25$, the optimal control $v(t)$ will have a small terminal arc with $v(t) = 0$.

---

4.3.2. Optimal control solution with fixed final time $T=16$ and delays $h_1=10.6$ in $u$ and $h=1.84$ in $v$

Here we choose the objective

$J(u) = p(T) + 0.5\, q(T).$

The discretization approach with $N=10000$ grid points yields the more complicated control structures

$u(t) = u_{\rm max} \, |\, 0\, |\, u_{\rm max} \, , \;\;\;\; v(t) = 0\, |\, v_{\rm max} \, |\, 0 \, .$

The control $u(t)$ switches at $t_1 = 5.4$ and $t_3 = 11.4$, while $v(t)$ switches at $t_2 = 9.16$ and $t_4 = 14.16$. We obtain

$p(T) = 921.85, \;\; q(T) = 508.45.$

The numerical results allow to verify that the switching functions $\phi_u$ and $\phi_v$ defined in (40) satisfy the control law (41). Note that the switching function $\phi_u(t)$ has a jump at $t_1 = T-h_1 = 5.4$ while $\phi_v(t)$ has a jump at $t_4 = T-h = 14.16$.

As it has been expected by increasing the value of $\alpha$ in the objective functional we have reached two effects (in spite of delays in control variables): the terminal value of $q$ is reduced and there are two switches of control variables. Moreover the control $v(t)$ is a shift of the optimal non-delayed control $v(t)$ to the left by the value of time delay $h=1.84$. It is however worth noticing that in this case we are not able (because of the delays) to verify sufficient optimality conditions. The terminal values of state variables are less than in non-delayed case. However their comparison is not fair. In the model with delays there are two antiangiogenic agents while in the model without delays we have only one such agent.

---

5. Discussion and conclusions

In our opinion, models of combined therapy with multiple delays in control have not been discussed before. In this paper we propose to describe the effects of combined therapy by a two-compartmental model with two control variables with multiple delays which represent the differences in pharmacokinetics of different agents and different goals of the therapy. While the primary goal is related to eradication of a tumor or at least survival benefits, the secondary one is to normalize the tumor vasculature thereby facilitating chemotoxic drug delivery. This leads to a complex multi-control problem, the complete solution of which is much more complicated than in the single control case.

We have discussed two aspects of this problem, one of which is related to the question of attainability of an arbitrary final state of the system using admissible treatment protocols. This question has been answered in this paper, at least partially, by using sufficient conditions of relative constrained controllability for semilinear systems with multiple delays in control. We have found that such conditions are satisfied in the proposed models. The interesting finding is that the results are not structurally sensitive in the sense that they do not depend qualitatively on the structure of the model (within a class of models discussed). Since the conditions are only local, the next step of our study was devoted to optimal control synthesis to ensure driving the dynamical system to a neighborhood of the required final state. We have used necessary conditions of optimality for systems with delays in control and constrains imposed on control and state variables [18]. The optimal control was found numerically using a two-stage computational algorithm. The first stage is based on large scale non-linear programming for a discretized version of the optimal control problem and the second one is related to optimization of switching times by the arc-parametrization method. We have analyzed how sensitive the solutions are with respect to the type of growth function in the tumor dynamics, introduction of the time delays in control variable and the form of the objective functional. The Gompertz-type growth function, which has been used most often in modeling tumor growth, is not mandatory to describe the unperturbed tumor growth slowdown observed in clinical and experimental data [36]. Its drawback is that for small ratios of tumor volume and vascular carrying capacity the relative tumor growth capacity is unbounded. This feature is absent in the case when logistic type growth is used. Yet another advantage of using this model is the absence of singular arcs in optimal protocols of anti-angiogenic treatment which are present, when the Gompertz-type growth function is used. In our study we found that this property is true also in the case when two control variables (representing two anticancer modalities) with multiple delays are considered in the model. In this sense the optimal control problem is structurally sensitive, since the use of Gompertz-type growth leads to optimal controls with singular intervals (which are practically unrealizable), whereas the logistic-type growth yields pure bang-bang control. Moreover, their numerical computation, especially in the case of the fixed terminal time, is near to impossible (because of a strong chattering effect). The fixed time of treatment seems to be much more consistent with clinical practice. Especially, if we include time delays in control the choice of terminal time greater than maximal time delay is the only reasonable solution. In the case of logistic growth we have analyzed both cases with free and fixed terminal time. For the models without delays, we are able to verify sufficient optimality conditions. On the other hand our controllability conditions are independent of the choice of the tumor growth function.

In literature usually the objective functional takes into account only the final size of the tumor (exceptions are our papers [38,39,40,41]). We have decided to compare results for such functional with the case when the performance index is a linear combination of final values of both state variables. Although, qualitatively, the solutions are similar, our results show that the choice of weights in such function allows for better control of the final size of the supporting vascularity network, which is especially important in view of different goals of antiangiogenic therapy discussed in section 1. It should be mentioned that although we have considered only positive weights the negative value of $\alpha$ in the objective functional may be also reasonable. Such choice allows to keep the angiogenic vasculature large enough to facilitate chemotoxic drug delivery which may be one of the goals of the combined therapy (see section 1).

Introduction of multiple delays in control variables in the models has led to some changes in understanding and testing conditions of controllability and optimality and their numerical computation. Qualitatively different machinery should be used for models with delays in state variables, as proposed in [10] and analyzed in [35]. Other notions of controllability should be applied and the optimality conditions, although based on the same version of the Maximum Principle used by us, lead to more complex mathematical formulas.

---

Acknowledgments

The authors wish to thank Prof. Ronald Hancock and dr Roman Jaksik for their assistance in preparation of the final version of the manuscript.

---