Computational and Mathematical Methods in Medicine

Volume 2014 (2014), Article ID 982978, 15 pages

http://dx.doi.org/10.1155/2014/982978

Research Article

Delay Differential Model for Tumour-Immune Response with Chemoimmunotherapy and Optimal Control

F. A. Rihan,1,2 D. H. Abdelrahman,¹ F. Al-Maskari,³ F. Ibrahim,4,5 and M. A. Abdeen²

¹Department of Mathematical Sciences, College of Science, UAE University, P.O. Box 15551, Al-Ain, UAE

²Department of Mathematics, Faculty of Science, Helwan University, Cairo 11795, Egypt

³Zayed bin Sultan Al Nahyan Center for Health Sciences, College of Medicine and Health Sciences, UAE University, P.O. Box 17666, Al-Ain, UAE

⁴Institut für Angewandte Mathematik, LS III, TU Dortmund, Vogelpothsweg 87, 44227 Dortmund, Germany

⁵Department of Mathematics, Faculty of Science, South Valley University, Qena 83523, Egypt

Received 13 May 2014; Revised 1 July 2014; Accepted 3 July 2014; Published 14 August 2014

Academic Editor: Loredana G. Marcu

Copyright © 2014 F. A. Rihan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We present a delay differential model with optimal control that describes the interactions of the tumour cells and immune response cells with external therapy. The intracellular delay is incorporated into the model to justify the time required to stimulate the effector cells. The optimal control variables are incorporated to identify the best treatment strategy with minimum side effects by blocking the production of new tumour cells and keeping the number of normal cells above 75% of its carrying capacity. Existence of the optimal control pair and optimality system are established. Pontryagin’s maximum principle is applicable to characterize the optimal controls. The model displays a tumour-free steady state and up to three coexisting steady states. The numerical results show that the optimal treatment strategies reduce the tumour cells load and increase the effector cells after a few days of therapy. The performance of combination therapy protocol of immunochemotherapy is better than the standard protocol of chemotherapy alone.

1. Introduction

Cancer is considered one of the most complicated diseases to be treated clinically and one of the main causes of death. Accordingly, a great research effort is being devoted to understand the interaction between the tumour cells and the immune system. The treatment of cancer is then one of the most challenging problems of modern medicine. The cancer treatment should kill cancer cells in the entire body and in the meantime keep the healthy cells above the minimum level. Chemotherapy is one of the most prominent cancer treatment modalities. However, it is not always a comprehensive solution for tumor regression. Other treatment options, including surgery, immunotherapy, and radiation, are often able to force the cancer into remission, but better and suitable treatments are needed to fulfil the requirements [1–3].

Recently, chemotherapy is used in combination with immunotherapy to protect the patient from opportunistic infection, as well as fighting the cancer itself [4, 5]. This is due to the fact that the chemotherapy treatment kills both cancerous and healthy cells and therefore it depletes the patients immune system, making the patient prone to dangerous infections. For this and other reasons, it is desirable to strengthen the immune system after an immune-depleting course of chemotherapy. Additionally, however, the ability to recruit the body’s own defenses to fight the cancer can be a powerful treatment strategy. Therefore, maintaining a strong immune system, by combining immunotherapy during chemotherapy, may be essential to successfully fight the cancer. However, the query now is how to most effectively combine cancer immunotherapy and chemotherapy?

Mathematical models, using ordinary, partial, and delay differential equations, play an important role in understanding the dynamics and tracking tumour and immune populations over time (see, e.g., [6–14]). Kuznetsov et al. [3] model the interactions of cytotoxic T lymphocyte (CTL) response and the growth of an immunogenic tumour. In the contributions of [15–17], the authors take into account the penetration of the tumour cells by the effector cells, which simultaneously causes the inactivation of effector cells. Recently, in [18], the authors adopted a predator-prey formulation of the tumour immunity problem as a battle between immune cells and tumour cells (predators and prey, resp.). Many research papers have also been done on formulations of the mathematical models describing the interaction between tumour cells and immune cells alone, between tumour cells and normal cells alone, and between tumour cells and chemotherapy treatment alone [19, 20]. We should mention here that the application of the optimal control theory requires the boundedness of the solutions of the model populations; see also [21, 22].

