Optimal control problems for differential equations applied to tumor growth: state of the art

In this paper, we shall focus on topics where methods from optimal control theory will be applied to biomedical problems arising from cancer treatments. Specifically, we will focus on optimal control problems for systems of differential equations. Optimal control for partial differential equations, stochastic differential equations, discrete difference equations and so on are not considered in this document.

There exists a vast literature in which optimal control methods are applied to cancer treatment. The early literature on applications of optimal control to cancer treatment is almost exclusively on cancer chemotherapy (see references [14, 32, 45]). Much of this work is on cell-cycle nonspecific mathematical models and focusing on different growth models [32]. The monograph [31] by Martin and Teo is entirely devoted to such models, but including increasingly more comprehensive medical conditions such as single and multi-drug resistance.

Later, a greater emphasis in research has been put on cell-cycle specific models. In reference [10], the authors considered phase-specific administration of cytotoxic agents to increase the selectivity of therapy. Reference [46] models the proliferation cycle of leukemia. Eisen [14] also gives some introductory discussion of cell-cycle specific chemotherapy and compartmental models (distinguishing between proliferating and quiescent cells), but mostly in the context of pharmacokinetics and clearly the emphasis is more on modeling than on analysis.

After that, a great effort was done to study the development and analysis of compartmental models which describe and analyze the actions of drugs on specific compartments from a control theory point of view [49, 50, 51]. In these models, the emphasis is on the cancerous cells and side effects are only measured indirectly through the total dose of drugs administered.

Thus, from the mathematical point of view, the modeling of chemotherapy for cancer using the optimal control theory has 40 years of history. Different frameworks have been considered in the literature, depending on the study approach, from macroscopic levels [45] to microscopic levels [25,26], including systems of differential equations (v. [9,11,16,40] and their references) and equations in partial derivatives (v. [45] and their references). Normally, the control variable represents the dose of the drug and the state variables represent the tumor cells. In general, the main objective is to determine the optimal application of the drugs that minimizes the number of tumor cells at the end of treatment, while maintaining the toxic effects of the drug at an acceptable level. To reduce the toxicity of the anticancer drug, which also affects normal cells, it is customary to add some control restrictions, isoperimetric or state. There are different types of restrictions depending on the assumptions of the model, mainly punctual or integral. For example, a maximum rate is prescribed at which the chemotherapeutic drug is administered at each time or a limit is imposed on the total dose of the drug in the therapy interval (see [9, 12, 27]).

On the other hand, mathematical modeling of tumor growth in general is very important when it comes to having a model that fits the biology of the problem [35]. There are many mathematical models that deal with the different aspects of the growth of a tumor. They use ordinary differential equations [2, 6, 13, 54], partial differential equations [1, 3, 34], stochastic processes [38], game theory [52], etc.

In this review, we focus on 4 applications of the control theory to the growth of tumors:

Section 2 shows several compartmental models, which have been useful in the understanding of the application of optimal control theory to the growth of tumors.

Another topic of great importance is the study of the administration of chemotherapy in brain tumors [6, 20, 37]. Thus, in this review we study a model developed in the reference [40] based on a work by Ribba [37]. This revision can be found in Section 3.

Next topic is the drug administration in presence of multi-drug resistance cells. There exists a vast literature on optimal administration of chemotherapeutical drugs in the presence of drug resistance [7, 8, 18, 26, 33]. A model of sensitive and resistant cells to the treatment can be found in Section 4.

Optimal control is useful as a methodology in response to mathematical models for novel cancer treatments as anti-angiogenic treatments [27, 39]. A mathematical model with anti-angiogenic and chemotherapeutical treatment can be found in Section 5.

Finally, the review ends with some Conclusions in Section 6.

Fig. 1

In this section, we analyze a set of differential equations for a class of cell-cycle specific compartmental models for cancer chemotherapy.

Each cell goes through a sequence of phases from birth to cell division: G₁→G₀→S→G₂→M. While the simplest mathematical models for the treatment of cancer with chemotherapy, applying the theory of optimal control, treat the complete cycle as a compartment, the multicompartmental models group the phases of the cell cycle with the purpose of representing the action of the different types of chemotherapeutic agents used. Cytotoxic agents generally act in the phase G₂/M because the cell walls become thin and porous in mitosis M and, therefore, the cell is more vulnerable to attack. The cytostatic agents aim to synchronize the transitions of cells through the cell cycle causing a brief and invisible inhibition of the DNA synthesis in the S phase and the retention of the cells in the G₁ phase. Recruitment agents cause the cells to leave the inactive G₀ compartment where they are not vulnerable to attacks by any chemotherapeutic drug. This classification is general, being able to find agents that have more than one effect depending on their concentration [43]. Depending on the specific type of tumor that is considered, and the drugs that are administered, an appropriate grouping of the phases of the cell cycle is carried out.

