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

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.


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][2][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 immunedepleting 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][7][8][9][10][11][12][13][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][16][17], the authors take into account 2 Computational and Mathematical Methods in Medicine 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 predatorprey 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.

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 tot ( ) = ( ) + ( ). The rate of binding of ECs to TCs and the rate of separation of ECs from TCs without damaging them are denoted by 1 and −1 , respectively. The rate at which EC-TC integrations program for lysis is denoted by 2 while the rate at which EC-TC interaction inactivate ECs is denoted by 3 . 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), F( ( ), ( )) = F( ( ), and ( )) > 0 (when ( ) > 0) describes the accumulation of ECs in the tumour side due to the presence of the tumour. However, 1 , 2 , and 3 are the coefficients of the processes of destruction and migration of , * , and * , respectively. The maximal growth of tumour is represented by the coefficient , and −1 is the environment capacity. If we assume that ( )/ ≈ 0, therefore ( ) ≈ ( ) ( ) where = 1 /( −1 + 2 + 3 ), 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 (0) = 0 and (0) = 0 . 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 F( , ) may take the form F( , ) = ( ) ( )/( + ( )) 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 Computational and Mathematical Methods in Medicine 3 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 ( ) = 1 ( ) and ( ) = 2 ( ), for all ∈ [− ,0]. 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 ≥ 0 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]).  Therefore, we arrive at the following proposition. )] , 2 = max ( 1 , (0)) .

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

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

Tumour cells T(t) Effector cells E(t)
Normal cells N(t) Chemotherapy u(t)

Attack
We assume that the drug kills all types of cells but that the killing rate differs for each type of cells; ( ) = (1 − − ) is the fraction cell kill for a given amount of drug, ( ), at the tumour site. We denote by 1 , 2 , and 3 the three different response coefficients. V( ) represents the amount of dose that is injected into the system, while 1 is the decay rate of the drug once it is injected. In this case, the quantity we will control directly is not ( ) but V( ). 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: = 2 identifies the parameter associated with tumour, and = 3 identifies the one associated with the normal tissue; and 2 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 − 1 and − 2 . Let C = C([− , 0], R 4 ) be the Banach space of continuous functions mapping the interval [− , 0] into R 4 with the topology of uniform convergence. It is easy to show that there exists a unique solution ( ( ), ( ), ( ), ( )) of system (8) with initial data ( 0 , 0 , 0 , 0 ) ∈ C. For biological reasons, we assume that the initial data of system (8) satisfy 0 ≥ 0, 0 ≥ 0, 0 ≥ 0, and 0 ≥ 0. For = 0, 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 2 = 1.
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, ( ) = 0, for all .

Drug-Free Steady States and Their Stability
In the absence of chemotherapy ( ( ) = 0), model (8) and * is a nonnegative solution for * + ( Computational and Mathematical Methods in Medicine 5 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.

Stability of Tumour-Free Steady State.
In this subsection, we investigate the parameter ranges for which the tumourfree steady state E 0 is locally asymptotically stable. The Jacobian matrix of the linearized system at the tumour-free steady state is given by with the eigenvalues 1 = − < 0, 2 = 2 − / − 1 , and 3 = − 3 < 0. Hence, the tumour-free steady state E 0 is locally stable if 2 < 0 if and only if This relates 2 (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.

Stability of Lethal Steady States.
The same analysis done above shows that the deadly steady state ( / , 0, 0) has the eigenvalues 1 = − < 0, 2 = 2 − / , and 3 = 3 > 0 and hence it is unstable saddle point, while the other deadly steady state ( ( * ), * , 0) 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.

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: = 0.2, = 0.3, = 0.003611, 2 = 1.03, 3 = 1, = 2 × 10 −3 , = 1, 1 = 0.00003, and 2 = 3 × 10 −9 . 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 = 1.4 and = 0.1. It shows that, for = 0.8, the steady state is asymptotically stable (Figure 3(a)), while, for = 1.2, 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 ( ( ) > 0) in the underlying model and establish the existence of an optimal control for the model and provide necessary conditions for the optimal control.

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: 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 V( ) 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 ad (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 ad , 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 V( ) 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 max V∈ ad which subject to the underlying DDEs and state constraint Here, V 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 0 ≤ V( ) ≤ V max < ∞, where V( ) = V max represents the maximum chemotherapy, while V( ) = 0 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).

Existence of an Optimal Solution
To prove the existence of the optimal solution of (15a)-(15g), we use the results of Fleming    (1) The set of admissible state is nonempty.
(2) The admissible set ad 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.
, of the objective functional is a concave on ad .
(5) There exist constants ℎ 2 , ℎ 3 > 0, and > 1 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 ( − )/( + ( − )) < , V max < ∞ and, by neglecting the negative terms in the model, we have System (18) 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 ad . System (15b)-(15e) is bilinear in the control variable V and can be rewritten as where ℎ 1 depends on the coefficients of the system, and We also note that the integrand of (V) is concave in ad . Finally, where ℎ 2 depends on the upper bounds of ( ) and ( ), and ℎ 3 = V /2. This completes the proof.
We also conclude that there exists an optimal control variable V * .

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 H ( , , , , , , V, ) To minimize the Hamiltonian functional, the Pontryagin's minimum principle [31] is used. Thus, we arrive at the the following theorem.

Immunochemotherapy
Model (8) is extended to include external source of immunotherapy treatment of the effector cells such as ACI.
We then add the term ( ) 1 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 ad is defined by ad = { (V, ) : (V, ) piecewise continuous, such that Similarly, the optimal solution of the optimization problem (33a)-(33g) satisfies the state equations * ( ) = + * ( − ) * ( − ) The adjoint state equations are

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 stepsize ℎ, such that = ℎ and − 0 = ℎ, where ( , ) ∈ N 2 . We define the state, adjoint, and control variables at the mesh points. An initial guess is given for the controls V 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 = 0.5, = 0.01, and = 1.2 (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 ( = 0.2, = 0.2, and = 1.5).  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 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.
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.

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.