The objective of this paper is to adopt a delay differential model and analyze and provide computationally an optimal way to combine chemotherapy and immunotherapy treatment strategies that identify the best treatment strategy that can minimize the tumor load while maximizing the strength of the immune system. We formulate and analyze a delay differential model describing immune response and tumour cells under the influence of chemotherapy alone and the combinations of chemotherapy and immunotherapy. The outline of the paper is as follows. In Section 2, we describe the model. In Section 3, we study the qualitative behaviour of the model without external therapy. In Section 4, we describe the optimal control problem governed by DDEs with only chemotherapy control variable. Existence of the solution and optimality conditions are discussed in Sections 5 and 6. In Section 7, we extend the control problem to include a combination of chemotherapy and immunotherapy treatments with two-contrail variable. Numerical simulations and conclusion are given in Sections 8 and 9.

2. Description of the Model

Let us recall Kuznetsov et al.’s model [3] that describes the dynamics of tumour immune system interactions by incorporating a system of five ordinary differential equations (ODEs) model; then we reduce it into two equations but with time delays. The model describes the response of the effector cells (ECs) to the growth of tumour cells (TCs). The penetration of TCs by ECs has been taken into account, which simultaneously causes the inactivation of ECs. It is assumed that interactions between the ECs and TCs are in vitro such that , , , , and denote the local concentrations of ECs, TCs, EC-TC conjugates, inactivated ECs, and “lethally hit” TCs, respectively. The total population of unattacked TCs is . The rate of binding of ECs to TCs and the rate of separation of ECs from TCs without damaging them are denoted by and , respectively. The rate at which EC-TC integrations program for lysis is denoted by while the rate at which EC-TC interaction inactivate ECs is denoted by . The model takes the form

Here, the parameter represents the normal rate (not increased by the presence of the tumour) of the flow of adult ECs into the tumour side (region), , , and (when ) describes the accumulation of ECs in the tumour side due to the presence of the tumour. However, , , and are the coefficients of the processes of destruction and migration of , , and , respectively. The maximal growth of tumour is represented by the coefficient , and is the environment capacity. If we assume that , therefore where , and the model can be reduced to two equations which describe the behavior of ECs and TCs only [2, 3]. That leads to the fact that the rate of stimulated accumulation has some maximum value as get large.

Thus, the reduced model that describes the interactions between the two populations, tumour cells and activated effector cells (such as cytotoxic T-cells or natural killer cells), is of the form

with the initial conditions and . The interaction between the effector cells and tumour cells can reduce the size of both populations with different rates. This is expressed as and to illustrate the interaction between the tumour cells and effector cells. As a result of this interaction, the immune effector cells decrease the population of tumour cells at the rate . Additionally, tumour cells infect some of the effector cells and, therefore, the population of uninfected effector cells decreases at the rate .

If one considers as prey and as predator, then may take the form which is Michaelis-Menton form. In this term, is the maximum immune response rate and is the steepness of immune response. The presence of the tumour cells virtually initiates the proliferation of tumour-specific effector cells to reach a saturation level parallel to the increase in the tumour populations. Hence, the recruitment function should be zero in the absence of the tumour cells, whereas it should increase monotonically towards a horizontal asymptote [23]. We also incorporate a discrete time-delay into the model, to describe the time needed by the immune system to develop a suitable response after recognizing the tumour cells. Accordingly, the model with discrete time delay takes the form

with the initial functions and , for all . Here represents the normal rate (not increased by the presence of the tumour) of the flow of adult effector cells into the tumour side (region). The source of the immune cells is considered to be outside of the system, so it is reasonable to assume a constant influx rate . Furthermore, in the absence of any tumour, the cells will die at the rate . The presence of tumour cells stimulates the immune response, represented by the positive nonlinear growth term for the immune cells where and are positive constants and is the time delay that presents the time needed by the immune system to develop a suitable response after recognizing the tumour cells. The saturation term (Michaelis-Menton form) with the compartment and logistic term with the compartment are consoled. The presence of the tumour cells virtually initiates the proliferation of tumour-specific effector cells to reach a saturation level parallel to the increase in the tumour populations. Let us first prove the nonnegativity and boundedness solutions of the underlying DDEs model (3) (see [24]).

