Numerical Simulations of a Delay Model for Immune System-Tumor Interaction

In this paper we consider a system of delay differential equations as a model for the dynamics of tumorimmune system interaction. We carry out a stability analysis of the proposed model. In particular, we show that the system can have up to two steady states: the tumor free steady state, which always exist, and the tumor persistent steady state, which exists only when the relative rate of increase of the tumor cells exceeds the ratio between the natural proliferation rate and the relative death rate of the effector cells. We also determine an upper bound for the delay, such that stability is preserved. Numerical simulations of the system under different parameter values are performed.


Introduction
he immune system (IS) is a very complex one.It is composed of a complex network of different cells which operate collectively by communicating information through signalling.The duty of these cells is to keep our system natural, i.e., free of foreign entities.When a foreign entity enters our system, particular immune system cells raise an alarm and send signals through the network to other immune system cells calling for an attack on the foreign entity.Once the foreign substance is located and identified as non-self, the immune system plans an attack strategy.Among the most deadly foreign entities are cancer cells.Cancer cells are difficult to deal with because of their ability to multiply so fast that the IS cannot keep up.Another property of cancer cells is that they are able to camouflage the antigens (the substance that triggers the IS) such that the IS cannot recognise them as foreign cells.
Once cancer cells are identified as foreign, the IS starts an attack to destroy them.The interactions between tumour cells and the immune system are very complex and require sophisticated models to describe them.Mathematical models, based on ordinary differential equations, delay differential equations or partial differential equations, have proven to be useful tools in analysing and understanding the IS-tumor interactions.Several mathematical models have been suggested to describe the interactions between tumour and immune system [1-9].

T
Many of these papers consider the tumor-IS interaction with other factors such as treatment.In this paper, we consider a mathematical model of IS-tumor interactions based on a system of delay differential equations.The delay is introduced to reflect the non-instantaneous outcome of the interaction between the immune system cells and the tumor cells.We carry out a stability analysis of the system and we simulate the system under different parameter values to explore the different asymptotic behaviors.
Many models that have appeared in the literature are based of the well-known Kuznetsov and Taylor model [1] which we describe in the next section.

Kuznetsov and Taylor model
The Kuznetsov and Taylor mathematical model to describe the tumor-immune system interaction is based on the following: 1.The growth of tumor cells population T (in the absence of the immune system cells) follows a logistic model: T and E cells which are bound to die.Kuznetsov and Taylor's complete model is given by the following set of differential equations where, again, are the concentrations at the tumor site of the immune system cells (effector cells), the tumor cells, the effector-tumor cells complex, the inactived effector cells, and the inactived tumor cells, respectively.The above model was simplified based on the fact that the T E  bond, which is formed at a rate of 1 k , lasts a relatively short period of time [10]  and/or by adding additional considerations (see [3]- [8], and references therein).
In [1], the authors analyzed model (3) with a Holling Type II form for F .In non-dimensionl variables, their model takes the form  [6], the authors considered Kuznetsov and Taylor's model (3) with Holling Type I response function , with and without delay, where the delay was introduced to reflect the delay in the response of the immune system before proliferating effector cells.
In this work, we consider (3) with a Holling Type I response function with delay introduced in both equations as described in the next section.

