Travelling Waves of a Diffusive Epidemic Model with Latency and Relapse

Abstract

An SEIR epidemic model with relapse and spatial diffusion is studied. By analyzing the corresponding characteristic equations, the local stability of each of the feasible steady states to this model is discussed. The existence of a travelling wave solution is established by using the technique of upper and lower solutions and Schauder's fixed point theorem. Numerical simulations are carried out to illustrate the main results.

1. Introduction

In most epidemic models, individuals are often divided into several classes such as susceptible, infective, and recovered classes. Some diseases such as human or bovine tuberculosis have a latency (see, e.g., [1, 2] and references therein): infected individuals do not infect other susceptible individuals until sometime later. Hence, it is reasonable to include a latent class in such disease models. Also, for some diseases, recovered individuals may relapse with reactivation of latent infection and revert back to the infective class. This recurrence of disease is an important feature of some animal and human diseases. For human tuberculosis, incomplete treatment can lead to relapse, but relapse can also occur in patients who took a full course of treatment and were declared cured [3–5]. Most tuberculosis in human adults (caused by Mycobacterium tuberculosis) in the USA results from reactivation of latent infection [3].

In order for a model to be more realistic, taking into consideration a general exposed distribution and the possibility of relapse, van den Driessche et al. [6] formulated and studied the following model: d𝑆(𝑡)(𝑡)𝐼(𝑡)d𝑡=𝜇𝑁−𝛽𝑆𝑁(𝑡),𝐸(𝑡)=−𝜇𝑆_𝑡0(𝜉)𝐼(𝜉)𝛽𝑆𝑁𝑒^{−𝜇(𝑡−𝜉)}𝑃()(𝑡)𝑡−𝜉d𝜉,d𝑅(𝑡)−𝑅(𝑡),𝐼(𝑡)(𝑡)(𝑡)(𝑡),d𝑡=𝛾𝐼𝜇+𝛿=𝑁−𝑆−𝐸−𝑅 (1) where 𝑁 is the size of the population; 𝑆,𝐸,𝐼, and 𝑅 are the population sizes of susceptible, latent, infective, and recovered classes, respectively; 𝛿>0 is a constant rate at which an individual in the recovered class reverts to the infective class; 𝜇>0 is the recruitment rate and the death rate (this guarantees that 𝑁 can be assumed as a constant); 𝛽>0 is the transmission coefficient; 𝛾⩾0 is the recovery rate; 𝑃(𝑡) denotes the probability that a latent individual still remains in the latent class 𝑡 time units after entering the latent class.

As in [6], when 𝑃(𝑡) is a step function, that is, 𝑃(𝑡)=1,0<𝑡<𝜏,0,𝑡>𝜏, (2) by rescaling 𝑆(𝑡)𝑁(𝑡),𝐸(𝑡)⟶𝑆𝑁(𝑡),𝐼(𝑡)⟶𝐸𝑁(𝑡),𝑅(𝑡)⟶𝐼𝑁⟶𝑅(𝑡), (3) model (1) reduces to the following delay differential equation system: 𝑆^′(𝑡)(𝑡)(𝑡)𝐼(𝑡),𝐼=𝜇−𝜇𝑆−𝛽𝑆_′(𝑡)=𝛽𝑒^−𝜇𝜏𝑆()𝐼()(𝑡)−𝐼(𝑡),𝑅𝑡−𝜏𝑡−𝜏+𝛿𝑅𝜇+𝛾_′(𝑡)=𝛾𝐼(𝑡)−𝜇+𝛿𝑅(𝑡), (4) where 𝜏>0 can be explained as the average latent period of the disease. The global stability of the steady states for system (1) and (4) was analyzed in [6, 7].

Note that it is implicitly assumed that the population is well mixed, and the spatial mobility of individuals has been ignored in model (1) as well as most other epidemic models. In reality, the environment in which an individual lives is often heterogeneous making it necessary to distinguish the locations, and due to the large mobility of individuals within an area or even worldwide, spatially uniform models are not sufficient to give a realistic picture of a disease's transmission. For this reason, the spatial effects cannot be neglected in studying the spread of epidemics. If the environment is spatially continuous, random diffusion is often used to describe the mobility of the population, leading to models in the form of reaction-diffusion equations (see, e.g., [8–10]). Noble [9] applied reaction-diffusion theory to describe the spread of plague through Europe in the mid-14th century. Maidana and Yang [10] proposed a spatial model to analyze the West Nile Virus propagation across the USA, from east to west.

Motivated by the work mentioned above, the main purpose of this paper is to study the effect of the spatial diffusion of the individuals on the dynamics of an epidemic model with latency and relapse. If we further consider the spatial diffusion, then the epidemic model corresponding to (4) should be described by the following delayed reaction diffusion system: 𝜕𝑆(𝑥,𝑡)𝜕𝑡=𝑑₁𝜕²𝑆(𝑥,𝑡)𝜕𝑥²()()𝐼(),()+𝜇−𝜇𝑆𝑥,𝑡−𝛽𝑆𝑥,𝑡𝑥,𝑡𝜕𝐼𝑥,𝑡𝜕𝑡=𝑑₂𝜕²𝐼(𝑥,𝑡)𝜕𝑥²+𝛽𝑒^−𝜇𝜏𝑆()𝐼()()−𝐼(),()𝑥,𝑡−𝜏𝑥,𝑡−𝜏+𝛿𝑅𝑥,𝑡𝜇+𝛾𝑥,𝑡𝜕𝑅𝑥,𝑡𝜕𝑡=𝑑₃𝜕²𝑅(𝑥,𝑡)𝜕𝑥²+𝛾𝐼(𝑥,𝑡)−𝜇+𝛿𝑅(𝑥,𝑡), (5) with initial conditions 𝑆(𝑥,𝑡)=𝜌₁(𝑥,𝑡),𝐼(𝑥,𝑡)=𝜌₂(𝑥,𝑡),𝑅(𝑥,𝑡)=𝜌₃(),[]𝑥,𝑡𝑡∈−𝜏,0,𝑥∈Ω. (6) In problem (5) and (6), 𝑑₁,𝑑₂, and 𝑑₃ are the spatial diffusion rates of susceptible, infective, and recovered individuals, respectively. The functions 𝜌_𝑖(𝑥,𝑡)∈𝒞=𝐶([−𝜏,0],𝑋)(𝑖=1,2,3) and 𝑋 is defined by 𝑋=𝑆,𝐼,𝑅∈𝑊^2,2(Ω)∶𝜕𝑆(𝑥,𝑡)=()𝜕𝑛𝜕𝐼𝑥,𝑡=()𝜕𝑛𝜕𝑅𝑥,𝑡,𝜕𝑛=0,𝑥∈𝜕Ω (7) with the inner product ⟨⋅,⋅⟩, where 𝜕/𝜕𝑛 denotes the outward normal derivative on 𝜕Ω.

The organization of this paper is as follows. In the next section, by analyzing the corresponding characteristic equations, the local stability of an endemic steady state and a disease-free steady state of system (5) is discussed. In Section 3, by constructing a pair of upper-lower solutions, we use the cross-iteration method and the Schauder's fixed point theorem to prove the existence of travelling wave solutions to system (5). In Section 4, numerical simulations are carried out to illustrate the main theoretical results, and by combining with numeric simulations, we discuss the effects of some parameters on the dynamics of system (5). A brief conclusion is given in the last section.

2. Local Stability

In this section, by analyzing the corresponding characteristic equations, we investigate the local stability of steady states to system (5) with the initial conditions (6) and the homogeneous Neumann boundary conditions 𝜕𝑆(𝑥,𝑡)=()𝜕𝑛𝜕𝐼𝑥,𝑡=()𝜕𝑛𝜕𝑅𝑥,𝑡𝜕𝑛=0,𝑡⩾0,𝑥∈𝜕Ω, (8) respectively. The homogeneous Neumann boundary conditions imply that the populations do not move across the boundary 𝜕Ω.

