Ultimate dynamics and optimal control of a multi-compartment model of tumor resistance to chemotherapy

In this paper, a non-trivial generalization of a mathematical model put forward in [ 35 ] to account for the development of resistance by tumors to chemotherapy is presented. A study of the existence and local stability of the solutions, as well as the ultimate dynamics of the model, is addressed. An analysis of different chemotherapeutical protocols using discretization and optimization methods is carried out. A number of objective functionals are considered and the necessary optimality conditions are provided. Since the control variable appears linearly in the associated problem, optimal controls are concatenations of bang-bang and singular arcs. A formula of the singular control in terms of state and adjoint variables is derived analytically. Bang-bang and singular controls from the numerical simulations are obtained where, in particular, singular controls illustrate the metronomic chemotherapy.

1. Introduction. Nowadays, cancer is a major health problem. Together with heart diseases, is one of the main causes of death in the industrialized world [47]. The annual global economic burden of cancer has been estimated to be above two trillion dollars and it is expected to increase in the next coming years. In addition, in 2015, 107 billion dollars were spent in cancer medicines, according to a report from the IMS Institute for Healthcare Informatics. While new drugs are being developed and tested in clinical trials, it is also important to optimize treatments protocols for the existing drugs, which could result in better outcomes and in a more effective use of current resources.
Chemotherapy is quite often an essential, if not the only, component of the therapeutical approaches administered to patients depending on the disease' state. In most common malignancies occurring worldwide, such as lung, breast and colon cancers, chemotherapy is used in 7 − 57%, 38 − 66% and 14 − 66%, respectively. However, its effect on 5-year survival in adult malignancies contributes only to about a 2 percent increase [34]. In the overwhelming majority of cases this is due to the resistance that tumors develop to treatments. For chemotherapy, this causes that subsequent large panels of drugs become ineffective to control the progression of malignant cells. Biological explanations to this phenomenon goes from Darwinian selection [19] to non-Darwinian (e.g. Lamarckian) dynamics, where a combination of both being the currently most accepted framework [38].
Therefore, when designing effective therapies, it is essential not only to target the total tumor mass but to also take into account the emergence of resistance and its dynamics [52]. This requires to incorporate heterogeneity in the tumor cell populations. Hence, mathematical models aimed at dealing with the interplay of cancer cells and chemotherapy must include these ingredients in their foundation. In Ref. [54] a simple mathematical model consisting of ordinary differential equations (ODEs) for two tumor cell populations, sensitive and resistant cells to the drug was proposed. Later, in [12] and [3], in addition to proliferation and resistance to chemotherapy, more than two cell populations were included in the analysis. In [50,30,5], other processes were incorporated into the equations, such as inhibition and/or diffusion. In Refs. [43,1], both ordinary and partial differential equation models were proposed to consider specific processes of resistance transfer mediated by direct contact, via microtubes, and nonlocal microvesicle dissemination.
The mathematical modelling of cancer chemotherapy by using optimal control theory has many years of history. There exist different areas in the field of Applied Mathematics where optimal control theory has been used. Without being exhaustive, we can cite references for systems of ODEs (see, for instance, [26,13,14,15,17,46,42] and references therein) and partial differential equations [9,10,51,40].
In the last years, optimal control problems for drug resistance tumor cells have attracted considerable attention [27,2,28,29]. The competition between sensitive and resistant cells was included in a ODE model in reference [11]. In this model, resistant cells were selected in the presence of drugs according to a Darwinian interpretation and optimal control analysis suggested alternative ways to schedule the chemotherapy doses that could result in a better control of the total cancer population. More complex models considering growth, resistance and chemotherapy using partial differential equations and with its corresponding optimal control analysis were developed in [40].
Since tumor resistant cells do not respond to the administered cytotoxic agent but quite often display a significantly lower proliferation rate in comparison with sensitive tumor cells, one alternative therapeutical pathway is to find a balance between sensitive and resistance cells in order to keep the total amount of tumor cells below a certain threshold. This is the idea of the adaptative therapy [18] in which it is explained that maintenance of the cancer as an isolated organ might be a preferable strategy over attempting to entirely kill it, as this will most likely produce an aggressive recurrence. Therefore, a continuous maintenance of this kind of strategy might be preferable over a maximum tolerated dose therapy approach. The final goal of this rationale is to make cancer a manageable chronic disease.
In Ref. [35], an increasingly complex number of mathematical models were developed for a specific brain tumor, low-grade glioma treated with the chemotherapeutic agent temozolomide. These models incorporated genetic and non-genetic acquired resistance together with the drug dynamics. Then the authors classified patients in groups according to their acquired resistance mechanisms. Here, our aim is to study one of the models that was validated with clinical data in [35] and provide new treatment protocols that could benefit those patients. Even if this model has been specifically developed for a low-grade glioma, our study is sufficiently general and suggestive that it could be useful for other cancer types by making the corresponding modifications.
The paper is organized as follows. In Section 2, we introduce the mathematical model given in [35] and propose a nontrivial generalization into a system of ODEs. Next, in Section 3, we study the dynamics and the stability of the model system using the qualitative theory of differential equations and analyze the ultimate dynamics of our extended model. In Section 4, we formulate optimal control problems with several objective functionals. We discuss the necessary optimality conditions of Pontryagin's Maximum (Minimum) Principle [39]. Since the control variable appears linearly in the dynamical system and objectives, optimal controls are concatenations of bang-bang and singular arcs. We derive a formula for the singular control in terms of state and adjoint variables. A further necessary condition is the generalized Legendre-Clebsch condition [21,45] which we evaluate explicitly. Section 5 presents optimal control solutions for different objectives and short/large treatment periods. The computations are based on discretization techniques and nonlinear programming methods. Finally, in Section 6, we discuss and summarize our conclusions.
2. The mathematical model. In Ref. [35], the authors developed various quantitative approaches to address the resistance to chemotherapy of low-grade glioma cells. One of the presented models was validated using clinical data from 121 lowgrade glioma patients and it assumes that the tumor is composed of sensitive S, damaged D and resistance R cells.
The chemotherapeutic treatment u(t) acts in the following way. The sensitive cells S are affected by u(t) and subsequently transformed into damaged cells D with a (positive) rate constant α s . These damaged cells may die, via apoptosis, with a (positive) rate constant µ d or, else, repair their DNA and/or alter their membrane receptors/transporters becoming resistance cells to the administered drug. This last process of phenotypic transformation, which is irreversible, takes place with a (positive) rate constant µ r . We also include the term −γy in the second equation and +γy in the third equation indicating that when there is no drug in the system, the damaged cells can repair their DNA with a rate γ and they become resistant if they get it, going to the third compartment. Moreover, we assume that, as more drug u(t) is administered, an increase in the probability that tumor cells become resistant is produced, i.e., the drug induces the emergence of a resistant cell subpopulation. For simplicity, here we consider that the dosage, the concentration and the effect of the drug are equal, i.e., the pharmacokinetics and pharmacodynamics are not explicitly modelled. Thus, the system of differential equations governing the dynamics of the sensitive, damaged and resistance cells and the effect of chemotherapy on those populations is given by with initial conditions x(0) = x 0 , y(0) = 0 and z(0) = 0, where K is the carrying capacity of the system. Therefore, since no damaged and resistant cells are assumed to be present in the tumor prior to the administration of the drug, no Darwinian selection takes place. Our scenario thus considers Lamarckian induction of a resistant phenotype by the presence of therapy as the sole trigger of the phenotypic transformation. Moreover, we assume that once sensitive cells have acquired a resistant phenotype they irreversibly remain in that state, and thus do not return to a sensitive state after any subsequent time elapsed. This assumption is in agreement with several works dealing with resistances in tumours (see for instance [38] and [41]).
The above system (1) can be considered as a particular case of the following more general system of differential equations x, y and z refer to tumor mean diameters and are measured in mm, (although they could also represent an average number of cells per unit of length) and where 1. The function f : [0, +∞) → R is a C 1 class and a decreasing function. 2. f (0) = 1 and f (K) = 0. A similar model was studied in reference [7].
3. Existence and local stability of solutions and ultimate dynamics of the model. Let R 3 +,0 = {x; y; z ≥ 0}. We state the following result of existence of solutions: Theorem 3.1. For any initial data (x 0 , y 0 , z 0 ) such that x 0 , y 0 , z 0 ≥ 0 and u ∈ C 1 such that u ∈ [0, u M ], the solutions to the system (4) exist for all t > 0, are nonnegative and unique.
Since the right hand side of Eqs. (4a), (4b) and (4c) is a continuous function in (x, y, z), existence of solutions follows from Cauchy-Peano theorem [31]. As the partial derivatives of the right hand side of the system are continuous and bounded, uniqueness follows from the Picard-Lindelof theorem. Now, we perform an asymptotic analysis for the system (4) for u(t) = u = constant. We find that the equilibrium points of the system being biologically meaningful are: P 1 = (0, 0, 0), P 2 = (0, 0, K). In order to check the stability of the fixed points, the Jacobian matrix J is calculated for system (4):