A delay model
As in [6], we consider ) , ( T E F of Holling Type I, but we consider the delay  according to the following biologically possible scenario.We assume that there is delay in both the proliferation of effector cells and in the disposal of tumor cells.In dimensionless form, the model we consider is where ) (t x and ) (t y are the nondimensional density at the tumor site of the effector cells and of the tumor cells, respectively.All parameters in (5) are nonnegative real numbers and are approximated (see [1] and [6]) by: The delay factor in model ( 5) models the natural delays in both the proliferation of effector cells and in the disposal of tumor cells.This is biologically meaningful as the IS takes time to react and migrate effector cells to the tumor site.Also, in the second equation of model (5), it models the delay in completely eradicating the tumor cells as they will not die immediately upon interacting with effector cells.
Our aim in studying the above model ( 5) is to investigate the effect, if any, of the time delay,  , in the disposal of the tumor cells.We carry out the stability analysis of the equilibrium states of the model ( 5) and determine necessary conditions for the local asymptotic stability.Also, we determine an upper bound for the delay parameter  such that stability is preserved.Numerical simulations of the model are performed to investigate the asymptotic behavior of the system.

Steady states
The steady states, ) , ( y x , of ( 5) are the solutions of the nonlinear homogeneous system Then depending on the parameter values, there may exist up to two more steady states ) , (

y x E
, where 1 x and 1 y are given by (8).
The value  / relative to  is very important as it determines how many steady states the system can asymptotically assume.The value  / measures the rate of increase of the density of the effector cells and  measures the rate of increase of the tumor cells.Thus if , this biologically means that the effector cells will take over the tumor cells, which as a result, are bound to die in the long run.On the other hand, if , the tumor cells will always exist.Now, we turn to the stability analysis of the steady states 0 E and 1 E .

Stability analysis
We know that a steady state is locally asymptotically stable if all the eigenvalues of the Jacobian matrix of the system at the steady state have negative real parts.
Before we continue, we recall an important stability criterion known the Mikhailov criterion [9]: Mikhailov criterion: Let () Ps and () Qs be two polynomials with deg( has no roots on the imaginary axis, then all roots of () We have the following proposition.Proposition 1 Before we proceed, we state and prove the following Lemma about the positivity of 0 C , which will be used later.
, we can write A as x and 1 y (Eq. ( 8)) we have .= 2

 
We have the following proposition for the stability of 1 , where 0  is given by (24).

Proof. Let
, and 1 E would be locally asymptotically stable, provided that ) ( W has no roots on the imaginary axis.Now, we check for pure imaginary roots of The solutions of (31) are formally written as .In this case, the positive solution of (32) is has only two pure imaginary roots 0 js  , where 0 s is given by (33).With 0 s as in (33), the first positive value of  such that 0 = ) (s has no imaginary roots, and 1 E is locally asymptotically stable.□

Numerical simulations
In this section, we perform a number of numerical simulations to reveal the dynamics of system (5).We consider the two cases (i)


which was proved to be asymptotically stable up to an upper bound for the delay  (see Proposition 1 and Eq. ( 21)).The  is calculated to be 0 = 0.0880838  .To confirm these results, we simulated the system for = 0, 0.01, 0.05, 0.08, 0.09  and various initial conditions.The results are displayed in Figures 6 and 7, where Figure 6 displays the phase portraits for = 0, 0.01, 0.05, 0.08, 0.09  and Figure 7 displays the solutions of the system for = 0.05, 0.08, 0.09  and initial condition (0) = 1, (0) = 30 xy .It is to be noted that for 0 = 0, 0.01, 0.05, 0.08 <  , the solutions converge to 1 E , confirming the local asymptotic stability of 1 E .However, in the case 0 = 0.9 >  , we notice the formation of the closed orbit around 1 E .This also can be seen in the right- most plot in Figure 7, where the effector and tumor cells oscillate around 1 E .This confirms the existence of a bifurcation with respect to parameter  as it crosses 0 = 0.0880838  .
It should be also mentioned that for certain initial conditions the solutions of the system become negative in finite time.

Conclusion
A mathematical model for the interaction between the immune system cells and tumor cells was proposed and studied.The model was derived from Kuznetsov and Taylor's model [1] by introducing a delay term which reflects the slow response of the immune system cells to the presence of tumor cells.The importance in studying this delay model is in determining the effect of the delay term on the dynamics of the system.A stability analysis of the equilibrium states was carried out.Numerical simulations of the system were done for a number of different initial conditions.The simulation results confirm the theoretical stability results.Numerically, it was also observed that as the delay term  crosses the upper bound 0  , the steady state 1 E loses its stability, resulting in the formation of closed orbits around 1 .E This is an important result in the sense that as the IS gets weaker (slower response, i.e., larger delay), the size of the tumor cells persists in an oscillating fashion (closed orbits).
first equation of(7) we have the quadratic equation for y 7) has a positive root, since the functions roots, then all its roots have negative real parts if and only if


Squaring both sides of (29) and (30) and adding, we get It is clear that a sufficient condition for the right hand side of (32) to be positive is

. 3 Case 2 :
and 0.3 , respectively.We remark that these upper bounds are valid only when the solutions()   xt and () yt remain non-negative.In fact, simulations reveal that the solutions become negative for much smaller values of  , The outcomes of the simulations for this case are displayed in Figures2, 4, and 5.In Figure2, we display the phase portraits of system (5to E .In Figures4 and 5, the solutions of system (5) are displayed with = 0.1  and = 0., respectively, each with = 0.05  .Both figures confirm the convergence of the solutions to the steady state 0 E .It is important to mention that for > 0.05  , under certain initial conditions, the solutions become negative at some time * t at which the system becomes biologically invalid.For the case all parameters as in (6) and simulate the system with different delay values  .In this case, we have two steady states,