Let ℛ₀=𝛽𝑒^−𝜇𝜏𝜇+𝛿𝜇.𝜇+𝛿+𝛾 (9)ℛ₀ is called the basic reproductive number of system (5), which describes the average number of newly infected individuals at the beginning of the infectious process. It is easy to show that system (5) always has a disease-free uniform steady state 𝐸₀(1,0,0). Further, if ℛ₀>1, system (5) has a unique endemic steady state 𝐸^∗(𝑆^∗,𝐼^∗,𝑅^∗), where 𝑆^∗=𝜇𝜇+𝛿+𝛾𝑒^𝜇𝜏𝛽,𝐼𝜇+𝛿_∗=𝛽𝑒^−𝜇𝜏𝜇+𝛿−𝜇𝜇+𝛿+𝛾𝛽,𝑅𝜇+𝛿+𝛾_∗=𝛾𝐼𝜇+𝛿^∗. (10)

Let 𝐸(𝑆⁰,𝐼⁰,𝑅⁰) represent any feasible uniform steady state of system (5). If we set 𝑆=𝑆−𝑆⁰, 𝐼=𝐼−𝐼⁰, 𝑅=𝑅−𝑅⁰, and drop the bars for simplicity of notations, then (5) can be transformed into the following equivalent system: 𝜕𝑆(𝑥,𝑡)𝜕𝑡=𝑑₁𝜕²𝑆(𝑥,𝑡)𝜕𝑥²−𝜇𝑆(𝑥,𝑡)−𝛽𝑆(𝑥,𝑡)𝐼(𝑥,𝑡)+𝐼⁰−𝛽𝑆⁰𝐼(),()𝑥,𝑡𝜕𝐼𝑥,𝑡𝜕𝑡=𝑑₂𝜕²𝐼(𝑥,𝑡)𝜕𝑥²+𝛽𝑒^−𝜇𝜏𝑆(𝑥,𝑡−𝜏)𝐼(𝑥,𝑡−𝜏)+𝐼⁰+𝑆⁰𝐼(()−𝐼(),()𝑥,𝑡−𝜏)+𝛿𝑅𝑥,𝑡𝜇+𝛾𝑥,𝑡𝜕𝑅𝑥,𝑡𝜕𝑡=𝑑₃𝜕²𝑅(𝑥,𝑡)𝜕𝑥²+𝛾𝐼(𝑥,𝑡)−𝜇+𝛿𝑅(𝑥,𝑡). (11) Denote 𝑈=(𝑆(⋅,𝑡),𝐼(⋅,𝑡),𝑅(⋅,𝑡))^𝑇. Then the linearization of (11) can be rewritten as an abstract differential equation in the phase space 𝒞=𝐶([−𝜏,0],𝑋), 𝑈̇(𝑡)=𝐷Δ𝑈(𝑡)+𝐿𝑈_𝑡, (12) where 𝐷=diag(𝑑₁,𝑑₂,𝑑₃), and 𝐿 is defined by 𝐿𝜙=−𝜇−𝛽𝐼⁰−𝛽𝑆⁰0𝛿𝜙0−𝜇+𝛾0𝛾−𝜇+𝛿¹(0)𝜙₂(0)𝜙₃(0)+00𝛽𝑒^−𝜇𝜏𝐼⁰𝛽𝑒^−𝜇𝜏𝑆⁰0𝜙000¹()𝜙−𝜏₂()𝜙−𝜏₃(−𝜏) (13) for 𝜙=(𝜙₁,𝜙₂,𝜙₃)^𝑇∈𝒞, and its characteristic equation is 𝜆𝑦−𝐷Δ𝑦−𝐿𝑒^𝜆𝑦=0, (14) where 𝑦∈dom(Δ) and 𝑦≠0, dom(Δ)⊂𝑋.

From the properties of the Laplacian operator defined on the bounded domain, the operator Δ on 𝑋 has the eigenvalues −𝑘² with the relative eigenfunctions 𝛽¹_𝑘=𝛾_𝑘,0,0^𝑇,𝛽²_𝑘=0,𝛾_𝑘,0^𝑇,𝛽³_𝑘=0,0,𝛾_𝑘^𝑇, (15) where 𝛾_𝑘=cos(𝑘𝑥), 𝑘∈𝑁₀={0,1,2,…}. Clearly, (𝛽¹_𝑘,𝛽²_𝑘,𝛽³_𝑘)^∞_𝑘=0 construct a basis of the phase space 𝑋, and therefore, any element 𝑦 in 𝑋 can be expanded as Fourier series in the following form: 𝑦=^∞_𝑘=0𝑌^𝑇_𝐾𝛽¹_𝑘,𝛽²_𝑘,𝛽³_𝑘^𝑇,𝑌^𝑇_𝐾=𝑦,𝛽¹_𝑘,𝑦,𝛽²_𝑘,𝑦,𝛽³_𝑘. (16) Some simple computations show that 𝐿𝜙^𝑇𝛽¹_𝑘,𝛽²_𝑘,𝛽³_𝑘^𝑇=𝐿𝜙^𝑇𝛽¹_𝑘,𝛽²_𝑘,𝛽³_𝑘^𝑇,𝑘∈𝑁₀. (17) From (16) and (17), (14) is equivalent to ^∞_𝑘=0𝑌^𝑇_𝐾𝜆𝐼₂+𝐷𝑘²−𝑀𝛽¹_𝑘,𝛽²_𝑘,𝛽³_𝑘^𝑇=0, (18) where 𝑀=−𝜇−𝛽𝐼⁰−𝛽𝑆⁰0𝛽𝑒^{−(𝜇+𝜆)𝜏}𝐼⁰𝛽𝑒^{−(𝜇+𝜆)𝜏}𝑆⁰−𝜇+𝛾𝛿0𝛾−𝜇+𝛿. (19) If 𝐸⁰=𝐸₀, (18) is equivalent to 𝜆+𝑘²𝑑₁×+𝜇𝜆+𝑘²𝑑₂+𝜇+𝛾𝜆+𝑘²𝑑₃−+𝜇+𝛿−𝛿𝛾𝜆+𝑘²𝑑₃+𝜇+𝛿𝛽𝑒^−𝜇𝜏𝑒^−𝜆𝜏=0. (20) Clearly, for any 𝑖⩾1, (20) has a negative root −𝑘²𝑑₁−𝜇. Its other roots satisfy 𝑔₁(𝜆)≜𝜆+𝑘²𝑑₂+𝜇+𝛾𝜆+𝑘²𝑑₃+𝜇+𝛿−𝛿𝛾−𝜆+𝑘²𝑑₃+𝜇+𝛿𝛽𝑒^−𝜇𝜏𝑒^−𝜆𝜏=0. (21)

If ℛ₀>1, then for 𝑖=1, we get that, for 𝜆 real, 𝑔₁(0)=𝜇𝜇+𝛿+𝛾−𝜇+𝛿𝛽𝑒^−𝜇𝜏<0,lim_𝜆→+∞𝑔₁(𝜆)=∞. (22) Hence, when 𝑖=1, (21) has a positive real root, and accordingly, 𝐸₀ is unstable.

If ℛ₀<1 and 𝜏=0, we have 𝑔₁(𝜆)=𝜆²+𝑘²𝑑₂+𝑘²𝑑₃𝜆+𝑘+2𝜇+𝛾+𝛿−𝛽₂𝑑₂+𝜇+𝛾−𝛽𝑘²𝑑₃+𝜇+𝛿−𝛿𝛾>0, (23) which means that the roots of (21) are negative constants, then 𝐸₀ is locally asymptotically stable. For 𝜏>0, if 𝑖𝜎(𝜎>0) is a solution of (21), separating the real and imaginary parts, we derive that −𝜎²+𝐴₀=𝐴₁cos𝜎𝜏+𝐴₂𝐴𝜎sin𝜎𝜏,₃𝜎=𝐴₂𝜎cos𝜎𝜏−𝐴₁sin𝜎𝜏, (24) where 𝐴₀=𝑘²𝑑₂+𝜇+𝛾𝑘²𝑑₃𝐴+𝜇+𝛿−𝛿𝛾,₁=𝑘²𝑑₃+𝜇+𝛿𝛽𝑒^−𝜇𝜏,𝐴₂=𝛽𝑒^−𝜇𝜏,𝐴₃=𝑘²𝑑₂+𝑑₃+2𝜇+𝛾+𝛿. (25) Squaring and adding the two equations of (24), it follows that 𝜎⁴+𝐴²₃−2𝐴₀−𝐴²₂𝜎²+𝐴²₀−𝐴²₁=0. (26)

