Travelling waves of a delayed SIRS epidemic model with spatial diffusion

Nonlinear Analysis: Real World Applications 12 (2011) 52–68

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa

Travelling waves of a delayed SIRS epidemic model with spatial diffusion✩ Qintao Gan ∗ , Rui Xu, Pinghua Yang Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China

article

info

Article history: Received 23 December 2009 Accepted 12 May 2010 Keywords: SIRS Reaction diffusion Travelling waves Upper-lower solutions Partial monotonicity

abstract This paper is concerned with the existence of travelling waves to an SIRS epidemic model with bilinear incidence rate, spatial diffusion and time delay. By analysing the corresponding characteristic equations, the local stability of a disease-free steady state and an endemic steady state to this system under homogeneous Neumann boundary conditions is discussed. By using the cross iteration method and the Schauder’s fixed point theorem, we reduce the existence of travelling waves to the existence of a pair of upper-lower solutions. By constructing a pair of upper-lower solutions, we derive the existence of a travelling wave solution connecting the disease-free steady state and the endemic steady state. Numerical simulations are carried out to illustrate the main results. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Epidemic models study the transmission dynamics of infectious diseases in a host population. Since Kermack and Mckendrick [1] constructed a mathematical model of an ODE to study epidemiology in 1927, the concept of ‘‘compartment modelling’’ has been used until now. Letting S (t ) represent the number of individuals who are susceptible to the disease, that is, who are not yet infected at time t ; I (t ) represents the number of infected individuals who are infectious and are able to spread the disease by contact with susceptible individuals, and R(t ) represents the number of individuals who have been infected and then removed from the possibility of being infected again or of spreading at time t, [2] considered the following SIRS epidemic model

 S˙ (t ) = A − dS (t ) − β S (t )I (t ) + δ R(t ), ˙I (t ) = β S (t )I (t ) − (γ + α + d)I (t ), ˙ R(t ) = γ I (t ) − (δ + d)R(t ),

(1.1)

where parameters A, d, β, δ, γ , α are positive constants in which A is the recruitment rate of the population, d is the natural death rate of the population, γ is the recovery rate of the infective individuals, δ is the rate at which recovered individuals lose immunity and return to the susceptible class, β is the transmission rate, α is the death rate due to disease. The threshold was found in [2] to determine whether the disease dies out or approaches an endemic equilibrium. The SIRS model assumes that the recovered individuals have only temporary immunity, namely, susceptible individuals become infectious, then removed with immunity after recovery from infection and then susceptible again when the temporary immunity fades

✩ This work was supported by the National Natural Science Foundation of China (No. 10671209) and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry. ∗ Corresponding address: Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, 97 Heping West Road, Shijiazhuang 050003, Hebei Province, PR China. Tel.: +86 311 87994035. E-mail addresses: [email protected] (Q. Gan), [email protected] (R. Xu), [email protected] (P. Yang).

1468-1218/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2010.05.035

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

53

away. The assumption has been proved to be reasonable in the study of the dynamics of communicable diseases such as influenza [3] and the sexually transmitted diseases such as gonorrhea [4]. Many epidemic models in the literature represent the dynamics of diseases by systems of ordinary differential equations without delays. However, the inclusion of temporal delays in such models makes them more realistic by allowing to describe the effects of the disease’s latency or immunity (see, for example, [5–10]). Cooke [11] introduced a time delay effect into an epidemic model by assuming that the force of infection at time t is given by β S (t )I (t − τ ), where τ > 0 is a fixed time during which the infectious agents develop in the vector and it is only after that time that the infected vector can infect a susceptible human. Following the work of [11,12] formulated a delayed SIRS epidemic model spread by vectors which have an incubation time to become infectious. By constructing a Lyapunov functional, they studied the global stability of endemic equilibrium of the model. Up to now, many authors have studied different kinds of SIRS epidemic models with time delays and a significant body of work has been carried out (see, for example, [13–15] and the references cited therein). We note that the spatial content of the environment has been ignored in the models aforementioned. The models have been traditionally formulated in relation to the time evolution of uniform population distributions in the habitat and are as such governed by ordinary differential equations. However, due to the large mobility of people within a country or even worldwide, spatially uniform models are not sufficient to give a realistic picture of a disease’s diffusion. For this reason, the spatial effects cannot be neglected in studying the spread of epidemics. Noble [16] applied reaction–diffusion theory to describe the spread of plague through Europe in the mid-14th century. By using the linear theory of semigroups, Saccomandi [17] investigated the existence and uniqueness of the solution for an SIR model with spatial inhomogeneity, nonlocal interactions, and an open population. In recent times, many investigators have introduced population movements into related equations for epidemiological modeling and simulations in efforts to understand the most basic features of spatially distributed interactions (see, for example, [18–22]). Motivated by the work of [2,11,16], in the present paper, we are concerned with the following delayed SIRS epidemic model with spatial diffusion

 ∂ 2S ∂S   = D + A − dS (x, t ) − β S (x, t )I (x, t − τ ) + δ R(x, t ),  S   ∂t ∂ x2   2 ∂I ∂ I = DI 2 + β S (x, t )I (x, t − τ ) − (γ + α + d)I (x, t ),  ∂t ∂x     ∂ 2R  ∂R  = DR 2 + γ I (x, t ) − (δ + d)R(x, t ), ∂t ∂x

(1.2)

with initial conditions S (x, t ) = ρ1 (x, t ), I (x, t ) = ρ2 (x, t ), R(x, t ) = ρ3 (x, t ),

¯. t ∈ [−τ , 0], x ∈ Ω

(1.3)

In problem (1.2)–(1.3), the positive constants DS , DI and DR denote the corresponding diffusion rates for the susceptible, infected and removed populations, respectively; Ω is a bounded domain in Rn with a smooth boundary ∂ Ω . The functions ¯ . In this paper, we assume ρi (x, t ) (i = 1, 2, 3) are nonnegative and Hölder continuous and satisfy ∂ρi /∂ x = 0 in [−τ , 0]× Ω that DS = DI = DR = D. In the biological context, it is important to analyse the epidemic waves which are described by travelling wave solutions propagating with a certain speed. In this paper, we are interested in the existence of travelling wave solutions to the SIRS epidemic model (1.2). The main tool to study the existence of travelling wave solutions for the reaction–diffusion equations with delays is the sub- and supersolution technique due to [23]. Wu and Zou [24,25] studied the existence of travelling wavefronts for delayed reaction–diffusion systems with reaction terms satisfying the so-called quasi-monotonicity or exponential quasi-monotonicity conditions, where the well-known monotone iteration techniques for elliptic systems with advanced arguments [26,27] were used. However, we note that the nonlinear reaction terms of system (1.2) do not satisfy either the quasi-monotonicity or the exponential quasi-monotonicity conditions. Therefore, the method of upper–lower solutions and its associated monotone iteration scheme developed by Wu and Zou [24,25] cannot be used to study the existence of travelling wave solutions to system (1.2). Recently, by constructing a pair of desirable upper–lower solutions, [28] got a subset, and employed the Schauder’s fixed point theorem in this subset to investigate the existence of travelling wave solutions of a class of delayed reaction diffusion systems with two equations in which the nonlinear reaction terms satisfy the partial quasimonotonicity and partial exponential quasi-monotonicity, respectively. Li et al. [29] investigated the existence of travelling wave solutions of a class of delayed reaction diffusion systems with two equations in which the reaction terms satisfy weak quasi-monotonicity and weak exponential quasi-monotonicity conditions, respectively. Clearly, all the results above cannot directly be applied to a system with more than three equations. Therefore, it remains important and challenging to study the existence of travelling wave solutions for delayed reaction diffusion systems with more than three equations in which the nonlinear reaction terms do not satisfy either the quasi-monotonicity or the exponential quasi-monotonicity conditions. The organisation of this paper is as follows. In the next section, we give some preliminaries. In Section 3, by analysing the corresponding characteristic equations, we analyse the local stability of steady states to system (1.2). In Section 4, we employ the cross iteration method and the Schauder’s fixed point theorem in a profile set to establish the existence of travelling wave solutions to a generalised system with three equations in which the nonlinear reaction terms do not satisfy either the quasi-monotonicity or the exponential quasi-monotonicity conditions. In Section 5, by constructing a pair of

54

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

upper–lower solutions, we use the results derived in Section 4 to prove the existence of travelling wave solutions to system (1.2). In Section 6, numerical simulations are carried out to illustrate the main results. A brief discussion is given in Section 7 to end this work. 2. Preliminaries Throughout this paper, we adopt the usual notations for the standard ordering in R3 . Thus, for u = (u1 , u2 , u3 )T and v = (v1 , v2 , v3 )T , we denote u ≤ v if ui ≤ vi , i = 1, 2, 3; u < v if u ≤ v but u ̸= v ; and u ≪ v if u ≤ v but ui ̸= vi , i = 1, 2, 3. If u ≤ v , we also denote (u, v] = {w ∈ R3 : u < w ≤ v}, [u, v) = {w ∈ R3 : u ≤ w < v}, and [u, v] = {w ∈ R3 : u ≤ w ≤ v}. We use | · | to denote the Euclidean norm in R3 and ‖ · ‖ to denote the supremum norm in C ([−τ , 0], R3 ). Throughout this paper, we study the existence of travelling wave solutions in x ∈ R. Before proving the existence of travelling wave solutions to system (1.2), we first investigate the following general delayed reaction–diffusion system:

 ∂ 2u ∂u   = D 2 + f1 (ut (x), vt (x), wt (x)),    ∂t ∂x   ∂v ∂ 2v = D 2 + f2 (ut (x), vt (x), wt (x)),  ∂t ∂x     ∂w ∂ 2w   = D 2 + f3 (ut (x), vt (x), wt (x)). ∂t ∂x