2.1. Nonnegativity and Boundedness Solutions of Model (3)

To show that the solutions of model (3) are bounded and remain nonnegative in the domain of its application for sufficiently large values of time , we recall the following lemma.

Lemma 1 (Gronwall’s Lemma [25, page 9]). Let , , and be real continuous functions defined in , for . One supposes that on one has the inequality . Then in .

Therefore, we arrive at the following proposition.

Proposition 2. Let be a solution of the DDEs model (3); then and for all sufficiently large time , where

Proof. Let denote the solution of model (3). From the second equation of system (3), we have . Thus, may be compared with the solution of

This proves that . From the first equation of system (3), we obtain

To show that is bounded, we use the generalized Gronwall Lemma. Since and , we have

The generalized Gronwall Lemma gives where is uniformly bounded. It follows that if is a solution of (3), then for all . This shows that the solutions of model (3) are uniformly bounded. This completes the proof.

From (1) and the solution , we arrive at the following result.

Corollary 3. If , then the solutions for model (3) are nonnegative for any nonnegative initial condition. However, if , then there exist nonnegative initial conditions such that becomes negative in a finite time interval.

2.2. Model with Chemotherapy

We extend model (3) to consider extra two variables, namely, amount of chemotherapy, , and normal cells, (see Figure 1). We also assume a homogeneity of the tumour cells. The modified model is

We assume that the drug kills all types of cells but that the killing rate differs for each type of cells; is the fraction cell kill for a given amount of drug, , at the tumour site. We denote by , , and the three different response coefficients. represents the amount of dose that is injected into the system, while is the decay rate of the drug once it is injected. In this case, the quantity we will control directly is not but . The tumour cells and normal cells are modelled by a logistic growth law, with the parameters representing the growth rate of two types of cells: identifies the parameter associated with tumour, and identifies the one associated with the normal tissue; and are the reciprocal carrying capacities of tumour cells and host cells, respectively. In addition, there are two terms representing the competition between the tumour and host cells and .

Figure 1: The interaction of tumour cells, immune cells, and normal cells in the presence of chemotherapy drug.

Let be the Banach space of continuous functions mapping the interval into with the topology of uniform convergence. It is easy to show that there exists a unique solution of system (8) with initial data . For biological reasons, we assume that the initial data of system (8) satisfy , , , and . For , the model is reduced to ODEs model developed by De Pillis and Radunskaya in [26].

Remark 4. We consider that the units of the model cells are normalized, so that .

The main objective in developing chemotherapy treatment, in system (8), is to reach either tumour-free steady state or coexisting steady state in which the tumour cells’ size is small, while the normal cells’ size is closed to its normalized carrying capacity. We next start the analysis with studying the stability of the system when there is no drug (treatment) input; that is, , for all .

3. Drug-Free Steady States and Their Stability

In the absence of chemotherapy (), model (8) has the following types of steady states:

(a)tumour-free steady state, where the tumour cells population is zero, while the normal cells survive. This steady state has the form ;(b)dead (lethal) steady state, where the normal cells population is zero, which has the following forms:

(i) in which the normal and tumour cell populations have died off,(ii) where the normal cells alone have died off and the tumour cells have survived, where

and is a nonnegative solution for

(c)coexisting steady state, where both normal and tumour cells coexist with nonzero populations. The steady state is given by where , and is a nonnegative solution of

The number of coexisting steady states mainly depends on the parameter values. There could be zero, one, two, or three of these steady states (see Figure 2). We next study the stability of the previously mentioned steady states.

(a)

(b)

Figure 2: (a) shows the regions of existence of coexisting equilibria of model (3) in -plane with the parameter values , , , , , , , , and . The dark blue region (0) represents the case where there is no equilibrium, the light blue regions (1a, 1b) represent the case where there is a unique equilibrium, the orange region (2) represents the case where there are two steady states, and the brown region (3) represents the case where there is three equilibria. (b) shows the nullclines of the model which has up to four steady states: tumour-free steady state “”; tumour dormancy steady state “”; medium concentration tumour steady state “”; and escape tumour steady state “.”