Note that if 𝑅₀<1, we obtain 𝐴²₃−2𝐴₀−𝐴²₂=𝑘²𝑑₂+𝜇+𝛾²+𝑘²𝑑₃+𝜇+𝛾²+2𝛿𝛾−𝛽²𝑒^−2𝜇𝜏⩾𝜇+𝛾²+𝜇+𝛿²+2𝛿𝛾−𝛽²𝑒^−2𝜇𝜏⩾𝜇²𝐴>0,₀−𝐴₁⩾𝜇+𝛿𝜇+𝛾−𝛽𝑒^−𝜇𝜏𝑒−𝛿𝛾=𝜇𝜇+𝛿+𝛾−𝛽𝜇+𝛿^−𝜇𝜏>0. (27) Hence, (26) has no positive roots for all 𝑖⩾1. Therefore, if ℛ₀<1, 𝐸₀ is locally asymptotically stable for all 𝜏>0.

If 𝐸⁰=𝐸^∗, (18) is equivalent to 𝜆³+𝑝₂𝜆²+𝑝₁𝜆+𝑝₀+𝑞₂𝜆²+𝑞₁𝜆+𝑞₀𝑒^−𝜆𝜏=0, (28) where 𝑝₀=𝑘²𝑑₁+𝜇+𝛽𝐼^∗𝑘²𝑑₂+𝜇+𝛾𝑘²𝑑₃,𝑝+𝜇+𝛿−𝛿𝛾₁=𝑘²𝑑₂+𝜇+𝛾𝑘²𝑑₃+𝑘+𝜇+𝛿−𝛿𝛾₂𝑑₁+𝜇+𝛽𝐼^∗𝑘²𝑑₂+𝑘²𝑑₃,𝑝+2𝜇+𝛿+𝛾₂=𝑘²𝑑₁+𝜇+𝛽𝐼^∗+𝑘²𝑑₂+𝜇+𝛾+𝑘²𝑑₃,𝑞+𝜇+𝛿₀=−𝛽𝑆^∗𝑒^−𝜇𝜏𝑘²𝑑₁+𝜇𝑘²𝑑₃,𝑞+𝜇+𝛿₁=−𝛽𝑆^∗𝑒^−𝜇𝜏𝑘²𝑑₁+𝑘²𝑑₃,𝑞+2𝜇+𝛿₂=−𝛽𝑆^∗𝑒^−𝜇𝜏. (29) When 𝜏=0, (28) becomes 𝜆³+𝑝₂+𝑞₂𝜆²+𝑝₁+𝑞₁𝜆+𝑝₀+𝑞₀=0. (30) It is not difficult to verify that 𝑝₂+𝑞₂>𝜇+𝛾−𝛽𝑆^∗=𝛿𝛾𝑝𝜇+𝛿>0,₀+𝑞₀>𝑘²𝑑₁+𝜇𝜇+𝛾𝑘²𝑑₃+𝜇+𝛿−𝛿𝛾−𝛽𝑆^∗𝑘²𝑑₁+𝜇𝑘²𝑑₃,>𝑘+𝜇+𝛿𝜇+𝛿²𝑑₁+𝜇𝛿𝛾𝑘𝜇+𝛿−𝛿𝛾²𝑑₁𝑝+𝜇=0,₁+𝑞₁𝑝₂+𝑞₂−𝑝₀+𝑞₀>𝑘²𝑑₂+𝜇+𝛾𝑘²𝑑₃+𝑘+𝜇+𝛿−𝛿𝛾₂𝑑₁+𝜇+𝛽𝐼^∗𝑘²𝑑₂+𝑘²𝑑₃+2𝜇+𝛿+𝛾−𝛽𝑆^∗𝑘²𝑑₁+𝑘²𝑑₃×𝑘+2𝜇+𝛿₂𝑑₂+𝑘²𝑑₃+2𝜇+𝛿+𝛾−𝛽𝑆^∗>0, (31) which indicates that all roots 𝜆_𝑖,1,𝜆_𝑖,2, and 𝜆_𝑖,3 of (30) have negative real parts.

Next we can conclude that there exists a positive constant 𝜉 such that Re𝜆_𝑖,1,Re𝜆_𝑖,2,Re𝜆_𝑖,3⩽−𝜉,𝑖⩾1. (32)

In fact, let 𝜆=𝜇_𝑖𝑣, then 𝜆³+𝑝₂+𝑞₂𝜆²+𝑝₁+𝑞₁𝜆+𝑝₀+𝑞₀=𝜇³_𝑖𝑣³+𝑝₂+𝑞₂𝜇²_𝑖𝑣²+𝑝₁+𝑞₁𝜇_𝑖𝑣+𝑝₀+𝑞₀≜ℎ_𝑖(𝑣). (33) Since 𝜇_𝑖→∞ as 𝑖→∞, it yields that lim_𝑖→∞ℎ_𝑖(𝑣)𝜇³_𝑖=𝑣³+𝑑₁+𝑑₂+𝑑₃𝑣²+𝑑₁𝑑₂+𝑑₁𝑑₃+𝑑₂𝑑₃𝑣+𝑑₁𝑑₂𝑑₃≜ℎ(𝑣). (34) It is easy to see that these three roots 𝑣₁,𝑣₂, and 𝑣₃ of ℎ(𝑣) have negative real parts from the Routh-Hurwitz criterion again, and thus, there exists a positive constant 𝜉 such that Re𝑣₁,Re𝑣₂,Re𝑣₃⩽−𝜉. (35)

By continuity, we see that there exists 𝑖₀ such that the three roots 𝑣_𝑖,1,𝑣_𝑖,2, and 𝑣_𝑖,3 of ℎ_𝑖(𝑣) satisfy Re𝑣_𝑖,1,Re𝑣_𝑖,2,Re𝑣_𝑖,3⩽−𝜉2,𝑖⩾𝑖₀, (36) and thus, Re𝜆_𝑖,1,Re𝜆_𝑖,2,Re𝜆_𝑖,3𝜇⩽−^𝑖𝜉2𝜇⩽−^𝑖0𝜉2,𝑖⩾𝑖₀. (37) Let −𝜉^∼=max_{1⩽𝑖⩽𝑖0}Re𝜆_𝑖,1,Re𝜆_𝑖,2,Re𝜆_𝑖,3, (38) then 𝜉^∼, and (32) holds for 𝜉∈(𝜉^∼,𝜇_𝑖0𝜉/2).

Consequently, the spectrum of ℒ, consisting only of eigenvalues, lies in Re𝜆⩽−𝜉. Hence, the infected steady state is locally asymptotically stable when 𝜏=0.