(2.1)

Corresponding to (2.1), we make the following hypotheses: (A1) There exist constants k1 , k2 , k3 > 0 such that fi (0, 0, 0) = fi (k1 , k2 , k3 ) = 0, i = 1, 2, 3. (A2) There exist three positive constants Li > 0 (i = 1, 2, 3) such that

|f1 (φ1 , ϕ1 , ψ1 ) − f1 (φ2 , ϕ2 , ψ2 )| ≤ L1 ‖Φ − Ψ ‖, |f2 (φ1 , ϕ1 , ψ1 ) − f2 (φ2 , ϕ2 , ψ2 )| ≤ L2 ‖Φ − Ψ ‖, |f3 (φ1 , ϕ1 , ψ1 ) − f3 (φ2 , ϕ2 , ψ2 )| ≤ L3 ‖Φ − Ψ ‖, for Φ = (φ1 , ϕ1 , ψ1 ), Ψ = (φ2 , ϕ2 , ψ2 ) ∈ C ([−τ , 0], R3 ) with 0 ≤ φi (s) ≤ M1 , 0 ≤ ϕi (s) ≤ M2 , 0 ≤ ψi (s) ≤ M3 , i = 1, 2. On substituting u(t , x) = φ(x + ct ), v(t , x) = ϕ(x + ct ), w(t , x) = ψ(x + ct ) and denote the travelling wave coordinate x + ct still by t, we derive from (2.1) that

 Dφ ′′ (t ) − c φ ′ (t ) + fc1 (φt , ϕt , ψt ) = 0, Dϕ ′′ (t ) − c ϕ ′ (t ) + fc2 (φt , ϕt , ψt ) = 0, Dψ ′′ (t ) − c ψ ′ (t ) + f (φ , ϕ , ψ ) = 0 c3 t t t

(2.2)

satisfying the following partial quasi-monotonicity conditions (PQM): (PQM) There exist three positive constants β1 , β2 , β3 > 0 such that fc1 (φ1 , ϕ1 , ψ1 ) − fc1 (φ2 , ϕ2 , ψ2 ) + β1 [φ1 (0) − φ2 (0)] ≥ 0, fc2 (φ1 , ϕ1 , ψ1 ) − fc2 (φ1 , ϕ2 , ψ1 ) + β2 [ϕ1 (0) − ϕ2 (0)] ≥ 0, fc2 (φ1 , ϕ1 , ψ1 ) − fc2 (φ2 , ϕ1 , ψ1 ) ≤ 0, fc2 (φ1 , ϕ1 , ψ1 ) − fc2 (φ1 , ϕ1 , ψ2 ) ≤ 0, fc3 (φ1 , ϕ1 , ψ1 ) − fc3 (φ2 , ϕ2 , ψ2 ) + β3 [ψ1 (0) − ψ2 (0)] ≥ 0,

(2.3)

where φi , ϕi , ψi ∈ C ([−τ , 0], R), i = 1, 2, with 0 ≤ φ2 (s) ≤ φ1 (s) ≤ M1 , 0 ≤ ϕ2 (s) ≤ ϕ1 (s) ≤ M2 , 0 ≤ ψ2 (s) ≤ ψ1 (s) ≤ M3 , s ∈ [−τ , 0], Mj > kj (j = 1, 2, 3) are positive constants, φt (ζ ) = φ(ζ + t ), ϕt (ζ ) = ϕ(ζ + t ), ψt (ζ ) = ψ(ζ + t ), and the functions fci : Xc = C ([−c τ , 0]; R3 ) → R3 (i = 1, 2, 3) are defined by fci (φ, ϕ, ψ) = fi (φ c , ϕ c , ψ c ),

ϕ (s) = ϕ(cs), c

φ c (s) = φ(cs), ψ (s) = ψ(cs), s ∈ [−τ , 0], c

i = 1, 2, 3.

If, for some c > 0, system (2.1) has a solution defined on R3 satisfying lim φ(t ) = φ− ,

lim ϕ(t ) = ϕ− ,

lim ψ(t ) = ψ− ,

t →−∞

t →−∞

t →−∞

t →+∞

t →+∞

t →+∞

lim φ(t ) = φ+ ,

lim ϕ(t ) = ϕ+ ,

lim ψ(t ) = ψ+ ,

(2.4)

where (φ− , ϕ− , ψ− ) and (φ+ , ϕ+ , ψ+ ) are steady states of (2.1). Then u(t , x) = φ(x + ct ), v(t , x) = ϕ(x + ct ), w(t , x) = ψ(x + ct ) is called a travelling wave solution of system (2.1) with speed c. Without loss of generality, we can assume (φ− , ϕ− , ψ− )=(0,0,0) and (φ+ , ϕ+ , ψ+ ) = (k1 , k2 , k3 ), and seek for a travelling wave solution connecting these two steady states.

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

55

In the next section, we will apply Schauder’s fixed point theorem, which requires the continuity of the operator under consideration. For this purpose, we need to introduce a topology in C (R, R3 ). Let µ > 0 and equip C (R, R3 ) with the exponential decay norm defined by

|Φ |µ = sup e−µ|t | |Φ (t )|R3 . t ∈R

Define Bµ (R, R3 ) = {Φ ∈ C (R, R3 ) : |Φ |µ < ∞}. Then it is easy to check that (Bµ (R, R3 ), | · |µ ) is a Banach space. We look for travelling wave solutions to system (2.1) in the following profile set:

  Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)) = φ(t ) ≤ φ(t ) ≤ φ(t ), ϕ(t ) ≤ ϕ(t ) ≤ ϕ(t ), ψ(t ) ≤ ψ(t ) ≤ ψ(t ) . Obviously, Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)) is non-empty, convex, closed and bounded. We also need the following definition of upper and lower solutions to system (2.1). Definition 2.1. A pair of continuous functions Φ = (φ, ϕ, ψ) and Φ = (φ, ϕ, ψ) is called a pair of upper–lower solutions of system (2.1) if Φ and Φ are twice differentiable almost everywhere in R and they are essentially bounded on R, and there hold ′′

′

Dφ (t ) − c φ (t ) + fc1 (φ t , ϕ t , ψ t ) ≤ 0,

a.e. in R,

(2.5)

Dϕ (t ) − c ϕ (t ) + fc2 (φ , ϕ t , ψ ) ≤ 0,

a.e. in R,

(2.6)

′′

′

t

′′

t

′

Dψ (t ) − c ψ (t ) + fc3 (φ t , ϕ t , ψ t ) ≤ 0,

a.e. in R,

(2.7)

and Dφ ′′ (t ) − c φ ′ (t ) + fc1 (φ , ϕ , ψ ) ≥ 0,

a.e. in R,

(2.8)

Dϕ (t ) − c ϕ (t ) + fc2 (φ t , ϕ , ψ t ) ≥ 0,

a.e. in R,

(2.9)

t

t

′′

t

′

t

Dψ ′′ (t ) − c ψ ′ (t ) + fc3 (φ , ϕ , ψ ) ≥ 0, t

t

t

a.e. in R.

(2.10)

Unlike the standard upper and lower solutions defined in [24], fc2 is evaluated in a cross iteration scheme given in (2.6) and (2.9). 3. Local stability In this section, by analysing the corresponding characteristic equations, we discuss the local stability of a disease-free steady state and an endemic steady state to system (1.2) with the initial conditions (1.3) and the homogeneous Neumann boundary conditions

∂ S (x, t ) ∂ I (x, t ) ∂ R(x, t ) = = = 0, t ≥ 0, x ∈ ∂ Ω , (3.1) ∂n ∂n ∂n respectively, where ∂/∂ n denotes the outward normal derivative on ∂ Ω , the homogeneous Neumann boundary conditions imply that the populations do not move across the boundary ∂ Ω . System (1.2) always has a disease-free steady state E1 (A/d, 0, 0). Further, if Aβ > d(γ + α + d), then system (1.2) has a unique endemic steady state E ∗ (S ∗ , I ∗ , R∗ ), where γ +α+d S∗ = , β (δ + d)[Aβ − d(γ + α + d)] I∗ = , β[(δ + d)(γ + α + d) − δγ ] γ [Aβ − d(γ + α + d)] R∗ = . β[(δ + d)(γ + α + d) − δγ ] Let

R0 =

Aβ d(γ + α + d)

.

R0 is called the basic reproductive number (sometimes called basic reproductive rate or basic reproductive ratio) of system (1.2), which describes the average number of newly infected cells generated from one infected cell at the beginning of the

56

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

infectious process. This quantity determines the thresholds for disease transmissions. It is easy to show that if R0 > 1, the endemic steady state E ∗ exists; if R0 < 1, E ∗ is not feasible. Let 0 = µ1 < µ2 < · · · be the eigenvalues of the operator −∆ on Ω with the homogeneous Neumann boundary conditions, and E (µi ) be the eigenspace corresponding to µi in C 1 (Ω ). Let X = [C 1 (Ω )]3 , {φij ; j = 1, . . . , dim E (µi )} be an orthonormal basis of E (µi ), and Xij = {c φij |c ∈ R3 }. Then [30–32]

X=

∞ 

dimE (µi )

Xi

and



Xi =

i =0

Xij .

j =1

Let D = diag(DS , DI , DR ), Z = (S , I , R), LZ = D ∆Z + G(Eˆ )Z , where

−d − β I 0 G(Eˆ )Z =  β I 0 

0

δ

0

−(γ + α + d) γ





× 0 −(δ + d)

