Hopf Bifurcation of an Epidemic Model with Delay.

A spatiotemporal epidemic model with nonlinear incidence rate and Neumann boundary conditions is investigated. On the basis of the analysis of eigenvalues of the eigenpolynomial, we derive the conditions of the existence of Hopf bifurcation in one dimension space. By utilizing the normal form theory and the center manifold theorem of partial functional differential equations (PFDs), the properties of bifurcating periodic solutions are analyzed. Moreover, according to numerical simulations, it is found that the periodic solutions can emerge in delayed epidemic model with spatial diffusion, which is consistent with our theoretical results. The obtained results may provide a new viewpoint for the recurrent outbreak of disease.

Introduction

Currently, new infectious diseases continuously emerge, and existing diseases recurrently outbreak [1-9]. Ebola virus was firstly discovered in 1976, which began to outbreak in Guinea in February 2014, then spread to West Africa. It caused serious death and social panic. After the outbreak of 2014, Ebola once again emerged in Guinea in March 2016 [10-12]. These diseases have brought a great threat to the public health. In order to provide some suggestions for the prevention and control of the disease, it is necessary to establish rational mathematics model based on infectious mechanism of disease, the route of transmission, and the symptoms of the infected individuals. In particular, the incidence rate describes the number of new infections per unit time, which largely reflects the transmission mechanism of the disease [13-17]. For example, Capasso et al. proposed saturated incidence rate βSI/(1 + kI) to model the cholera epidemics in Bari in 1973, which reflects the psychological effect or the inhibition effect [18]. By taking appropriate preventive measures, May and Anderson gave nonlinear incidence rate β(SI/(1 + αS)) [19]. Therefore, some reasonable suggestions can be provided for the prevention and effective control of infectious diseases.

It takes an individual a period of time to show the corresponding symptoms based on the infectious mechanism of disease, after an individual is infected disease, such as, dengue, rabies, cholera and so on [20-27]. Therefore, time delay describing the incubation period of disease is a significant quantity. In fact, these potentially asymptomatic individuals (incubation individuals) may promote the wide spread of disease [28, 29]. Thus, it is necessary for us to introduce time delay in the epidemic models.

Because of all the species living in the space environment, and they could diffuse the surrounding area. The individual diffusion in space has an effect on the disease contagion. For example, Zhang et al. indicated that dog movement led to the traveling wave of dog and human rabies and had a large influence on the minimal wave speed [30]. However, previous works on epidemic models did not account for the spatial diffusion factors. McCluskey proved that the endemic equilibrium was globally asymptotically stable whenever it existed for an SIR epidemic model with delay and nonlinear incidence rate [31]. A delayed predator-prey system with disease in the prey was investigated by Han et al., they considered the existence of Hopf bifurcation with time delay in terms of degree 2 [32]. Hence, it is more suitable for us to consider time delay and spatial factor in epidemic model.

This paper is organized as below. In Section II, the eigenpolynomials of spatiotemporal epidemic model with nonlinear incidence rate are given, we further analyze the existence of Hopf bifurcation for two cases. In Section III, by using the normal form theory and the center manifold theorem, some properties of Hopf bifurcation are showed. In Section IV, on the basis of numerical simulations, we show that the epidemics will display recurrent behavior if time delay exceeds a critical point. Finally, some conclusions are obtained.

Materials and Methods

Existence of Hopf bifurcation

We consider a SI epidemic model with nonlinear incidence rate βSI with p > 0, q > 0. This form of nonlinear incidence rate was firstly proposed by Liu et al., and exhibited qualitatively different dynamical behaviors [33, 34]. Therefore, it is helpful to interpret some complex epidemic phenomena. In this paper, let p = 1 and q = 2. Since the time delay describing incubation period of transmission process widely exists in most epidemiological models [35-37], thus we need to introduce the time delay into the infected population. Furthermore, we consider Neumann boundary conditions. Consequently, the following system with Neumann boundary conditions is given: where S(x, t) represents the number of the susceptible at location x and time t, I(x, t) the number of the infectious at location x and time t, A represents the recruitment rate of the susceptible, d and μ are natural death rate and the disease-related death rate due to the infected, respectively. d1 and d2 are diffusion coefficients. x represents the one dimensional space, and denotes the usual Laplacian operator.

Assuming ϕ = (ϕ1, ϕ2) ∈ ℘ = C([−τ,0], X), τ > 0 and X is defined as with the inner product 〈⋅, ⋅〉.

System (1) without diffusion and delay corresponds to the following system:

