STABLE PERIODIC OSCILLATIONS IN A TWO-STAGE CANCER MODEL OF TUMOR AND IMMUNE SYSTEM INTERACTIONS

This paper presents qualitative and bifurcation analysis near the degenerate equilibrium in a two-stage cancer model of interactions between lymphocyte cells and solid tumor and contributes to a better understanding of the dynamics of tumor and immune system interactions. We first establish the existence of Hopf bifurcation in the 3-dimensional cancer model and rule out the occurrence of the degenerate Hopf bifurcation. Then a general Hopf bifurcation formula is applied to determine the stability of the limit cycle bifurcated from the interior equilibrium. Sufficient conditions on the existence of stable periodic oscillations of tumor levels are obtained for the two-stage cancer model. Numerical simulations are presented to illustrate the existence of stable periodic oscillations with reasonable parameters and demonstrate the phenomenon of long-term tumor relapse in the model.

1. Introduction.According to the study in Boyle et al. [3], every year millions of people suffer with cancer and die from this disease throughout the world.In the new century cancer still remains one of the most dangerous killers of humankind.In order to study the progress of cancer and seek effective treatment strategies, researchers have proposed comprehensive approaches, including biological, computational and mathematical methods, to study the disease.An important aspect in studying cancer is to understand the interactions between the immune system and tumor cells.In order to simulate the host's own immune response to destroy and eliminate tumor cells, various types of mathematical models have been proposed, see Albert et al. [1], de Pillis et al. [6], d'Onofrio [7], Kirschner and Panetta [15], Kirschner and Tsyvintsev [16], Kuznetsov et al. [17], Lejeune et al. [19], Liu et al. [22], and Swan [30].We refer to a recent survey by Eftimie at al. [9] and the references cited 348 DAN LIU, SHIGUI RUAN AND DEMING ZHU therein for more references on modeling tumor and immune system interactions, in particular using ordinary differential equations.
In this work we are interested in studying mathematical models of carcinogenesis in which the tumor cells elicit an immune response.Based on some reasonable hypotheses, the following model was proposed by DeLisi and Rescigno [5] to qualitatively estimate the function of the immune surveillance (Lefever and Garay [18]), where L(t) and C(t) denote respectively the number of free lymphocytes and the total number of tumor cells at time t, C and C f are the total number of free tumor cells and the number of free cells on a tumor surface (i.e., not bounded by lymphocytes), respectively.The values of λ 1 and α 1 may be dependent on environmental stimuli and the type of tissues where the tumor is growing.The last factor in the first equation 1 − L (b) The death rate of the stage-2 lymphocytes is a first-order process with constant α 3 .(c) Stage-1 lymphocytes are transformed into stage 2 at a constant rate λ 1 .(d) Stage-1 lymphocytes are produced at a fixed rate λ 1 L 0 (in the absence of tumor cells) plus a rate proportional to the product of the number of free tumor cells and stage-2 lymphocytes.(e) In the absence of lymphocytes, the reproduction rate of the tumor cells is at a constant rate λ 2 .(f) The killing rate of the tumor cells by lymphocytes is proportional to the frequency of interaction with stage-2 lymphocytes with constant α 2 .(g) There is an equilibrium relation between free tumor cells and stage-2 lymphocytes.(h) The progeny of L 2 must first pass through stage one before becoming mature lymphocytes to approximate some biological delay in the formation of L 2 by tumor cell stimulation.Now the interactions between the tumor cells and lymphocytes are described a system of three ordinary differential equations (see [26]): where α 1 exp(−L 2 /L c ) represents the saturation of the system which restricts the maximum number L c of lymphocytes, (C − C f )/C f L 2 is a constant K under the assumption that the relation between stage-2 lymphocytes and free tumor cells is equilibrium controlled.Therefore, C f = C/(1 + KL 2 ) and ( 2) is written as Some potential oscillation phenomena were observed by Rescigno and DeLisi [26] in the development of cancer and lymphocyte cells.In this paper we study the nonlinear dynamics of system (3) and prove that there is a stable periodic solution bifurcated from the interior degenerate equilibrium under certain assumptions of the parameters, which means that all the trajectories in the neighborhood of this equilibrium spiral towards the bifurcated limit cycle as time increases.More precisely, by applying the Hopf bifurcation theorem (Zhang [32]) for 3-dimensional systems, we first discuss the Hopf bifurcation in model ( 3) and obtain the nonexistence of degenerate Hopf bifurcation.The noticeable point lies in the determination of conditions on the stability of the bifurcated periodic orbit so that this cancer model exhibits stable periodic oscillations between the solid tumor cells and lymphocyte cells.Therefore, the tumor levels will oscillate inescapably once the initial values of tumor cells and lymphocyte cells are close enough to the degenerate equilibrium.This oscillatory phenomenon of tumor levels has been observed clinically and is known as Jeff's Phenomenon (Thomlinson [31]).The theoretical results in this work will be helpful for the further study of the dynamical behaviors in the mathematical models of carcinogenesis and for the better understanding of the development of cancer when the tumor cells elicit an immune response.
The paper is organized as follows.Section 2 provides the qualitative analysis, the existence of Hopf bifurcation and the nonexistence of degenerate Hopf bifurcation.In section 3, an applicable transformation of a 3×3 Jacobian matrix is given.Section 4 discusses the stability of the bifurcated periodic orbit and presents the sufficient conditions on the stable case.Numerical simulations supporting the theoretical results of section 4 are presented in section 5. Section 6 provides a brief discussion on this work.