This model deals with the application of a single cytotoxic agent that is activated in the G₂/M phase of the cell cycle. Taking into account the sensitivity of the drug to the phase, the cell cycle is divided into two compartments, one combining the second growth phase G₂ and the mitosis M and the other formed by the rest of the phases. The state of the system is described by a two-dimensional vector where N₁(t) represents the mean number of cancer cells in the first compartment at time t (composed of the phases G₀, G₁ and S) and N₂(t) the average number of carcinogenic cells in the second compartment at time t (composed of G₂ and M). The equation for the second compartment has the form

(1)N˙<!-- ˙ -->2(t)=−<!-- − -->a2N2(t)+a1N1(t)$$\begin{equation}\dot{N}_2(t)=-a_2N_2(t)+a_1N_1(t) \end{equation}$$

where a_i > 0 is the average inverse transit time through the i-th compartment. The equation for the first compartment becomes

(2)N˙<!-- ˙ -->1(t)=−<!-- − -->a1N1(t)+2a2N2(t).$$\begin{equation}\label{eq1_c1} \dot{N}_1(t)=-a_1N_1(t)+2a_2N_2(t). \end{equation}$$

In matrix form, the system given by the equations (1) and (2) can be expressed as N = AN where N = (N₁,N₂)^T and A 2 ℝ^2×2 with

A=−<!-- − -->a12a2a1−<!-- − -->a2.$$A=\begin{pmatrix} -a_1 & 2a_2\\ a_1 & -a_2 \end{pmatrix}.$$

Now, we consider a cytotoxic agent whose concentration in blood is represented by the function u acting in the G₂/M phase. It is assumed that all cells are sensitive to treatment (homogeneous population). According to the log-kill hypothesis, it is assumed that the concentration of the drug u(t) kills a fraction of the output flow a₂N₂(t) of cells of the G₂/M compartment and, therefore, the number of cells dead is given by fu(t)a₂N₂(t), ϕ being a parameter of constant chemotherapeutic death. The control set is a compact interval [0,u_max], with u_max the maximum dose rate or the maximum concentration. In the model, the u control always appears together with the constant ϕ and so, in order to keep the number of free parameters as low as possible, it is combined with the maximum dose rate in an amount represented as u_max under the assumption that u_max≤ 1. If the concentration is high enough, then really u_max = 1 is realistic: almost all the cancer cells in that compartment can be killed. The cells that die in G₂/ M leave this compartment, that is, they are counted as outputs from the second compartment, but they no longer enter the first compartment. In this sense, it is considered to kill the cell even if apoptosis is not induced. The remaining fraction (1-u)a₂N₂ suffers cell division and therefore the controlled mathematical model becomes

(3)N˙<!-- ˙ -->1(t)=−<!-- − -->a1N1(t)+2(1−<!-- − -->u)a2N2(t),N1(0)=N10,$$\begin{eqnarray}\dot{N}_1(t)&=&-a_1N_1(t)+2(1-u)a_2N_2(t),\quad N_1(0)=N_{10},\\ \end{eqnarray}$$

(4)N˙<!-- ˙ -->2(t)=a1N1(t)−<!-- − -->a2N2(t),N2(0)=N20.$$\begin{eqnarray}\dot{N}_2(t)&=&a_1N_1(t)-a_2N_2(t),\quad N_{2}(0)=N_{20}. \end{eqnarray}$$

The matrix expression for the system has the form [41]

(5)N˙<!-- ˙ -->(t)=(A+uB)N(t),N(0)=N0,$$\begin{equation}\dot{N}(t)=(A+uB)N(t),\quad N(0)=N_0, \end{equation}$$

where A and B are

(6)A=−<!-- − -->a12a2a1−<!-- − -->a2,B=0−<!-- − -->2a200.$$\begin{equation}A=\begin{pmatrix} -a_1 & 2a_2\\ a_1 & -a_2 \end{pmatrix}, \quad B=\begin{pmatrix} 0 & -2a_2\\ 0 & 0 \end{pmatrix}. \end{equation}$$

Consider now as a functional objective

(7)J=rN(T)+∫<!-- ∫ -->0TqN(t)+su(t)dt→<!-- → -->min$$\begin{equation}J=rN(T)+\int_{0}^{T}qN(t)+su(t) dt\rightarrow \text{min} \end{equation}$$

where T is the final fixed time, r = (r₁,r₂) is a positive vector of weights and q = (q₁,q₂) is a vector row of weights not negative. The penalty term rN(T) = r₁N₁(T)+r₂N₂(T) represents an average weight of the total number of cancer cells at the end of treatment and the term qN(t) = q₁N₁(t)+q₂N₂(t) measures the volume of the tumor during treatment. Side effects of the treatment are only indirectly included in the model through the minimization of the total dose, ∫<!-- ∫ -->0Tu(t)dt$\int^{T}_{0}u(t)dt$. The positive coefficient s in this integral is redundant since it could be absorbed by the other weights. Recall that the number of cancer cells that do not undergo cell division at time t and that are consider eliminated is given by u(t)a₂N₂, that is, u(t) is proportional to the fraction of inefficient divisions. Since the drug kills healthy cells at a similar rate, the integral ∫<!-- ∫ -->0Tu(t)dt$\int^{T}_{0}u(t)dt$ represents the cumulative negative effects of the treatment on the normal tissue or its toxicity. The aim is to minimize tumor volume in both compartments at the end of treatment, but taking into account the tumor volume in both compartments during treatment, while minimizing side effects. Thus, the following optimal control problem arises:

Theorem 2

If (N_*,u_*) is an optimal controlled trajectory for problem 1, then there does not exist an interval on which the control u_* is singular.

We show some examples of extremal bang-bang controls that have been computed using the values in Table 1, which are shown in Figure 2.

Fig. 2

Table 1

Numerical values for the coefficients and parameters used in numerical computations for the optimal control problem 1, extracted from [41].

Coefficient	Interpretation	Numerical value
a1	Inverse transit time through G0/G1+ S	0.197
a2	Inverse transit time through G2/M	0.356
N10	Initial condition for N1	0.7012
N20	Initial condition for N2	0.2988
umax	Maximum dose rate/concentration effectiveness of the drug	0.90
s	Penalty/weight for the total dose of cytotoxic agent	0.50
q1	Penalty/weight in the objective for the average number of cancer cells in G0/G1+ S during therapy	0.10
q2	Penalty/weight in the objective for the average number of cancer cells in G2/M during therapy	0.10
r1	Penalty/weight in the objective for the average number of cancer cells in G0/G1+ S at the end of therapy	3
r2	Penalty/weight in the objective for the average number of cancer cells in G2/M at the end of therapy	3
T	Therapy horizon	T = 7,21,60

This model was originally formulated by Swierniak et al. [47], where a cytotoxic agent, which can act in S or G₂/M phases, is combined with a recruitment agent, which attracts inactive cells in the G₀ compartment to re-enter the active cell cycle and thus be attacked by the cytotoxic agent. One of the biggest problems in the chemotherapeutic treatment of some leukemias is a large residue of inactive cells in the G₀ phase, which are not sensitive to most of the cytotoxic agents. This problem has also been reported in breast and ovarian cancer. In this model, three compartments are considered: G₀, G₁ and S +G₂/M. The state space is thus the first orthant P in R₃. The variables N₀, N₁ and N₂ represent the states. The newborn cells initiate the process of cell division, but they can remain dormant and remain in the quiescent stage G₀. Let p₀ and p₁ be positive numbers, with p₀+ p₁ =1, representing the probabilities that the daughter cells will pass through the first growth phase without obstacles (p₁) or fall into the state of quiescence (p₀). A recruiting agent is applied to reduce the average transit time through the G₀ compartment, resulting in a higher outflow from G₀. If w represents the concentration of the recruitment agent, it is assumed that the outflow is increased by a factor of 1 +w, 0 ≤ w ≤ w_max, where w = 0 represents that no drug is being applied and w_max corresponds to a full dose treatment. A cytotoxic agent u, 0 ≤ u ≤ u_max, is applied in the third compartment, where u_max corresponds to the maximum dose (see Figure 3 ). The combination of these drugs with the standard characteristics of cell division and cell death (and under the log-kill hypothesis), results in a bilinear control system, of the form

Fig. 3

(18)N˙<!-- ˙ -->(t)=A+uB1+wB2N(t),N(0)=N0,$$\begin{equation}\label{sistema_general} \dot{N}(t)=\left(A+uB_1+wB_2\right)N(t), \quad N(0)=N_0, \end{equation}$$

where the matrices A, B₁, B₂ are

(19)A=−<!-- − -->a002p0a2a0−<!-- − -->a12p1a20a1−<!-- − -->a2,B1=00−<!-- − -->2p0a200−<!-- − -->2p1a2000,B2=−<!-- − -->a000a000000.$$\begin{equation}A= \begin{pmatrix} -a_0 & 0 & 2p_0a_2\\ a_0 & -a_1 & 2p_1a_2\\ 0 & a_1 & -a_2 \end{pmatrix}, \quad B_1= \begin{pmatrix} 0 & 0 & -2p_0a_2\\ 0 & 0 & -2p_1a_2\\ 0 & 0 & 0 \end{pmatrix}, \quad B_2= \begin{pmatrix} -a_0 & 0 & 0\\ a_0 & 0 & 0\\ 0 & 0 & 0 \end{pmatrix}. \end{equation}$$

As in the 2– compartment model considered earlier, the a_i are positive coefficients related to the inverse transit times of the cancer cells through the compartments.

Thus, we considere the following optimal control problem:

Theorem 4

Let (N_*,u_*,w_*) be an optimal controlled trajectory for problem 2. If the weight s₂in the Lagrangian for the recruiting agent w is zero, s₂ = 0 then there does not exist an interval on which the control w_* is singular. However, if this coefficient is positive, s₂> 0, then the Legendre-Clebsch condition for optimality of a singular control w is satisfied.