For 𝜏>0, if 𝑖𝜔(𝜔>0) is a solution of (28), separating the real and imaginary parts, we derive that 𝜔³−𝑝₁𝜔=𝑞₁𝜔cos𝜔𝜏+𝑞₂𝜔²−𝑞₀𝑝sin𝜔𝜏,₂𝜔²−𝑝₀=𝑞₁𝜔sin𝜔𝜏−𝑞₂𝜔²−𝑞₀cos𝜔𝜏. (39) Squaring and adding the two equations of (39), it follows that 𝜔⁶+𝑝²₂−2𝑝₁−𝑞²₂𝜔⁴+𝑝²₁−2𝑝₀𝑝₂−𝑞²₁+2𝑞₀𝑞₂𝜔²+𝑝²₀−𝑞²₀=0. (40) By calculation, it is not difficult to verify that, for all 𝑖⩾1, 𝑝²₂−2𝑝₁−𝑞²₂=𝑘²𝑑₁+𝜇+𝛽𝐼^∗2+𝑘²𝑑₂+𝜇+𝛾²+𝑘²𝑑₃+𝜇+𝛿²+2𝛿𝛾−𝛽𝑆^∗𝑒^−𝜇𝜏2>𝜇+𝛾²−𝛽²𝑆^∗2𝑝>0,₂₁−2𝑝₀𝑝₂−𝑞²₁+2𝑞₀𝑞₂>𝑘²𝑑₂+𝜇+𝛾𝑘²𝑑₃+𝜇+𝛿−𝛿𝛾²−𝛽𝑆^∗𝑒^−𝜇𝜏2𝑘²𝑑₃+𝜇+𝛿²>𝜇²𝜇+𝛿+𝛾²−𝛽𝑆^∗𝑒^−𝜇𝜏𝜇+𝛿²𝑝=0,₂₀−𝑞²₀=𝑝₀−𝑞₀𝑝₀+𝑞₀>0. (41)

Hence, (40) has no positive roots. Therefore, if ℛ₀>1, 𝐸^∗ is locally asymptotically stable for all 𝜏>0.

From the discussion above, we can obtain the following results.

Theorem 1. For system (5) with initial conditions (6) and homogeneous Neumann conditions (8), if ℛ₀<1, 𝐸₀ is locally asymptotically stable; if ℛ₀>1, 𝐸₀ becomes unstable, and 𝐸^∗ is locally asymptotically stable.

3. Existence of Travelling Waves

In this section, we investigate the existence of travelling wave solutions to system (5). The technique of the proofs is to use the Schauder's fixed point theorem, the method of upper-lower solutions, and its associated cross iteration scheme.

Let 𝑆=1−𝑆, then system (5) is transformed into (omitting the hat on 𝑆 for simplicity) 𝜕𝑆(𝑥,𝑡)𝜕𝑡=𝑑₁𝜕²𝑆(𝑥,𝑡)𝜕𝑥²()((,()−𝜇𝑆𝑥,𝑡+𝛽𝐼𝑥,𝑡)(1−𝑆𝑥,𝑡))𝜕𝐼𝑥,𝑡𝜕𝑡=𝑑₂𝜕²𝐼(𝑥,𝑡)𝜕𝑥²+𝛽𝑒^−𝜇𝜏𝐼((()−𝐼(),()𝑥,𝑡−𝜏)(1−𝑆𝑥,𝑡−𝜏))+𝛿𝑅𝑥,𝑡𝜇+𝛾𝑥,𝑡𝜕𝑅𝑥,𝑡𝜕𝑡=𝑑₃𝜕²𝑅(𝑥,𝑡)𝜕𝑥²+𝛾𝐼(𝑥,𝑡)−𝜇+𝛿𝑅(𝑥,𝑡). (42)

A travelling wave solution of (42) is a solution (𝑆(𝑥,𝑡),𝐼(𝑥,𝑡),𝑅(𝑥,𝑡)) of the form 𝑆(𝑥,𝑡)=𝜙(𝑥+𝑐𝑡), 𝐼(𝑥,𝑡)=𝜑(𝑥+𝑐𝑡), 𝑅(𝑥,𝑡)=𝜓(𝑥+𝑐𝑡) where 𝜙,𝜑,𝜓∈𝐶²(ℝ,ℝ) and 𝑐>0 is a constant corresponding to the wave speed. On substituting 𝑆(𝑥,𝑡)=𝜙(𝑥+𝑐𝑡), 𝐼(𝑥,𝑡)=𝜑(𝑥+𝑐𝑡), and 𝑅(𝑥,𝑡)=𝜓(𝑥+𝑐𝑡) denoting the travelling wave coordinate 𝑥+𝑐𝑡 still by 𝑡, we derive from (42) that 𝑑₁𝜙^′′(𝑡)−𝑐𝜙^′(𝑡)+𝑓_𝑐1𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑=0,₂𝜑^′′(𝑡)−𝑐𝜑^′(𝑡)+𝑓_𝑐2𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑=0,₃𝜓^′′(𝑡)−𝑐𝜓^′(𝑡)+𝑓_𝑐3𝜙_𝑡,𝜑_𝑡,𝜓_𝑡=0, (43) where 𝑓_𝑐1𝜙_𝑡,𝜑_𝑡,𝜓_𝑡(𝑡)(𝑡)(𝑡),𝑓=−𝜇𝜙+𝛽𝜑1−𝜙_𝑐2𝜙_𝑡,𝜑_𝑡,𝜓_𝑡=𝛽𝑒^−𝜇𝜏𝜑()()(𝑡)−𝜑(𝑡),𝑓𝑡−𝑐𝜏1−𝜙𝑡−𝑐𝜏+𝛿𝜓𝜇+𝛾_𝑐3𝜙_𝑡,𝜑_𝑡,𝜓_𝑡=𝛾𝜑(𝑡)−𝜇+𝛿𝜓(𝑡). (44) Equation (43) will be solved subject to the following conditions: Lim_{𝑡→−∞}𝜙(𝑡),𝜑(𝑡),𝜓(𝑡)=(0,0,0),lim_𝑡→+∞𝜙(𝑡),𝜑(𝑡),𝜓(𝑡)=𝑘₁,𝑘₂,𝑘₃≜1−𝑆^∗,𝐼^∗,𝑅^∗. (45) We define the upper and lower solutions of system (43) as follows.

Definition 2. A pair of continuous functions Φ=(𝜙,𝜑,𝜓) and Φ=(𝜙,𝜑,𝜓) are called a pair of upper-lower solutions of system (43) if Φ and Φ are twice differential almost everywhere in ℝ and they are essentially bounded on ℝ, and there hold 𝑑₁𝜙^′′(𝑡)−𝑐𝜙^′(𝑡)+𝑓_𝑐1𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑⩽0,a.e.inℝ,₂𝜑^′′(𝑡)−𝑐𝜑^′(𝑡)+𝑓_𝑐2𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑⩽0,a.e.inℝ,₃𝜓^′′(𝑡)−𝑐𝜓^′(𝑡)+𝑓_𝑐3𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑⩽0,a.e.inℝ,₁𝜙^′′(𝑡)−𝑐𝜙^′(𝑡)+𝑓_𝑐1𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑⩾0,a.e.inℝ,₂𝜑^′′(𝑡)−𝑐𝜑^′(𝑡)+𝑓_𝑐2𝜙_𝑡,𝜑_𝑡,𝜓_𝑡𝑑⩾0,a.e.inℝ,₃𝜓^′′(𝑡)−𝑐𝜓^′(𝑡)+𝑓_𝑐3𝜙_𝑡,𝜑_𝑡,𝜓_𝑡⩾0,a.e.inℝ. (46)

As to 𝑓_𝑐1,𝑓_𝑐2, and 𝑓_𝑐3, we have the following.

Lemma 3. 𝑓_𝑐1,𝑓_𝑐2, and 𝑓_𝑐3 satisfy the following partial quasimonotonicity conditions.(C1) There exist three positive constants 𝛽₁,𝛽₂,𝛽₃ such that 𝑓_𝑐1𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐1𝜙₂,𝜑₂,𝜓₂+𝛽₁𝜙₁(0)−𝜙₂(0𝑓)⩾0,_𝑐2𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐2𝜙₁,𝜑₂,𝜓₂+𝛽₂𝜑₁(0)−𝜑₂(0𝑓)⩾0,_𝑐2𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐2𝜙₂,𝜑₁,𝜓₁𝑓⩽0,_𝑐3𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐3𝜙₂,𝜑₂,𝜓₂+𝛽₃𝜓₁(0)−𝜓₂(0)⩾0, (47) where 𝜙_𝑖,𝜑_𝑖,𝜓_𝑖∈𝐶([−𝜏,0],𝑅), 𝑖=1,2 with (0,0,0)⩽(𝜙₂,𝜑₂,𝜓₂)⩽ (𝜙₁,𝜑₁,𝜓₁)⩽ (𝑀₁, 𝑀₂, 𝑀₃), 𝑀_𝑗>𝑘_𝑗(𝑗=1,2,3) are positive constants.

