Analysis of tumour-immune evasion with chemo-immuno therapeutic treatment with quadratic optimal control

ABSTRACT A simple mathematical model for the growth of tumour with discrete time delay in the immune system is considered. The dynamical behaviour of our system by analysing the existence and stability of our system at various equilibria is discussed elaborately. We set up an optimal control problem relative to the model so as to minimize the number of tumour cells and the chemo-immunotherapeutic drug administration. Sensitivity analysis of tumour model reveals that parameter value has a major impact on the model dynamics. We numerically illustrate how does these delay can change the stability region of the immune-control equilibrium and display the different impacts to the control of tumour. Finally, epidemiological implications of our analytical findings are addressed critically.


Introduction
Many mathematical models have been developed to describe the immunological response to infection with different types of biological models, for example human immunodeficiency virus (HIV), H1N1 influenza-A, tumour models and so on [16][17][18][19]21]. Modelling tumour-immune interaction has attracted much attention in the last decades. This interaction is very complex and mathematical models can help to shape our understanding of dynamics this biological phenomenon [1,3,4,[6][7][8]15,22,23,27,29]. Much research has focussed on how to enhance the anti-tumour activity, by stimulating the immune system with vaccines or by direct injection of T cells or cytokines [24,25]. For instance, the mathematical model of Krischner and Panetta [15], which also focuses on the tumour-immune interaction, indicates that the dynamics among tumour cells, immune cells and the cytokine interleukin(IL-2) can explain both short-term oscillations in tumour size as well as long-term tumour elapse. Of course, the development of powerful cancer immunotherapies requires first an understanding of the mechanisms governing the dynamics of tumour growth. Accordingly, a great research effort is being devoted to understand the interaction between the tumour cells and the immune system.
The interaction between cellular populations of tumour and an immune system is, from an ecological point of view, competitive, and involves many events and molecules. Thus CONTACT P. Krishnapriya priyaprithu1205@gmail.com it has a high degree of complexity implying that the immune system is not able to eliminate a neoplasm in all cases. For this reason, it is desirable to strengthen the immune system after an immune-depleting course of chemotherapy. Chemotherapy is one of the most prominent cancer treatment modalities. However, it is not always a comprehensive solution for tumour regression. Chemotherapy depletes a patient's immune system, making the patient prone to dangerous infections. The goal of this paper is to model mathematically, analyse and explore computationally potentially optimal ways to combine chemo-immunotherapy treatment strategies that can minimize a tumour while maximizing the strength of the immune system, with minimal toxicity to the patient. We formulate and analyse a delay differential model describing immune response and tumour cells under the influence of chemotherapy alone and the combinations of chemo-immunotherapy.
In this study, Wilson et al. [26] examined the both vaccine and TGF-β inhibitors were given, the model predicts that the tumour size will reach its peak on day 5 and tumour eradication will occur on day 21. The conclusion of the experimental study in Wilson et al., was that TGF-β with vaccine treatment can able to eradicate the tumour level. Our mathematical model highlights just one possible way of combining tumour treatments to promote tumour eradication through an immune response. Within the framework of modelling interacting populations by systems of ordinary differential equations in [26], as follows: with initial conditions Although the addition of TGF-β to the model indicates qualitatively parallel behaviour with the original model in [15], several important quantitative differences also occur. For using the treatment with IL-2 can be clear the tumour, and the immune system grows without bounds causing a side effect. Hence, by adding the time delay effect 'τ ' in effector cells and immune system,the uncontrolled growth of the immune system situation is under control [16]. Hence our modified mathematical model as follows: where T(t), B(t), R(t), E(t) and I(t) represent the populations of tumour size, TGF-β concentration, regulatory T cells, effector cells and IL-2 (Interleukin-2), respectively. The model presented is a stiff system of differential equations and an appropriate non-dimensional scaling is essential for numerical accuracy. The equations are nondimensionalized as follows: After dropping the over-bar notation for convenience, the system described as follows: These scaling need to be chosen to help adjust for the fact that this is a numerically 'stiff' system. That is, without scaling, or inappropriate scaling, the numerical routines used to solve these equations will fail. This is due to very large changes in some of the variables over very short ranges of time. The parameter values are described in Table 1. The organization of this paper is as follows: After this introduction, we develop the model and to show that the non-negativity and boundedness of solutions and also local stability, analysing the existence of local Hopf bifurcation through the tumour-free steady state for tumour-immune system in Section 2. We study the optimal control problem governed by delay differential equations (DDEs) with only control variable. Existence of the solution and optimality condition are also discussed in Section 3. We investigate the sensitivity analysis in Section 4. We examine the stability results and mathematical findings for the dynamical behaviour of the tumour-immune model with optimal control are also numerically verified in Section 5 through MAPLE and MATLAB packages. Finally we give short conclusion in Section 6.