The proof of this theorem can be found again in [47].

For small positive weights s₂, and this is the important case for this model, bang-bang controls are optimal. Figure 4 shows the graphs of such controls u and w and their corresponding trajectories when the weight s₂ in the Lagrangian is chosen as s₂ = 0 and s₂ = 0:1 for parameter values specified in Table 2.

Fig. 4

Table 2

Numerical values for the coefficients and parameters used in numerical computations for the optimal control problem 2 with cytotoxic and recruitment agents, extracted from [41].

Coefficient	Interpretation	Numerical value
a0	Inverse transit time through G0	0.05
a1	Inverse transit time through G1	0.5
a2	Inverse transit time through S +G2/M	1
p0	Probability that cells enter G0	0.9
p1 = 1– p0	Probability that cells enter G1	0.1
N0(0)	Initial condition for N0	0.8589
N1(0)	Initial condition for N1	0.0954
N2(0)	Initial condition for N2	0.0456
umax	Maximum dose rate/concentration effectiveness	0.95
wmax	Maximum dose rate/concentration effectiveness	6
s1	Penalty/weight at the cytotoxic agent	1
s2	Penalty/weight at the recruiting agent	0, 0.1
q0	Penalty/weight in the objective for the average number of cancer cells in G0 during therapy	3
q1	Penalty/weight in the objective for the average number of cancer cells in G1 during therapy	1
q2	Penalty/weight in the objective for the average number of cancer cells in S +G2/M during therapy	1
r0	Penalty/weight in the objective for the average number of cancer cells in G0 at the end of therapy	3
r1	Penalty/weight in the objective for the average number of cancer cells in G1 at the end of therapy	1
r2	Penalty/weight in the objective for the average number of cancer cells in S +G2/M at the end of therapy	1
T	Therapy horizon	T = 21

In this model the combination of two drugs is considered: a cytotoxic agent u, which is responsible for killing tumor cells and acts in the G₂/M phase, and a blocking agent v (cytostatic), which slows the development of tumor cells in the S phase, preventing them from passing into the cell division phase [48]. Therefore, three compartments are considered: the first formed by the phases G₀ and G₁, the second by the S phase and the third is formed by the phases G₂ and M. The space of states is the first orthant P in ℝ³, being N_i(t), i = 1,2,3, the number of cancer cells in the i-th compartment at time t.

As a result, the flow of tumor cells from the second compartment to the third is reduced by a factor given by v(t), from a value of a₂N₂(t) to (1-v(t))a₂N₂(t), where 0 ≤ v(t) ≤ v_max< 1. The value v(t) = 0 represents that the treatment is not being applied, while v_max represents a complete dose. The control u represents the dose of cytotoxic agent administered with u = 0 indicating the absence of treatment and u = 1 the maximum dose (see Figure 5). Under all assumed conditions for the previous compartmental models (homogenous tumor population, constituted by cells sensitive to therapy and under the log-kill hypothesis), we get back to a bilinear system again of the form (18), but now the matrices involved have the form

Fig. 5

A=−<!-- − -->a102a3a1−<!-- − -->a200a2−<!-- − -->a3,B1=00−<!-- − -->2a3000000,B2=0000a200−<!-- − -->a20.$$A= \begin{pmatrix} -a_1 & 0 & 2a_3\\ a_1 & -a_2 & 0\\ 0 & a_2 & -a_3 \end{pmatrix}, \quad B_1= \begin{pmatrix} 0 & 0 & -2a_3\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{pmatrix}, \quad B_2= \begin{pmatrix} 0 & 0 & 0\\ 0 & a_2 & 0\\ 0 & -a_2& 0 \end{pmatrix}.$$

We thus have the following optimal control problem

Theorem 6

Suppose that

(25)a1+a2(1−<!-- − -->vmax)−<!-- − -->2a3≥<!-- ≥ -->0.$$\begin{equation}\label{c3} a_1+a_2(1-v_{\text{max}})-2a_3\geq 0. \end{equation}$$

If (N_*,u_*,v_*) is an optimal controlled trajectory for problem 3, there also does not exist an interval on which the control u_* is singular.

We close this section with some numerical experiments of locally optimal controls for problem 3. The parameter values are given in Table 3 and condition (25) is satisfied. Figure 6 shows two different scenarios.

Fig. 6

Table 3

Numerical values for the coefficients and parameters used in computations for the optimal control problem 3 with cytostatic and cytotoxic agents, extracted from [41].