S (x, t ) I (x, t ) R(x, t )





0 + 0 0

−β S 0 βS0 0

0 S (x, t − τ ) 0  I ( x, t − τ ) , R(x, t − τ ) 0





and Eˆ (S 0 , I 0 , R0 ) represents any feasible uniform steady state of system (1.2). The linearization of system (1.2) at Eˆ is of the form Zt = LZ . For each i ≥ 1, Xi is invariant under the operator L, and λ is an eigenvalue of L if and only if it is an eigenvalue of −µi D + G(Eˆ ) for some i ≥ 1, in which case, there is an eigenvector in Xi . The characteristic equation of −µi D + G(E1 ) takes the form

  Aβ −λτ (λ + µi D + d)(λ + µi D + δ + d) λ + µi D + γ + α + d − e = 0. d

(3.2)

Clearly, for any i ≥ 1, Eq. (3.2) always has two negative real roots −(µi D + d) and −(µi D + δ + d). Its other roots are determined by the following equation

λ + µi D + γ + α + d −

Aβ −λτ e = 0. d

(3.3)

Let g1 (λ) = λ + µi D + γ + α + d −

Aβ −λτ e . d

If R0 > 1, then for i = 1, we have that, for λ real, g1 (0) = γ + α + d −

Aβ d

< 0,

lim g1 (λ) = ∞.

λ→∞

Hence, when i = 1, (3.3) has a positive real root. Therefore, there is a characteristic root λ with a positive real part in the spectrum of L. Accordingly, if R0 > 1, E1 (A/d, 0, 0) is unstable. If R0 < 1, it is easy to see that the roots of (3.3) are negative constants when τ = 0, then E1 (A/d, 0, 0) is locally asymptotically stable when τ = 0. If iσ (σ > 0) is a solution of (3.3), separating the real and imaginary parts, we derive that d(µi D + γ + α + d) = Aβ cos σ τ , dσ = −Aβ sin σ τ .



(3.4)

Squaring and adding the two equations of (3.4), it follows that d2 σ 2 + d2 (µi D + γ + α + d)2 − A2 β 2 = 0.

(3.5)

Note that if R0 < 1, d(µi D + γ + α + d) − Aβ > 0. Hence, Eq. (3.5) has no positive roots for all i ≥ 1. Therefore, if R0 < 1, then E1 (A/d, 0, 0) is locally asymptotically stable for all τ > 0. The characteristic equation of −µi D + G(E ∗ ) is of the form

λ3 + p2 λ2 + p1 λ + p0 + (q2 λ2 + q1 λ + q0 )e−λτ = 0, where p0 = (µi D + d + β I ∗ )(µi D + γ + α + d)(µi D + δ + d) − βγ δ I ∗ , p1 = (µi D + d + β I ∗ )(µi D + γ + α + d) + (µi D + d + β I ∗ )(µi D + δ + d) + (µi D + γ + α + d)(µi D + δ + d), p2 = (µi D + d + β I ∗ ) + (µi D + γ + α + d) + (µi D + δ + d), q0 = −dβ S ∗ (δ + d), q1 = −β S ∗ (δ + 2d), q2 = −β S ∗ .

(3.6)

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

57

When τ = 0, Eq. (3.6) becomes

λ3 + (p2 + q2 )λ2 + (p1 + q1 )λ + p0 + q0 = 0. (3.7) It is easy to verify that for all i ≥ 1 p2 + q2 = (µi D + d + β I ∗ ) + µi D + (µi D + δ + d) > 0, p1 + q1 ≥ (µi D + d + β I ∗ )(µi D + δ + d) + (2µi D + β I ∗ )(µi D + γ + α + d) > 0, p0 + q0 ≥ β(γ + α + d)(δ + d)I ∗ − βγ δ I ∗ > 0, (p1 + q1 )(p2 + q2 ) − (p0 + q0 ) > 0. Then it follows from the Routh–Hurwitz criterion that all roots λi,1 , λi,2 and λi,3 of Eq. (3.7) have negative parts. Next we can conclude that there exists a positive constant ξ such that Re{λi,1 }, Re{λi,2 }, Re{λi,3 } ≤ −ξ , i ≥ 1. (3.8) In fact, let λ = µi ν , then λ3 + (p2 + q2 )λ2 + (p1 + q1 )λ + p0 + q0 = µ3i ν 3 + (p2 + q2 )µ2i ν 2 + (p1 + q1 )µi ν + p0 + q0 , hi (ν). Since µi → ∞ as i → ∞, it yields that hi (ν) = ν 3 + 3Dν 2 + 3D2 ν + D3 , h(ν). lim i→∞ µ3 i It is easy to see that these three roots ν1 , ν2 and ν3 of h(ν) have negative real parts from the Routh–Hurwitz criterion again, and thus there exists a positive constant ξ¯ such that Re{ν1 }, Re{ν2 }, Re{ν3 } ≤ −ξ¯ . By continuity, we see that there exists i0 such that the three roots νi,1 , νi,2 and νi,3 of hi (ν) satisfy ξ¯ Re{νi,1 }, Re{νi,2 }, Re{νi,3 } ≤ − , i ≥ i0 , 2

and thus Re{λi,1 }, Re{λi,2 }, Re{λi,3 } ≤ −

µi ξ¯ 2

≤−

µi0 ξ¯ 2

,

i ≥ i0 .

Let

  −ξ˜ = max Re{λi,1 }, Re{λi,2 }, Re{λi,3 } , 1≤i≤i0

then ξ˜ , and (3.8) holds for ξ ∈ (ξ˜ , µi0 ξ¯ /2).

Consequently, the spectrum of L, consisting only of eigenvalues, lies in Reλ ≤ −ξ . Hence, the endemic steady state E ∗ is locally asymptotically stable when τ = 0. If iω(ω > 0) is a solution of (3.6), separating real and imaginary parts, we derive that



ω3 − p1 ω = q1 ω cos ωτ + (q2 ω2 − q0 ) sin ωτ , p2 ω2 − p0 = q1 ω sin ωτ − (q2 ω2 − q0 ) cos ωτ .

(3.9)

Squaring and adding the two equations of (3.9), it follows that

ω6 + (p22 − 2p1 − q22 )ω4 + (p21 − 2p0 p2 − q21 + 2q0 q2 )ω2 + p20 − q20 = 0.

(3.10)

Let z = ω , then Eq. (3.10) becomes 2

z 3 + (p22 − 2p1 − q22 )z 2 + (p21 − 2p0 p2 − q21 + 2q0 q2 )z + p20 − q20 = 0.

(3.11)

In the following, we need to seek the conditions under which (3.11) has at least one positive root. Denote h(z ) = z 3 + (p22 − 2p1 − q22 )z 2 + (p21 − 2p0 p2 − q21 + 2q0 q2 )z + p20 − q20 .

(3.12)

By calculation it follows that for all i ≥ 1 p22 − 2p1 − q22 = (µi D + d + β I ∗ )2 + µ2i D2 + 2(γ + α + d)µi D + (µi D + δ + d)2 > 0, p21 − 2p0 p2 − q21 + 2q0 q2 ≥ (µi D + d + β I ∗ )2 (µi D + γ + α + d)2 + (µi D + d + β I ∗ )2 (µi D + δ + d)2 > 0, p20 − q20 = (p0 − q0 )(p0 + q0 ) > 0. Hence, we know that Eq. (3.12) has no positive roots. Therefore, if R0 > 1, E ∗ is locally asymptotically stable for all τ > 0. We therefore obtain the following results. Theorem 3.1. For system (1.2) with initial conditions (1.3) and homogeneous Neumann boundary conditions (3.1), we have (i) If R0 > 1, the disease-free steady state E1 (A/d, 0, 0) is unstable; if R0 < 1, E1 is locally asymptotically stable. (ii) Let R0 > 1, the endemic steady state E ∗ (S ∗ , I ∗ , R∗ ) is asymptotically stable for all τ > 0.

58

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

4. Existence of travelling waves for system (2.1) In this section, we study the existence of travelling wave solutions to system (2.1) with nonlinear reaction terms satisfying (PQM). We assume that a pair of upper–lower solutions (φ(t ), ϕ(t ), ψ(t )) and (φ(t ), ϕ(t ), ψ(t )) are given such that (P1) (0, 0, 0) ≤ (φ(t ), ϕ(t ), ψ(t )) ≤ (φ(t ), ϕ(t ), ψ(t )) ≤ (M1 , M2 , M3 ), t ∈ R, (P2) limt →−∞ (φ(t ), ϕ(t ), ψ(t )) = (0, 0, 0), limt →+∞ (φ(t ), ϕ(t ), ψ(t )) = (k1 , k2 , k3 ), ′

′

(P3) sups≤t φ(t ) ≤ φ(t ); φ (t +) ≤ φ (t −) and φ ′ (t +) ≥ φ ′ (t −) for t ∈ R. For the constants β1 , β2 , β3 > 0 in (PQM), define H : C (R, R3 ) → C (R, R3 ) by H1 (φ, ϕ, ψ)(t ) = fc1 (φt , ϕt , ψt ) + β1 φ(t ),

φ, ϕ, ψ ∈ C (R, R),

(4.1)

H2 (φ, ϕ, ψ)(t ) = fc2 (φt , ϕt , ψt ) + β2 ϕ(t ),

φ, ϕ, ψ ∈ C (R, R),

(4.2)

H3 (φ, ϕ, ψ)(t ) = fc3 (φt , ϕt , ψt ) + β3 ψ(t ),

φ, ϕ, ψ ∈ C (R, R).

(4.3)

