STABILITY OF LIMIT CYCLE IN A DELAYED MODEL FOR TUMOR IMMUNE SYSTEM COMPETITION WITH NEGATIVE IMMUNE RESPONSE

This paper is devoted to the study of the stability of limit cycles of a system of nonlinear delay differential equations with a discrete delay. The system arises from a model of population dynamics describing the competition between tumor and immune system with negative immune response. We study the local asymptotic stability of the unique nontrivial equilibrium of the delay equation and we show that its stability can be lost through a Hopf bifurcation. We establish an explicit algorithm for determining the direction of the Hopf bifurcation and the stability or instability of the bifurcating branch of periodic solutions, using the methods presented by Diekmann et al.


Introduction
We consider in this paper the model which provides a description of tumor cells in competition with the immune system.This description is described by many authors, using ordinary and delayed differential equations to model the competition between immune system and tumor, in particular see [20,24,25].Other similar models provide a description of the modelling, analysis, and control of tumor immune system interaction, see, [17,26,29].
Other authors use kinetic equations to model the competition between immune system and tumor.Although they give a complex description in comparison with other simpler models, they are, for example, needed to model the differences of virulence between viruses, see [1][2][3][4][5]10].Several other fields of biology use kinetic equations, for instance [12,13] give a kinetic approach to describe population dynamics, [2] deals with the development of suitable general mathematical structures including a large variety of Boltzmann-type models.
The reader interested in a more complete bibliography about the evolution of a cell, and the pertinent role that has cellular phenomena to direct the body towards the recovery 1. Kinetic scheme describing interactions between ECs and TCs (see [20]).
or towards the illness, is addressed to [6,21].A detailed description of virus, antivirus, body dynamics can be found in the references [8,14,27,28].
The mathematical model with which we are dealing was proposed in a recent paper by Gałach [20].In this paper, the author developed a new simple model with one delay of tumor immune system competition, this idea is inspired from [25] and he recalled some numerical results in [25] in order to compare them with those obtained in his paper [20].

Mathematical model
The model proposed in [25] describes the response of effector cells (ECs) to the growth of tumor cells (TCs).This model differs from others because it takes into account the penetration of TCs by ECs, which simultaneously causes the inactivation of ECs.It is assumed that interactions between ECs and TCs in vitro can be described by the kinetic scheme shown in Figure 2.1, where E, T, C, E * , and T * are the local concentrations of ECs, TCs, EC-TC complexes, inactivated ECs, and "lethally hit" TCs, respectively, k 1 and k −1 denote the rates of bindings of ECs to TCs and the detachment of ECs from TCs without damaging them, k 2 is the rate at which EC-TC interactions program TCs for lysis, and k 3 is the rate at which EC-TC interactions inactivate ECs.
Kuznetsov and Taylor model is as follows: where s is the normal (i.e., not increased by the presence of the tumor) rate of the flow of adult ECs into the tumor site, F(C,T) describes the accumulation of ECs in the tumor site, d 1 , d 2 , and d 3 are the coefficients of the processes of destruction and migration for E, E * , and T * , respectively, a is the coefficient of the maximal growth of tumor, and b is the environment capacity.
In [25], it is claimed that experimental observations motivate the approximation dC/dt ≈ 0. Therefore, it is assumed that C ≈ KET, where K = k 1 /(k 2 + k 3 + k −1 ), and the model can be reduced to two equations which describe the behavior of ECs and TCs only.