Coefficient	Interpretation	Numerical value
a1	Inverse transit time through G1/G0	0.197
a2	Inverse transit time through S	0.395
a3	Inverse transit time through G2/M	0.107
N1(0)	Initial condition for N1	0.3866
N2(0)	Initial condition for N2	0.1722
N3(0)	Initial condition for N3	0.4412
umax	Maximum dosage/concentration/effectiveness of cytotoxic agent	0.95
vmax	Maximum blocking effect of cytostatic agent	0.30
s1	Penalty/weight at the cytotoxic agent	1
s2	Penalty/weight at the cytostatic agent	0.01
q1	Penalty/weight in the objective for the average number of cancer cells in G1/G0 during therapy	1, resp. 0.1
q2	Penalty/weight in the objective for the average number of cancer cells in S during therapy	1, resp. 0.1
q3	Penalty/weight in the objective for the average number of cancer cells in G2/M during therapy	1, resp. 0.1
r1	Penalty/weight in the objective for the average number of cancer cells in G1/G0 at the end of therapy	1, resp. 8.25
r2	Penalty/weight in the objective for the average number of cancer cells in S at the end of therapy	1, resp. 8.25
r3	Penalty/weight in the objective for the average number of cancer cells in G2/M at the end of therapy	1, resp. 8.25
T	Therapy horizon	T = 21

In this section, we consider a low grade glioma (LGG), which is a term denoting a set of primary progressive brain tumor of astrocytic or oligodendroglial origin characterized by slow and continuous growth preceding the anaplastic transformation [30]. Such model was considered by Ribba and colaborators [37] and assumes that a LGG is composed of proliferative P and nonproliferative quiescent tissue Q. The chemotherapy treatment u(t) eliminates proliferative cells and the quiescent cells are also affected by the treatment and become damaged quiescent cells Q_p. Damaged quiescent cells, when re-entering the cell cycle, be able to repair their DNA and become proliferative once again or can die because of unrepaired damages. For simplicity, here it is assumed that the dosage, the concentration and the effect of the drug are equal. Thus, it is assumed that the drug kills a proportion of u of proliferating cells, γuP, and damages a proportion of u of quiescent cells, γuQ. Therefore, the system of differential equations that describes the proliferative, non-proliferative and damaged cells and the effect of chemotherapy on these populations is given by [40]

(26)dPdt=cpP1−<!-- − -->NK+cQPQp−<!-- − -->cPQP−<!-- − -->γ<!-- γ -->Pu,P(0)=P0,$$\begin{eqnarray}\label{sist_eq0} \frac{dP}{dt}&=& c_p P\left(1-\frac{N}{K}\right)+c_{QP}Q_p-c_{PQ}P-\gamma Pu,\quad P(0)=P_0, \end{eqnarray}$$

(27)dQdt=cPQP−<!-- − -->γ<!-- γ -->Qu,Q(0)=Q0,$$\begin{eqnarray} \frac{dQ}{dt}&=&c_{PQ}P-\gamma Qu,\quad Q(0)=Q_0, \end{eqnarray}$$

(28)dQpdt=γ<!-- γ -->Qu−<!-- − -->cQPQp−<!-- − -->δ<!-- δ -->QPQp,Qp(0)=0,$$\begin{eqnarray} \frac{dQ_p}{dt}&=&\gamma Qu-c_{QP}Q_p-\delta_{QP}Q_p,\quad Q_{p}(0)=0, \end{eqnarray}$$

where N = P+Q+Qp is the total number of cells and P,Q,Q_p ≥ 0. Consider now the following functional objective

(29)Jw(u)=(1−<!-- − -->w)(P(T)+Q(T))+w∫<!-- ∫ -->0T(P(t)+Q(t))dt,w∈<!-- ∈ -->[0,1],$$\begin{equation}\label{funcional1} J_{w}(u)=(1-w)(P(T)+Q(T))+w\int_{0}^{T}(P(t)+Q(t))dt, \quad w\in[0,1], \end{equation}$$

where w is a weight parameter. We obtain two different functionals with biological relevance for w = 0 and w = 1. Thus, for w = 0 we obtain a functional in the form of Mayer that focuses on the tumor cells at the end of the treatment,

J0(u)=P(T)+Q(T),$$J_0(u)=P(T)+Q(T),$$

while for w =1 we obtain a functional in the form of Lagrange that takes into account the tumor cells throughout the dynamic process:

J1(u)=∫<!-- ∫ -->0TP(t)+Q(t)dt.$$J_1(u)=\int_{0}^{T}P(t)+Q(t)dt.$$

In addition, we introduce a restriction for additional control by imposing a limit on the total dose of the drug:

(30)∫<!-- ∫ -->0Tu(t)dt≤<!-- ≤ -->M.$$\begin{equation}\label{restr1} \int_{0}^{T}u(t)dt\leq M. \end{equation}$$

Thus, we have the following optimal control problem:

We solve the optimal control problem 4 using AMPL software [19], which is linked to the IPOPT Interior Point optimization code. We use a high number N 2 [3000,10000] of grid points and implement the trapezoidal integration scheme [5].

The objective of this section is to calculate a standard chemotherapy protocol for a LGG, where the total amount of TMZ dose is M = 5 and the treatment time is T = 1 (30 days). We use the parameters of Table 4 and the initial conditions

P(0)=0.924,Q(0)=42.3,Qp(0)=0.$$P(0)=0.924, \quad Q(0)=42.3, \quad Q_p(0)=0.$$

Fig. 7