The operators H1 , H2 and H3 admit the following properties: Lemma 4.1. Assume that (A1) and (PQM) hold, then H2 (φ1 , ϕ1 , ψ1 )(t ) ≤ H2 (φ2 , ϕ1 , ψ1 )(t ), H2 (φ1 , ϕ2 , ψ1 )(t ) ≤ H2 (φ1 , ϕ1 , ψ1 )(t ), H2 (φ1 , ϕ1 , ψ1 )(t ) ≤ H2 (φ1 , ϕ1 , ψ2 )(t ), for t ∈ R with 0 ≤ φ2 (t ) ≤ φ1 (t ) ≤ M1 , 0 ≤ ϕ2 (t ) ≤ ϕ1 (t ) ≤ M2 , 0 ≤ ψ2 (t ) ≤ ψ1 (t ) ≤ M3 . Proof. By (PQM), direct calculation shows that H2 (φ1 , ϕ1 , ψ1 )(t ) − H2 (φ2 , ϕ1 , ψ1 )(t ) = fc2 (φ1t , ϕ1t , ψ1t ) − fc2 (φ2t , ϕ1t , ψ1t ) ≤ 0, H2 (φ1 , ϕ1 , ψ1 )(t ) − H2 (φ1 , ϕ2 , ψ1 )(t ) = fc2 (φ1t , ϕ1t , ψ1t ) − fc2 (φ1t , ϕ2t , ψ1t ) + β2 [ϕ1 (t ) − ϕ2 (t )] ≥ 0, H2 (φ1 , ϕ1 , ψ1 )(t ) − H2 (φ1 , ϕ1 , ψ2 )(t ) = fc2 (φ1t , ϕ1t , ψ1t ) − fc2 (φ1t , ϕ1t , ψ2t ) ≤ 0. This completes the proof.



Lemma 4.2 ([24]). Assume that (A1) and (PQM) hold. Then for any (0, 0, 0) ≤ (φ, ϕ, ψ) ≤ (M1 , M2 , M3 ), we have (i) H1 (φ, ϕ, ψ)(t ), H3 (φ, ϕ, ψ)(t ) ≥ 0, t ∈ R. (ii) H1 (φ2 , ϕ2 , ψ2 )(t ) ≤ H1 (φ1 , ϕ1 , ψ1 )(t ) and H3 (φ2 , ϕ2 , ψ2 )(t ) ≤ H3 (φ1 , ϕ1 , ψ1 )(t ) for t ∈ R with 0 ≤ φ2 (t ) ≤ φ1 (t ) ≤ M1 , 0 ≤ ϕ2 (t ) ≤ ϕ1 (t ) ≤ M2 , 0 ≤ ψ2 (t ) ≤ ψ1 (t ) ≤ M3 . In terms of H1 , H2 and H3 , system (5.2) can be rewritten as

 Dφ ′′ (t ) − c φ ′ (t ) − β1 φ(t ) + H1 (φ, ϕ, ψ)(t ) = 0, Dϕ ′′ (t ) − c ϕ ′ (t ) − β2 ϕ(t ) + H2 (φ, ϕ, ψ)(t ) = 0, Dψ ′′ (t ) − c ψ ′ (t ) − β ψ(t ) + H (φ, ϕ, ψ)(t ) = 0. 3 3

(4.4)

Define

λ1 = λ3 = λ5 =

c−



c 2 + 4β1 D

 2D c − c 2 + 4β2 D  2D c − c 2 + 4β3 D 2D

,

λ2 =

,

λ4 =

,

λ6 =

c+



c 2 + 4β1 D

 2D c + c 2 + 4β2 D  2D c + c 2 + 4β3 D 2D

, , .

Let CK (R, R3 ) = {(φ, ϕ, ψ) ∈ C (R, R3 ) : (0, 0, 0) ≤ (φ, ϕ, ψ) ≤ (M1 , M2 , M3 )}, and define F = (F1 , F2 , F3 ) : CK (R, R3 ) → C (R, R3 ) by F1 (φ, ϕ, ψ)(t ) = F2 (φ, ϕ, ψ)(t ) = F3 (φ, ϕ, ψ)(t ) =

1 D(λ2 − λ1 ) 1 D(λ4 − λ3 ) 1 D(λ6 − λ5 )

[∫

t

eλ1 (t −s) H1 (φ, ϕ, ψ)(s)ds +

∫

eλ3 (t −s) H2 (φ, ϕ, ψ)(s)ds +

∫

eλ5 (t −s) H3 (φ, ϕ, ψ)(s)ds +

∫

−∞

[∫

t

t

−∞

]

eλ2 (t −s) H1 (φ, ϕ, ψ)(s)ds ,

t

−∞

[∫

∞

∞

]

eλ4 (t −s) H2 (φ, ϕ, ψ)(s)ds ,

t

t

∞

]

eλ6 (t −s) H3 (φ, ϕ, ψ)(s)ds ,

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

59

for (φ, ϕ, ψ) ∈ CK (R, R3 ). It is easy to see that F1 (φ, ϕ, ψ), F2 (φ, ϕ, ψ) and F3 (φ, ϕ, ψ) satisfy

 DF1′′ (φ, ϕ, ψ) − cF1′ (φ, ϕ, ψ) − β1 F1 (φ, ϕ, ψ) + H1 (φ, ϕ, ψ) = 0, DF ′′ (φ, ϕ, ψ) − cF ′ (φ, ϕ, ψ) − β2 F2 (φ, ϕ, ψ) + H2 (φ, ϕ, ψ) = 0, DF2′′ (φ, ϕ, ψ) − cF2′ (φ, ϕ, ψ) − β F (φ, ϕ, ψ) + H (φ, ϕ, ψ) = 0. 3 3 3 3 3

(4.5)

Corresponding to Lemmas 4.1 and 4.2, we have the following result for F . Lemma 4.3. Assume that (A1) and (PQM) hold. For any (0, 0, 0) ≤ (φ, ϕ, ψ) ≤ (M1 , M2 , M3 ), we have F1 (φ2 , ϕ2 , ψ2 ) ≤ F1 (φ1 , ϕ1 , ψ1 ), F2 (φ1 , ϕ1 , ψ1 ) ≤ F2 (φ2 , ϕ1 , ψ1 ), F2 (φ1 , ϕ2 , ψ1 ) ≤ F2 (φ1 , ϕ1 , ψ1 ), F2 (φ1 , ϕ1 , ψ1 ) ≤ F2 (φ1 , ϕ1 , ψ2 ), F3 (φ2 , ϕ2 , ψ2 )(t ) ≤ F3 (φ1 , ϕ1 , ψ1 )(t ) for t ∈ R with 0 ≤ φ2 (t ) ≤ φ1 (t ) ≤ M1 , 0 ≤ ϕ2 (t ) ≤ ϕ1 (t ) ≤ M2 , 0 ≤ ψ2 (t ) ≤ ψ1 (t ) ≤ M3 . We next verify the continuity of F . Lemma 4.4. Assume that (A2) holds, then F = (F1 , F2 , F3 ) is continuous with respective to the norm | · |µ in Bµ (R, R3 ). Proof. For any fixed ϵ > 0, let σ < ϵ/(L1 eµc τ +β1 ), then for Φ = (φ1 , ϕ1 , ψ1 ), Ψ = (φ2 , ϕ2 , ψ2 ) ∈ Bµ (R, R3 ) with

|Φ − Ψ |µ = sup |Φ (t ) − Ψ (t )|e−µ|t | < δ, t ∈R

direct calculation shows that

|H1 (φ1 , ϕ1 , ψ1 ) − H1 (φ2 , ϕ2 , ψ2 )|e−µ|t | ≤ |f1 (φ1t , ϕ1t , ψ1t ) − f1 (φ2t , ϕ2t , ψ2t )|e−µ|t | + β1 |φ1 − φ2 |µ ≤ L1 ‖Φt − Ψt ‖Xc e−µ|t | + β1 |φ1 − φ2 |µ = L1 sup |Φ (s + t ) − Ψ (s + t )|e−µ|t | + β1 |φ1 − φ2 |µ s∈[−c τ ,0]

≤ L1 sup |Φ (s + t ) − Ψ (s + t )|e−µ|t +s| sup eµ|t +s| e−µ|t | + β1 |φ1 − φ2 |µ s∈[−c τ ,0]

s∈[−τ ,0]

−µ|t | µ|t | µc τ

≤ L1 |Φ − Ψ |µ e e e + β1 |Φ − Ψ |µ ≤ (L1 eµc τ + β1 )|Φ − Ψ |µ ≤ ϵ. For t > 0, we see that

|F1 (φ1 , ϕ1 , ψ1 ) − F1 (φ2 , ϕ2 , ψ2 )|e−µ|t | [∫ t 1 eλ1 (t −s) |H1 (φ1 , ϕ1 , ψ1 )(s) − H1 (φ2 , ϕ2 , ψ2 )(s)|ds = D(λ2 − λ1 ) −∞ ] ∫ ∞ + eλ2 (t −s) |H1 (φ1 , ϕ1 , ψ1 )(s) − H1 (φ2 , ϕ2 , ψ2 )(s)|ds e−µt t

[ ] ∫ 0 ∫ t ∫ ∞ ϵ eλ1 t e−(λ1 +µ)s ds + eλ1 t e(µ−λ1 )s ds + eλ2 t e(µ−λ2 )s ds e−µt D(λ2 − λ1 ) −∞ 0 t [ ] ϵ λ2 − λ1 2µ (λ1 −µ)t = e + D(λ2 − λ1 ) λ21 − µ2 (µ − λ1 )(λ2 − µ) [ ] ϵ λ2 − λ1 2µ ≤ + . D(λ2 − λ1 ) λ21 − µ2 (µ − λ1 )(λ2 − µ) ≤