3.1. Stability of Tumour-Free Steady State

In this subsection, we investigate the parameter ranges for which the tumour-free steady state is locally asymptotically stable. The Jacobian matrix of the linearized system at the tumour-free steady state is given by

with the eigenvalues , , and . Hence, the tumour-free steady state is locally stable if if and only if

This relates (the growth rate of the tumour cells) to the (the resistance coefficient), which measures how efficiently the immune system competes with the tumour cells. If this tumour-free steady state is unstable, then no amount of chemotherapy will be able to completely eradicate the tumour cells.

3.2. Stability of Lethal Steady States

The same analysis done above shows that the deadly steady state has the eigenvalues , , and and hence it is unstable saddle point, while the other deadly steady state can be either stable or unstable depending on the model parameters and the value of the time-delay . This can be shown by using Routh Harwatz test and Rouche’s theorem as shown in detail in the previous chapters. Since the stability of this steady state is not needed for the developing treatment therapy, we will not introduce more details in this part.

3.3. Stability of Coexisting Steady States

To study the stability of the coexisting steady states, we vary the two parameters (the immune cells growth rate) and (the normal flow rate of immune cells), with fixing the other parameters: , , , , , , , , and . Table 1 summarizes the existence and stability results of the coexisting steady states as present in different regions of Figure 2. It shows that the light blue region (1a) represents the “escape” case where there is a unique stable node steady state with high tumour size, while the light blue region (1b) represents the case where there is a unique steady state with low tumour size. It is stable spiral for , while at the limit cycle occurs due to Hopf bifurcation. Furthermore, the orange region (2) represents the case where there are two steady states; one is stable node and the other is unstable sacdle node. To this end, the brown region (3) represents the case where there are three steady states: one is stable node, one is unstable node, and the last steady state is spiral stable for , while the limit cycle occurs at . Of interest are the existence and stability of steady states where a small tumour population size might coexist with a large number of normal cells. Figure 3 presents the phase space for the cell populations in the case where and . It shows that, for , the steady state is asymptotically stable (Figure 3(a)), while, for , a limit cycle is born around the steady state (Figure 3(b)). We utilize MIDDE code [27] to solve the DDEs model, which is suitable to simulate stiff and nonstiff problems, using monoimplicit RK methods [28]. We next consider the chemotherapy treatment () in the underlying model and establish the existence of an optimal control for the model and provide necessary conditions for the optimal control.

Table 1: The stability results for the coexisting steady states by using the functions and , while fixing the rest of the parameters as mentioned in the text.

Region

in Figure 2

Steady state

Eigenvalues  

, ,
Stability  

for

Light blue

(1a)
0.1
0.05
(0.0269, 486.9244, 0.9999)

Stable node

Light blue

(1b)
1.4
0.1
(1.0299, 0.0238, 0.9999)
,  

Stable spiral for ,   

stable limit cycles at

Orange

(2)
0.2
0.23
(0.1648, 419.9672, 0.9999)

Stable node

(0.8656, 79.7712, 0.99)

Unstable saddle node

Brown

(3)
0.6
0.1
(0.0789, 461.688, 0.99)

Stable node

(0.7236, 148.7078, 0.99)

Unstable saddle node

(1.0298, 0.0623, 0.99)
,  

Stable spiral for ,  

stable limit cycles at

(a)

(b)

Figure 3: The phase space for the cell populations in the case where and . (a) shows that, for , the steady state is asymptotically stable. (b) shows that, for , a limit cycle is born around the steady state.

4. Optimal Control Problem Governed by DDEs

Once a suitable model of interacting cell populations is constructed, we then focus on the design of an efficient treatment protocol, where we employ the tools of optimal control theory.

Consider the optimal control problem with pure state constraints and control bounds, as follows:

subject to the DDEs

with the control constraint

and state constraint

