Stability and Hopf Bifurcation Analysis of an Epidemic Model with Time Delay.

Epidemic models are normally used to describe the spread of infectious diseases. In this paper, we will discuss an epidemic model with time delay. Firstly, the existence of the positive fixed point is proven; and then, the stability and Hopf bifurcation are investigated by analyzing the distribution of the roots of the associated characteristic equations. Thirdly, the theory of normal form and manifold is used to drive an explicit algorithm for determining the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions. Finally, some simulation results are carried out to validate our theoretic analysis.

1. Introduction

Today, the serious epidemics, such as SARS and H1N1, are still threatening the life of people continually. Plenty of mathematical models have been proposed to analyze the spread and the control of these diseases [1-7].

However, many infectious diseases, for instance, gonorrhea and syphilis, occur and spread amongst the mature, while some epidemics, for example, chickenpox and FMD, only result in infection and death in immature. For this reason, stage structure should be taken into consideration in models. Aliello and Freedman [8] proposed a stage-structured model described bywhere x(t) is the immature population density and y(t) represents the density of the mature population. α, γ, τ, and β are all positive constants. α is the birth rate, and γ is the natural death rate; τ is the time from birth to maturity; β is the death rate of the mature because of the competition with each other.

And then, many infectious diseases with sage structure have been built and investigated [9-14]. Xiao and Chen [15] improved (1) by separating the population into mature and immature and supposing that only the immature were susceptible to the infection.

Based on the model in [15], supposing that only the mature were susceptible, Jia and Li [16] built a new one as follows:

where x(t), y(t), and R(t) are the susceptible, infectious, and recovered mature population densities, respectively; z(t) denotes the immature population density. All the parameters are positive constants. α, β, γ, and τ are the same as those in (1); m is the transmission coefficient describing the infection between the susceptible and the infectious; c is the death rate because of the epidemic; g is the recovery rate; αe−x(t − τ) denotes the population who were born at t − τ and survive at t.

In systems (1) and (2), the time delay was also taken into consideration. Indeed, time delay plays an important role in the epidemic system, making the models more accurate. In recent years, delays have been introduced in more and more epidemic and predator-prey systems [17-19].

In this paper, on the basis of (2), we further assume that

Both the susceptible and the infectious have fertility, while in (2), only the susceptible is fertile

For the infectious, there is competition with all the susceptible and the infectious, while for the susceptible, there is only competition between generations

Meanwhile, all the death of the susceptible, the same as that in (1), is only due to the competition. To simplify model (2), we denote γ + c + g = d, and let b(x(t) + y(t)) present the transmission from immature to mature.

As a consequence, the new epidemic model could be described as follows:where w is the death rate of the mature because of the competition.

We can notice that z(t) depends on x(t) and y(t) and R(t) depends on y(t); however, x(t) and y(t) have nothing to do with z(t) and R(t). According to Qu and Wei [20], we will mainly focus on x(t) and y(t), that is,

The rest of the paper is organized as follows. In Section 2, we calculate the steady states of system (4) and prove the existence and uniqueness of the positive equilibrium in particular. And then, the stability of the two nonzero equilibria and the existence of the Hopf bifurcation are investigated in Sections 3 and 4, respectively. In Section 5, the direction and stability of the Hopf bifurcation at the positive equilibrium are studied by using the center manifold theorem and the normal form theory [21]. And in the last section, some numerical simulations are carried out to validate the theoretical analysis.

2. The Existence and Uniqueness of the Positive Equilibrium of the Model

In this section, we discuss the existence of the equilibria of (4) and the positive one in particular.

The equilibria are the solutions of the equations (5),

Clearly, E1(0, 0) and E2(b/w, 0) are two equilibria of (4).

In the following, we will focus on the existence of the positive equilibrium.

Theorem 1 .

If b(m − w) > dw, (4) has one positive equilibrium E3(x∗, y∗), where

Proof

Positive equilibrium is the positive solution of the equations (7),

From the second equation of (7), we have

Taking (8) into the first equation of (7), we can obtain

which leads to

where

Together with

we can know that both of the two solutions of (9) are positive,

