The Dynamical Behaviors for a Class of Immunogenic Tumor Model with Delay.

This paper aims at studying the model proposed by Kuznetsov and Taylor in 1994. Inspired by Mayer et al., time delay is introduced in the general model. The dynamic behaviors of this model are studied, which include the existence and stability of the equilibria and Hopf bifurcation of the model with discrete delays. The properties of the bifurcated periodic solutions are studied by using the normal form on the center manifold. Numerical examples and simulations are given to illustrate the bifurcation analysis and the obtained results.

1. Introduction

For the longest time, malignant tumors have posed a threat to human lives. Effective strategies based on the immune system have been championed to complement traditional methods of cancer treatment, such as surgery, radiation, and chemotherapy. Cancer immunotherapy is the use of the immune system to treat cancer. Immunotherapy is used to provoke the immune system into attacking the tumor cells by using these cancer antigens as targets. Cell types that can be used in this way are natural killer cells, lymphokine-activated killer cells, cytotoxic T cells, and dendritic cells. There has been much interest in mathematical models describing the interaction between tumor cells and effector cells of the immune system; see [1-6]. An ideal model can provide insight into the dynamics of interactions of the immune response with the tumor and may play an important role in understanding of cancer and developing effective drug therapies. Developing ideal models for such complex processes is not easy. Simple models which display some useful immunological phenomena have been proposed and studied extensively. See Bi and Xiao [7], Galach [4], and Yafia [8-10], and the references cited therein. In 1994, Kuznetsov et al. [11] introduced a model, which describes competition between the tumor and immune cells. It also describes the response of effector cells (ECs) to the growth of tumor cells (TCs). It is assumed that the tumor cells are “immunogenic” and thus subject to immune attack by cytotoxic (CT) or natural killer (NK) cells. This model studies the infiltration of TCs by ECs and also the possibility of ECs inactivation. It is assumed that interactions between ECs and TCs in vitro can be described by the following kinetic scheme describing interactions between ECs and TCs:where E, T, C, T, E are the local concentrations of ECs, TCs, EC-TC complexes, inactivated ECs, and TCs, respectively. The parameters k1, k−1, k2 and k3 are nonnegative constants, which describe the conversion rates of differential cells. Then Kuznetsov and Taylor's model is as follows:where C is the normal rate of the flow of adult ECs into the tumor site, F(C, T) describes the accumulation of effector cells in the tumor cells localization due to the presence of the tumor. F(C, T) = fC/(g + T), (f, g are constants). If Ttot ≈ T  (Ttot = T + C), it is reasonable to assume dC/dt ≈ 0, that is, C ≈ kET  (k = k1/(k−1 + k2 + k3)), then F(C, T) = F(E, T). Then we only need to analyze the dynamical behavior of the first two equations.

In 2003, Galach [4] suggested that the function F is in the Michaelis-Menten form F(E, T) = θET, and thus (2) can be simplified aswhere x is the local concentrations of ECs, and y is the local concentrations of TCs, m1 = kk2, m2 = kk3, and all the parameters are positive. The properties of the model (3) is studied in [7-9]. The dynamical behaviors and the bifurcation of the model with delay are also studied by [7, 10].

In this paper, we analyze the dynamics of an immune response function with Michaelis-Menten form F(E, T) = pET/(g + T), where p and g are positive constants, that is,In order to simplify the original model, we nondimensionalize (4) by choosing scale for E and T cell population, respectively. Let τ = nT0t, and replace τ with t; thus the model can be written as follows:where y = T/g, x = E/g, τ = tng, σ = s/ng2, α = p/ng, μ = m/n, δ = d/ng, f = a/ng, β = bg and all the above parameters are positive.

This paper is organized as follows: In Section 2, the model without delay is considered, and the conditions for existence and stability of equilibria are given. In Section 3, the model with delay is studied. The existence of Hopf bifurcation is obtained; the direction and stability of bifurcated periodic solutions are also given with the help of center manifold and bifurcation theories. At the end of this paper, numerical results are given to illustrate the main results of this paper.

2. Existence and Stability of the Equilibria