is called objective functional and the integrand is called the Lagrangian of objective functional. Furthermore, and are called the lower and upper bounds. The function is called an admissible control if and only if it fulfills inequality constraints (14d). The set of all admissible controls is called the admissible set and we referred to it by (where “ad” stands for the admissible). The state enters with a delay as in the system of the state equations (14b) while it is evaluated at the time as in the objective functional (14a). The set of all admissible states , which satisfy the state equations and the state constraint, is called the set of admissible state.

The goal of chemotherapy is to eradicate the tumour cells, while maintaining adequate amounts of healthy tissue. From a mathematical point of view, adequate destruction to tumour cells might mean forcing the system out of the basin of an unhealthy spiral node, out of a limit cycle, and into a basin of attraction of a stable tumour-free equilibrium. Alternatively, if the therapy pushes the system into a limit cycle in which the size of the tumour is small for a long period of time (as long as the life of the patient, for example), this could also be considered a “cure.”

Optimality in treatment might be defined in a variety of ways. Some studies have been done in which the total amount of drug administered is minimized, or the number of tumour cells is minimized. The general goal is to keep the patient healthy while killing the tumour. Since our model takes into account the toxicity of the drug to all types of cells, our control problem consists of determining the function that will maximize the amount of effector cells and minimize the number of tumour cells and the cost of the control with the constraint that we do not kill too many normal cells. If the units of cells are normalized, so that the carrying capacity of normal cells is 1, we then require that the number of normal cells stays above three-fourths of the carrying capacity. Therefore, our main objective is to optimize the functional

which subject to the underlying DDEs

with control constraint

and state constraint

Here, is a weight factor that describes the patient’s acceptance level of chemotherapy. We choose the control parameter as a class of piecewise continuous functions defined for all such that , where represents the maximum chemotherapy, while represents the case where there is no chemotherapy. Thus, we depict the class of admissible controls as

We next prove the existence of the optimal solution of the underlying system (15a)–(15g).

5. Existence of an Optimal Solution

To prove the existence of the optimal solution of (15a)–(15g), we use the results of Fleming and Rishel [29, Theorem 4.1, pages 68-69] and Lukes [30, Theorem 9.2.1, page 182].

Theorem 5. There exists an optimal solution for the optimal control problem (15a)–(15g) such that

where if the following conditions are satisfied.

(1)The set of admissible state is nonempty.(2)The admissible set is nonempty, convex, and closed.(3)The right-hand side of the state system is bounded by a linear combination of the state and control variables.(4)The integrand, , of the objective functional is a concave on .(5)There exist constants , , and such that .

Proof. In order to verify the above conditions, we should first prove the existence of the solution for system of the state equations (15b)–(15e). Since , and, by neglecting the negative terms in the model, we have

System (18) can be rewritten in the form

where . This system is linear in finite time with bounded coefficients. Then the solutions of this linear system are uniformly bounded. Therefore, the solution of the nonlinear system (15b)–(15e) is bounded and exists [30]. Hence, condition one is satisfied.

Clearly, the second condition is satisfied by the definition of . System (15b)–(15e) is bilinear in the control variable and can be rewritten as

where , , and and are the vector valued functions of and , respectively. Using the fact that the solutions are bounded, we have

where depends on the coefficients of the system, and

We also note that the integrand of is concave in . Finally,

where depends on the upper bounds of and , and . This completes the proof.

We also conclude that there exists an optimal control variable .

6. Optimality Conditions

In this section, we establish the necessary conditions for the optimal solution of the optimization problem (15a)–(15g); we use Pontryagin’s minimum (maximum) principle derived by Göllmann et al. [31] for the retarded optimal control problem with mixed control-state constraints. To this end, we define the augmented Hamiltonian function involving the inequality constraints by

where

and () are the adjoint variables that satisfy

Here denotes the indicator function of the interval and defined by

To minimize the Hamiltonian functional, the Pontryagin’s minimum principle [31] is used. Thus, we arrive at the the following theorem.

Theorem 6. Let be the optimal solutions of (15a)–(15g), where . Then, there exists an adjoint state , defined by (26), such that the triple satisfies the state equation