The system (2) has three equilibria, , the saddle and the stable node E*(S*, I*), where Based on the biological meaning, E1 and E* are satisfied the following conditions [38]: Let , , then system (1) can be transformed into: One can define and for i, j, l = 0, 1, 2, …, let

In the phase space ℘ = C([−τ,0];X), the abstract differential equation of the system (3) is where , We set U(θ) = U(t+θ),ϕ = (ϕ1, ϕ2) ∈ ℘, ϕ(θ) = U(θ), and θ ∈ [−τ, 0].

Let L: ℘ → X and F: ℘ → X are given by where and where , , , , , .

The linearized part of system (4) is given by then we set U(t) = ye, and y = (y1, y2), hence the characteristic equation is where and .

On the basis of the Laplacian operator in the bound domain, on X have eigenvalues −k2 with the corresponding eigenfunctions , , k ∈ N0 = {0, 1, 2…}, namely, a basis of the phase space X is . Thus, for ∀ y ∈ X, y can be expanded as Fourier series in the following form:

Furthermore, through simple computations, we have

From the above Eqs (7) and (8), Eq (6) can be written as then the eigenpolynomial associated with λ of system (1) is given by: where a11 = −β(I*)2−d, a13 = −2βS* I*, a21 = β(I*)2, a22 = −(d+μ), a23 = 2βS* I*, then a11 a23−a13 a21 < 0 can be derived.

Considering k = 0, the eigenpolynomial Eq (10) becomes

By replacing λ with iw(w > 0) in Eq (11), then

Through separating the real and imaginary parts of above equations, the following equations are obtained:

Further, by squaring and adding the two parts of Eq (12), we get where . Thus we give if the formula B < 0, and B2 − 4C > 0 then w2 > 0 and can hold simultaneously. Moreover,the corresponding condition is therefore, Eq (11) has two groups of simple imaginary roots ±iw0, .

In the following part, we take into account the imaginary roots ±iw0, the other one is similar.

From Eq (12), we can obtain Moreover, some simple derivations show that

λ(τ) = α(τ) + iw(τ) is the root of Eq (11) near , which satisfies and , where j = 0, 1, 2, ….

Next, taking the derivative with respect to τ on two sides of Eq (11), then we derive By the above expression, one can derive

So the transversality condition is deduced.

Theorem 2.1 If (A1), (A2) and (A3) are all satisfied, system (1) without diffusion experiences a spatially homogeneous Hopf bifurcation at equilibrium E* = (S*, I*) when , and period solution will appear.

Lemma 2.1 (S1) If there is a certain k0 ∈ N = {1, 2, …} such that then Eq (10) has a pair purely imaginary roots ±iw, and

Proof: If we assume k = k0 ∈ N, and λ = iw(w > 0) be a root of Eq (10). By inserting iw(w > 0) into Eq (10) and using the same method as before, then Eq (10) can be translated into: where , C = [d1 d2 k4−d2 a11 k2−d1 a23 k2+a11 a23]2−[d1 a22 k2−(a11 a22−a12 a21)]2, for ∀k ∈ N. Besides, we set C = C1 × C2, where

Further, we can deduce

It is clear that C1(k) > 0 for ∀k ∈ N, according to assumption (S1), C(k0) < 0 is given, then we get ,

Lemma 2.1 shows that the critical value of bifurcation parameter τ can be found. Similar to the method for the case of k = 0, then we get

Let λ(τ) = α(τ) + iw(τ) be the root of (10) near which satisfies and , where j = 0, 1, 2, ….

Lemma 2.2 If condition (S1) is established, then the transversality condition is derived.

The proof can be found in S1 File.

Theorem 2.2 In the presence of space, if the conditions (A1) and (S1) are satisfied, then system (1) undergoes a Hopf bifurcation at E* = (S*, I*) when , and period solution will emerge.

Results

The properties of Hopf bifurcating period solutions

The above section gives the conditions of the existence of Hopf bifurcation for two cases. In this section, we investigate properties of these bifurcating periodic solutions from the positive constant steady state E*(S*, I*) of system (1) by employing the normal form theory and the center manifold theorem of partial functional differential equations (PFDEs) [39-42], these properties include the direction, stability and period. It’s simple for mathematical calculation to mark .

Let , , then system (3) can be expressed as

In the space ℘ = C([−1, 0], X), let τ = τ + α (α ∈ R), , then system (18) can be rewritten as:

Let ϕ = (ϕ1, ϕ2) ∈ ℘, U(θ) = U(t+θ), and ϕ(θ) = U(θ) for θ ∈ [−1, 0]. Defining L(b)(⋅): R × ℘ → X (b is τ or α) and F: ℘ × R → X as and where