It is easy to see that (5) has a tumor-free equilibrium P0 = (σ/δ, 0). In order to find the positive equilibria P1 = (x, y), we need to solve the following equations:Equation (6) is reduced to the following cubic polynomial:where a1 = fβμ, a2 = f(−αβ + βδ − μ + βμ), a3 = fα − fδ + fβδ − fμ + σ, a4 = −fδ + σ. Obviously, system (5) in R2+ is as “well behaved” just as in biology. We have the following lemma by qualitative analysis.

Lemma 1 .

If μ > α, then the solutions of (5) are invariable in D:

Proof

If x = 0 and 0 < y < 2/β, then Similarly, if y = 0 and 0 < x < 2σ/δ, one has Also, if y = 2/β and 0 < x < 2σ/δ, we have and, then, there exists a positive number x, which is given by Such that dx/dt < 0 when x > x. Similarly, dx/dt > 0 when x < x. On the other hand, one has with the help of μ > α. That is 0 < x < 2σ/δ.

When x = 2σ/δ and 0 < y < 2/β, we can prove the result in a similar way. This finishes the proof.

It is easy to see that (5) has one tumor-free equilibrium P0 = (σ/δ, 0). With the distribution of the eigenvalues, we can easily obtain the following results.

Theorem 2 .

If f > σ/δ, then P0 is unstable. If f = σ/δ, P0 is stable. If f < σ/δ, P0 is asymptotically stable.

Theorem 3 .

If >fδ, δβ > μ > α, then P0 is globally asymptotically stable.

When >δf, βδ > μ > α, then it is easy to see that This implies that P0 is the only positive equilibrium.

Noting σ > df, we know that P0 is locally stable with the help of Theorem 2. In the remaining part of the proof, we only need to prove the global stability. By Lemma 1, we know that the following domain D is invariable: Then it is easy to prove that dx/dt < 0 as x > 2σ/δ. In addition, when y0 > 2/β, This shows that the vector fields are moving towards D as t increases.

Let Dulac function B = 1/y. Then Hence there are no closed trajectory surrounds P0 in field D. That is, the result follows.

Using the original parameters, we can give the results as follows.

Theorem 4 .

If s > da/n, db > m > p/g, system (5) has only one critical equilibrium P0 = (σ/δ, 0). Furthermore, P0 is globally asymptotically stable.

Remark 5 .

Theorem 4 is an instructive results to kill the tumor cells. The tumor cells will be killed out by the immune cells sooner or later under the above conditions; then we only need to take the necessary measures to control the parameters to satisfy the inequity in Theorem 4.

In the following, we will study the existence and stability of the tumor-present equilibrium P1(x, y). Similar to the proof of Lemma  2.4 in [12], we can easily obtain the following results.

Theorem 6 .

For the number of positive equilibria, we can get the following results:

If >δ > μ/β > μ, σ < fδ, and Δ < 0, then (5) has three distinct positive roots.

If one of the following conditions is satisfied, then (5) has two distinct positive roots:

α > δ > μ/β > μ, σ > fδ, and Δ < 0.

α > δ > μ/β > μ, σ > fδ and Δ = 0.

Assume one of the following conditions is satisfied, then (5) has one positive root:

α > δ > μ/β > μ, σ < fδ, and Δ < 0.

α > δ > μ/β > μ, σ < fδ, and A = B = 0.

δ > α, β > 1, α + βδ < μ, σ < δf, and Δ > 0.

α > δ > μ/β > μ, σ < fδ, and Δ = 0.

The corresponding Jacobian matrix at P1(x, y) isThus, we can give the following results.

Theorem 7 .

If σβ > μf, then P1 is stable. And no Hopf bifurcation appears around the equilibrium P1.

It is easy to see Tr(J(P1)) = −(σ/x + fβy) < 0. Noting σβ > μf, it is easy to see that Thus the results are proved.

Theorem 8 .

Let (x, y) be the coordinate of the positive equilibrium P1. Then the following results hold:

If fσβ ≥ μx, then P1 is stable.

If fσβ < μx, then P1 is stable for and unstable for .

It is easy to see where a = 1/x(σfβ − μx).

Let Then the sign of the det⁡(J(P1)) is the same as that of h(y). Let Δ = 4α(μ − σfβ/x). If fσβ = μx for any y > 0, then a = Δ = 0. Hence, Then the first result can be obtained easily.