Similarly, for t < 0, we have

|F1 (φ1 , ϕ1 , ψ1 ) − F1 (φ2 , ϕ2 , ψ2 )|e−µ|t | ≤

[ ] ϵ λ2 − λ1 2µ − , D(λ2 − λ1 ) λ22 − µ2 (λ1 + µ)(λ2 + µ)

which implies that F1 : Bµ (R, R3 ) → Bµ (R, R3 ) is continuous with respect to the norm | · |µ in Bµ (R, R3 ). By using a similar argument as above, it can be shown that F2 , F3 : Bµ (R, R3 ) → Bµ (R, R3 ) are continuous. Thus, we obtain that F = (F1 , F2 , F3 ) is continuous with respect to the norm | · |µ in Bµ (R, R3 ). This completes the proof.  Lemma 4.5. Assume that (A1) and (PQM) hold, then F (Γ ((φ, ϕ, ψ), (φ, ϕ, ψ))) ⊂ Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)).

60

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

Proof. For any (φ, ϕ, ψ) with (φ, ϕ, ψ) ≤ (φ, ϕ, ψ) ≤ (φ, ϕ, ψ), it follows from Lemma 4.3 that F1 (φ, ϕ, ψ) ≤ F1 (φ, ϕ, ψ) ≤ F1 (φ, ϕ, ψ), F2 (φ, ϕ, ψ) ≤ F2 (φ, ϕ, ψ) ≤ F2 (φ, ϕ, ψ),

(4.6)

F3 (φ, ϕ, ψ) ≤ F3 (φ, ϕ, ψ) ≤ F3 (φ, ϕ, ψ). By the definition of upper–lower solutions, we have ′′

′

Dφ (t ) − c φ (t ) − β1 φ(t ) + H1 (φ, ϕ, ψ)(t ) ≤ 0.

(4.7)

Choosing (φ, ϕ, ψ) = (φ, ϕ, ψ) in the first equation of (4.5), and denoting φ 1 (t ) = F1 (φ, ϕ, ψ)(t ), we get ′′

′

Dφ 1 (t ) − c φ 1 (t ) − β1 φ 1 (t ) + H1 (φ, ϕ, ψ)(t ) = 0.

(4.8)

Setting ω(t ) = φ 1 (t ) − φ(t ) and combining (4.7) and (4.8), we have Dω′′ (t ) − c ω′ (t ) − β1 ω(t ) ≥ 0.

(4.9)

Repeating the proof of Lemma 3.3 in [24] shows that ω(t ) ≤ 0, which implies that F1 (φ, ϕ, ψ) ≤ φ . By a similar argument, we know that F1 (φ, ϕ, ψ) ≥ φ, F2 (φ, ϕ, ψ) ≥ ϕ, F2 (φ, ϕ, ψ) ≤ ϕ, F3 (φ, ϕ, ψ) ≥ ψ, F3 (φ, ϕ, ψ) ≤ ψ , then F (Γ ((φ, ϕ, ψ), (φ, ϕ, ψ))) ⊂ Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)). This completes the proof.



Lemma 4.6. Assume that (PQM) holds, then F : Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)) → Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)) is compact. Proof. Noting that F1′ (φ, ϕ, ψ)(t ) =

λ1 eλ1 t D(λ2 − λ1 )

t

∫

e−λ1 s H1 (φ, ϕ, ψ)(s)ds + −∞

λ2 eλ2 t D(λ2 − λ1 )

∞

∫

e−λ2 s H1 (φ, ϕ, ψ)(s)ds, t

it follows from in Lemma 4.3 that F1′ (φ, ϕ, ψ)(t ) ≥ 0. By (i) in Lemma 4.2 and the fact that λ1 < 0 < λ2 , we have

∫ ∞ λ2 eλ1 t e−λ2 s H1 (φ, ϕ, ψ)(s)ds D(λ2 − λ1 ) t ∫ ∞ H1 (φ, ϕ, ψ)(t ) e−λ2 s ds

0 ≤ F1′ (φ, ϕ, ψ)(t ) ≤

≤ ≤

λ2 eλ2 t D(λ2 − λ1 ) 1 D(λ2 − λ1 )

t

H (φ, ϕ, ψ)(t ).

Hence, (P1) implies that there exists a constant N1 such that |F1′ (φ, ϕ, ψ)(t )|µ ≤ N1 . For any (φ, ϕ, ψ) ∈ Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)), F2′ (φ, ϕ, ψ)(t ) =

λ3 eλ3 t D(λ4 − λ3 )

∫

t

e−λ3 s H2 (φ, ϕ, ψ)(s)ds + −∞

λ4 eλ4 t D(λ4 − λ3 )

∞

∫

e−λ4 s H2 (φ, ϕ, ψ)(s)ds. t

Thus, we have

 ∫ t  λ3 eλ3 t |F2′ (φ, ϕ, ψ)(t )|µ = sup  e−λ3 s H2 (φ, ϕ, ψ)(s)ds t ∈R D(λ4 − λ3 ) −∞  ∫ ∞  λ4 eλ4 t + e−λ4 s H1 (φ, ϕ, ψ)(s)ds e−µ|t | D(λ4 − λ3 ) t ∫ t −λ3 ≤ sup eλ3 t −µ|t | e−λ3 s eµ|s| e−µ|s| H2 (φ, ϕ, ψ)(s)ds D(λ4 − λ3 ) t ∈R −∞ ∫ ∞ λ4 λ4 t −µ|t | + sup e e−λ4 s eµ|s| e−µ|s| H2 (φ, ϕ, ψ)(s)ds D(λ4 − λ3 ) t ∈R t ∫ t −λ3 ≤ |H2 (φ, ϕ, ψ)(t )|µ sup eλ3 t −µ|t | e−λ3 s eµ|s| ds D(λ4 − λ3 ) t ∈R −∞ ∫ ∞ λ4 λ4 t −µ|t | + |H2 (φ, ϕ, ψ)(t )|µ sup e e−λ4 s eµ|s| ds. D(λ4 − λ3 ) t ∈R t

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

61

If t > 0, we get

] [∫ 0 ∫ t −λ3 (µ−λ3 )s −(λ3 +µ)s (λ3 −µ)t e ds |H2 (φ, ϕ, ψ)(t )|µ e ds + |F2 (φ, ϕ, ψ)(t )|µ ≤ sup e D(λ4 − λ3 ) t ∈R 0 −∞ ∫ ∞ λ4 e(µ−λ4 )s ds|H2 (φ, ϕ, ψ)(t )|µ + sup e(λ4 −µ)t D(λ4 − λ3 ) t ∈R t λ4 λ3 |H2 (φ, ϕ, ψ)(t )|µ + |H2 (φ, ϕ, ψ)(t )|µ ≤ D(µ + λ3 )(λ4 − λ3 ) D(λ4 − µ)(λ4 − λ3 ) [ ] λ3 λ4 1 + |H2 (φ, ϕ, ψ)(t )|µ . = D(λ4 − λ3 ) µ + λ3 λ4 − µ ′

If t < 0, we obtain

−λ3 sup e(λ3 +µ)t |F2 (φ, ϕ, ψ)(t )|µ ≤ D(λ4 − λ3 ) t ∈R ′

∫

t

e−(λ3 +µ)s ds|H2 (φ, ϕ, ψ)(t )|µ −∞

] [∫ 0 ∫ ∞ λ4 (µ−λ4 )s −(µ+λ4 )s (λ4 +µ)t e ds |H2 (φ, ϕ, ψ)(t )|µ e ds + + sup e D(λ4 − λ3 ) t ∈R 0 t λ3 λ4 ≤ |H2 (φ, ϕ, ψ)(t )|µ + |H2 (φ, ϕ, ψ)(t )|µ D(µ + λ3 )(λ4 − λ3 ) D(λ4 − µ)(λ4 − λ3 ) [ ] 1 λ3 λ4 = + |H2 (φ, ϕ, ψ)(t )|µ . D(λ4 − λ3 ) µ + λ3 λ4 − µ Noting that (0, 0, 0) ≤ (φ, ϕ, ψ) ≤ (M1 , M2 , M3 ), it follows from Lemma 4.1 that |H2 (φ, ϕ, ψ)(t )|µ is bounded by a positive number. Therefore, there exists a constant N2 such that |F2′ (φ, ϕ, ψ)(t )|µ ≤ N2 . Similar to the proof of F1 , we have that there exists a constant N3 such that |F3′ (φ, ϕ, ψ)(t )|µ ≤ N3 . The above estimate for F ′ shows that F (Γ ((φ, ϕ, ψ), (φ, ϕ, ψ))) is equicontinuous. It follows from the proof of Lemma 4.5 that F (Γ ((φ, ϕ, ψ), (φ, ϕ, ψ))) is uniformly bounded. Next, we define F (φ, ϕ, ψ)(t ), F (φ, ϕ, ψ)(n), F (φ, ϕ, ψ)(−n),

 F (φ, ϕ, ψ)(t ) = n

t ∈ [−n, n], t ∈ (n, ∞), t ∈ (−∞, −n).

Then, for each n ≥ 1, F n (Γ ((φ, ϕ, ψ), (φ, ϕ, ψ))) is also equicontinuous and uniformly bounded. Now, in the interval

[−n, n], it follows from the Ascoli–Arzela Theorem that F n is compact. On the other hand, F n → F in Bµ (R, R3 ) as n → ∞, since sup |F n (φ, ϕ, ψ)(t ) − F (φ, ϕ, ψ)(t )|e−µ|t | = t ∈R

sup

t ∈(−∞,−n)

(n,∞)

|F n (φ, ϕ, ψ)(t ) − F (φ, ϕ, ψ)(t )|e−µ|t |

≤ 2K e−µn → 0, n → ∞. By Proposition 2.12 in [33], we have that F : Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)) → Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)) is compact. This completes the proof.  Theorem 4.1. Assume that (A1), (A2) and (PQM) hold. Suppose there is a pair of upper–lower solutions Φ = (φ, ϕ, ψ) and Ψ = (φ, ϕ, ψ) for (2.1) satisfying (P1), (P2) and (P3). Then, system (2.1) has a travelling wave solution. Proof. Combining Lemmas 4.1–4.6 with the Schauder’s fixed point theorem, we know that there exists a fixed point (φ ∗ (t ), ϕ ∗ (t ), ψ ∗ (t )) of F in Γ ((φ, ϕ, ψ), (φ, ϕ, ψ)), which gives a solution to (2.1). In order to prove that the solution is a travelling wave solution, we need to verify the asymptotic boundary conditions (2.4). By (P2) and the fact that