Thus
• Substituting the equilibrium point P 1 in the matrix, using the hypothesis on f and calculating the eigenvalues, we obtain λ 1 = ρ r , λ 2 = −(µ d + µ r u + γ) and λ 3 = ρ s −α s u. Thus, P 1 is a saddle point and therefore an unstable point. • Substituting P 2 in J and calculating the eigenvalues, we get λ 1 = −α s u, λ 2 = −(µ d + µ r u + γ) and λ 3 = ρK ∂f ∂p (0, 0, K). Since, for hypothesis, f is a decreasing function, the three eigenvalues are negative, and therefore P 2 is a stable node. Now, we are going to study the ultimate dynamics of our model. The ultimate global dynamics of various cancer models with (or without) therapy has attracted the attention of many researchers in recent years, see, for example [20,44,48,49].
The reason for this attention is related to the interest in the patient's health forecast, depending on the tumor type and the parameters of the treatment methods used. Thus, in this section we are going to deal with this aspect for system (4).
We introduce notations for the following sets: . and introduce some preliminary definition: We denote L F h 1 as the Lie derivative with respect to the function F , where F denotes the vector field corresponding to the system (4) and h 1 is a C ∞ class function.
Moreover, we establish the following theorem, which can be found in reference [22], and will be useful in the proof of lemma 3.3.  (4). Examples of convergent trajectories to P 2 for x 0 + y 0 + z 0 ≤ K. Parameters used to calculate the phase portrait are given in Table 1.
Let h be a differential function which is bounded below on the set D ⊂ R n . Then the solution ϕ(x, t) of the systemẋ = v(x), with initial condition Next, we present the following lemma: Proof. Let us apply the function h 1 = x + y + z with respect to D 1 (ε), with ε > 0. It follows from the properties of the function f that max w≥K+ε f (w) < 0.
By −δ(ε) we denote the value of this maximum and let Then we have So by Theorem 3.2, for arbitrarily small ε > 0 any trajectory in D 1 (ε) leaves D 1 (ε) in a finite time. Therefore either this trajectory leaves D 1 (0) in a finite time or this trajectory is contained in D 1 (0) and in virtue of the LaSalle theorem tends to the maximal invariant set M 1 in D 3 ∪ D 4 . In order to study M 1 we should consider 1. The equation for D 3 together with the equation obtained by time differentiation: f (x + y + z)(ρ s x + ρ r z) − µ d y = 0, with y > 0 2. Equations for D 3 and D 4 together with equations obtained by time differentiation: Thus solving these equations we derive that M 1 = P 2 . It means that P 2 is the maximal compact invariant set in D 3 ∪D 4 which is locally stable. Hence, a trajectory through a point in D 3 ∪ D 4 either go inside D 2 or tends to P 2 . Now, we are in disposition to prove our next theorem: Theorem 3.4. Any trajectory in R 3 +,0 tends to one of the equilibrium points P 1 or P 2 .
Proof. In virtue of Lemma 3.3 it is sufficient to examine trajectories only within D 2 . Let us take the function h 2 = −z with respect to D 2 . Then and it is equal to zero with Then, by the LaSalle theorem (see for instance [23,37]) each trajectory in D 2 approaches to the maximal invariant set M 2 in considered on D 2 are satisfied only at P 1 we have that M 2 = P 1 ∪ P 2 . Hence, we arrive to the desired assertion.
Remark 1. Let M st (P 1 ) be the stable manifold of P 1 and u = 1. Then each point in R 3 +,0 which is not contained in M st (P 1 ) tends to P 2 . This fact corresponds to the worst health prognosis for a patient. (see Fig. 1). 4. Optimal control problem. In standard chemotherapy protocols, drugs are administered in cycles at maximum tolerated doses (MTD) alternating with periods of rest. The underlying rationale is that the more anti-cancer drug is given to the patient, the more tumor cells are killed. Hence, the chemotherapeutic schedule traditionally used in the clinical practice consists of delivering the maximal dose of a chemotherapeutic (or a combination) that can be tolerated by the patient.
Alternatively, the concept of metronomic therapy appears, which is the almost continuous administration of chemotherapeutic drug at significantly lower doses rather than MTD therapy. The hope is that, in absence of side effects, it is possible to give chemotherapy over prolonged time intervals, and, therefore, to extend the time horizon of treatment and to obtain an improvement of the treatment in comparison with repeated short MTD doses.
On this basis, the following sections study these protocols.