with the initial conditions

the adjoint state equations

with transversality conditions

and the optimal control

Proof. The optimal control can be solved from the optimality condition ; that is, . By using the handedness of the control set , it is easy to obtain in the form of (32).

7. Immunochemotherapy

Model (8) is extended to include external source of immunotherapy treatment of the effector cells such as ACI. We then add the term to represent the input rate of externally administered antitumour effector cells, where is the control parameter. Our goal is to maximize an objective functional subject to the new model with a combination of chemotherapy and ACI and constraints on the control and the state:

subject to DDEs

the control constraints

and the state constraint

where is a weight factor that describes a patient’s acceptance level of immunotherapy and the set of all admissible controls is defined by

Similarly, the optimal solution of the optimization problem (33a)–(33g) satisfies the state equations

The adjoint state equations are

with the transversality conditions , , and the minimum condition

When (without immunotherapy), system (35)–(37) reduces to system (28)–(32).

Remark 7. In the case of immunotherapy alone (), the objective functional becomes

8. Numerical Simulations of the Optimal Control System

Numerical simulations leading to the approximation of the optimal controls (35)–(37) are carried out using the forward Euler method for the state system and backward difference approximation for the adjoint system. We assume the step-size , such that and , where . We define the state, adjoint, and control variables at the mesh points. An initial guess is given for the controls and , which are then updated continuously until the objective functional satisfies the conditions. However, there are several major problems to be overcome when solving delay differential equations. These include stability, stiffness, and discontinuities in the right-hand side of the equation. Stability and stiffness can be handled by the correct choice of implicit solvers [27]. The delay terms can create a whole suite of discontinuities; see [32, 33].

We choose a different set of parameter values (in stable and unstable regions). In the current simulations, we vary the three parameters , , and , and fix the other parameters:

We solve the optimality system to determine the optimal control situation (i.e., the drug strategy) and predict the evolution of the tumour cells, effector cells, and normal cells of each control strategy in 30 days.

Figure 4 shows the numerical simulations of the state system before and after chemotherapy treatment using optimality system (28)–(32) when , , and (in the stable region). We note that, in the presence of chemotherapy with optimal control, the effector cells population grows up significantly, while the tumour cells population decreases and is totally eradicated after 20 days. In the meantime, the normal cells population remains over 75%. Yet, Figure 5 shows the impact of chemotherapy treatments (with optimal control) when we choose the parameter values in an unstable region (, , and ). The tumour and effector cells populations are oscillating over time in the absence of chemotherapy, while the presence of treatment helps the immune system to keep the growth of the tumour cells under its control.

Figure 4: Simulations of system (28)–(32), in the stable region, before and after the treatments with control with the initial conditions , , and and the parameter values are given in the text.

Figure 5: Simulations of the system (28)–(32), in an unstable region, before and after the chemotherapy treatment with the control and initial conditions , , and and the parameter values are given in the text.

Figure 6 presents the evolution of system (35)–(37) in the case of combination of chemotherapy and ACI. The parameters values are chosen in the stable region. We notice that the tumour cells population can be eradicated after day 12 which is faster compared to the results of Figure 4 when we used the chemotherapy alone. In other words, the numerical results show that using the combination immunochemotherapy is more effective than using chemotherapy treatment alone.

Figure 6: Simulations of system (35)–(37), in the stable region, before and after the immunochemotherapy treatments with controls. It shows that the tumour cells population can be eradicated in day 12.

However, Figure 7 shows evolution of the system with only immunotherapy (i.e., without chemotherapy). We may notice from the figure that this case reflects the best therapeutic strategies for treatment of tumour, where the recovery becomes faster with high dosage of immunotherapy where can reach the value of 3.5 level compared with the combination it was in level 2.

Figure 7: Simulations of the tumour cells population of system (35)–(37), before and after immunotherapy with control. It shows that the tumour cells can be eradicated at day 7 with high immunotherapy dosage where the control value reaches the value of 3.5.

9. Concluding Remarks