(0, 0, 0) ≤ (φ(t ), ϕ(t ), ψ(t )) ≤ (φ ∗ (t ), ϕ ∗ (t ), ψ ∗ (t )) ≤ (φ(t ), ϕ(t ), ψ(t )) ≤ (M1 , M2 , M3 ), we get that lim (φ ∗ (t ), ϕ ∗ (t ), ψ ∗ (t )) = (0, 0, 0),

t →−∞

and lim (φ ∗ (t ), ϕ ∗ (t ), ψ ∗ (t )) = (k1 , k2 , k3 ).

t →∞

62

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

Therefore, the fixed point (φ ∗ (t ), ϕ ∗ (t ), ψ ∗ (t )) satisfies the asymptotic boundary conditions (2.4). This completes the proof.  5. Existence of travelling waves for system (1.2) In this section, we use the results developed in Section 4 to study the existence of travelling wave solutions to system (1.2). Denoting N = S + I + R, then system (1.2) is equivalent to the following system

 ∂ 2N ∂N   = D 2 + A − α I (x, t ) − dN (x, t ),    ∂t ∂x   ∂I ∂ 2I = D 2 + β(N (x, t ) − I (x, t ) − R(x, t ))I (x, t − τ ) − (γ + α + d)I (x, t ),  ∂t ∂x     ∂ R ∂ 2R   = D 2 + γ I (x, t ) − (δ + d)R(x, t ). ∂t ∂x By making a change of variables N˜ = A/d − N , ˜I = I, R˜ = R and dropping the tildes, system (5.1) becomes  ∂ 2N ∂N   = D 2 + α I (x, t ) − dN (x, t ),    ∂t ∂x   ∂I ∂ 2I A = D 2 + β( − N (x, t ) − I (x, t ) − R(x, t ))I (x, t − τ ) − (γ + α + d)I (x, t ),  ∂ t ∂ x d    ∂R 2  ∂ R   = D 2 + γ I (x, t ) − (δ + d)R(x, t ). ∂t ∂x It is easy to show that if R0 > 1, system (5.2) has two steady states (0, 0, 0) and (k1 , k2 , k3 ), where k1 = A/d − S ∗ − I ∗ − R∗ ,

k2 = I ∗ ,

(5.1)

(5.2)

k3 = R∗ .

On substituting N (x, t ) = φ(x + ct ), I (x, t ) = ϕ(x + ct ), R(x, t ) = ψ(x + ct ) and denote the travelling wave coordinate x + ct still by t, we derive from (5.2) that

 Dφ ′′ (t ) − c φ ′ (t ) + αϕ(t ) − dφ(t ) = 0, Dϕ ′′ (t ) − c ϕ ′ (t ) + β(A/d − φ(t ) − ϕ(t ) − ψ(t ))ϕ(t − c τ ) − (γ + α + d)ϕ(t ) = 0, Dψ ′′ (t ) − c ψ ′ (t ) + γ ϕ(t ) − (δ + d)ψ(t ) = 0,

(5.3)

satisfying the following asymptotic boundary conditions lim (φ(t ), ϕ(t ), ψ(t )) = (0, 0, 0),

t →−∞

lim (φ(t ), ϕ(t ), ψ(t )) = (k1 , k2 , k3 ).

t →+∞

Lemma 5.1. The nonlinear reaction terms of system (5.2) satisfy (PQM). Proof. For any φi , ϕi , ψi ∈ C ([−τ , 0], R), i = 1, 2, with 0 ≤ φ2 (s) ≤ φ1 (s) ≤ M1 , 0 ≤ ϕ2 (s) ≤ ϕ1 (s) ≤ M2 , 0 ≤ ψ2 (s) ≤ ψ1 (s) ≤ M3 , s ∈ [−τ , 0], we have fc1 (φ1t , ϕ1t , ψ1t ) − fc1 (φ2t , ϕ2t , ψ2t ) = αϕ1 (0) − dφ1 (0) − αϕ2 (0) + dφ2 (0)

≥ −d(φ1 (0) − φ2 (0)). Let β1 = d > 0. We derive that fc1 (φ1 , ϕ1 , ψ1 ) − fc1 (φ2 , ϕ2 , ψ2 ) + β1 [φ1 (0) − φ2 (0)] ≥ 0. For fc2 (Nt (x), It (x), Rt (x)), it follows that fc2 (φ1t , ϕ1t , ψ1t ) − fc2 (φ1t , ϕ2t , ψ1t ) = β(A/d − φ1 (0) − ϕ1 (0) − ψ1 (0))ϕ1 (−τ ) − (γ + α + d)ϕ1 (0)

− β(A/d − φ1 (0) − ϕ2 (0) − ψ1 (0))ϕ2 (−τ ) + (γ + α + d)ϕ2 (0) ≥ −(β M2 + γ + α + d)(ϕ1 (0) − ϕ2 (0)). Let β2 = β M2 + γ + α + d > 0. We know that fc2 (φ1 , ϕ1 , ψ1 ) − fc2 (φ1 , ϕ2 , ψ1 ) + β2 [ϕ1 (0) − ϕ2 (0)] ≥ 0. fc2 (φ1t , ϕ1t , ψ1t ) − fc2 (φ2t , ϕ1t , ψ1t ) = β(A/d − φ1 (0) − ϕ1 (0) − ψ1 (0))ϕ1 (−τ ) − (γ + α + d)ϕ1 (0) − β(A/d − φ2 (0) − ϕ1 (0) − ψ1 (0))ϕ1 (−τ ) + (γ + α + d)ϕ1 (0) = β(φ2 (0) − φ1 (0))ϕ1 (−τ ) ≤ 0. fc2 (φ1t , ϕ1t , ψ1t ) − fc2 (φ1t , ϕ1t , ψ2t ) = β(A/d − φ1 (0) − ϕ1 (0) − ψ1 (0))ϕ1 (−τ ) − (γ + α + d)ϕ1 (0) − β(A/d − φ1 (0) − ϕ1 (0) − ψ2 (0))ϕ1 (−τ ) + (γ + α + d)ϕ1 (0) = β(ψ2 (0) − ψ1 (0))ϕ1 (−τ ) ≤ 0.

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

63

For fc3 (Nt (x), It (x), Rt (x)), we have fc3 (φ1t , ϕ1t , ψ1t ) − fc3 (φ2t , ϕ2t , ψ2t ) = γ ϕ1 (0) − (δ + d)ψ1 (0) − γ ϕ2 (0) + (δ + d)ψ2 (0)

≥ −(δ + d)(ψ1 (0) − ψ2 (0)). Let β3 = δ + d > 0. We obtain fc3 (φ1 , ϕ1 , ψ1 ) − fc3 (φ2 , ϕ2 , ψ2 ) + β3 [ψ1 (0) − ψ2 (0)] ≥ 0. This completes the proof.

√



Let c > c = 2 D(β A/d − (γ + α + d)). There exist ηi > 0 (i = 1, 2, 3, 4) such that ∗

 2 Dη1 − c η1 − d + α M2 /M1 = 0,    Dη2 − c η2 + β A/d − (γ + α + d) = 0, 2

Dη32 − c η3 + M2 γ /M3 − (δ + d) = 0,  2   Dη4 − c η4 − (γ + α + d) = 0, η1 , η3 ≤ η2 .

We find that there exist εi > 0 (i = 0, 1, . . . , 6) satisfying

d(k + ε ) − α M > ε , 1 1 2 0   (γ + α + d)(k2 + ε2 ) − β(A/d − k1 + ε4 − k2 − ε2 − k3 + ε6 )M2 > ε0 ,   (δ + d)(k3 + ε3 ) − γ M2 > ε0 , α(k2 − ε5 ) − d(k1 − ε4 ) > ε0 ,     β(A/d − M1 − k2 + ε5 − M3 ) − (γ + α + d) > ε0 , γ (k2 − ε5 ) − (δ + d)(k3 − ε6 ) > ε0 .

(5.4)

For the above constants, suitable constants ti (i = 1, 2, . . . , 7) and M0 = k2 eη4 t5 satisfying t2 ≥ max(t1 , t3 ), min(t4 , t7 ) ≥ t6 , t2 ≥ max(t4 , t7 ) and t6 − c τ > t5 , we define the continuous functions Φ (t ) = (φ1 (t ), ϕ1 (t ), ψ1 (t )) and Ψ (t ) = (φ2 (t ), ϕ2 (t ), ψ2 (t )) as follows

φ1 (t ) = ϕ1 (t ) = ψ1 (t ) =

k1 eη1 t , k1 + ε1 e−ηt ,

t ≤ t1 , t > t1 ,