Non-negativity and boundedness
We denote by C the Banach space of continuous functions ϕ : The initial condition of system (4) is given as To show that the solutions of model (4), are bounded and remain nonnegative in the domain of its application for sufficiently large values of time 't', we recall the following lemma. Using the above Lemma 2.1, we derived the following propositions for solving the boundedness of solution (4).
The details of the proofs are found in Appendix 1.

Stability analysis
The first equilibrium is the trivial state where all the populations are zero, namely I 0 = (0, 0, 0, 0, 0). The eigenvalues of the Jacobian matrix for I 0 is 0. Therefore, I 0 is always unstable saddle point. Consider the another equilibrium is tumour-free state, namely and as follows: We obtained the characteristic equation of the above systems as follows: Now, where and The details of our analysis are given in Appendix 2.

Theorem 2.3: Consider a coefficient parameterized quasi polynomial
The details of the proofs are found in Appendix 2. From the above analyses, we conclude that the stability of equilibria is important from physiological view point. Similar analyses can be performed using any of the system parameters in order to determine conditions for the appearance or disappearance of equilibria and to determine equilibrium stability. Hence our proposed model (4) with the scaled values, the tumour-free equilibrium is always locally stable only when τ < τ * , otherwise it is unstable. Then, only in order to better determine under what circumstances the tumour can be eliminated, thus we implement optimal control problem.

Epidemic model with control
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 cell is minimized. The general goal is to keep the patient healthy while killing the tumour. In the context of mathematical modelling in cancer growth with chemo-immunotherapy, it is essential to frame an optimal control problem so that the total amount of drug used is minimized. While, we minimize drug doses because we assume that toxic side-effects are a concern, and that the smaller the dose, the better. We minimize an objective functional of a form that reflects the trade-off we require in minimizing both tumour size and drug-doses: (10) J, which involves a 'quadratic control' because it is quadratic in the treatment terms, must be minimized which subject to system, where ρ(t) and η(t) are control constraint Here, B ρ and B η are weight factor. The function ρ(t) is control describing the percentage of adoptive cellular immunotherapy given. B ρ is a weight factor that describes a patients acceptance level of immunotherapy and B η is a weight factor that describes a patient's acceptance level of chemotherapy. We choose as our control class piecewise continuous functions defined for all 't' such that 0 ≤ ρ ≤ 1 where ρ(t) = 1 represents maximal immunotherapy and ρ(t) = 0 represents no immunotherapy and η(t) = 1 represents maximal chemotherapy and η(t) = 0 represents no chemotherapy.

Analysis of the super solutions
To prove the existence of the optimal solution of (10)- (11) if the following conditions are met: (1) The set of all admissible state is non-empty.
(2) The admissible set U is non-empty, convex and closed.
(3) The right-hand side of the state system is bounded by a linear combination of the state and control variables.
The details of the proofs are given in Appendix 3.
Since we have the existence of an optimal control triple, we next determine the necessary conditions associated with it via Pontryagin's Maximum Principle.