Next, the linear part of the system (19) is given by

From the conclusions of section II, an equilibrium of the system (20) is the origin, the corresponding characteristic equation of the system (20) at origin has two pairs of purely imaginary eigenvalues ±iw0 τ, for k = 0, and only a pair of purely imaginary eigenvalues ±iw τ for k ∈ N. We account for purely imaginary eigenvalues ±iw0 τ for the case k = 0, and set Λ0 = {iw τ, −iw τ}, (k = 0, k0).

Considering the ordinary functional differential equation:

For ϕ ∈ C([−1, 0], X), according to the Riesz representation theorem, there is a 2 × 2 matrix function η(θ, τ)(−1 ≤ θ ≤ 0), then we have [40] where

For ϕ ∈ C([−1, 0], X), defining semigroup induced by the solution of the linear eq (20), and the infinitesimal generator A(τ) of the semigroup is

For ψ ∈ C([0, 1], X), the formal adjoint operators of A(τ) is A*(τ) which denotes [43]

Here, the bilinear pairing form associated A(τ) with A*(τ) is

On the basis of the discussion of section II, A(τ) has a pair purely imaginary eigenvalues ±iw τ, which are also eigenvalues of A*(τ). Furthermore, the generalized eigenspaces of A(τ) and A*(τ) associated with Λ0 are the center subspaces P and P*, respectively. P* is the adjoint space of P and dimP = dimP* = 2 [42].

By some computations, the following Lemma is directly given:

Lemma 2.3 A basis of P with Λ0 is and a basis of P* with Λ0 is where

Φ = (Φ1, Φ2) and are obtained by separating the real and imaginary parts of q1(θ) and , respectively. Obviously, Φ is the basis of P, Φ* is the basis of P*, and

According to the bilinear pairing form Eq (26), we can compute:

Next, we construct a new basis Ψ for P*, where Ψ = (Ψ1, Ψ2) = (Φ*, Φ)−1 Φ* and . (Ψ, Φ) = I2 needs to be satisfied. In addition, , where , . For c = (c1, c2) ∈ C([−1, 0], X), we define .

On the basis of the theory of decomposition of the phase space, we have ℘ = P℘ + P℘, where P℘ is the center subspace of linear Eq (20), and P℘ is the complement subspace of P℘.

Since the infinitesimal generator A(τ) is induced by the solution of Eq (20), then Eq (18) can be translated into: where

According to the phase space decomposition ℘ = P℘ + P℘ and Eq (28), the solution of Eq (19) is written as where , and h(x1, x2, α) ∈ P℘, h(0,0,0) = 0, Dh(0,0,0) = 0. Moreover, the solution of Eq (19) on center manifold is

Let z = x1 − ix2, Ψ(0) = (Ψ1(0), Ψ2(0)), and q1 = Φ1 + iΦ2, thus

By using the previous variable substitution, Eq (31) can be transformed into: where , and setting

According to the conclusions of Ref. [42], z satisfies where and setting

From f(ϕ, a) and Eq (31), it is easy to compute where and i, j = 0, 1, 2, …, m = 1, 2.

Let (ψ1, ψ2) = Ψ1−iΨ2, then we can obtain

Since the expression of g21 containing W20(θ) and W11(θ) for θ ∈ [−1, 0], it is necessary to compute them. From Eq (34), we can derive

Meanwhile, from the conclusion of literature [42], where with H ∈ P*, i, j = 0, 1, 2….

Therefore, from Eqs (35) and (37)–(41), the following form can be given by:

Because A(τ) has only two eigenvalues ±iw τ, Eq (42) has unique solution W in the following form:

From Eq (41), for −1 ≤ θ < 0,

Thus, for −1 ≤ θ < 0,

For θ = 0, , we have

Based on the definition of infinitesimal generator A(τ), then Eq (42) is transformed into and −1 ≤ θ < 0.

From q1(θ) = q1(0)ei, − 1 ≤ θ ≤ 0, we have further we obtain

For k = 0, θ = 0, in the light of the definition of A(τ) and Eq (49), the first Eq of Eq (42) becomes

So we can derive

From Eq (51), the formula of C1 can be derived where

For −1 ≤ θ < 0, similar to the above case, W11(θ) can be obtained further we derive where

Through the above calculations of W20(θ) and W11(θ), we obtain the expression of g21. Consequently, in order to determine the properties of Hopf bifurcating period solutions at the critical value τ, we can compute the following values:

μ2 > 0 (μ2 < 0) determines the direction of the Hopf bifurcation is supercritical (τ > τ) (subcritical (τ < τ)); if β2 < 0 (β2 > 0) indicates that the bifurcating period solutions on center manifold are asymptotically stable (unstable); furthermore, T2 can determine the period of the bifurcating period solutions, namely, T2 < 0 (T2 > 0) represents the decrease (increase) of the period.

Numerical results

Compared with the theoretical analyses, we perform a series of extensive numerical simulations of the spatiotemporal epidemic model with nonlinear incidence rate in one-dimensional space, and investigate the incubation period how to affect the spread of epidemics. We solve the numerical solutions of system (1) by using Matlab. The reaction-diffusion system is solved in a discrete domain with N × N lattice sites. The Laplacian describing diffusion is approximated by using finite differences, and we also discretize the time evolution.

In case k = 0, we set d1 = 6, d2 = 1, A = 1, β = 32, μ = 1.8, d = 1, the equilibrium is E* = (S*, I*) = (0.43, 0.20). By some calculations, , c1(0) = −9.81 + 22.15i are obtained. Through the formulae of properties of Hopf bifurcating period solutions in section III, we get μ2 > 0, β2 < 0 and T2 > 0. These parameter values shows E* is asymptotically stable for 0 ≤ τ < τ. With the increase of τ, E* loses its stability and Hopf bifurcation occurs at critical point τ, these bifurcating period solutions are stable, the direction of bifurcation is forward and the period increases, which are presented in Fig 1.

Fig 1

Hopf bifurcation with k = 0.

(a) The constant steady state E* is asymptotically stable for τ = 1.2 with initial conditions S(x, t) = 0.42, I(x, t) = 0.20, (x, t) ∈ [0, π] × [−1.2, 0]; (b) The bifurcating periodic solutions are asymptotically stable for τ = 1.6 with initial conditions S(x, t) = 0.42, I(x, t) = 0.20, (x, t) ∈ [0, π] × [−1.6, 0].

In case k = 1, setting d1 = 6, d2 = 1, A = 1, β = 32, μ = 1.8, d = 1, then the equilibrium is E* = (S*, I*) = (0.43, 0.20). Furthermore, by using the formulae derived in section III, we compute , c1(0) = 2.30 × 102 − 7.55 × 101i. By computing the formulae (58), μ2 > 0, β2 < 0 and T2 < 0 are obtained, which indicates that these bifurcating period solutions are stable, the direction of Hopf bifurcation is forward, and the period decreases. These phenomena are showed in Fig 2.

Fig 2

Hopf bifurcation with k = 1.

(a) When τ = 0.3, the constant steady state E* is asymptotically stable with initial conditions S(x, t) = 0.42, I(x, t) = 0.20, (x, t) ∈ [0, π] × [−0.3, 0]; (b) When τ = 1.5, the bifurcating periodic solutions are asymptotically stable with initial conditions S(x, t) = 0.42, I(x, t) = 0.20, (x, t) ∈ [0, π] × [−1.5, 0].

Discussion

In this study, the characteristic equation at the positive constant steady state E*(S*, I*) is derived. In order to study the influence of incubation period on epidemic transmission, we choose time delay τ as a bifurcation parameter. Moreover, we get the two classes conditions of the existence of Hopf bifurcation: one is the absence of diffusion k = 0, the other is the presence of diffusion k = k0 ∈ N. With increasing of parameter τ, the stability of positive constant steady state E*(S*, I*) will change, and Hopf bifurcation will concurrently occur in system (1) at the critical point τ( or ). In the following, we obtain the properties of bifurcating period solutions including direction, stability and period by utilizing the normal formal theory and the center manifold theorem of partial functional differential equations (PFDs).

It should be noted that spatial pattern may be found in epidemic model (1). Based on pattern dynamics of model (1), one can obtain the pattern structures in different parameters space [44, 45]. In this case, we can reveal the distributions of disease with high density or low density and thus provide useful control measures to eliminate the disease.

Conclusion

The numerical results validate our theoretical findings, which show that the length of the incubation period have significant impacts on epidemic transmission. The biennial outbreaks of measles is the signature of an endemic infectious disease, which becomes non-endemic if there were a minor increase in infectivity or a decrease in the length of the incubation period [15]. Based on this paper, we provide a possible mechanism to explain the recurrent outbreak of disease.

(EPS)

Click here for additional data file.

(EPS)

Click here for additional data file.

Transversality condition.

The relationship between real part of eigenvalues and time delay.

(PDF)

Click here for additional data file.