k2 eη2 t , k2 + ε2 e−ηt ,

t ≤ t2 , t > t2 ,

ϕ2 (t ) =

k3 eη3 t , k3 + ε3 e−ηt ,

t ≤ t3 , t > t3 ,

ψ2 (t ) =









0, k1 − ε4 e−ηt ,

t ≤ t4 , t > t4 ,

 0,

t ≤ t5 , t5 < t ≤ t6 , t > t6 ,

φ2 (t ) =

k2 eη4 t − M0 , k − ε e−ηt , 2 5



0, k3 − ε6 e−ηt ,

t ≤ t7 , t > t7 ,

where M0 = k2 eη4 t6 − k2 + ε5 e−ηt6 , η > 0 is a constant to be chosen later. It is easy to know that M1 = supt ∈R φ1 > k1 , M2 = supt ∈R ϕ1 > k2 , M3 = supt ∈R ψ1 > k3 , Φ (t ) = (φ1 (t ), ϕ1 (t ), ψ1 (t )) and Ψ (t ) = (φ2 (t ), ϕ2 (t ), ψ2 (t )) satisfy (5.4), (P1), (P2) and (P3). Lemma 5.2. Φ (t ) = (φ1 (t ), ϕ1 (t ), ψ1 (t )) is an upper solution of system (5.3). Proof. If t ≤ t1 , φ1 (t ) = k1 eη1 t and ϕ1 (t ) = k2 eη2 t . It follows that Dφ1′′ (t ) − c φ1′ (t ) + αϕ1 (t ) − dφ1 (t ) ≤ (Dη12 − c η1 − d + α M2 /M1 )k1 eη1 t = 0. If t > t1 , φ1 (t ) = k1 + ε1 e−ηt . We have Dφ1′′ (t ) − c φ1′ (t ) + αϕ1 (t ) − dφ1 (t ) ≤ I1 (η), where I1 (η) = (Dε1 η2 + c ε1 η)e−ηt − d(k1 + ε1 e−ηt ) + α M2 . It follows from (5.4) that I1 (0) < 0 and there exists η1∗ > 0 such that Dφ1′′ (t ) − c φ1′ (t ) + αϕ1 (t ) − dφ1 (t ) < 0 for all η ∈ (0, η1∗ ). If t ≤ t2 , ϕ1 (t ) = k2 eη2 t . It follows that Dϕ1′′ (t ) − c ϕ1′ (t ) + β(A/d − φ2 (t ) − ϕ1 (t ) − ψ2 (t ))ϕ1 (t − c τ ) − (γ + α + d)ϕ1 (t )

≤ (Dη22 − c η2 + β A/d − (γ + α + d))k2 eη2 t = 0. If t > t2 , ϕ1 (t ) = k2 + ε2 e−ηt . We get Dϕ1′′ (t ) − c ϕ1′ (t ) + β(A/d − φ2 (t ) − ϕ1 (t ) − ψ2 (t ))ϕ1 (t − c τ ) − (γ + α + d)ϕ1 (t ) ≤ I2 (η), where I2 (η) = (Dε2 η2 + c ε2 η)e−ηt + β(A/d − k1 + ε4 e−ηt − k2 − ε2 e−ηt − k3 + ε6 e−ηt )M2

− (γ + α + d)(k2 + ε2 e−ηt ).

64

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

It follows from (5.4) that I2 (0) < 0 and there exists η2∗ > 0 such that Dϕ1′′ (t )− c ϕ1′ (t )+β(A/d −φ2 (t )−ϕ1 (t )−ψ2 (t ))ϕ1 (t − c τ ) − (γ + α + d)ϕ1 (t ) < 0 for all η ∈ (0, η2∗ ). If t ≤ t3 , ϕ1 (t ) = k2 eη2 t and ψ1 (t ) = k3 eη3 t . It follows that Dψ1′′ (t ) − c ψ1′ (t ) + γ ϕ1 (t ) − (δ + d)ψ1 (t ) ≤ (Dη32 − c η3 + M2 γ /M3 − (δ + d))k3 eη3 t = 0. If t > t3 , ψ1 (t ) = k3 + ε3 e−ηt . We get Dψ1′′ (t ) − c ψ1′ (t ) + γ ϕ1 (t ) − (δ + d)ψ1 (t ) ≤ I3 (η), where I3 (η) = (Dε3 η2 + c ε3 η)e−ηt + γ M2 − (δ + d)(k3 + ε3 e−ηt ). It follows from (5.4) that I3 (0) < 0 and there exists η3∗ > 0 such that Dψ1′′ (t ) − c ψ1′ (t ) + γ ϕ1 (t ) − (δ + d)ψ1 (t ) < 0 for all η ∈ (0, η3∗ ). Taking η ∈ (0, min(η1∗ , η2∗ , η3∗ )), we see that the conclusion is true. This completes the proof.  Lemma 5.3. Ψ (t ) = (φ2 (t ), ϕ2 (t ), ψ2 (t )) is a lower solution of system (5.3). Proof. If t ≤ t4 , φ2 (t ) = 0. It follows that Dφ2′′ (t ) − c φ2′ (t ) + αϕ2 (t ) − dφ2 (t ) = αϕ2 (t ) ≥ 0. If t > t4 , φ2 (t ) = k1 − ε4 e−ηt and ϕ2 (t ) = k2 − ε5 e−ηt . We have Dφ2′′ (t ) − c φ2′ (t ) + αϕ2 (t ) − dφ2 (t ) = I4 (η), where I4 (η) = −(Dε4 η2 + c ε4 η)e−ηt − d(k1 − ε4 e−ηt ) + α(k2 − ε5 e−ηt ). It follows from (5.4) that I4 (0) > 0 and there exists η4∗ > 0 such that Dφ2′′ (t ) − c φ2′ (t ) + αϕ2 (t ) − dφ2 (t ) > 0 for all η ∈ (0, η4∗ ). If t ≤ t5 , ϕ2 (t ) = 0. It follows that Dϕ2′′ (t ) − c ϕ2′ (t ) + β(A/d − φ1 (t ) − ϕ2 (t ) − ψ1 (t ))ϕ2 (t − c τ ) − (γ + α + d)ϕ2 (t ) = 0. If t5 < t ≤ t6 , ϕ2 (t ) = k2 eη4 t − M0 . We get Dϕ2′′ (t ) − c ϕ2′ (t ) + β(A/d − φ1 (t ) − ϕ2 (t ) − ψ1 (t ))ϕ2 (t − c τ ) − (γ + α + d)ϕ2 (t )

≥ (Dη42 − c η4 − (γ + α + d))k2 eη4 t = 0. If t > t6 , ϕ2 (t ) = k2 − ε5 e−ηt . We know Dϕ2′′ (t ) − c ϕ2′ (t ) + β(A/d − φ1 (t ) − ϕ2 (t ) − ψ1 (t ))ϕ2 (t − c τ ) − (γ + α + d)ϕ2 (t ) ≥ I5 (η), where I5 (η) = −(Dε5 η2 + c ε5 η)e−ηt + β(A/d − M1 − k2 + ε5 e−ηt − M3 )(k2 − ε5 e−η(t −c τ ) ) − (γ + α + d)(k2 − ε5 e−ηt ). It follows from (5.4) that I5 (0) > 0 and there exists η5∗ > 0 such that Dϕ2′′ (t )− c ϕ2′ (t )+β(A/d −φ1 (t )−ϕ2 (t )−ψ1 (t ))ϕ2 (t − c τ ) − (γ + α + d)ϕ2 (t ) > 0 for all η ∈ (0, η5∗ ). If t ≤ t7 , ψ2 (t ) = 0. It follows that Dψ2′′ (t ) − c ψ2′ (t ) + γ ϕ2 (t ) − (δ + d)ψ2 (t ) = γ ϕ2 (t ) ≥ 0. If t > t7 , ψ2 (t ) = k3 − ε6 e−ηt . We get Dψ2′′ (t ) − c ψ2′ (t ) + γ ϕ2 (t ) − (δ + d)ψ2 (t ) ≥ I6 (η), where I6 (η) = −(Dε3 η2 + c ε3 η)e−ηt + γ (k2 − ε5 e−ηt ) − (δ + d)(k3 − ε6 e−ηt ). It follows from (5.4) that I6 (0) > 0 and there exists η6∗ > 0 such that Dψ2′′ (t ) − c ψ2′ (t ) + γ ϕ2 (t ) − (δ + d)ψ2 (t ) > 0 for all η ∈ (0, η6∗ ). Letting η ∈ (0, min(η4∗ , η5∗ , η6∗ )), we see that our conclusion is true. This completes the proof.  Applying Lemmas 5.1–5.3, we know that if R0 > 1, system (5.2) has a travelling wave solution with speed c > c ∗ connecting the steady states (0, 0, 0) and (k1 , k2 , k3 ). Accordingly, we have the following conclusion.

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

65

Fig. 1. The temporal solution found by numerical integration of problem (1.2) under conditions (1.3) and (3.1) with D = 1, A = 0.5, β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 1 and ρ1 (x, t ) = 0.1, ρ1 (x, t ) = 0.001, ρ1 (x, t ) = 0.0001, x ∈ [0, 1], t ∈ [−1, 0].