J0(u∗<!-- ∗ -->)=37.0623,P(T)=0.8396,Q(T)=36.7801,Qp(T)=5.5286.$$J_0(u^*)=37.0623,\quad P(T)=0.8396,\quad Q(T)=36.7801,\quad Q_p(T)=5.5286.$$

w =1. Figure 8 shows the optimal control and the state variables for the weight w =1 in the objective (29). The optimal treatment is based on DMT therapy, but now, in contrast to the solution for w = 0, the initial dose is the maximum dose until the total dose M = 5 is reached. Subsequently, no dose is administered. We get the following final value for the objective and the state variables

Fig. 8

J1(u∗<!-- ∗ -->)=38.0623,P(T)=0.8607,Q(T)=36.7804,Qp(T)=5.4084.$$J_1(u^*)=38.0623,\quad P(T)=0.8607,\quad Q(T)=36.7804,\quad Q_p(T)=5.4084.$$

w = w_s = 0:004067. For this value, the optimal control is given by

(38)0≤<!-- ≤ -->w<ws:u∗<!-- ∗ -->(t)=0for0≤<!-- ≤ -->t≤<!-- ≤ -->251for25<t≤<!-- ≤ -->30$$\begin{eqnarray}0\leq w<w_s&:&u^*(t)= \begin{cases} 0 & \text{for}\quad 0\leq t\leq 25\\ 1 & \text{for} \quad 25<t\leq 30 \end{cases} \end{eqnarray}$$

(39)ws<w≤<!-- ≤ -->1:u∗<!-- ∗ -->(t)=1for0≤<!-- ≤ -->t≤<!-- ≤ -->50for5<t≤<!-- ≤ -->30$$\begin{eqnarray}w_s<w\leq 1&:&u^*(t)= \begin{cases} 1 & \text{for}\quad 0\leq t\leq 5\\ 0 & \text{for} \quad 5<t\leq 30 \end{cases} \end{eqnarray}$$

Therefore, for w ≠ w_s all controls are bang-bang with a single switching that occurs at the time when the total amount of dose, M =5, has been administered. But the structure of the control "switch" in the weight w_s = 0:004067 being, for this value, a totally singular optimal control.

Figure 9 show this fact.

Fig. 9

Now, we solve the optimal control problem (29) for T = 30 months.

The singular control u_sing is obtained in (36). On the other hand, we can obtain the change times t₁ = 20:1 and t₂ = 25:2. The control and state variables are shown in Figure 10. We get the following final value for the objective and the state variables

Fig. 10

J0(u∗<!-- ∗ -->)=0.9059,P(T)=0.2545,Q(T)=0.6514,Qp(T)=38.2335.$$J_0(u^*)=0.9059,\quad P(T)=0.2545,\quad Q(T)=0.6514,\quad Q_p(T)=38.2335.$$

Here, the treatment starts with zero doses, then jumps to a single dose with fairly small values followed by a total dose at the end of the treatment.

2. w = 1. Here, we consider the functional

J1(u)=∫<!-- ∫ -->0T(P(t)+Q(t))dt,$$J_1(u)=\int_{0}^{T}(P(t)+Q(t))dt,$$

and we find again the control structure bang-singular-bang

(41)u∗<!-- ∗ -->(t)=1for0≤<!-- ≤ -->t<t1,using(x(t),λ<!-- λ -->(t))fort1≤<!-- ≤ -->t≤<!-- ≤ -->t2,0fort2<t≤<!-- ≤ -->30.$$\begin{equation}u^*(t)= \begin{cases} 1 & \text{for}\quad 0\leq t<t_1,\\ u_{\text{sing}}(x(t),\lambda(t)) & \text{for}\quad t_1\leq t\leq t_2,\\ 0 & \text{for}\quad t_2< t\leq 30. \end{cases} \end{equation}$$

The switching times for are t₁ =3:00 and t₂ =18:4 showing a fairly long singular arc. We get the following final value for the objective and the state variables

J1(u∗<!-- ∗ -->)=139.31,P(T)=4.7982,Q(T)=1.6042,Qp(T)=20.0689.$$J_1(u^*)=139.31,\quad P(T)=4.7982,\quad Q(T)=1.6042,\quad Q_p(T)=20.0689.$$

Figure 11 shows the state variable and the optimal control.

Fig. 11

Resistance to drugs is one of the most important problems in the treatment of cancer, causing a high number of deaths. There are several probabilistic models for the development of drug resistance [8, 53]. The amplification of a gene is an increase in the number of copies of that gene present in the cell after cell division, the depletion corresponds to a decrease in its number of copies. The hypothesis of amplification of a gene to a copy establishes that in cell division at least one of the two daughter cells will be an exact copy of the mother cell while the second, with some positive probability, will undergo genetic amplifications. Next, we present a model studied by Ledzewicz et al. in [26].