where

If ,

then,

So, is dropped.

If ,

then,

w > 0, d > 0, and , so ,

Then, taking into (8), we have

Therefore, if b(m − w) > dw, (4) has the unique positive equilibrium E3(x∗, y∗).☐

3. Stability Analysis of the Equilibrium E2(b/w, 0)

In this section, we analyze the stability of the equilibrium E2(b/w,  0).

For convenience, the new variables u(t) = x(t) − b/w and v(t) = y(t) are introduced, and then, around E2(b/w,  0), the system (4) could be linearized as (18):

whose characteristic equation is given by

from which, we can get that

or

Obviously, if b(m − w) > dw,

then

which implies that the equilibrium E2 is unstable.

If b(m − w) < dw,

then

As a consequence, we will discuss the other roots of (19), that is, the roots of (21), under the condition b(m − w) < dw.

For τ = 0, equation (21) becomes

whose root is

For τ > 0, if iω(ω > 0) is a root of (21), then

(27) can be obtained by separating the real and the imaginary parts,

which leads to

from which, we can get the unique positive root

Let

Then, when τ = τ, (21) has a pair of purely imaginary roots ±iω0.

Suppose

which is the root of (21) such that

To investigate the distribution of λ(τ), we will discuss the trend of λ(τ) at τ = τ.

Substituting λ(τ) into (21) and taking the derivative with respect to τ, we can get

which yields,

Together with (27), we have

which means that when undergoes τ = τ, λ(τ) will add a pair of roots with positive real parts. That is, with the increase of τ, the number of roots with positive real part is increasing, leading to the change of the stability of the system (4).

Therefore, the distribution of the roots of (21) could be obtained.

Lemma 1 .

Let ω0 and  τ ( j = 0, 1, 2,  ⋯) be defined by (29) and (30), respectively.

If b(m − w) > dw, then (19) has at least one positive root

If b(m − w) < dw, and τ = 0, then both roots of (19) are negative

If b(m − w) < dw, and τ > 0, then (19) has a pair of simple imaginary roots ±iω0 at τ = τ; furthermore, if τ < τ0, then all the roots of (19) have negative real parts; if τ ∈ (τ, τ), (19) has 2(j + 1) roots with positive real parts

Together with condition (35), the Hopf bifurcation theorem [21], and Lemma 2, the following theorem could be stated.

Theorem 2 .

Let  τ (j = 0, 1, 2,  ⋯) be defined by (30), then we have

If b(m − w) > dw, then the equilibrium E2(b/w, 0) of (4) is unstable

If b(m − w) < dw, then the equilibrium E2(b/w, 0) is asymptotically stable when τ ∈ [0, τ0), and it is unstable when τ > τ0

If b(m − w) < dw, then system (4) undergoes a Hopf bifurcation at the equilibrium E2(b/w, 0) for τ = τ

4. Stability Analysis of Positive Equilibrium E3(x∗, y∗)

In this section, we analyze the stability of the positive equilibrium E3(x∗, y∗).

For convenience, the new variables u(t) = x(t) − x∗ and v(t) = y(t) − y∗ are introduced, and then, around E3(x∗, y∗), the system (4) could be linearized as (36):

whose characteristic equation is given by

where

For τ = 0, equation (37) becomes

Firstly, computing λ1λ2, we have

where Δ is the same as that in (10).

So,

which implies that the real parts of  λ1 and λ2 have the same signs.

Then, (λ1 + λ2) is calculated:

Together with (41), we can get that both the real parts of the two roots of (39) are negative.

For τ > 0, equation (37) can be rewritten as

where

If iω(ω > 0) is a root of (37), then

Separating the real and the imaginary parts, we havewhich leads to

Let z = ω2 > 0, and then, (47) can be rewritten as

Firstly, computing z1z2, we getwhere

has been proved in (41).

Then, we will calculate

If H(4-1): ,

then

and then,

which implies that (48) has one unique positive solution z0 = ω02,

where

If H(4-2): ,

then

and then,

Let

If ∆1 < 0, then (48) has no real roots.