Hopf bifurcation.
For system (3), we make the following transformations: The above system has two possible fixed points: α1yB exp( yB yc )).The Jacobian matrix of the linearized system of (4) is From the biological point of view, the state space of system (4) is the non-negative octant R 3 + = {X = (x, y, z) T ∈ R 3 : x 0, y 0, z 0}, where T means the transpose of a matrix or a vector throughout this paper.Denote the positive open subset of R 3 + as R 3 p = {X = (x, y, z) T ∈ R 3 : x > 0, y > 0, z > 0}.If we set R 0 = α3λ2 α2λ1x0 , then the following conclusion can be obtained directly from [26].
Lemma 2.1.If R 0 < 1, then for any solution (x(t), y(t), z(t)) of (4), we have that and the boundary equilibrium point A is the unique and globally stable equilibrium of (4) in R 3 + .From Lemma 2.1 we know that for any choice of initial conditions, the trajectory of system (4) always approaches to A when R 0 < 1.In other words, cancer will be controlled and tumor cells will be annihilated eventually as long as the increasing rate of lymphocyte cells is large enough.
Point B exists as a nonnegative equilibrium only if R 0 1.In this case, point A turns to be a saddle with a one-dimensional unstable manifold from a stable node, which implies that for each orbit starting in R 3 p the number of cancer cells will not tend to zero.A bifurcation occurs at R 0 = 1.
Remark 1.The role of R 0 is similar to the basic reproduction number in epidemic models.Now we intend to study the qualitative properties of (4) near the equilibrium B. The Jacobian matrix of (4) at B is and the corresponding characteristic equation is where 1+yB > 0. Applying Routh-Hurwitz criterion to the Jacobian matrix J B , we achieve that the sufficient condition for the stability of B is Nevertheless, equation (7) always has at least one negative real root γ no matter what the sign of a 1 a 2 − a 3 is.Also γ is the principle eigenvalue, and the absolute value of γ is greater than that of the real parts of the other eigenvalues.Under the assumption a 1 a 2 − a 3 = 0, it is clear to see that J B has one negative eigenvalue −(α 3 + λ 1 ) and two pure imaginary eigenvalues ±iw, where w > 0 satisfies Although the stability of B is not determined, Hopf bifurcation may occur at this equilibrium.We resort to the following Hopf bifurcation theorem in three-dimensional differential systems to determine the occurrence of Hopf bifurcation which can be found in Zhang [32].Lemma 2.2.Let Ω ⊆ R 3 be an open set containing O(x 1 , y 1 , z 1 ) and let S ⊆ R be an open set with 0 ∈ S. Let f : Ω × S → R 3 be an analytic function such that f (O, µ) = 0 for any µ ∈ S. Assume the variational matrix Df (O, µ) of f has one real eigenvalue γ(µ) and two conjugate imaginary eigenvalues α(µ) ± iβ(µ) with γ(0) < 0, α(0) = 0, β(0) > 0. Furthermore, suppose that the eigenvalues cross the imaginary axis with nonzero speed, that is, dα dµ (0) = 0. Then the following differential system undergoes Hopf bifurcation near the equilibrium point O at µ = 0.
Remark 2. If α (0) > 0 and O is a stable but not asymptotically stable equilibrium when µ = 0, then the solutions of system (8) on some surface in the neighborhood of O are all periodic orbits.
In order to clarify how system (4) undergoes Hopf bifurcation at point B, we choose the death rate α 3 of stage-2 lymphocytes as a perturbed parameter.The equation a 1 a 2 − a 3 = 0 brings out where we need 0 We also need to determine the sign of the real part of dλ/dµ at µ = 0 when the above equation is valid.Differentiating both sides of ( 7) with respect to µ, we have where µ in x B is omitted for simplicity.Therefore, by ( 9), we have where x = α3λ2 α2λ1 .Obviously, the degeneracy Re( dλ dµ )| µ=0 = 0 occurs at when h(x 0 ) vanishes (see Figure 1).
According to Lemma 2.2 and the above calculations, we obtain the following conclusion.
Theorem 2.3.For any given x 0 ∈ [0, x] and x 0 = x * 0 , system (4) undergoes nondegenerate Hopf bifurcation at the nonhyperbolic equilibrium B when the parameters satisfy equation (9).When x 0 = x * 0 , the real part of dλ dµ becomes zero, then the transversality condition fails.The occurrence of degeneracy may induce degenerate Hopf bifurcation or no Hopf bifurcation.We will continue to differentiate (10) associating to µ to determine the existence of degenerate Hopf bifurcation.Using (10), it can be easily deduced after tedious and mechanical calculations that (12) On the assumption of (9) and the condition x 0 = x * 0 , it follows that Re( dλ dµ | µ=0 ) = 0 and sign{Re( )w+ where the last equation comes from h(x 0 ) = 0 in (11), which implies that λ 1 y B (α Three cases of α(µ) when |µ| > 0 is small enough.
Theorem 2.4.For parameters satisfying equation (9), there is no Hopf bifurcation for system (4) at the nonhyperbolic equilibrium B when x 0 = x * 0 .From Figure 2 we can see that if x 0 = x * 0 , then when µ passes through the origin from any direction, the interior equilibrium is always locally stable.In terms of the notations in Golubitsky and Langford [12], α (0) = 0 is equivalent to a µ (0, 0) = 0. Since at the interior equilibrium of model (4) a µ (0, 0) = 0 holds but α (0) < 0, the normal form of a(x 2 , µ) corresponds to the cases (2) or (7) or (10) in Theorem 3.19 of [12] for codimension less than 3.However, by the stability analysis, the corresponding bifurcation diagram should be displayed as in Fig. 4.2 of that paper.
3. Normalized form of the Jacobian matrix.In this section we provide an applicable transformation to normalize the Jacobian matrix.Consider the following where 3 ) for any µ ∈ S and the variational matrix of f = (f 1 , f 2 , f 3 ) T at point X * is denoted as A = (a ij ) 3×3 .Assume A has one negative real eigenvalue γ and a pair of purely imaginary eigenvalues ±iw (w > 0) as µ = 0. Then the following lemma gives the normalized form of A. Lemma 3.1.([4, pp.106])Let A be a real n × n matrix.Then there exists a real nonsingular matrix P such that Ã = P −1 AP has the real canonical form consisting of real square matrices A 1 , . . ., A k , B 1 , . . ., B m down the main diagonal.Each A j has the form where 0 2 is the 2 × 2 zero matrix, I 2 is the 2 × 2 identity matrix, and Then B j has the form .
We can validate this lemma by considering the canonical form of the abovementioned 3 × 3 matrix A for µ = 0 in explicit transformations.In fact, set  Consider the first two equations of (17), we can solve one of the nonzero solutions as The following choice of a series of linear transformations is motivated by the methods in Lu and Luo [23].Choose then we obtain that Choose and set T = T 1 T 2 T 3 , then we obtain The lemma is proved.
4. Stability of periodic orbits.In order to determine the stability of periodic orbits bifurcated from the interior equilibrium of system (4), we apply the general formula in Hsü and Kazarinoff [14], which turns out to be more applicable than the original bifurcation formulas given by Hopf [13] and Friedrichs [11], and the center manifold theorem in Perko [24].Firstly we need to outline this method and some notations.Consider the following class of autonomous differential systems: where (20) In virtue of Lemma 3.1, there exists a nonsingular matrix P which can transform A into the normal form with the Jordan block as A j and B j .For the sake of simplicity, in the following discussion we suppose A has the normal form as Hopf [13] obtained the bifurcation theorem on the existence of a periodic analytic solution p(t, ) with period T ( ) of system (26) under the above assumptions, where the parameter is related to µ by an analytic function µ = µ( ) such that µ(0) = 0, p(t, 0) = 0, T (0) = 2π/w and p(t, ) = 0 for all sufficiently small | | > 0. Furthermore, these periodic solutions exist for exactly one of three cases: µ > 0, µ < 0, or µ = 0.By Hopf's theorem, the characteristic exponents of this bifurcated periodic solution p(t, ) are the eigenvalues of the eigenvalue problem where V (t) has the same period T ( ) and L(t, ) = A(µ( )) + F X (p(t, ), µ( )).
The interested readers are referred to Hopf [13] for more details.
Based on the expressions ( 24), ( 30) and an algebraic computation, Hsü and Kazarinoff [14] derived a formula for the direction of the bifurcation of periodic orbits of system (19).Lemma 4.1.( [14]) If µ = 0, the direction of the bifurcation of (4) is determined by the following equation, The stability of the periodic orbit bifurcated from the degenerate equilibrium can be determined by the following lemma.Lemma 4.2.( [14]) Under the assumptions in Hopf 's theorem, if A(0) has exactly two purely imaginary eigenvalues and the other n − 2 eigenvalues have negative real part and if µ 2 Re(α (0)) > 0, then a bifurcating periodic solution whose existence is asserted by Hopf 's theorem is asymptotically orbitally stable with asymptotic phase; however, if µ 2 Re(α (0)) < 0, these small periodic solutions are unstable.Moreover, if any one of the other n − 2 eigenvalues has positive real part, then the bifurcating periodic solution is orbitally unstable.Now let us return to determine the stability of limit cycles bifurcated from the equilibrium B of system (4).It is necessary first to transform the Jacobian matrix J B expressed in (6) into the normal form.For simplicity, let us set where b and c are both positive, then J B can be written as Suppose system (4) undergoes Hopf bifurcation at point B, then the parameters in J B satisfy bc = (α 3 − a)(λ 1 + α 3 ) (33) and J B has one negative eigenvalue −δ = −(λ 1 + α 3 and a pair of purely imaginary eigenvalues ±wi, where It is easy to check that the inverse matrix of P turns out to be and L(0) P −1 J B P becomes the normal form as the right side in (18), where = w(w 2 + δ 2 )/b.
In the following a series of linear transformations will be introduced for system (4).Let then the positive equilibrium B is translated into the origin of the system: where Q i 1 (i = 1, 3) are C ∞ power series of y 1 and z 1 with power higher than 3, f yy , f yz , f yyy , f yyz and h yy , h yz , h yyy , h yyz are the corresponding partial derivatives at point B.
If we denote and system (34) is changed into ) are C ∞ power series of x 2 , y 2 , z 2 with power higher than 3. Figure 5 shows three solution orbits near the equilibrium point B on different initial conditions in R 3 as µ = 0.0051, where initial values are chosen as (a) (0.1, 1.3, 0.3) (outside the limit cycle), (b) (0.194, 1.184, 0.27) (on the limit cycle), and (c) (0.09, 1.1, 0.1) (inside the limit cycle), respectively.For the sake of convenience, we denote the limit cycle by L. It can be observed that the degenerate equilibrium point B becomes a stable equilibrium point in this case.Figure 6 presents the projections of the periodic solution curve on the coordinate planes, where (a), (b) and (c) represent the projections on the xy-plane, yz-plane, xz-plane, respectively.Figure 7 corresponds to diagrams between the coordinate components and time t of the solution curve in Figure 6, where (a), (b) and (c) represent the curves of x(t), y(t), z(t), respectively.6. Discussion.Periodic oscillations with various periods have been observed in some cancer data (Fortin and Mackey [10]) and oscillations are common in the immune system (Stark et al. [29]).So it is reasonable to expect that the oscillatory phenomenon occurs in tumor and immune system interaction models.Rescigno and DeLisi [26] proposed a three equation model (3) to describe two different stages of lymphocytes in the interactions of the solid tumor and immune system.The oscillations existing in the signs of the derivatives of L 1 , L 2 and C subject to the time t were observed.
In this paper we established the existence and stability of periodic oscillations in the two-stage cancer model (3) by choosing the parameters suitably.After making a series of transformations of the variables, the qualitative analysis and some bifurcation results near the degenerate equilibrium were given.We have shown that the system could undergo Hopf bifurcation in the small neighborhood of the interior degenerate equilibrium.It is valuable to find out that any degenerate Hopf bifurcations cannot occur in this model, which excludes more complex dynamical behaviors.
The conditions on the stability of the bifurcated periodic orbits guarantee the appearance of stable periodic oscillations of the tumor levels.When a stable periodic orbit exists, it can be seen as "safe" in the neighborhood of this closed orbit since the trajectories originating from there will overwhelmingly go towards it.Therefore,  8 in terms of the time t when a 1 a 2 − a 3 < 0.
the solid tumor and the immune system can coexist for a long term although the cancer is not eliminated.The oscillatory dynamics in the tumor and immune system interaction models demonstrate the phenomenon of long-term tumor relapse and have been observed in some related tumor and immune system models (d'Onofrio [7], Eftimie at al. [9], Kirschner and Panetta [15], Kuznetsov et al. [17], and Lejeune et al. [19]).We can interpret this situation from the biological point of view that while the immune system fights with cancer in the host, there exists a balance between them because of the periodic changes in the internal tissues and the external circumstances such that they coexist in a bounded region.It is necessary to mention at the end that our discussions on the equilibrium point B are based on the hypothesis (9).If we choose parameters in (4) such that the left hand side term a 1 a 2 − a 3 in ( 9) is positive, then B turns to be an asymptotically stable equilibrium point, which is similar to the first case of the numerical simulations in the last section.If a 1 a 2 − a 3 < 0, then we find that B is an unstable nondegenerate equilibrium by numerical simulations (see Figures 8-9).In this case, the cancer cells will increase uncontrollably because their development are not controlled by the immune system any longer.
It will be very interesting to include immunotherapy in the tumor and immune system interaction model (Smyth et al. [28]), study the effects of adoptive cellular immunotherapy on the model, and explore conditions under which the tumor can be eliminated (de Pillis et al. [6], Kirschner and Panetta [15]).We leave this for future consideration.
Finally, we would like to mention that in model (3) the development from the immature stage to the mature stage of lymphocytes was modeled by two ordinary differential equations in L 1 (t) and L 2 (t) explicitly.As pointed out by the reviewer, this process can be described by a time delay τ > 0. Thus, we would like to propose the following two-stage cancer model described by two differential equations with a time delay: where the time delay τ > 0 is the time for immature lymphocytes to develop into mature ones and L(t) represents the density of the mature lymphocytes at time t.Recently, delayed models of tumor and immune response interactions have been studied extensively, we refer to d'Onofrio et al. [8], Bi and Ruan [2] and the references cited therein.In particular, Bi and Ruan [2] have shown that various bifurcations, including Hopf bifurcation, Bautin bifurcation and Hopf-Hopf bifurcation, can occur in such models.It would be interesting to consider the nonlinear dynamics of the delayed model (45).

Figure 6 .Figure 7 .
The projections of the periodic orbit L on the (a) xy-, (b) yz-, and (c) xz-planes, respectively when µ = 0.0051 and x 0 = x * 0 .The plots of the components (a) x(t), (b) y(t), and (c) z(t), respectively, of the periodic orbit in terms of the time t when µ = 0.0051 and x 0 = x * 0 .

Figure 8 .Figure 9 .
Figure 8. Trajectories spiral away from the unstable equilibrium point B when a 1 a 2 − a 3 < 0.