Cancer is fundamentally a genetic disease: it arises as a consequence of pathological changes in the information contained in the DNA. That is, the uncontrolled growth of the cells is the result of alterations or mutations in the genetic material. These mutations lead to a process of altering the cell division and the inability to undergo programmed cell death (or apoptosis). The growth of a single cancer (or primary cancer) does not usually lead to the death of the organism. However, some cancer cells can acquire the ability to enter the bloodstream, travel to another part of the body and begin to grow in a different organ. This process is known as metastasis. Generally, metastatic growth is what kills the body. When a cancer cell is created, it can multiply and lead to a tumor of small size without the need for blood supply. At this stage, the growth of an avascular tumor stops. Subsequently, to continue growing, since the tumor needs oxygen and nutrients to grow, provided by the blood, the tumor sends out chemical signals called promoters. This induces the blood vessels to grow towards the tumor and leads to the complete vascularization of the tumor, which allows it to increase in size. The formation of new blood vessels is called angiogenesis. The understanding of the role of angiogenesis in the development of cancer has advanced significantly due to the studies conducted by Folkman [17]. We now know that without angiogenesis, the size of most primary tumors does not exceed 3 mm. One type of very important promoters are the so-called Vascular Endothelial Growth Factors (VEGF). These factors stimulate cellular responses joining the tyrosine kinase receptors (a type of enzyme) on the surface of the endothelial cells and stimulate their proliferation and migration to the tumor. Other pro-angiogenic factors are PDGF, PIGF or FGF. Therefore, anti-angiogenic therapies are important in the treatment of certain cancers. But, by itself, the antiangiogenic treatment only prevents the tumor from developing its support of blood vessels and does not directly destroy the cancer cells. In fact, the benefits of using anti-angiogenic therapies alone are, at best, transient and are followed by the restoration of tumor growth and progression [4]. As a consequence, it is common not to limit the treatment to a body, but to combine different treatments in the hope of achieving a synergistic effect. In this way, antiangiogenic treatment should be combined with other treatments, such as radiotherapy or chemotherapy. In this section, we combine antiangiogenic therapy with chemotherapy.

The mathematical model studied in [39] is

(52)dpdt(t)=−<!-- − -->ξ<!-- ξ -->p(t)logp(t)s(t)−<!-- − -->φ<!-- φ -->v(t)p(t),$$\begin{eqnarray}\frac{dp}{dt}(t)&=&-\xi p(t)\log\left(\frac{p(t)}{s(t)}\right)-\varphi v(t) p(t), \end{eqnarray}$$

(53)dsdt(t)=χ<!-- χ -->q(t)−<!-- − -->μ<!-- μ -->s(t),$$\begin{eqnarray}\label{eq1_P5} \frac{ds}{dt}(t)&=&\chi q(t)-\mu s(t), \end{eqnarray}$$

(54)dqdt(t)=−<!-- − -->χ<!-- χ -->q(t)+bp(t)−<!-- − -->dp(t)2/3q(t)−<!-- − -->γ<!-- γ -->u(t)q(t),$$\begin{eqnarray}\frac{dq}{dt}(t)&=&-\chi q(t)+bp(t)-dp(t)^{2/3}q(t)-\gamma u(t)q(t), \end{eqnarray}$$

where p(t) is the size of the tumor, s(t) is the number of stable vessels, q(t) is the number of unstable vessels and u(t) and v(t) are the antiangiogenic and cytotoxic agents respectively. Table 5 shows the values of the parameters of the model.

Thus, the problem that we want to solve is to minimize the size of the tumor p(t):

(55)J(u)=p(T)$$\begin{equation}\label{funcional_Prob5} J(u)=p(T) \end{equation}$$

for a free final time T and subject to the constraints

(56)∫<!-- ∫ -->0Tu(t)dt≤<!-- ≤ -->M,∫<!-- ∫ -->0Tv(t)dt≤<!-- ≤ -->N,$$\begin{equation}\int_{0}^{T}u(t)dt\leq M,\quad \int_{0}^{T}v(t)dt\leq N, \end{equation}$$

where the set of admissible controls is given by

u∈<!-- ∈ -->U={u(t)|0≤<!-- ≤ -->u(t)≤<!-- ≤ -->uM,∀<!-- ∀ -->t∈<!-- ∈ -->[0,T]},v∈<!-- ∈ -->V={v(t)|0≤<!-- ≤ -->v(t)≤<!-- ≤ -->vM,∀<!-- ∀ -->t∈<!-- ∈ -->[0,T]}.$$\begin{eqnarray*}u\in U&=&\{u(t) | 0\leq u(t)\leq u_{M}, \forall t\in[0,T]\},\\ v\in V&=&\{v(t) | 0\leq v(t)\leq v_{M}, \forall t\in[0,T]\}. \end{eqnarray*}$$

Thus, we have the following problem:

Proposition 9

If the optimal control u^* is singular over an interval I =(a,b), then λ₂,λ₃,λ₄,λ₅>0. In addition, the control u^* ends in an interval (t₁,T) where u^* = 0.

Then, we suppose the control u is singular and let’s calculate that control. After some tedious calculations, we get the following singular optimal control