If fσβ > μx, that is, a > 0, Δ < 0, then h(y) > 0 for any y > 0. That is to say, P1 is stable as fσβ > μx. On the other hand, one knows a < 0 as fσβ < μx. Noting α > μ, then Then h(y) = 0 has roots and . From y > 0, we have h(y) < 0 as and h(y) > 0 as . That is to say, the second result holds.

In the following, we give some simulation results of the above results. We consider the system (5) and the parameters suggested by V. A. Kuznetsov et al. Then system (5) becomesObviously, system (24) has three positive equilibria: E1(1.6348,0.366269), E2(0.6538,300.188), and E3(0.1906,441.757). By simple computation, it is easy to know that the eigenvalues of E1 are −0.0367189 ± 0.599721i, and E2 has eigenvalues −1.45811 and 0.295245, and E3 has eigenvalues −1.68973 and −0.375447. These results are represented in Figures 1 and 2.

Figure 1

(a) Stable equilibrium E1. (b) Saddle-node equilibrium E2.

Figure 2

(a) Stable focus equilibrium E3. (b) Saddle node (E2) and stable focus (E3).

3. Dynamical Behaviors of the Model with Delay

In this section, we shall consider model (5) with delayThe existence of the equilibria is the same as those of (5). The dynamical behaviors of the trivial equilibria and semitrivial equilibria are not difficult; we only study the positive equilibrium P1 = (x, y) here. The characteristic equation of linearized system for (25) at P1 takes the following form:where Similar to the proof in [7], then the following results can be obtained easily.

Theorem 9 .

(1) If then all roots of (26) have negative real parts for all τ ≥ 0.

(2) If C2 − E2 < 0 or D2 − B2 + 2C > 0 and (D2 − B2 + 2C)2 = 4(C2 − E2), then (26) has a pair of purely imaginary roots ±iω+ at τ = τ+.

(3) If D2 − B2 + 2C > 0, C2 − E2 > 0 and (D2 − B2 + 2C)2 > 4(C2 − E2), then (26) has a pair of purely imaginary roots ±iω+  (±iω−, resp.) at τ = τ+  (τ = τ−, resp.), where

Theorem 10 .

(1) If the conditions of (1) in Theorem 9 hold, then P1 is asymptotically stable for any τ > 0.

(2) If C2 < E2 or the conditions of (3) in Theorem 9 hold, then E undergoes Hopf bifurcation as τ = τ0.

(3) If D2 − B2 + 2C > 0, (D2 − B2 + 2C)2 = 4(C2 − E2) and B2 + 2DτE + τ2 > D2 + E2, then E undergoes a Hopf bifurcation as τ = τ0.

From the analysis of the above, we only need to compute

In third case, d(Reλ)/dτ| = 0, then we need the sign of the second derivative of Reλ(τ) of (26) at the point where λ is passing through λ0.

We can easily obtainwith

Since sgn⁡Re(dλ/dτ)−1| = 0, sgn⁡Im(dλ/dτ)−1| ≠ 0. Thus it is easy to see that E2 + D2ω2(2ω2 − 2C − 2DτE − τ2) − 2E2A − A2 < 0. This finishes the proof.

In the above, we obtained the conditions under which a family of periodic solutions bifurcated from the positive equilibrium at τ = τ0. In the following, we derive the explicit formulae for determining the properties of the Hopf bifurcated solution by using the normal form and the center manifold theory. Throughout this section, we always assume that the system (25) undergoes Hopf bifurcations at the positive equilibrium P1(x, y) for τ = τ0 and ±iω are the corresponding pure imaginary roots.

Let t = τt′, x(t) = x(τt′). Then (25) becomesSet τ = τ0 + p, p ∈ R. Then p = 0 is a Hopf bifurcation value for (33). Equation (5) can be written aswhere x(t) = (x1(t), x2(t)) ∈ R2, and L : C → R, F : R × C → R are given, respectively, bywith ϕ(0) = (ϕ1(0), ϕ2(0)) ∈ C and

By Riesz representation theorem, there exists a function η(θ, p) of bounded variation for θ ∈ [−1,0] such thatFor ϕ ∈ C1([−1,0], R3), defineThen system (34) is equivalent towhere x(θ) = x(t + θ) for θ ∈ [−1,0]. For Ψ ∈ C1([0,1], (R3)), defineand a bilinear inner productwhere η(θ) = η(θ, 0). Then A(0) and A are joint operators.