If ∆1 > 0 and z1 + z2 = 2p − r2 + s2 < 0, then (48) has two negative roots, and there no positive ω for (47);

If ∆1 > 0 and z1 + z2 = 2p − r2 + s2 > 0,

then (48) has two positive roots, and there are two positive ω for (47), which are

Lemma 2 .

If , then (47) has one positive root ω0

When , if ∆1 < 0, or z1 + z2 = 2p − r2 + s2 < 0, then (47) has no positive roots

When , if ∆1 > 0 and z1 + z2 = 2p − r2 + s2 > 0, then (47) has two positive roots

In the following, we will discuss the expression of τ.

Lemma 3 .

If

where a > 0, t > 0, then

If f(a) > 0, g(a) > 0, then

(2) If f(a) < 0, g(a) > 0, then

(3) If f(a) > 0, g(a) < 0, then

(4) If f(a) < 0, g(a) < 0, then

In conclusion,where j = 0, 1, 2, ⋯

If g(a) > 0, then at = arccos(f(a)) + 2jπ

If g(a) < 0, then at = 2π − arccos(f(a)) + 2jπ,

According to (46), we have

If , substituting ω0 defined in (54) into (64), we can get f(ω0) and g(ω0). Together with Lemma 5, the expression of τ could be obtained.

If

then

If

then

That is, when τ = τ, the characteristic equation (37) has a pair of purely imaginary roots ±iω0.

Suppose λ(τ) = α(τ) + iω(τ) is the root of (37), and then, we have

To investigate the distribution of the λ(τ), we will discuss the trend of λ(τ) at τ = τ.

Substituting λ(τ) into (37) and taking the derivative with respect to τ, we can get

which leads to

Together with (46) and (54), we have

which means that when undergoes τ = τ, λ(τ) will add a pair of roots with positive real parts. That is, with the increase of τ, the number of roots with positive real part is increasing, leading to the change of the stability of the system (4).

If , ∆1 > 0 and z1 + z2 = 2p − r2 + s2 > 0, then f(ω0) and g(ω0) could be calculated by substituting ω± defined in (58) into (64). According to Lemma 5, we can get the expression of τ±:

If

then

If

then

That is, when τ = τ±, the characteristic equation (37) has a pair of purely imaginary roots ±iω±.

Let λ(τ) = α(τ) + iω(τ) be the root of (37), satisfying

To investigate the distribution of the λ(τ), we will discuss the trend of λ(τ) at τ = τ±.

Using the same method, we have

This implies that

and

which means that when undergoes τ = τ+, λ(τ) will add a pair of roots with positive real parts, while undergoes τ = τ−, λ(τ) will lose a pair of roots with positive real parts; if τ− > τ+, then the characteristic equation (37) must have roots with positive real parts for τ > τ+.

In conclusion, the distribution of the roots of (37) could be obtained.

Lemma 4 .

Let ω0 be defined by (54), and  τ ( j = 0, 1, 2,  ⋯) be defined by (66) or (68), and ω± be defined by (58), and τ± ( j = 0, 1, 2,  ⋯) be defined by (74) or (76), respectively.

When , if ∆1 < 0, or z1 + z2 = 2p − r2 + s2 < 0, then all the roots of (37) are with negative real parts

When , (37) has a pair of simple imaginary roots ±iω0 at τ = τ; furthermore, if τ ∈ [0, τ0), then all the roots of (37) are with negative real parts; if τ ∈ (τ, τ), then (37) has 2(j + 1) roots with positive real parts

When , if ∆1 > 0, and z1 + z2 = 2p − r2 + s2 > 0, then (37) has a pair of simple imaginary roots ±iω± at τ = τ±; if τ ∈ [0, τ0+) or τ ∈ (τ−, τ+), then all the roots of (37) are with negative real parts; if τ ∈ (τ+, τ−), then (37) has a pair of roots with positive real parts; if τ− > τ+, for τ > τ+, (37) must have roots with positive real parts

Together with conditions (72) and (78), the Hopf bifurcation theorem [21], and Lemma 6, the following theorem could be stated.

Theorem 3 .