Necessary conditions for optimality
In this section, we establish the necessary conditions for the optimal solution of the optimization problem (10) and (11), we use Pontrygian's minimum (maximum) principle is derived by and λ i , i = (1, 2, 3, 4, 5, 6) are the adjoint variables that satisfy Here χ [0,t f −τ ] denotes the indicator function of the interval [0, t f − τ ] and defined by To minimize the Hamiltonian functional, the Pontryagian's minimum principle [12] is used. Thus, we arrive at the following theorem.

Sensitivity to parameter changes
In this section, we show the sensitivity analysis with respect to the parameter is considered. We would like to consider how a small shift in the parameters would change the stability of the tumour-free equilibrium for this model. It is quite usual for a model to display high sensitivity to small variations in some parameters, while displaying robustness to variations in other parameters. In a more recent report [2], Baker and Rihan formally derive sensitivity equations for DDE models, as well as the equations for the sensitivity of parameter estimates with respect to observations. Now, we consider a linearized system (6) of parameter dependent DDEs with vector parameter w i = [δ 0 , a 0 , d, p 1 , g 1 , r, q 1 , q 2 , μ 2 , δ 1 ] T ∈ R 5 , for i = 1, 2, . . . , 10, given by The corresponding sensitivity of system (6), with respect to the parameter 'δ 0 ' is as follows: The corresponding sensitivity of system (6), with respect to the parameter 'a 0 ' is as follows: The corresponding sensitivity of system (6), with respect to the parameter 'p 1 ' is as follows: The corresponding sensitivity of system (6), with respect to the parameter 'r' is as follows: r (t, r), The corresponding sensitivity of system (6), with respect to the parameter 'δ 1 ' is as follows: We may observe that a small change in a 0 , δ 0 , δ 1 , r and p 1 can produce the significant changes in the level of tumour cells. The parameter of model (6) is perturbed both positive and negative of their base case values to determine the effect on the output solutions. Figure 1(a-e) shows that the sensitivity of tumour cell population u 1 , due to small perturbation in a 0 , δ 0 , δ 1 , r and p 1 . We notice that from Figure 1(b), is insensitive with increasing the value of a 0 into the earlier interval, after some time it become highly sensitive. But the other parameters except that a 0 , are highly very sensitive with increasing their level.