4.1.
A general formulation of the optimal control problem. In this subsection, the optimal control problem for Eqs. (4) is studied. The control variable u represents the drug dose. In general, the main objective is to determine the optimal drug schedule which minimizes, as much as possible, the number of tumor cells during all the treatment and/or at the final time of the treatment period using a maximum quantity of drug. The set of admissible controls for the dynamical system is given by Notice that u is now a measurable function. Here, u M represents the MTD, which cannot be exceeded during the course of the therapy. In addition, we introduce another control constraint by imposing a bound on the total amount of drug dose administered: From a clinical point of view, the objective function should be designed in such a way that the tumor cell populations x and z are minimized by an admissible control u satisfying the dynamical system (4) and the constraint (6). This can be done by using the following objective functional: where s 1 and s 2 are a weight parameters. By introducing the additional equationẇ = u(t), w(0) = 0, the constraint (6) is equivalent to the terminal constraint w(T ) ≤ M . Then the control system can be written as:ẋ The state variable vector of the control problem can be cast as v = (x, y, z, w) T ∈ R 4 Since the control variable u appears linearly in the dynamical system, we can write the equation (8) 4.2. Pontryagin's Minimum Principle: Necessary optimality conditions. In this subsection, we are going to evaluate the necessary optimality conditions of Pontryagin's Minimum Principle [39,45] for the optimal control problem with the objective functional J s (u) given in (7). Defining the adjoint variable λ = (λ 1 , λ 2 , λ 3 , λ 4 ), the Hamiltonian is The adjoint equations becomė where a = (s 2 , 0, s 2 , 0) T , b = (s 1 , 0, s 1 ) T and Dh and Dg denote the matrices of the partial derivatives of the vector fields h and g, respectively: Thus, we can write the set of differential equations for the adjoint variables λ i , i = 1..4: On the other hand, the optimal control u * minimizes the Hamiltonian H with respect to u ∈ [0, u M ] and therefore, we obtain the switching function φ(v, λ) as φ(v, λ) = λ, g(v) = −λ 1 α s x + λ 2 (α s x − µ r y) + λ 3 µ r y + λ 4 (11) and writing φ(t) = φ(v(t), λ(t)), the optimal control satisfies the switching condition with singular controls possible if the switching function vanishes over an open interval I s ⊂ [0, T ]. Thus, the control is not determined a priori by the minimum condition at times when φ(t) = 0 and we need to analyze the derivatives of φ(t).