Let  τ ( j = 0, 1, 2,  ⋯) be defined by (66) or (68), and τ± ( j = 0, 1, 2,  ⋯) be defined by (74) or (76), respectively, and then we have

When , if ∆1 < 0, or z1 + z2 = 2p − r2 + s2 < 0, the equilibrium E3(x∗, y∗) of (4) is asymptotically stable

When , then the equilibrium E3(x∗, y∗) is asymptotically stable when τ ∈ [0, τ0), and it is unstable when τ > τ0. System (4) undergoes a Hopf bifurcation at the equilibrium E3(x∗, y∗) for τ = τ

When , if ∆1 > 0, and z1 + z2 = 2p − r2 + s2 > 0, then there is a positive h, such that when τ ∈ [0, τ0+) ∪ (τ0−, τ1+) ∪ (τ1−, τ2+) ∪ ⋯∪(τ−, τ+), the equilibrium E3(x∗, y∗) is asymptotically stable, and it is unstable when τ ∈ (τ0+, τ0−) ∪ (τ1+, τ1−) ∪ (τ2+, τ2−) ∪ ⋯∪(τ+, τ−) ∪ (τ+, +∞). We call that system (4) undergoes (h + 1) switches

System (4) undergoes a Hopf bifurcation at the equilibrium E3(x∗, y∗) for τ = τ±.

5. The Direction and Stability of Hopf Bifurcation at E3(x∗, y∗)

In the previous section, we have already gotten some conditions making that the system (4) undergoes a Hopf bifurcation at the positive equilibrium E3(x∗, y∗) when τ = τ±, j = 0, 1, 2,  ⋯. In this section, under the conditions in Theorem 7, the direction of Hopf bifurcation and stability of the periodic solutions from E3 will be investigated by using the center manifold and normal form theories [21].

Without loss of the generality, let be the critical value of τ = τ+(τ−), j = 0, 1, 2,  ⋯.

For convenience, let , ρ ∈ R, u(t) = x(τt) − x∗, v(t) = y(τt) − y∗, and then system (4) undergoes Hopf bifurcation at ρ = 0; with τ normalized by the time scaling t⟶t/τ, (4) could be rewritten as

Choose the space as C = C([−1, 0], R2); for any Φ = (Φ1, Φ2) ∈ C, let

and

where s and n are the same as those in (37).

According to Riesz's representation theorem, there is a 2 × 2 matrix η(θ, ρ) (θ ∈ [−1, 0]), whose components are bounded variation functions such that

In fact, η(θ, ρ) could be chosen as

where

For any Φ ∈ C1([−1, 0], R2), we define

and

And then, system (4) could be translated into

where

For Ψ ∈ C1([−1, 0], (R2)∗), define

and a bilinear form

where

According to (87) and (91), we can get that A(0) and A∗ are adjoint operators.

From the analysis in Section 4, we know that are a pair of eigenvalues of A(0) and also eigenvalues of A∗, where ω is ω+ or ω− defined in (58).

It is easy to verify that vectors q(θ) and q∗(s) are the eigenvalues of A(0) and A∗ corresponding to the eigenvalues and , where

For convenience, let

then

We choose

and then,

Using the same method as that in [21], the center manifold C0 at ρ = 0 is first computed. Suppose that x is the solution of (81) when ρ = 0, and define

Then, on the center manifold C0, we can get that

where

In the direction q∗ and , z and are local coordinates for center manifold C0. It is not difficult to note that when x is real, W is real. And the real solutions are considered only.

For any solution x ∈ C0, since ρ = 0, we can get that

where

We rewrite

with

Together with (83), we have

Comparing with (105), we get

All the numbers in the expressions of g20,  g11, and  g02 are known; however, W20 and W11 in g21 are unknown, which will be computed in the following.

From (89) and (99),

where

When θ ∈ [−1, 0],

Comparing the coefficients in the expansions of (109) and (110), we have

Taking (5-18) into (5-17), we obtain

By (111) and (112), the following can be gotten

Because

by integration, we have

and

where E1 and E2 are unknown.

From (112) and the definition of A, we obtain

and

Then, together with (108), we get

Substituting (115) and (119) into (117), we have

That is