In this paper, we provided a delay differential model with control variables that describe the interactions of immune cells, tumour cells, normal cells, and immunochemotherapy treatment with control variables. A pontryagin-type maximum principle is derived, for retarded optimal control problems with delays in the state variable when the control system is subject to a mixed controlstate constraint, in order to minimize the cost of treatment, reduce the tumour cells load, and keep the number of normal cells above 75% of its carrying capacity. We presented an efficient numerical technique, based on forward difference approximation to the state system and backward difference scheme to the adjoint system, to solve the optimal control problem and identify the best treatment strategy when we adopt the chemotherapy treatment alone or a combination of chemoimmunotherapy, with minimum side effects. The numerical results show and confirm that the optimal treatment strategies reduce the tumour cells load and increase the effector cells after few days of therapy. The performance of combination therapy protocol was better than the standard protocol of chemotherapy alone. The numerical simulations show the rationality of the model presented, which in some degree meets the natural facts.

This work can be extended to more sophisticated problems with delays in both state and control variables, when the control system is subject to a mixed controlstate constraints.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by UAE University. The authors would like to thank the referees and Professor Loredana G. Marcu for their valuable and contractive comments that improved this paper.

References

B. Joshi, X. Wang, S. Banerjee, H. Tian, A. Matzavinos, and M. A. J. Chaplain, “On immunotherapies and cancer vaccination protocols: a mathematical modelling approach,” Journal of Theoretical Biology, vol. 259, no. 4, pp. 820–827, 2009.

D. Kirschner and J. C. Panetta, “Modeling immunotherapy of the tumor—immune interaction,” Journal of Mathematical Biology, vol. 37, no. 3, pp. 235–252, 1998.

V. A. Kuznetsov, I. A. Makalkin, M. A. Taylor, and A. S. Perelson, “Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis,” Bulletin of Mathematical Biology, vol. 56, no. 2, pp. 295–321, 1994.

L. G. de Pillis, W. Gu, and A. E. Radunskaya, “Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations,” Journal of Theoretical Biology, vol. 238, no. 4, pp. 841–862, 2006.

L. G. de Pillis, K. R. Fister, W. Gu et al., “Optimal control of mixed immunotherapy and chemotherapy of tumors,” Journal of Biological Systems, vol. 16, no. 1, pp. 51–80, 2008.

R. P. Araujo and D. L. S. McElwain, “A history of the study of solid tumour growth: the contribution of mathematical modelling,” Bulletin of Mathematical Biology, vol. 66, no. 5, pp. 1039–1091, 2004.

N. Bellomo, N. K. Li, and P. K. Maini, “On the foundations of cancer modelling: selected topics, speculations, and perspectives,” Mathematical Models and Methods in Applied Sciences, vol. 18, no. 4, pp. 593–646, 2008.

H. M. Byrne, T. Alarcon, M. R. Owen, S. D. Webb, and P. K. Maini, “Modelling aspects of cancer dynamics: a review,” Philosophical Transactions of the Royal Society of London A, vol. 364, no. 1843, pp. 1563–1578, 2006.

M. A. J. Chaplain, “Modelling aspects of cancer growth: insight from mathematical and numerical analysis and computational simulation,” in Multiscale Problems in the Life Sciences, vol. 1940 of Lecture Notes in Mathematics, pp. 147–200, Springer, Berlin, Germany, 2008.

M. L. Martins, S. C. Ferreira Jr., and M. J. Vilela, “Multiscale models for the growth of avascular tumors,” Physics of Life Reviews, vol. 4, no. 2, pp. 128–156, 2007.

J. D. Nagy, “The ecology and evolutionary biology of cancer: a review of mathematical models of necrosis and tumor cell diversity,” Mathematical Biosciences and Engineering, vol. 2, no. 2, pp. 381–418, 2005.

T. Roose, S. J. Chapman, and P. K. Maini, “Mathematical models of avascular tumor growth,” SIAM Review, vol. 49, no. 2, pp. 179–208, 2007.

R. Yafia, “Dynamics analysis and limit cycle in a delayed model for tumor growth with quiescence,” Nonlinear Analysis: Modelling and Control, vol. 11, no. 1, pp. 95–110, 2006.