Proof. It is not difficult to verify that 𝑓_𝑐1𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐1𝜙₂,𝜑₂,𝜓₂=𝛽1−𝜙₁(0)𝜑₁(0)−𝜇𝜙₁(0)+𝜇𝜙₂(0)+𝛽1−𝜙₂(0)𝜑₂(0)⩾−𝜇+𝛽𝜙₂(0)𝜙₁(0)−𝜙₂(0),𝑓_𝑐2𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐2𝜙₁,𝜑₂,𝜓₂=𝛽𝑒^−𝜇𝜏1−𝜙₁(−𝜏)𝜑₁(−𝜏)−𝜑₂()𝜓−𝜏+𝛿₁(0)−𝜓₂(0)−𝜇+𝛾𝜑₁(0)−𝜑₂(0𝜑)⩾−𝜇+𝛾₁(0)−𝜑₂(0),𝑓_𝑐2𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐2𝜙₂,𝜑₁,𝜓₁=𝛽𝑒^−𝜇𝜏𝜑₁(−𝜏)𝜙₂(−𝜏)−𝜙₁()𝑓−𝜏⩽0,_𝑐3𝜙₁,𝜑₁,𝜓₁−𝑓_𝑐3𝜙₂,𝜑₂,𝜓₂⩾−𝜇+𝛿𝜓₁(0)−𝜓₂(0)⩾0. (48) Let 𝛽₁=𝜇+𝛽𝜙₂(0),𝛽₂=𝜇+𝛾, and 𝛽₃=𝜇+𝛿, then the proof is completed.