Applications with numerical simulations
In this section, we have discussed the dynamical behaviour of the systems, we have analysed, graphically. Numerical simulations are carried out using MAPLE and MATLAB packages. We have also tried to show a comparative study between the systems with no therapy, with chemotherapy and with chemo-immunotherapy for tumour-immune evasion system. We present some numerical results of system (4), supporting the theoretical analysis. Using the value of Table 1, we consider the following system Clearly the positive tumour-free equilibrium is I 1 = (0, 0, 0, 0.000246207455, 0.005) T . From (7), we obtain which has only one positive real roots λ = 0.01945999975 and any other roots have negative real part. Thus, ω 0 = √ λ = 0.1394991030437. For different τ values, we plot the characteristic equation, which is illustrated by numerical simulation in Figure 2(a,b) with initial value (10 5 , 10 5 , 10 5 , 10 5 , 10 5 ) T .  (Table 1).  Figure 2(a,b) shows that the characteristic equation of (25) has negative real parts. Then the system (24) is always stable in tumour-free equilibrium.
The optimal system has been solved numerically and the results have been presented graphically. There are two systems of DDEs, the first system (11) being the state equations involving the control and the second (16) being the adjoint equations λ i 's (i = 1,2, . . . ,5). An initial guess was made for λ i 's (i = 1, 2, . . . , 5) gives an initial guess for the control. From here the state equations were solved using the initial condition. Our findings leading to the approximation of the optimal controls (11) are carried out using the forward Euler method for the state system and backward difference approximation for the adjoint system. We assume that the step size h, such that τ = mh and t f − t 0 = nh, where (m, b) ∈ N 2 . 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 overcome when solving DDEs. We choose a set of parameter values are taken from Table 2. We solve the optimality system to determine the optimal control situation (i.e., drug strategy), and predict the evolution of the system had taken each control strategy in 10 days. Figure 3(a-c) shows results of our simulations in the three treatment regimes along with the corresponding experimental data for Mouse data (Using Table 2). Figure 3(a) shows that tumour size level of no treatment therapy. In this case the tumour growth becomes immediately uncontrolled. Using chemotherapy with tumour-immune system, the tumour growth became to control for 6 days, after that the tumour growth immediately uncontrolled which can be shown from Figure 3(b). Finally, we show the both chemo-immuno therapy treatment of our model, the system becomes control within 8 days only. After, the tumour-immune system grows up quickly becomes uncontrolled. Figure 4(a-c) shows results of our simulations in the three treatment regimes along with the corresponding experimental data for human data (Using Table 2). Figure 4(a) shows that tumour size level of no treatment therapy. In this case the tumour growth becomes immediately uncontrolled. Even though, using chemotherapy with tumour-immune system, the tumour growth became to uncontrolled, which it can be shown from Figure 4(b). Finally, we show the both chemo-immuno therapy treatment for our model, the system becomes control within 7 days only. After, the tumour-immune system grows up quickly becomes uncontrolled.    Table 2 with B ρ = 1 and B η = 2 for mouse and human data.
To compare the tumour growth in all treatment regimes. It is clear that while monotherapy results in a slowing down of the tumour growth, the tumour is still able to escape immuno surveillance and grow uncontrolled. Only in the case of dual therapy is the immune system able to eradicate the tumour within 8 days.

Conclusion
We have examined a model incorporating interacting tumour and immune cell populations and their responses to chemo-immunotherapy treatment. The dynamics of the system without treatment reveal two equilibrium points for a specific parameter case: trivial equilibrium and tumour-free equilibrium. We presented the non-negativity and boundedness of solutions, existence of steady states of our model. The immune system inhibitory effects (such as blocking IL-2 production and inhibiting antigen-specific T-cell activation) and tumour-stimulating effects (such as increasing blood supply to tumour cells in order to enhance the tumour's ability to metastasize) of TGF-β provide explanation for enhanced tumour growth and failure of the host immune system. In order to counteract this occurrence, the model was extended to include a novel therapeutic strategy using the chemo-immunotherapy treatment without TGF-β cells. For using the both therapies level, from the very beginning the tumour controlled within 8 days. After some time, the tumour can be grow up uncontrolled.

Disclosure statement
No potential conflict of interest was reported by the authors.
If we assume that H 1 > 0 and H 2 > 0, then the above Equation (A12) has no positive root. In fact, it is observed that, has no positive real root by Descarte's rule of sign. Thus if, H 1 > 0 and H 2 > 0 then there is no ν such that iν is an eigenvalue of the characteristic equation (7), that is, λ will never be a purely imaginary root of Equation (7). Thus all the real parts of all eigenvalues of (7) are negative for all τ ≥ 0.

Proof of Theorem 2.3.:
Now if any one of H 1 and H 2 is negative then (y) = 0 and hence Equation (A12) has positive root ν 0 . This implies that the characteristic equation (7) has a pair of purely imaginary roots ±iν 0 . From (A10), we have Now, we determine sign (d Re(λ)/dτ )| τ =τ * where sign is the signum function and Re(λ) is a real part of λ. By using the following mathematical calculation we can say that the tumour-free steady state of model (6) remains stable for τ < τ * and Hopf bifurcation occurs when τ = τ * . Differentiating (7) with respect to τ , we get which implies, Therefore, We also note that the integrand of J(ρ, η) is concave in U. Hence where m 1 depends on the upper bounds of B(t), R(t), E(t) and I(t), then m 2 = (B ρ + B η )/2. This completes the proof.