4.3.
Chemotherapy as an optimal control problem and singular controls. As pointed out above, in order to study the optimal control one needs to analyze the derivatives of the switching function φ(t) and in these calculations the use of the Lie bracket of vector fields is extremely useful [45].
Definition 4.1. Given two continuously differentiable vector fields f and g defined on some open set G ⊂ R n , f, g : G → R n , their Lie bracket [f, g] is another vector field on G whose expression is .

Proposition 1.
Let v(t) be a solution of the dynamics (8) for the control u and let λ be a solution of the corresponding adjoint equation (10). For any continuously differentiable vector field m, we define the function Ψ(t) as and the derivative of this function is given bẏ where [h + ug, m] denotes the Lie bracket of the vector fields h + ug and m.
Using this proposition, we can calculate the first derivative of the switching function φ, equation (11) Since this formula does not depend on the control, we can differentiateφ once more to getφ and the control appears in the second derivative of the switching function if The order of a singular control u on an interval I s = [t 1 , t 2 ] is defined as the minimum integer q such that ∂ ∂u d 2q φ dt 2q = 0 Thus, if G 1 (v, λ) = 0 the order of a singular control in our problem is q = 1. Moreover, a further necessary optimality condition is the generalized Legendre-Clebsch condition [4,8,45]: . If this quantity is positive we say strict Legendre-Clebsch condition holds. Thus, for our problem, the strict Legendre-Clebsch condition is satisfied if Assuming this condition, the relationφ(t) = 0 yields a formula for the singular control in terms of the state variable v and adjoint variable λ: where To calculate the functions G 1 and G 2 in (17) we compute the Lie brackets: These previous results allow us to determine G 1 (v, λ): and therefore, to check if the strict Legendre-Clebsch condition (16): is satisfied on the singular interval I s . This is a necessary condition for existence of singular control. Finally, we can compute G 2 (v, λ), equation (18): and, with these expressions for G 1 (v, λ) and G 2 (v, λ), to obtain the singular control u sing , equation (17). Thus, we have the following result, which summarizes these general calculations:

b) Sensitive cells x(t) c) Damaged cells y(t) d) Resistant cells z(t).
The time horizon was fixed to T = 1 month. In this case, the singular control is given by

Numerical solutions.
In this section, we solve numerically the optimal control problem (4), for f (x + y + z) = 1 − (x + y + z)/K. Thus, we are solving the problem stated in reference [35]. To do it, we use the direct method, which consists of first discretize and then optimize. The optimal control problem is transformed into a large-scale nonlinear programming (NLP) problem using a suitable scheme of discretization [6,33]. The AMPL software was employed in order to obtain the solution of the optimal control problem. Moreover, we resorted to the Applied Model Interior Point optimization (IPOPT) method in our own codes which is very useful to solve efficiently the problem [53]. We will solve (4) for different objective functionals in order to reach a detailed study of the problem. Thus, a first objective functional corresponds to a functional

) Sensitive cells x(t) c) Damaged cells y(t) d) Resistant cells z(t).
The time horizon was fixed to T = 30 months.
in a Mayer form which focusses on the tumor cells for terminal time T : Another different objective corresponds to a functional in a Bolza form: which takes into account the tumor cells during the whole dynamical process. Both functionals, (20) and (21) are particular cases of the functional (7) previously studied, for s 1 = 1, s 2 = 0 and s 1 = 0, s 2 = 1, respectively.

5.1.
Optimal control for the functional J 0,1 (u): Short and long times. Here we aim to obtain the protocol of chemotherapy for a total amount of drug M = 1/6 and time of treatment T = 1 month for the objective J s1,s2 (u) for s 1 = 0, s 2 = 1. This value of M means that we have drug for a sixth part of the month. Recall that the parameters and the initial conditions are given in Table 1.
The optimal control together with the switching function and the state variables are shown in Figure 2. The results show that the optimal treatment corresponds to the MTD therapy, where the initial dose is maximal up to when the total dose M = 1/6 is reached. After that, no further dose is administered. Thus, we obtain a bang-bang control as in (12): where t 1 = 1/6.  To show the optimality of the control u * (t) in Figure 2 a), we apply the secondorder sufficient conditions (SSC) in Osmolovskii, Maurer [36]. Since there is only one switching point t 1 and one terminal constraint w(T ) = 1/6, the SSC degenerate into a first-order sufficient condition, which are fulfilled if the following strict bangbang property holds Figure 2 a) shows that the strict bang-bang property holds. Now, we solve the optimal control problem (1) for J 1 (u) and for T = 30 months. We increase the total amount of drug to M = 5. Thus, we have a total amount of drug M to be delivered in 30 months. For this case, we obtain a control u * with the structure bang-singular-bang: The optimal control together with the switching function and the state variables is shown in Figure 3 . The singular control u sing is given by (17). The switching times from bang to singular control and singular to bang control are respectively t 1 = 2.8 and t 2 = 24.1, exhibiting a rather long singular arc.
Here the treatment starts with a maximum total dose, then it jumps to a singular dose with small values and it finally ends with a zero dose. It is an example of a metronomic therapy. We can verify that the strict generalized Legendre-Clebsch