(62)using∗<!-- ∗ -->=λ<!-- λ -->1ψ<!-- ψ -->1+λ<!-- λ -->2ψ<!-- ψ -->2+λ<!-- λ -->3(φ<!-- φ -->bv+ψ<!-- ψ -->3)2bγ<!-- γ -->λ<!-- λ -->3,$$\begin{equation}u^{*}_{\text{sing}}=\frac{\lambda_1\psi_1+\lambda_2\psi_2+\lambda_3(\varphi b v+\psi_3)}{2b\gamma\lambda_3}, \end{equation}$$

where

ψ<!-- ψ -->1(p,q,s)=−<!-- − -->χ<!-- χ -->ξ<!-- ξ -->qs,ψ<!-- ψ -->2(p,q,s)=2χ<!-- χ -->b+χ<!-- χ -->qp(μ<!-- μ -->−<!-- − -->χ<!-- χ -->)−<!-- − -->dp−<!-- − -->1/3qχ<!-- χ -->,ψ<!-- ψ -->3(p,q,s)=bξ<!-- ξ -->logps−<!-- − -->χ<!-- χ -->−<!-- − -->dp2/3.$$\begin{eqnarray*}\psi_1(p,q,s)&=&-\chi\frac{\xi q}{s},\\ \psi_2(p,q,s)&=&2\chi b+\chi\frac{q}{p}(\mu-\chi)-dp^{-1/3}q\chi,\\ \psi_3(p,q,s)&=&b\left[\xi\log\left(\frac{p}{s}\right)-\chi-dp^{2/3}\right]. \end{eqnarray*}$$

Regarding to the control v, analytically it is complicated to characterize it for what the numerical methods are used. Thus, in all the numerical simulations, ϕ< 0 and thus the control v is of the bang-bang type.

Some numerical solutions are shown in Figure 12 and Figure 13.

Fig. 12

Fig. 13

In this survey, a series of problems of optimal control theory applied to cancer treatments have been addressed. First of all, a class of cell cycle specific compartmental models for cancer chemotherapy was analyzed. By considering the phases of the cell cycle separately, it is possible to appropriately model the different actions of various drugs involved. A common characteristic of the models analyzed in the second section is that it is implicitly assumed that the cancer population is homogeneous and consists of cells that are sensitive to the chemotherapy applied. If we then minimize a weighted average of the tumor population and the total dose of chemotherapy given over a fixed therapy interval, it will be seen that in this scenario it indeed is optimal to give chemotherapy in one full dose session upfront at the beginning of the therapy interval. These results are fully consistent with and confirm as optimal the classical MTD (maximum tolerated dose) regimen. Essentially, the underlying scenario in these models is one where cancer is growing rapidly at a critical level, but sensitive to chemotherapeutic agents. The intuitive optimal solution then is to hit it as hard as possible, as soon as possible.

In the third section, low grade glioma, a type of brain tumor, together with the action of a cytotoxic agent, as the temozolomide, is studied [39]. For this problem, the optimal control solutions depend on the treatment time. If the treatment time is long enough (T = 30 months), the optimal control will be unique. On the other hand, if the treatment time is T = 30 days, the control will be of the bang-bang type.

Many tumors consist of heterogeneous agglomerations of subpopulations of cells that show widely varying sensitivities toward the actions of a particular chemotherapeutic agent [15, 29]. Coupled with the fact that growing tumors also exhibit considerable evolutionary ability to enhance cell survival in an environment that is becoming hostile, this leads to multi-drug resistance of some strains of the cells. Naturally, it makes sense to combine drugs with different activation mechanisms to reach a larger population of the tumor cells-and this is what is being done-but the sad fact remains that some cells develop multi-drug resistance to a wide variety of even structurally unrelated drugs. There may even exist subpopulations of cells that are not sensitive to the treatment from the beginning (ab initio, intrinsic resistance). For certain types of cancer cells, there are simply no effective agents known. We have explored this aspect in section four, studying a 2– compartment model with sensitive and resistant populations.

Finally, we have considered tumor’s microenvironment. One of the most important component of a tumor’s microenvironment is its vasculature. In order to grow beyond a small size, a growing tumor needs to develop its own network of blood vessels and capillaries that will provide it with nutrients and oxygen. This process is called angiogenesis and was already pointed out as a therapeutic target by J. Folkman in the early 1970s [23,24]. Antiangiogenenic treatments aim at depriving the tumor of this needed vasculature by either disrupting the signaling process that the tumor uses to recruit surrounding, mature, host blood vessels or by directly inhibiting the growth of endothelial cells that form the lining of the newly developing blood vessels and capillaries. We have studied this kind of treatment together with the chemotherapy [39] which is based on this model by Hahnfeldt et al. [21].

We hope that this survey will be useful in understanding the optimal treatment protocols in tumors. Mathematical models need to be developed to aid in our understanding of how to implement this work and hopefully to show how to develop new optimal treatment strategies.