Radouane Yafia 3
Moreover, in [20], it is suggested that the function F should be in the following form: F(C,T) = F(E,T) = θET.Therefore, the model (2.1) takes the form where α 1 = θ − m, and a, b, s have the same meaning as in (2.1) All coefficients except α 1 are positive.The sign of α 1 depends on the relation between θ and m.If the stimulation coefficient of the immune system exceeds the neutralization coefficient of ECs in the process of the formation of EC-TC complexes, then α 1 > 0. We use the dimensionless form of model (2.2): where x denotes the dimensionless density of ECs, y stands for dimensionless density of the population of TCs, α = a/Kk 2 T 0 , β = bT 0 , δ = d/Kk 2 T 0 , σ = s/nE 0 T 0 , and ω = α 1 /n is the immune response to the appearance of the tumor cells, E 0 and T 0 are the initial conditions.In [20], the author studies the existence, uniqueness, and nonnegativity of solutions and he shows the nonexistence of nonnegative periodic solution of system (2.3).
The delayed mathematical model corresponding to (2.3) is given by the following system [20]: where the parameter τ is the time delay which the immune system needs to develop a suitable response after the recognition of nonself cells (see [20]).Time delays in connection with the tumor growth also appear in [6,7,9,18,19].
The existence and uniqueness of solutions of system (2.4) for every t > 0 are established in [20], and in the same paper it is shown that (1) if ω ≥ 0, these solutions are nonnegative for any nonnegative initial conditions (biologically realistic case), (2) if ω < 0, there exists nonnegative initial condition such that the solution becomes negative in a finite-time interval.
Our goal in this paper is to consider the case (2) when the immune response is negative (i.e., ω < 0) with the following conditions: αδ > σ and α 2 (βδ − ω) 2 + 4αβσω > 0. We give some local stability results for the unique nontrivial equilibrium P 2 and that it undergoes a Hopf bifurcation (see [32]).We establish a systematic criteria for determining the direction of Hopf bifurcation and the stability or instability of the bifurcating branch of periodic solutions.The case (1) when the immune response is positive (i.e., ω > 0) is treated in [30,31].
This paper is organized as follows.In Sections 3 and 4, we establish some results on the stability of the nontrivial steady states of the delayed system (2.4).The existence of a critical value of the delay in which the non-trivial steady state changes stability is investigated.Based on the Hopf bifurcation theorem, we show the occurrence of Hopf bifurcation as the delay crosses some critical value of the parameter delay (see [32]).The main result of this paper is given in Section 5, we establish a systematic criteria for determining the direction of Hopf bifurcation and the stability or instability of the bifurcating branch of periodic solutions.
Let u = x − x 2 and v = y − y 2 .By linearizing system (2.4) around the nontrivial equilibrium point P 2 , we obtain the following linear system: The characteristic equation of (3.2) has the form where p = δ + αβy 2 > 0, r = δαβy 2 > 0, s = −ωy 2 > 0, and q = αωy 2 (1 − 2βy 2 ).The stability of the equilibrium point P 2 is a result of the localization of the roots of the equation then we have the following theorem.
Theorem 3.1.Assume αδ > σ, α > 0, and β > 0 are close enough to 0. Then, there exists τ l > 0 such that P 2 is asymptotically stable for τ < τ l and unstable for τ > τ l , where For the proof of Theorem 3.1, we need the following lemma.
Lemma 3.2 [11].Consider the equation where p, r, q, and s are real numbers.
From the hypothesis αδ > σ, we deduce that q + r > 0. Therefore, the hypotheses (H 1 ), (H 2 ) of Lemma 3.2 are satisfied.Then all roots of the characteristic (3.3) have negative real parts for τ = 0 and the steady state P 2 is asymptotically stable for τ = 0.By Rouche's theorem, it follows that the roots of (3.3) have negative real parts for some critical value of the delay τ.
We want to determine if the real part of some root increases to reach zero and eventually becomes positive as τ varies.If iζ is a root of (3.3), then Separating the real and imaginary parts, we have The two roots of the above equation can be expressed as follows: ), the sign of r 2 − q 2 is deduced from the sign of (δβ − ω 2 (1 − 2βy 2 )) = (2αβδ − √ Δ)/α which is negative (because β is very small and α > 0).Therefore, r 2 − q 2 < 0 and the hypothesis (H 3 ) of Lemma 3.2 is satisfied.From Lemma 3.2, the unique solution of (3.10) has the following form: and there exists a unique critical value such that the equilibrium point P 2 is asymptotically stable for τ ∈ [0,τ l ) and unstable for τ > τ l .For τ = τ l , the characteristic (3.3) has a pair of purely imaginary roots ±iζ l .
In Section 4, we will study the occurrence of Hopf bifurcation when the delay passes through the critical value of the delay τ = τ l .

Hopf bifurcation occurrence
According to the Hopf bifurcation theorem [22], we come to the main result of this paper.

Direction of Hopf bifurcation
For determining the direction of Hopf bifurcation, there exist many formulas, we cite (1) the formulas using the theory of normal forms, see Hassard et al. [23], (2) Faria and Magalhães [16], (3) Diekmann [15].In this section, we follow methods presented in [15], where the direction and stability of the bifurcating branch are obtained by the Taylor expansion of the delay function τ that describes the parameter of bifurcation near the critical value τ l (see Sections 3 and 4).Namely, this direction and stability are determined by the sign of the first nonzero term of Taylor expansion, that is, and the sign of τ 2 determines either the bifurcation is supercritical (if τ 2 > 0) and periodic orbits exist for τ > τ l , or it is subcritical (if τ 2 < 0) and periodic orbits exist for τ < τ l .The term τ 2 may be calculated, see [15], using the formula where M is the characteristic matrix of (3.14) given by D 2 M(iζ l ,τ l ) denotes the derivative of M with respect to τ at the critical point (iζ l ,τ l ), the constant c is defined as follows: qD 3 1 f 0,τ l P 2 (θ),P(θ) + qD 2  1 f 0,τ l e 0. M −1 0,τ l D 2 1 f 0,τ l P(θ),P(θ) ,P(θ) + 1 2 qD 2 1 f 0,τ l e 2iζl.M −1 2iζ l ,τ l )D 2 1 f 0,τ l P(θ),P(θ) ,P(θ) , (5.4) where f is the nonlinear part of (3.14), D i 1 f , i = 2,3, denotes the ith derivative of f with respect to ϕ, P(θ) denotes the eigenvector of A, P(θ) denotes the conjugate eigenvector, and p and q are defined later.Now, we describe all the above operators and vectors precisely.Let denote the linear part of (3.14).Using the Riesz representation theorem, one obtains, see [22], where, δ(•) denotes the Dirac function.