For the constants 𝛽₁,𝛽₂,𝛽₃ above, define 𝐻₁𝜙,𝜑,𝜓(𝑡)=𝑓_𝑐1𝜙_𝑡,𝜑_𝑡,𝜓_𝑡+𝛽₁𝜙(𝑡),𝐻₂𝜙,𝜑,𝜓(𝑡)=𝑓_𝑐2𝜙_𝑡,𝜑_𝑡,𝜓_𝑡+𝛽₂𝜑(𝑡),𝐻₃𝜙,𝜑,𝜓(𝑡)=𝑓_𝑐3𝜙_𝑡,𝜑_𝑡,𝜓_𝑡+𝛽₃𝜓(𝑡),𝐹₁(𝑡)=1𝜙,𝜑,𝜓𝑑₁𝜆₂−𝜆₁×^𝑡_−∞𝑒^{𝜆1(𝑡−𝑠)}𝐻₁(𝑠)+𝜙,𝜑,𝜓d𝑠_∞𝑡𝑒^{𝜆2(𝑡−𝑠)}𝐻₁𝑠),𝐹𝜙,𝜑,𝜓(d𝑠₂𝑡)=1𝜙,𝜑,𝜓(𝑑₂𝜆₄−𝜆₃×^𝑡_−∞𝑒^{𝜆3(𝑡−𝑠)}𝐻₂(𝑠)+𝜙,𝜑,𝜓d𝑠_∞𝑡𝑒^{𝜆4(𝑡−𝑠)}𝐻₂(𝑠),𝐹𝜙,𝜑,𝜓d𝑠₃(𝑡)=1𝜙,𝜑,𝜓𝑑₃𝜆₆−𝜆₅×^𝑡_−∞𝑒^{𝜆5(𝑡−𝑠)}𝐻₃(𝑠)+𝜙,𝜑,𝜓d𝑠_∞𝑡𝑒^{𝜆6(𝑡−𝑠)}𝐻₃𝜙,𝜑,𝜓(𝑠)d𝑠, (49) where (0,0,0)⩽(𝜙,𝜑,𝜓)⩽(𝑘₁,𝑘₂,𝑘₃) and 𝜆₁=𝑐−𝑐²+4𝛽₁𝑑₁2𝑑₁,𝜆₂=𝑐+𝑐²+4𝛽₁𝑑₁2𝑑₁,𝜆₃=𝑐−𝑐²+4𝛽₂𝑑₂2𝑑₂,𝜆₄=𝑐+𝑐²+4𝛽₂𝑑₂2𝑑₂,𝜆₅=𝑐−𝑐²+4𝛽₃𝑑₃2𝑑₃,𝜆₆=𝑐+𝑐²+4𝛽₃𝑑₃2𝑑₃. (50) Define the following profile set: Γ𝜙,𝜑,𝜓,𝜙,𝜑,𝜓=𝜙,𝜑,𝜓∈𝐶ℝ,ℝ³∣𝜙,𝜑,𝜓⩽𝜙,𝜑,𝜓⩽𝜙,𝜑,𝜓. (51) From Lemmas 3.1–3.6 of [11], we have the following results.

Lemma 4. 𝐹(Γ((𝜙,𝜑,𝜓),(𝜙,𝜑,𝜓)))⊂Γ((𝜙,𝜑,𝜓),(𝜙,𝜑,𝜓)), where 𝐹=(𝐹₁,𝐹₂).

Lemma 5. 𝐹∶Γ((𝜙,𝜑,𝜓),(𝜙,𝜑,𝜓))→Γ((𝜙,𝜑,𝜓),(𝜙,𝜑,𝜓)) is compact.

From Lemmas 4 and 5, we see that the existence of travelling wave solutions for system (5) follows from the existence of a pair of upper and lower solutions (𝜙,𝜑,𝜓) and (𝜙,𝜑,𝜓) of (43) satisfying the following conditions.(C2)(0,0,0)⩽(𝜙(𝑡),𝜑(𝑡),𝜓(𝑡))⩽(𝜙(𝑡),𝜑(𝑡),𝜓(𝑡))⩽(𝑀₁,𝑀₂,𝑀₃),𝑡∈ℝ. (C3)lim_{𝑡→−∞}(𝜙(𝑡),𝜑(𝑡),𝜓(𝑡))=(0,0,0),lim_𝑡→+∞(𝜙(𝑡),𝜑(𝑡),𝜓(𝑡))=(𝑘₁,𝑘₂,𝑘₃).

Let 𝑑₁=𝑑₂=𝑑₃ and 𝑐>𝑐^∗≜max{2𝑑₁𝛽𝑘₂,2𝑑₂𝛽𝑘₁𝑘₂𝑒^−𝜇𝜏}, then there exists 𝜆₀ such that 𝑑₁𝜆²₀−𝑐𝜆₀+𝛽𝑘₂<0,𝑑₂𝜆²₀−𝑐𝜆₀+𝛽𝑘₁𝑘₂𝑒^−𝜇𝜏<0. (52) We can choose 𝜀_𝑖>0(𝑖=1,2,…,6) such that 𝛽𝜀₂1−𝑘₁−𝛽𝜀₁𝑘₂−𝜇𝜀₁<0,𝛽𝑒^−𝜇𝜏𝜀₂1−𝑘₁+𝜀₄+𝑘₂𝜀₄+𝛿𝜀₃−𝜇+𝛾𝜀₂<0,𝛾𝜀₂−𝜇+𝛿𝜀₃<0,𝛽𝜀₅1−𝑘₁+𝜀₄−𝛽𝑘₂+𝜇𝜀₄<0,𝛿𝜀₆−𝜇+𝛾𝜀₅<0,𝛾𝜀₅−𝜇+𝛿𝜀₆<0. (53)

For the above constants and suitable 𝑡_𝑖(𝑖=1,2,…,6) satisfying 𝑡₄>0 and 𝑡₅<𝑡₆, we define the continuous functions Φ(𝑡)=(𝜙₁(𝑡),𝜑₁(𝑡),𝜓₁(𝑡)) and Ψ(𝑡)=(𝜙₂(𝑡),𝜑₂(𝑡),𝜓₂(𝑡)) as follows: 𝜙₁(𝑡)=𝑘₁𝑒^𝜆0𝑡,𝑡⩽𝑡₁,𝑘₁+𝜀₁𝑒^−𝜆𝑡,𝑡>𝑡₁,𝜑₁(𝑡)=𝑘₂𝑒^𝜆0𝑡,𝑡⩽𝑡₂,𝑘₂+𝜀₂𝑒^−𝜆𝑡,𝑡>𝑡₂,𝜓₁(𝑡)=𝑘₃𝑒^𝜆0𝑡,𝑡⩽𝑡₃,𝑘₃+𝜀₃𝑒^−𝜆𝑡,𝑡>𝑡₃,𝜙₂(𝑡)=0,𝑡⩽𝑡₄,𝑘₁−𝜀₄𝑒^−𝜆𝑡,𝑡>𝑡₄,𝜑₂(𝑡)=0,𝑡⩽𝑡₅,𝑘₂−𝜀₅𝑒^−𝜆𝑡,𝑡>𝑡₅,𝜓₂(𝑡)=0,𝑡⩽𝑡₆,𝑘₃−𝜀₆𝑒^−𝜆𝑡,𝑡>𝑡₆. (54)

Lemma 6. Φ(𝑡)=(𝜙₁(𝑡),𝜑₁(𝑡),𝜓₁(𝑡)) is an upper solution of system (43).

Proof. Denote 𝑝₁(𝑡)=𝑑₁𝜙₁^′′(𝑡)−𝑐𝜙^′₁(𝑡)+𝛽𝜑₁(𝑡)1−𝜙₁(𝑡)−𝜇𝜙₁(𝑡),𝑝₂(𝑡)=𝑑₂𝜑₁^′′(𝑡)−𝑐𝜑^′₁(𝑡)+𝛽𝑒^−𝜇𝜏𝜑₁(𝑡−𝑐𝜏)1−𝜙₂(𝑡−𝑐𝜏)+𝛿𝜓₁(𝑡)−𝜇+𝛾𝜑₁(𝑡),𝑝₃(𝑡)=𝑑₃𝜓₁^′′(𝑡)−𝑐𝜓^′₁(𝑡)+𝛾𝜑₁(𝑡)−𝜇+𝛿𝜓₁(𝑡). (55)
If 𝑡⩽𝑡₁,  𝜙₁(𝑡)=𝑘₁𝑒^𝜆0𝑡,𝜑₁(𝑡)⩽𝑘₂𝑒^𝜆0𝑡, then 𝑝₁(𝑡)⩽𝑑₁𝑘₁𝜆²₀𝑒^𝜆0𝑡−𝑐𝑘₁𝜆₀𝑒^𝜆0𝑡+𝛽1−𝜙₁(𝑡)𝑘₂𝑒^𝜆0𝑡−𝜇𝑘₁𝑒^𝜆0𝑡⩽𝑘₁𝑒^𝜆0𝑡𝑑₁𝜆²₀−𝑐𝜆₀+𝛽𝑘₂<0. (56) If 𝑡>𝑡₁,  𝜙₁(𝑡)=𝑘₁+𝜀₁𝑒^−𝜆𝑡,𝜑₁(𝑡)⩽𝑘₂+𝜀₂𝑒^−𝜆𝑡. It follows that 𝑝₁(𝑡)⩽𝑒^−𝜆𝑡𝑑₁𝜀₁𝜆²+𝑐𝜀₁𝜆+𝛽𝜀₂1−𝑘₁−𝛽𝜀₁𝑘₂−𝜇𝜀₁. (57) Since 𝛽(1−𝑘₁)𝜀₂−𝛽𝜀₁𝑘₂−𝜇𝜀₁<0, then there exists a 𝜆^∗₁>0 such that 𝑝₁(𝑡)<0 for all 𝜆∈(0,𝜆^∗₁).
If 𝑡⩽𝑡₂,  𝜑₁(𝑡)=𝑘₂𝑒^𝜆0𝑡,𝜑₁(𝑡−𝑐𝜏)=𝑘₂𝑒^{𝜆0(𝑡−𝑐𝜏)},𝜓₁(𝑡)⩽𝑘₃𝑒^𝜆0𝑡, then 𝑝₂(𝑡)⩽𝑘₂𝑒^𝜆0𝑡𝑑₂𝜆²₀−𝑐𝜆₀+𝛽𝑒^{−(𝜇+𝑐𝜆0)𝜏}+𝛿𝑘₃𝑘₂−𝜇+𝛾⩽𝑘₂𝑒^𝜆0𝑡𝑑₂𝜆²₀−𝑐𝜆₀+𝛽𝑘₁𝑘₂𝑒^−𝜇𝜏<0. (58) If 𝑡>𝑡₂, we obtain that 𝑝₂(𝑡)⩽𝑒^−𝜆𝑡𝑑₂𝜀₂𝜆²+𝑐𝜀₂𝜆+𝛽𝑒^{−(𝜇−𝑐𝜆)𝜏}×𝜀₂1−𝑘₁+𝜀₄𝑒^{−𝜆(𝑡−𝑐𝜏)}+𝑘₂𝜀₄+𝛿𝜀₃−𝜇+𝛾𝜀₂⩽𝑒^−𝜆𝑡𝑑₂𝜀₂𝜆²+𝑐𝜀₂𝜆+𝐼₂(𝜆), (59) where 𝐼₁(𝜆)=𝛽𝑒^{−(𝜇−𝑐𝜆)𝜏}[𝜀₂(1−𝑘₁+𝜀₄)+𝑘₂𝜀₄]+𝛿𝜀₃−(𝜇+𝛾)𝜀₂, then 𝐼₁(0)<0 indicates that there exists a 𝜆^∗₂>0 such that 𝑝₂(𝑡)<0 for all 𝜆∈(0,𝜆^∗₂).
If 𝑡⩽𝑡₃,  𝜑₁(𝑡)⩽𝑘₂𝑒^𝜆0𝑡,  𝜓₁(𝑡)=𝑘₃𝑒^𝜆0𝑡, then 𝑝₃(𝑡)⩽𝑒^𝜆0𝑡𝑘₃𝑑₃𝜆²₀−𝑐𝑘₃𝜆₀+𝛾𝑘₂−𝜇+𝛿𝑘₃=𝑘₃𝜆₀𝑒^𝜆0𝑡𝑑₃𝜆₀−𝑐<0. (60) If 𝑡>𝑡₃,  𝜑₁(𝑡)⩽𝑘₂+𝜀₂𝑒^−𝜆𝑡,𝜓₁(𝑡)=𝑘₃+𝜀₃𝑒^−𝜆𝑡, then it follows 𝑝₃(𝑡)⩽𝑒^−𝜆𝑡𝑑₃𝜀₃𝜆²+𝑐𝜀₃𝜆+𝛾𝜀₂−𝜇+𝛿𝜀₃. (61) Note that 𝛾𝜀₂−(𝜇+𝛿)𝜀₃<0, then there exists a 𝜆^∗₃>0 such that 𝑝₃(𝑡)<0 for all 𝜆∈(0,𝜆^∗₃).
Clearly, for all 𝜆∈(0,min{𝜆^∗_𝑖,𝑖=1,2,3}),  𝑝_𝑖(𝑡)<0(𝑖=1,2,3). This completes the proof.

Lemma 7. Ψ(𝑡)=(𝜙₂(𝑡),𝜑₂(𝑡),𝜓₂(𝑡)) is a lower solution of system (43).

Proof. Denote 𝑞₁(𝑡)=𝑑₁𝜙₂^′′(𝑡)−𝑐𝜙^′₂(𝑡)+𝛽𝜑₂(𝑡)1−𝜙₂(𝑡)−𝜇𝜙₂(𝑡),𝑞₂(𝑡)=𝑑₂𝜑₂^′′(𝑡)−𝑐𝜑^′₂(𝑡)+𝛽𝑒^−𝜇𝜏𝜑₂(𝑡−𝑐𝜏)1−𝜙₁(𝑡−𝑐𝜏)+𝛿𝜓₂(𝑡)−𝜇+𝛾𝜑₂(𝑡),𝑞₃(𝑡)=𝑑₃𝜓₂^′′(𝑡)−𝑐𝜓^′₂(𝑡)+𝛾𝜑₂(𝑡)−𝜇+𝛿𝜓₂(𝑡). (62) If 𝑡⩽𝑡₄,  𝜙₂(𝑡)=0,𝑞₁(𝑡)=𝛽𝜑₂(𝑡)⩾0.
If 𝑡>𝑡₄,  𝜙₂(𝑡)=𝑘₁−𝜀₄𝑒^−𝜆𝑡,𝜑₂(𝑡)⩾𝑘₂−𝜀₅𝑒^−𝜆𝑡, then 𝑞₁(𝑡)⩾−𝜀₄𝑒^−𝜆𝑡𝑑₁𝜆²+𝑐𝜆+𝛽𝜀₅1−𝑘₁+𝜀₄𝜀₄−𝛽𝑘₂−𝜇. (63) Note that 𝛽𝜀₅(1−𝑘₁+𝜀₄)/𝜀₄−𝛽𝑘₂−𝜇<0, then there exists a 𝜆^∗₄>0 such that 𝑞₁(𝑡)>0 for all 𝜆∈(0,𝜆^∗₄).
If 𝑡⩽𝑡₅, 𝑞₂(𝑡)=0. If 𝑡>𝑡₅, 𝜑₂(𝑡)=𝑘₂−𝜀₅𝑒^−𝜆𝑡,  𝜑₂(𝑡−𝑐𝜏)⩾𝑘₂−𝜀₅𝑒^{−𝜆(𝑡−𝑐𝜏)},  𝜙₁(𝑡−𝑐𝜏)⩽𝑘₁+𝜀₁𝑒^{−𝜆(𝑡−𝑐𝜏)},   𝜓₂(𝑡)⩾𝑘₃−𝜀₆𝑒^−𝜆𝑡. It follows that 𝑞₂(𝑡)⩾−𝑒^−𝜆𝑡𝑑₂𝜀₅𝜆²+𝑐𝜀₅𝜆+𝛿𝜀₆−𝜇+𝛾𝜀₅. (64) Since 𝛿𝜀₆−(𝜇+𝛾)𝜀₅<0, then there exists a 𝜆^∗₅>0 such that 𝑞₂(𝑡)>0 for all 𝜆∈(0,𝜆^∗₅).
If 𝑡⩽𝑡₆, 𝜓₂(𝑡)=0,𝑞₃(𝑡)=𝛾𝜑₂(𝑡)⩾0. If 𝑡>𝑡₆, 𝜑₂(𝑡)=𝑘₂−𝜀₅𝑒^−𝜆𝑡,𝜓₂(𝑡)=𝑘₃−𝜀₆𝑒^−𝜆𝑡, then 𝑞₃(𝑡)⩾−𝜀₆𝑒^−𝜆𝑡𝑑₃𝜆²+𝑐𝜆+𝛾𝜀₅𝜀₆−𝜇+𝛿. (65) Note that 𝛾𝜀₅/𝜀₆−(𝜇+𝛿)<0, then there exists a 𝜆^∗₆>0 such that 𝑞₃(𝑡)>0 for all 𝜆∈(0,𝜆^∗₆).
Clearly, for all 𝜆∈(0,min{𝜆^∗_𝑖,𝑖=4,5,6}), 𝑞_𝑖(𝑡)>0(𝑖=1,2,3). This completes the proof.

Combining Lemmas 3–7 with the Schauder's fixed point theorem, we know that there exists a fixed point (𝜙^∗(𝑡),𝜑^∗(𝑡),𝜓^∗(𝑡)) of 𝐹 in Γ((𝜙,𝜑,𝜓),(𝜙,𝜑)), which gives a solution of (43). Furthermore, from (C3), we can obtain that Lim_{𝑡→−∞}𝜙^∗(𝑡),𝜑^∗(𝑡),𝜓^∗(𝑡)=(0,0,0),lim_𝑡→+∞𝜙^∗(𝑡),𝜑^∗(𝑡),𝜓^∗(𝑡)=𝑘₁,𝑘₂,𝑘₃, (66) which indicates that the fixed point satisfies the asymptotic boundary conditions (45). Therefore, there exists a travelling wave solution for system (42) connecting the steady state (0,0,0) and (𝑘₁,𝑘₂,𝑘₃), and accordingly, we have the following conclusion.

Theorem 8. Suppose ℛ₀>1 and 𝑑₁=𝑑₂=𝑑₃. For every 𝑐>𝑐^∗, system (5) always has a travelling wave solution with speed 𝑐 connecting the disease-free steady state 𝐸₀(1,0,0) and the endemic steady state 𝐸^∗(𝑆^∗,𝐼^∗,𝑅^∗).

4. Numerical Simulations

In this section, focusing on the travelling wave solutions of system (5), we perform some numerical simulations.

Firstly, we illustrate the results stated in Theorems 1 and 8. As in [6], we take parameter values motivated by bovine tuberculosis in a cattle herd, with time unit of one year. It is assumed that the exposed period is 9 months (i.e., 𝜏=0.75). Taking 𝜇=0.1, 𝛿=𝛾=0.5 (equal periods of infection and recovery before relapse), 𝛽=0.6, and 𝑑₁=𝑑₂=𝑑₃=0.01, system (5) with above coefficients has a unique endemic equilibrium 𝐸^∗(0.3294,0.3394,0.2828). For convenience, we truncate the spatial domain ℝ by [−10,10]. It is easy to show that ℛ₀=3.0363>1. If we use the following trivial functions as initial conditions: (𝑆(𝑥,𝑡),𝐼(𝑥,𝑡),𝑅(𝑥,𝑡))=(0.9,0.1,0),−10⩽𝑥⩽10,−𝜏⩽𝑡⩽0, (67) then the numerical simulations shown in Figures 1(a)–1(c) indicate that the endemic steady state 𝐸^∗ is stable. If we change the initial conditions, 𝐸^∗ keeps stable (see Figure 1(d); for simplicity, we show only the solutions of 𝐼(𝑥,𝑡)), which indicates that 𝐸^∗ is in fact globally asymptotically stable, although we fail to give a proof theoretically at present.

(a)

(b)

(c)

(d)

Figure 1 

The temporal solution found by numerical integration of system (5) with parameters 𝜇=0.1, 𝛿=𝛾=0.5, 𝛽=0.6, 𝜏=0.75,  𝑑1=𝑑2=𝑑3=0.01, and the initial conditions are (0.9, 0.1, 0) in (a)–(c), (0.6, 0.3, 0.1) in (A) of (d), and (0.99, 0.01, 0) in (B) of (d).

Now, to illustrate the existence of travelling wave solutions, we use the following piecewise functions as initial conditions: (𝑆()()(=()𝐸𝑥,𝑡,𝐼𝑥,𝑡,𝑅𝑥,𝑡))1,0,0,−10⩽𝑥<0,−𝜏⩽𝑡⩽0,_∗,0⩽𝑥⩽10,−𝜏⩽𝑡⩽0. (68) The numerical simulations shown in Figure 2 indicate that system (5) has a travelling wave solution with speed 𝑐≈0.07 connecting 𝐸₀ and 𝐸^∗.

(a)

(b)

(c)

(d)

Figure 2 

The traveling wave observed in system (5) with parameters: 𝜇=0.1, 𝛿=𝛾=0.5, 𝛽=0.6, 𝑑1=𝑑2=𝑑3=0.01, and 𝜏=0.75.

Secondly, we address the effects of parameters on the dynamics of system (5). For this part, the following parameter values will be fixed: 𝜇=0.1,𝛾=0.5, (69) and the initial conditions are the same as in (68).

We first investigate the impact of the diffusion rates. To this end, we fix 𝛿=0.5, 𝛽=0.6,  𝜏=0.75, and let the diffusivity 𝑑₁,𝑑₂, and 𝑑₃ vary. In this case, the values of ℛ₀ and 𝐸^∗ are the same as those in Figure 2 (Note that ℛ₀ and 𝐸^∗ are independent of the diffusion rates). Figure 3 is the simulation results of (5). In Figure 3(a), 𝑑₁=𝑑₂=𝑑₃=0.01, but in Figures 3(b)–3(d), we increase the values of 𝑑₁,𝑑₂, and 𝑑₃ to 0.04, respectively. By contrasting these figures, we can see that the travelling waves in Figures 3(c) and 3(d) are much fast than those in Figure 3(a), but the wave speed in Figure 3(b) is just the same as that in Figure 3(a), which indicates that both the diffusion rates of infective and recovered individuals can advance the time to arrive at the infection steady state, but the diffusion rate of susceptible individuals cannot. Even so, there is an obvious difference between the wave profiles in Figure 3(b) and the other three ones: the former has a hump in wave profile for 𝐼, but the latter does not seem to show that any more; that is, the diffusion rate of susceptible individuals could lead to nonmonotone travelling waves. In biological meaning, large diffusion of susceptible individuals may enhance the accumulative effect of infection.

(a)

(b)

(c)

(d)

Figure 3 

The traveling wave observed in system (5) with parameters: 𝜇=0.1, 𝛿=𝛾=0.5, 𝛽=0.6, 𝜏=0.75, and 𝑑1=𝑑2=𝑑3=0.01 in (a); 𝑑1=0.04, 𝑑2=𝑑3=0.01 in (b); 𝑑2=0.04, 𝑑1=𝑑3=0.01 in (c); 𝑑3=0.04, 𝑑1=𝑑3=0.01 in (d).

Next, we study the impact of 𝜏, that is, the average latent period of the disease. To this end, we fix 𝛿=0.5, 𝛽=0.6, 𝑑₁=0.02, and 𝑑₂=𝑑₃=0.01, and let 𝜏 vary. In this case, the values of ℛ₀ and 𝐸^∗ are different. (note that both ℛ₀ and 𝐸^∗ are decreasing in 𝜏). Figure 4 is the simulation results of (5). In Figure 4(a), 𝜏=0.75, but in Figure 4(b), we increase 𝜏 to 7.5. By contrasting the two figures, we can see that the travelling wave speed in Figure 4(b) is less than that in Figure 4(a), and the humps in Figure 4(a) disappear in Figure 4(b), which indicates that the average latent period of the disease can impact not only the wave speed but also the monotonicity of the wave profile, as well as the level of the positive steady state. In biological meaning, large latent period of the disease could lead to low and slow infection, and may weaken the accumulative effect of infection.

(a)

(b)

Figure 4 

The traveling wave observed in system (5) with parameters: 𝜇=0.1, 𝛿=𝛾=0.5,𝛽=0.6, 𝑑1=0.02, 𝑑2=𝑑3=0.01, and 𝜏=0.75 in (a); 𝜏=7.5 in (b).

Now we investigate the impact of transmission coefficient 𝛽. To this end, we fix 𝛿=0.5, 𝜏=0.75, 𝑑₁=𝑑₂=𝑑₃=0.01, and let 𝛽 vary. In this case, the values of ℛ₀ and 𝐸^∗ are different (note that both ℛ₀ and 𝐸^∗ are increasing in 𝛽). Figure 5 is the simulation results of (5). In Figure 5(a), 𝛽=0.3, but in Figure 5(b), we increase 𝛽 to 1.2. By contrasting the two figures, we can see that the travelling wave in Figure 5(b) is rather fast than that in Figure 5(a), but the travelling waves keep monotone as 𝜏 varies, which indicates that the transmission coefficient can obviously impact the wave speed, but it cannot change the monotonicity of the wave profile. In biological meaning, large transmission coefficient could lead to high and fast infection.

(a)

(b)

Figure 5 

The traveling wave observed in system (5) with parameters: 𝜇=0.1, 𝛿=𝛾=0.5, 𝜏=0.75, 𝑑1=𝑑2=𝑑3=0.01, and 𝛽=0.3 in (a); 𝛽=1.2 in (b).

At last, we study the impact of an important parameter 𝛿, that is, the rate at which an individual in the recovered class reverts to the infective class. In the following, for convenience, we call 𝛿 the relapse rate. Now we fix 𝛽=0.6, 𝜏=0.75, and 𝑑₁=𝑑₂=𝑑₃=0.01, and let 𝛿 vary. In this case, the values of ℛ₀ and 𝐸^∗ are also different (note that both ℛ₀ and 𝐸^∗ are increasing in 𝛿). Figure 6 is the simulation results of (5). In Figure 6(a), 𝛿=0.1, but in Figure 6(b), we increase 𝛿 to 2. By contrasting the two figures, we can see that the travelling wave in Figure 6(b) is rather fast than that in Figure 6(a), but the travelling waves keep monotone as 𝛿 varies, which indicates that the relapse rate can obviously impact the wave speed, but it cannot change the monotonicity of the wave profile. In biological meaning, large transmission coefficient could lead to high and fast infection, which is similar to the case of transmission coefficient.

(a)

(b)

Figure 6 

The traveling wave observed in system (5) with parameters: 𝜇=0.1, 𝛾=0.5, 𝛽=0.6, 𝜏=0.75, 𝑑1=𝑑2=𝑑3=0.01, and 𝛿=0.1 in (a); 𝛿=2 in (b).

To end this section, we point out that for the parameters given above, the following properties of the spread speed 𝑐 can be obtained from the simulations and discussions above: 𝑐 is increasing in 𝑑₂,𝑑₃,𝛽, and 𝛿 but decreasing in 𝜏. In addition, small 𝜏 and large 𝑑₁ may lead to nonmonotone travelling waves.

5. Conclusions

In this paper, we formulated an SEIR model with latent period, relapse and spatial diffusion. The dynamics of problem (5) and (6) was addressed. It was shown that the basic reproductive number of system (5) is given by ℛ₀=𝛽𝑒^−𝜇𝜏(𝜇+𝛿)/(𝜇(𝜇+𝛿+𝛾)), which describes the average number of newly infected individuals at the beginning of the infectious process. We have shown when the average number of newly infected individuals is less than unity, that is, ℛ₀<1; system (5) has a unique disease-free steady state 𝐸₀(1,0,0), which is asymptotically stable; when the average number of newly infected individuals is greater than unity, that is, ℛ₀>1, system (5) admits two steady states, 𝐸₀(1,0,0) and 𝐸^∗(𝑆^∗,𝐼^∗,𝑅^∗), where 𝑆^∗,𝐼^∗, and 𝑅^∗ are defined in (10). In this case, the endemic steady state 𝐸^∗ is always stable, while the disease-free steady state 𝐸₀ is unstable. Clearly, the spatial diffusion cannot destabilize the spatially homogenous steady state. Then by using the technique of upper and lower solutions and Schauder's fixed point theorem, we derived the existence of a travelling wave solution connecting the disease-free steady state and the endemic steady state. Furthermore, we gave some numeric simulations to illustrate the main results, and combining with numeric simulations, we discussed the effects of some parameters on the dynamics of system (5).

From the discussion in Section 4, we could conclude that the travelling wave speed is increasing in the diffusion rates of infective and recovered individuals, the transmission coefficient, and the relapse rate, but decreasing in the average latent period of the disease. In addition, small latent period or large diffusion rate of susceptible individuals may lead to nonmonotone travelling waves. In other words, to slow down or weaken the infection, we can reduce the diffusion rates, the transmission coefficient, or the relapse rate.

Acknowledgments

The authors are very grateful to the anonymous referees for careful reading and valuable comments which led to the improvement of their original paper. This work was supported by the National Natural Science Foundation of China (No. 11071254) and the Foundation of Shijiazhuang Mechanical Engineering College (Nos. JCB1201, YJJXM11004, and YJJXM12010).