Theorem 5.1. Let R0 > 1. For every c > c ∗ , regardless of the value of τ ≥ 0, system (1.2) always has a travelling wave solution with speed c connecting the uninfected steady state E1 (A/d, 0, 0) and the infected steady state E ∗ (S ∗ , I ∗ , R∗ ). Remark. The travelling wave solution in Theorem 5.1 may not be monotonic. The fact is illustrated by the following numerical simulations. 6. Numerical simulations In this section, by using the classical implicit format solving the partial differential equations and the method of steps for differential difference equations, we give the numerical simulations to illustrate the theoretical results above. Example 1. In system (1.2), let D = 1, A = 0.3, β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 1. System (1.2) with the above coefficients has a unique disease-free steady state E1 (1.5, 0, 0). It is easy to show that the basic reproduction number of system (1.2) is R0 = 0.8333 < 1. By Theorem 3.1 we see that E1 is locally stable. A numerical simulation illustrates our result (see Fig. 1). Example 2. In system (1.2), we set D = 1, A = 0.5, β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 1. System (1.2) with the above coefficients has a disease-free steady state E1 (2.5, 0, 0) and an endemic steady state E ∗ (1.8, 0.3333, 0.2). It is easy to show that the basic reproduction number of system (1.2) is R0 = 1.3889 > 1. Theorem 3.1 shows that E1 is unstable and E ∗ is locally stable. A numerical simulation illustrates our result (see Fig. 2). It follows from Theorem 5.1 that system (1.2) always has a travelling wave solution with speed c > c ∗ connecting E1 and E ∗ . The fact is illustrated by the numerical simulation in Fig. 3. For simplicity, we show only the travelling waves of I (x, t ) throughout this paper. Example 3. In system (1.2), we choose D = 1, A = 1, β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 1. System (1.2) with above coefficients has a disease-free steady state E1 (5, 0, 0) and an endemic steady state E ∗ (1.8, 1.5238, 0.9143). It is easy to show that the basic reproduction number of system (1.2) is R0 = 2.7778 > 1. By Theorem 3.1 we see that E1

66

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

Fig. 2. The temporal solution found by numerical integration of problem (1.2) under conditions (1.3) and (3.1) with D = 1, A = 0.5, β = 0.5, d = 0.2,

γ = 0.6, δ = 0.8, α = 0.1 and τ = 1 and ρ1 (x, t ) = 0.1, ρ1 (x, t ) = 0.001, ρ1 (x, t ) = 0.0001, x ∈ [0, 1], t ∈ [−1, 0].

Fig. 3. Solution of the travelling waves of system (1.2) starting from initial conditions as S (0) = 0.1, I (0) = 0.001, R(0) = 0.0001 with D = 1, A = 0.5,

β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 1.

is unstable and E ∗ is locally stable. It follows from Theorem 5.1 that system (1.2) always has a travelling wave solution with speed c > c ∗ connecting E1 and E ∗ . The fact is illustrated by the numerical simulation in Fig. 4. Clearly, Fig. 4 shows that the travelling wave solution does not possess monotonicity, and this seems to be due to the feature of the prey–predator interaction.

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

67

Fig. 4. Solution of the travelling waves of system (1.2) starting from initial conditions as S (0) = 0.1, I (0) = 0.001, R(0) = 0.0001 with D = 1, A = 1,

β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 1.

Fig. 5. Solution of the travelling waves of system (1.2) starting from initial conditions as S (0) = 0.1, I (0) = 0.001, R(0) = 0.0001 with D = 1, A = 0.5,

β = 0.5, d = 0.2, γ = 0.6, δ = 0.8, α = 0.1 and τ = 0.01.

7. Discussion In this paper, we have dealt with the existence of travelling wave solutions for an SIRS epidemic model with a bilinear incidence rate, spatial diffusion and time delay. At first, by analysing the corresponding characteristic equations, we discussed the local stability of a disease-free steady state and an endemic steady state to system (1.2) under homogeneous Neumann boundary conditions. We have shown in Theorem 3.1 that time delay and spatial diffusion are negligible for the local stability of the steady states to system (1.2). By using the cross iteration method and Schauder’s fixed point theorem, we reduced the existence of travelling wave solutions to the existence of a pair of upper–lower solutions. By constructing a pair of upper–lower solutions, we derived the existence of a travelling wave solution connecting the uninfected steady state E1 and the endemic steady state E ∗ . In fact, we also find that time delay τ can influence the monotone of the travelling wave solution connecting the diseasefree steady state E1 and the endemic steady state E ∗ . The effect of time delay on the travelling wave solution is illustrated by comparing the numerical simulations in Figs. 3 and 5. References [1] M. Kermack, A. Mckendrick, Contributions to the mathematical theory of epidemics, Proc. R. Soc. A 115 (1927) 700–721. [2] J. Mena-Lorca, H.W. Hethcote, Dynamic models of infectious diseases as regulators of population size, J. Math. Biol. 30 (1992) 693–716. [3] C. Castillo-Chavez, H.E. Hethcote, V. Andreasen, S.A. Levin, W.M. Liu, Epidemiological models with age structure, proportionate mixing, and crossimmunity, J. Math. Biol. 27 (1989) 233–258. [4] H.W. Hethcote, J.A. Yorke, Gonorrhea Transmission Dynamics and Control, in: Lecture Notes in Biomath., vol. 56, Springer, Berlin, 1984.

68

Q. Gan et al. / Nonlinear Analysis: Real World Applications 12 (2011) 52–68

[5] E. Beretta, Y. Takeuchi, Global stability of an SIR epidemic model with time delays, J. Math. Biol. 33 (1995) 250–260. [6] N.Y. Kyrychko, B.K. Blyuss, Global properties of a delayed SIR model with temporary immunity and nonlinear incidence rate, Nonlinear Anal. RWA 6 (2005) 495–507. [7] W. Ma, M. Song, Y. Takeuchi, Global stability of an SIR epidemic model with time delay, Appl. Math. Lett. 17 (2004) 1141–1145. [8] H.L. Smith, M. Williams, Stability in a cyclic epidemic models with delay, J. Math. Biol. 11 (1981) 95–103. [9] R. Xu, Z. Ma, Global stability of a SIR epidemic model with nonlinear incidence rate and time delay, Nonlinear Anal. RWA 10 (2009) 3175–3189. [10] R. Xu, Z. Ma, Stability of a delayed SIRS epidemic model with a nonlinear incidence rate, Chaos Solitons Fractals 41 (2009) 2319–2325. [11] K.L. Cooke, Stability analysis for a vector disease model, Rocky Mountain J. Math. 9 (1979) 31–42. [12] Z. Jin, Z. Ma, M. Han, Global stability of an SIRS epidemic model with delays, Acta Math. Sci. 26 (2006) 291–306. [13] T. Zhang, Z. Teng, Global behavior and permanence of SIRS epidemic model with time delay, Nonlinear Anal. RWA 9 (2008) 1409–1424. [14] T. Zhang, Z. Teng, Permanence and extinction for a nonautonomous SIRS epidemic model with time delay, Appl. Math. Model. 33 (2009) 1058–1072. [15] T. Zhang, J. Liu, Z. Teng, Stability of Hopf bifurcation of a delayed SIRS epidemic model with stage structure, Nonlinear Anal. RWA 11 (2010) 293–306. [16] J.V. Noble, Geographic and temporal development of plagues, Nature 250 (1974) 276–279. [17] G. Saccomandi, The spatial diffusion of diseases, Math. Comput. Modelling 25 (1997) 83–95. [18] M.A. Fuentes, M.N. Kuperman, Cellular automata and epidemiological models with spatial dependence, Physica A 267 (1999) 471–486. [19] K. Kim, Z. Lin, L. Zhang, Avian-human influenza epidemic model with diffusion, Nonlinear Anal. RWA 11 (2010) 313–322. [20] M.N. Kuperman, H.S. Wio, Front propagation in epidemiological models with spatial dependence, Physica A 272 (1999) 206–222. [21] N.A. Maidana, H. Yang, Spatial spreading of West Nile Virus described by traveling waves, J. Theoret. Biol. 258 (2009) 403–417. [22] R. Peng, S. Liu, Global stability of the steady states of an SIS epidemic reaction–diffusion model, Nonlinear Anal. 71 (2009) 239–247. [23] C. Atkinson, G.E. Reuter, Deterministic epidemic waves, Math. Proc. Cambridge Philos. Soc. 80 (1976) 315–330. [24] J. Wu, X. Zou, Travelling wave fronts of reaction diffusion systems with delay, J. Dynam. Differential Equations 13 (2001) 651–687. [25] X. Zou, J. Wu, Local existence and stability of periodic traveling wave of lattice functional differential equations, Can. Appl. Math. Q. 6 (1998) 397–418. [26] A.W. Leung, Systems of Nonlinear Partial Differential Equations with Applications to Biology and Engineering, Kluwer Academic Publishers, Dordrecht, 1989. [27] J.D. Murray, Mathematical Biology, Springer, New York, 1989. [28] J. Huang, X. Zou, Travelling wave solutions in delayed reaction diffusion systems with partial monotonicity, Acta Math. Appl. Sin. Engl. Ser. 22 (2) (2006) 243–256. [29] W. Li, G. Lin, S. Ruan, Existence of travelling wave solutions in delayed reaction–diffusion systems with applications to diffusion-competition systems, Nonlinearity 19 (2006) 1253–1273. [30] M. Wang, Stationary patterns for a prey–predator model with prey-dependent and ratio-dependent functional responses and diffusion, Physica D 196 (2004) 172–192. [31] R. Peng, M. Wang, Global stability of the equilibrium of a diffusive Holling-Tanner prey–predator model, Appl. Math. Lett. 20 (2007) 664–670. [32] M. Wang, Stability and Hopf bifurcation for a prey–predator model with prey-stage structure and diffusion, Math. Biosci. 212 (2008) 149–160. [33] E. Zeidler, Nonlinear Functional Analysis and its Applications, I, Fixed-point Theorems, Springer-Verlag, New York, 1986.