R. Yafia, “A study of differential equation modeling malignant tumor cells in competition with immune system,” International Journal of Biomathematics, vol. 4, no. 2, pp. 185–206, 2011.

R. Eftimie, J. L. Bramson, and D. J. D. Earn, “Interactions between the immune system and cancer: a brief review of non-spatial mathematical models,” Bulletin of Mathematical Biology, vol. 73, no. 1, pp. 2–32, 2011.

F. A. Rihan and D. H. A. Rahman, “Delay differential model for tumor-immune dynamics with HIV infection of CD4⁺ T-cells,” International Journal of Computer Mathematics, vol. 90, no. 3, pp. 594–614, 2013.

F. A. Rihan, M. Safan, M. A. Abdeen, and D. H. Abdel Rahman, “Qualitative and computational analysis of a mathematical model for tumor-immune interactions,” Journal of Applied Mathematics, vol. 2012, Article ID 475720, 19 pages, 2012.

F. A. Rihan, D. H. Abdel Rahman, S. Lakshmanan, and A. Alkhajeh, “A time delay model of tumour—immune system interactions: global dynamics, parameter estimation, sensitivity analysis,” Applied Mathematics and Computation, vol. 232, pp. 606–623, 2014.

L. G. de Pillis and A. Radunskaya, “The dynamics of an optimally controlled tumor model: a case study,” Mathematical and Computer Modelling, vol. 37, no. 11, pp. 1221–1244, 2003.

G. W. Swan, “Optimal control applications in the chemotherapy of multiple myeloma,” IMA Journal of Mathematics Applied in Medicine and Biology, vol. 2, no. 3, pp. 139–160, 1985.

F. Castiglione and B. Piccoli, “Cancer immunotherapy, mathematical modeling and optimal control,” Journal of Theoretical Biology, vol. 247, no. 4, pp. 723–732, 2007.

D. Kirschner and A. Tsygvintsev, “On the global dynamics of a model for tumor immunotherapy,” Mathematical Biosciences and Engineering, vol. 6, no. 3, pp. 573–583, 2009.

M. Villasana and A. Radunskaya, “A delay differential equation model for tumor growth,” Journal of Mathematical Biology, vol. 47, no. 3, pp. 270–294, 2003.

M. Bodnar, “The nonnegativity of solutions of delay differential equations,” Applied Mathematics Letters. An International Journal of Rapid Publication, vol. 13, no. 6, pp. 91–95, 2000.

A. Halanay, Differential Equations: Stability, Oscillations, Time Lags, Academic Press, New York, NY, USA, 1966.

L. G. De Pillis and A. Radunskaya, “A mathematical tumor model with immune resistance and drug therapy: an optimal control approach,” Journal of Theoretical Medicine, vol. 3, no. 2, pp. 79–100, 2001.

F. A. Rihan, E. H. Doha, M. I. Hassan, and N. Kamel, “Numerical treatments for Volterra delay integro-differential equations,” Computational Methods in Applied Mathematics, vol. 9, no. 3, pp. 292–308, 2009.

C. T. H. Baker, G. Bocharov, A. Filiz et al., “Numerical modelling by delay and volterra functional differential equations,” in Proceeding of the Topics in Computer Mathematics and its Applications, pp. 113–137, LEA, Athens, Greece, 1999.

W. H. Fleming and R. W. Rishel, Deterministic and Stocha stic Optimal Control, Springer, New York, NY, USA, 1975.

D. L. Lukes, Differential Equations: Classical to Controlled, Academic Press, New York, NY, USA, 1982.

L. Göllmann, D. Kern, and H. Maurer, “Optimal control problems with delays in state and control variables subject to mixed control-state constraints,” Optimal Control Applications & Methods, vol. 30, no. 4, pp. 341–365, 2009.

F. A. Rihan, Numerical treatment of delay differential equation [Ph.D. thesis], University of Manchester, Manchester, UK, 2000.

F. A. Rihan and B. F. Rihan, “Numerical modelling of biological systems with memory using delay differential equations,” Applied Mathematics & Information Sciences, vol. 9, no. 3, pp. 1–14, 2015.