from which, E1 can be determined.

By (115), (117), and (119), we have

That is

from which, E2 can be determined.

Substituting E1 and E2 into (115) and (116), W20 and W11 could be obtained; furthermore, g21 can be calculated.

Then, the important parameters can be obtained:

which determine the quantities of the bifurcation of system (4) at , where μ2 determines the direction of the bifurcation: Hopf bifurcation is subcritical (supercritical) if μ2 < 0 (μ2 > 0) and the bifurcating periodic solutions exist for (); β2 determines the stability of the bifurcating periodic solutions, which are stable (unstable) if β2 < 0 (β2 > 0); T2 determines the period of the periodic solutions, which decreases (increases) if T2 < 0 (T2 > 0).

From (78) in Section 4, we know that when τ = τ+, and when τ = τ−.

Theorem 4 .

If , then the Hopf bifurcation of system (4) at positive equilibrium E3(x∗, y∗) when τ = τ is supercritical (subcritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if Re[C1(0)] < 0 (Re[C1(0)] > 0).

Theorem 5 .

If , ∆1 > 0, and z1 + z2 = 2p − r2 + s2 > 0, then when τ = τ+, the Hopf bifurcation of system (4) at the positive equilibrium E3(x∗, y∗) is supercritical (subcritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if Re[C1(0)] < 0 (Re[C1(0)] > 0);

when τ = τ−, the Hopf bifurcation of system (4) at the positive equilibrium E3(x∗, y∗) is subcritical (supercritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if Re[C1(0)] < 0 (Re[C1(0)] > 0).

6. Numerical Simulations

In this section, some numerical simulations are carried out to support our theoretical analysis.

There are so many different cases that only the most particular one in (3) of Theorem 7 is considered in this section. The coefficients are chosen as follows: b = 0.4, d = 0.8, w = 1, and m = 6. Then, the conditions of (3) in Theorem 7 are satisfied, where

.

By direct calculation, we have

and

where τ2− > τ3+.

From Theorem 7, we know that the positive equilibrium E3(0.169906, 0.049528) should be asymptotically stable when τ ∈ [0, τ0+) ∪ (τ0−, τ1+) ∪ (τ1−, τ2+), and it is unstable when τ ∈ (τ0+, τ0−) ∪ (τ1+, τ1−) ∪ (τ2+, +∞). The system (4) undergoes 3 switches.

All the simulation results for τ in the six different intervals are in consistent with the theoretical analysis, which are shown in Figures 1–6.

Figure 1

The dynamic behaviors of system (4) with τ = 3.8 ∈ [0, τ0+). (a) and (b) are the waveform plot and phase for system (4), respectively. E3(x∗, y∗) is asymptotically stable and the initial condition is (0.16, 0.1).

Figure 2

The dynamic behaviors of system (4) with τ = 4.4 ∈ (τ0+, τ0−). (a) and (b) are the waveform plot and phase for system (4), respectively. E3(x∗, y∗) is unstable with the emergence of oscillation and the initial condition is (0.18, 0.12).

Figure 3

The dynamic behaviors of system (4) with τ = 11 ∈ (τ0−, τ1+). (a) and (b) are the waveform plot and phase for system (4), respectively. E3(x∗, y∗) is asymptotically stable and the initial condition is (0.1,0.06).

Figure 4

The dynamic behaviors of system (4) with τ = 23 ∈ (τ1+, τ1−). (a) and (b) are the waveform plot and phase for system (4), respectively. E3(x∗, y∗) is unstable with the emergence of oscillation and the initial condition is (0.2, 0.2).

Figure 5

The dynamic behaviors of system (4) with τ = 30 ∈ (τ1−, τ2+). (a) and (b) are the waveform plot and phase for system (4), respectively. E3(x∗, y∗) is asymptotically stable and the initial condition is (0.15, 0.02).

Figure 6

The dynamic behaviors of system (4) with τ = 36 ∈ (τ2+, +∞). (a) and (b) are the waveform plot and phase for system (4), respectively. E3(x∗, y∗) is unstable with the emergence of oscillation and the initial condition is (0.15, 0.14).