Suppose that q(θ) = (1, ρ)e is the eigenvector of A(θ) corresponding to iω0τ0, with ρ = −y/(fβy + iω0τ0). Then q(s) = D(1, ρ)e is the eigenvector of A with ρ = (1/y)[αy/(1 + y) − δ − μye − iω0τ0], such that 〈q(s), q(θ)〉 = 1, .

In the following, we use the ideas in Adam and Bellomo [1] to compute the coordinates describing center manifold C0 at p = 0. DefineOn the center manifold C0, we havewhere z and are local coordinates for the center manifold C0 in the direction of q and Note that W is real if x is real. For x ∈ C0 of (39), we havewhereFrom (43) and (44), we have and q(θ) = (1, ρ)e. Comparing the coefficients of the above equality with (45), we obtain In order to determine g21, we also need to compute W20(θ) and W11(θ) as follows:where

On the center manifold C0 near the origin,Comparing the coefficients of z, we obtainFrom (52) and the definition of A, one hasSolving the above ODE, thenwhere E1 = (E1(1), E1(2)) ∈ R2 is a constant vector. Similarly, we obtainwhere E2 = (E2(1), E2(2)) ∈ R2 is a constant vector.

From the definition of A and the above equations, we can obtain with In the same way, it is easy to get In order to determine E1 and E2, we can calculate W20(θ) and W11(θ) from (54) and (55). Furthermore, we can also determine g21. Thus we can be able to compute the following values:

Theorem 11 .

Equation (25) undergoes Hopf bifurcation as τ = τ0. If μ2 < 0 (μ2 > 0), then the bifurcated periodic solution is supercritical (subcritical); if β < 0 (β < 0), then the bifurcating solution is stable (unstable). The period of the bifurcated periodic solution is T2.

In the following, we present some numerical results of system (25) at different values of τ. We choose the parameters the same as before; the system (25) becomesObviously, system (60) has three positive equilibria: E1(1.6348,0.366269), E2(0.6538,300.188), and E3(0.1906,441.757). Equilibrium E1 is stable for all τ, while equilibrium E2 is unstable saddle node for all τ. When τ = 0, the positive equilibrium E3 is a stable focus, as τ is growing to τ0 = 0.707732, λ′(τ0) = 0.604929 − 0.176726i. When τ passes the critical value τ0, E3 loses its stability and a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcate from E1. This is clearly demonstrated by Figure 3.

Figure 3

(a) Stable equilibrium E1 as τ = 1. (b) Bifurcated solution around E3 as τ = 0.7.

4. Conclusions

In 1994, Kuznetsov et al. [11] introduced a model, which describes competition between the tumor and immune cells. This mode describes the response of effector cells to the growth of tumor cells. Later, the model is simplified to a two-dimensional differential equation with Lotka-Volterra form immune response function F(E, T) = θET; the rich dynamical behaviors are studied such as stability, qualitative behaviors, and bifurcation behaviors.

In this paper, we studied the model with an immune response function with Michaelis-Menten form F(E, T) = pET/(g + T). In this case, the properties of system (4) are more complicated. There are so many tedious computations even for the existence of the equilibria. The dynamical behaviors of the model are more rich.

We have studied the nonlinear dynamics of a two-dimensional general differential system. We first provided linear analysis of the system at the possible equilibria, namely, the semitrivial and positive equilibria, and discussed the existence of Hopf bifurcation at the positive equilibrium. Then we consider the system with delay; we investigated the Hopf bifurcation of the system. The existence and stability of periodic solutions were given. Numerical simulations were presented to illustrate the theoretical analysis and results. Cancer immunosurveillance functions are taken as an important defense to cancer; this is the elimination process. In fact, the existence and stability of the semitrivial equilibrium correspond to the elimination process. Our results on the existence and stability of the Hopf bifurcated periodic solutions describe the equilibrium process. If a stable periodic orbit exists, then the tumor and the immune system can coexist for a long time, although the cancer cannot be eliminated. Furthermore, the parameters are important in controlling the development and progression of the tumor, which is decided by the conditions of the existence of the bifurcations.

The existence of oscillatory solutions in the tummor and immune system interaction models demonstrates that the phenomenon has been observed in some related models [7, 10, 13–15]. The initial values and delay are also important in the oscillatory coexistence of the tumor cells and the effector cells. Numerical simulations indicate this information well.