b) Sensitive cells x(t) c) Damaged cells y(t) d) Resistant cells z(t).
The time horizon was fixed to T = 1 month.
condition (16) is satisfied on the singular interval [t 1 , t 2 ]. In this case, it is not possible to check the sufficient optimality conditions, since to the best of our knowledge numerically verifiable sufficient conditions are not available in the literature. Since this protocol can be difficult to implement in the clinical practice, we can compare this singular control with a suboptimal protocol of the bang-bang type that is the one that could be included in clinical practice. We have chosen for this case that the treatment starts at the beginning of the time horizon due to the similarity that this treatment has with the singular control obtained previously.
Thus, figure 4 shows such a suboptimal bang-bang protocol together with the state variables.
For the optimal protocol, the value of the functional objective and the terminal state variables are As can be seen, the differences between the optimal protocol and the suboptimal protocol are minimal and, due to the difficulty of implementing a singular arc in clinical practice, it could be appropriate to use the suboptimal protocol (bang-bang control) instead of the optimal protocol (singular control) for this case.

) Sensitive cells x(t) c) Damaged cells y(t) d) Resistant cells z(t).
The time horizon was fixed to T = 30 months.
state variables are shown in Figure 5 for T = 1 month. The results show that the optimal treatment corresponds to the MTD therapy. Thus, we obtain a bang-bang control. Initially, the dose level is zero and then jumps to the maximum tolerated dose (MTD) at t 1 = 5/6 (in units of months). Since the strict bang-bang property holds: the sufficient optimality condition is satisfied [36]. Thus, it proves the optimality of control u * (t). Now, we consider T = 30 months. We obtain a control u * with the structure bang-singular-bang: The optimal control together with the switching function and the states variables is shown in Figure 6. Again, u sing is given by (17). The switching times are computed as t 1 ≈ 0.2 and t 2 ≈ 26.7. Then, the treatment starts with a maximum total dose, then it jumps to a singular dose with a long singular arc and finally it ends with a maximum total dose again. We also verify that the strict generalized Legendre-Clebsch condition (16) is satisfied on the singular interval [t 1 , t 2 ]. Again, it is not possible to check the sufficient optimality conditions. Again, this protocol can be difficult to implement in the clinical practice, then we compare this singular control with two suboptimal protocols of the bang-bang type. Unlike the previous case, we have implemented two different protocols: 1. For the first protocol, we choose that the treatment ends in the final value of the time horizon since the objective function focuses on minimizing the tumor cells at the end of the treatment. Thus, figure 7 shows such a suboptimal bang-bang protocol together with the state variables.
Comparing the optimal protocol with this protocol, for the optimal protocol, the value of the functional objective and the terminal state variables are Notice that the number of resistance cells in the optimal protocol is less than the suboptimal protocol. The difference in percentage between the values of both functional objectives is approximately 15%. Moreover, this protocol allows a moderate growth of sensitive cells at the beginning of the schedule of treatment.
Thus, for this case, we think that the differences are relevant enough not to rule out the optimal protocol in favor of the suboptimal protocol, although the former is more difficult to implement in clinical practice than the latter. 2. For the second protocol, we chose one that approximates as much as possible the singular control obtained, but that can be implemented in the practice clinic. Thus, for this protocol, the treatment starts from the beginning until reaching 1 month, there is a long rest period, and the rest of the treatment is given the last 4 months. Figure 8 shows such a suboptimal bang-bang protocol together with the state variables.
Comparing the optimal protocol with this protocol, the values of the functional objective and the terminal state variables for the optimal protocol are In the comparison of both protocols, it is observed that the differences are not relevant, being able then to choose the suboptimal protocol instead of the optimal protocol in this case. 6. Conclusions and discussion. In this paper, we have put forward and studied a generalized dynamical model for the chemotherapeutic treatment of low grade gliomas, based on a previous work in [35], where the biological variables of the model are sensitive, damaged and resistant tumor cells.
We have presented a study on the local stability analysis and ultimate dynamics of the model followed by the introduction of an optimal control problem assuming different types of objective functionals. Since the control variable appears linearly both in the control system and in the objective functionals, the evaluation of the necessary optimality conditions shows that optimal controls are concatenations of bang-bang and singular controls. Using tools from geometric optimal control theory we determined the optimal protocols and the corresponding states associated to these optimal protocols. Moreover, an explicit formula for the singular control in terms of state and adjoint variables has been derived.
Our numerical simulations show two different therapeutical scenarios: For short treatment times we obtained simple bang-bang controls with only one switch which represents MTD therapies. This type of therapies are used in the clinic, where the chemotherapeutic schedule traditionally used consists of the administration of the maximal dose of the chemotherapeutical agent that can be tolerated by the patient. For longer treatment times we have found bang-singular-bang controls. These types of controls provide an example of metronomic therapies, which is an alternative to MTD therapy. We have also compared such metronomic therapies with optimal subprotocols which can be carried out more easily in clinical practice, obtaining in general good approximations to the optimal protocol Thus, these results offer a theoretical test benchmark to which realistic and easily realizable schemes can be compared and thus the efficiency of these protocols can be evaluated.
Finally, we have used other types of functions f = (1 − (x + y + z)/K) n , with n = 2, 3 in the system of equations, but the optimal controls obtained are, qualitatively speaking, similar to those obtained in this work. Therefore, they were not considered here.
In future works, investigations into more complex models will address the ability of combined therapies to control the resistant cells to the therapy which are the major cause of treatment failure in cancer patients.