Dynamical Properties of a Delay Prey-Predator Model with Disease in the Prey Species Only

Abstract

A three-dimensional ecoepidemiological model with delay is considered. We first investigate the existence and stability of the equilibria. We then study the effect of the time delay on the stability of the positive equilibrium. The existence of a Hopf bifurcation at the positive equilibrium is obtained through the study of an exponential polynomial equation with delay-dependent coefficients. Numerical simulation with a hypothetical set of data has been carried out to support the analytical findings.

1. Introduction

Recently, epidemiological models have become important tools in analyzing the spread and control of infectious diseases after the seminal model of Kermack-Mckendrick on SIRS (susceptible-infected-removed-susceptible) systems [1], and numerous mathematical models were developed to study a disease transmission, to evaluate the spread of epidemics, to understand the mechanisms of epidemics in order to prevent them or minimise the transmission of diseases and other measures. After the seminal models of Vito Volterra and Alfred James Lotka in the mid 1920s for predator-prey interactions [2, 3], predator-prey models have received a great deal of interest in different fields of science, from evolutionary and behavioural ecology to conservation and population biology and even economics [4–8].

It is well known that species do not exist alone in the natural world. It is of more biological significance to study the persistence-extinction threshold of each population in systems of two or more interacting species subjected to disease. Few investigators [9–12] have paid attention to study prey-predator model with infection.

In order to study the influence of disease on an environment where two or more interacting species are present. In this paper, we will put emphasis on such an ecoepidemiological system consisting of three species, namely, the sound prey (which is susceptible), the infected prey (which becomes infective by some viruses), and the predator population.

We have two populations:(1)the prey, whose total population density is denoted by 𝑁(𝑡),(2)the predator, whose population density is denoted by 𝑦(𝑡).

We make the following assumptions.(𝐴₁)In the absence of infection and predation, the prey population density grows logistically with carrying capacity 𝐾(𝐾>0) and an intrinsic birth rate constant 𝑟(𝑟>0) [13, 14], 𝑑𝑆(𝑡)(𝑡)𝑆(𝑡)𝑑𝑡=𝑟𝑆1−𝐾. (1.1)(𝐴₂)In the presence of disease, the total prey population 𝑁(𝑡) are divided into two distinct classes, namely, susceptible populations, 𝑆(𝑡), and infected populations, 𝐼(𝑡) [13, 14]. Therefore, at any time 𝑡, the total density of prey population is 𝑁(𝑡)=𝑆(𝑡)+𝐼(𝑡). (1.2)(𝐴₃)We assume that only susceptible prey 𝑆(𝑡) are capable of reproducing with logistic law (1.1); that is, the infected prey 𝐼(𝑡) are removed by death (say its death rate is a positive constant 𝜇), or by predation before having the possibility of reproducing. However, the infective population 𝐼(𝑡) still contributes with 𝑆(𝑡) to population growth toward the carrying capacity [13, 14].(𝐴₄)We assume that the force of infection at time 𝑡 is given by 𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼(𝑡−𝜏), where 𝛽 is the average number of contacts per infective per day, 𝜇₁ is disease-induced death rate and 𝜏>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 prey [15]. Hence, the 𝑆𝐼 model of the infected prey is 𝑑𝑆(𝑡)(𝑡)𝑆(𝑡)𝑑𝑡=𝑟𝑆1−𝐾−𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼(),(𝑡)𝑡−𝜏𝑑𝐼𝑑𝑡=𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼(𝑡−𝜏)−𝜇𝐼(𝑡). (1.3)(𝐴₅)Numerous field studies show that infected prey are more vulnerable to predation compared with their noninfected counterpart [16–18]. Lafferty and Morris [17] quantified that the predation rates on infected prey may be 31 times higher compared to that on susceptible prey. Thus, we consider the case when the predator mainly eats the infected prey. We assume that the predator eats only the infected prey with Leslie-Gower ratio-dependent schemes [19–23]. That is to say, the predator consumes the prey according to the ratio-dependent functional response [24–26] and the predator grows logistically with intrinsic growth rate 𝛿 and carrying capacity proportional to the prey populations size 𝐼(𝑡).

From the above assumptions we have the following model: 𝑑𝑆(𝑡)(𝑡)𝑆(𝑡)𝑑𝑡=𝑟𝑆1−𝐾−𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼(),(𝑡)𝑡−𝜏𝑑𝐼𝑑𝑡=𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼()−(𝑡)𝐼(𝑡)𝑡−𝜏𝑐𝑦(𝑡)(𝑡)(𝑡),(𝑡)𝑚𝑦+𝐼−𝜇𝐼𝑑𝑦(𝑡)(𝑡)𝑑𝑡=𝛿𝑦1−ℎ𝑦𝐼(𝑡). (1.4)

All the parameters of system (1.4) are positive. The initial conditions for system (1.4) defined in the Banach space 𝒞⁺=𝜑∈𝒞[−𝜏,0],ℝ³₊∣𝜑₁(𝜃)=𝑆(𝜃),𝜑₂(𝜃)=𝐼(𝜃),𝜑₃(𝜃)=𝑦(𝜃), (1.5) where ℝ³₊={(𝑆,𝐼,𝑦)∈ℝ³∣𝑆≥0,𝐼≥0,𝑦≥0}. We assume that 𝜑_𝑖(0)>0(𝑖=1,2,3) by the biological meaning.

It is well known by the fundamental theorem of functional differential equations that system (1.4) has a unique solution (𝑆(𝑡),𝐼(𝑡),𝑦(𝑡)) satisfying initial conditions (1.5).

In this paper, we will discuss the dynamical properties (e.g., the stability of equilibrium, the existence of Hopf bifurcation around the positive equilibrium) of the delay ecoepidemic model (1.4). This paper is organized as follows. Section 2 gives positivity and boundedness of solutions of system (1.4). Section 3 gives the analysis of existence and uniqueness of positive equilibrium as well as its stability, by studying the characteristic equation associated with system (1.4) which takes the form of a third-degree exponential polynomial with delay-dependent coefficients. In Section 4, some numerical simulations will be given to illustrate our theoretical results. Section 5 discusses the results.

2. Positivity and Boundedness of Solutions

It is important to show positivity and boundedness for the system (1.4) because the variables 𝑆, 𝐼 and 𝑦 represent populations. Positivity implies that the populations survive and boundedness may be interpreted as a natural restriction to growth as a consequence of limited resources.

Let us put system (1.4) in a vector form by setting 𝑋=col(𝑆,𝐼,𝑦)∈ℝ³, 𝐺(𝑋)=𝐺₁(𝑋)𝐺₂(𝑋)𝐺₃(𝑋)=𝑆𝑟𝑆1−𝐾−𝛽𝑒^−𝜇1𝜏𝑆𝐼(𝑡−𝜏)𝛽𝑒^−𝜇1𝜏𝑆𝐼(𝑡−𝜏)−𝑐𝑦𝐼𝑚𝑦+𝐼−𝜇𝐼𝛿𝑦1−ℎ𝑦𝐼, (2.1) where 𝐺∶𝒞⁺→ℝ³ and 𝐺∈𝒞^∞(ℝ³). Then system (1.4) becomes 𝑋̇=𝐺(𝑋), (2.2) with 𝑋(𝜃)=(𝜙₁(𝜃),𝜙₂(𝜃),𝜙₃(𝜃))∈𝒞⁺ and 𝜙_𝑖(0)>0(𝑖=1,2,3). It is easy to check in (2.1) that whenever choosing 𝑋(𝜃)∈𝒞⁺ such that 𝑥_𝑖=0, then 𝐺_𝑖(𝑋)|_{𝑥𝑖(𝑡)=0,𝑋∈ℝ}³₊≥0, (𝑖=1,2,3), where 𝑥₁=𝑆(𝑡), 𝑥₂=𝐼(𝑡) and 𝑥₃=𝑦(𝑡).

Due to lemma in [27] any solutions of (2.2) with 𝑋(𝜃)∈𝒞⁺, that is, 𝑋(𝑡)=𝑋(𝑡,𝑋(𝜃)), is such that 𝑋(𝜃)∈ℝ³ for all 𝑡>0.

Because the analysis of the boundedness of solutions of system (1.4) is similar to the one in Section 2 in [28], we omit it. Let 𝑀≥0, and we can get that all solutions of system (1.4) with initial condition enter the region 𝐵={(𝑆(𝑡),𝐼(𝑡),𝑦(𝑡))∣0≤𝑆(𝑡),𝐼(𝑡),𝑦(𝑡)≤𝑀}.

3. Stability Analysis and Hopf Bifurcation

In this section, we focus on investigating the stability of the equilibria and Hopf bifurcation around the positive equilibrium of the system (1.4). System (1.4) has the boundary equilibrium 𝐸₁(𝜇/𝛽𝑒^−𝜇1𝜏,(𝑟/𝛽𝑒^−𝜇1𝜏)(1−𝜇/𝛽𝑒^−𝜇1𝜏𝐾),0)≜(𝑆₁,𝐼₁,𝑦₁) and the positive equilibrium 𝐸₂(𝑆₂,𝐼₂,𝑦₂), where 𝑆₂=(𝜇ℎ+𝑐+𝜇𝑚)/𝛽𝑒^−𝜇1𝜏(𝑚+ℎ), 𝐼₂=𝑟(𝛽𝑒^−𝜇1𝜏𝑚𝐾+𝛽𝑒^−𝜇1𝜏ℎ𝐾−𝜇ℎ−𝑐−𝜇𝑚)/(𝛽𝑒^−𝜇1𝜏)²𝐾(𝑚+ℎ), 𝑦₂=𝑟(𝛽𝑒^−𝜇1𝜏𝑚𝐾+𝛽𝑒^−𝜇1𝜏ℎ𝐾−𝜇ℎ−𝑐−𝜇𝑚)/(𝛽𝑒^−𝜇1𝜏)²ℎ𝐾(𝑚+ℎ). Clearly, if 1−𝜇/𝛽𝑒^−𝜇1𝜏𝐾>0, then 𝐸₁ exists and remains positive, and 𝐸₂ exists and remains positive if 𝛽>(1/𝐾𝑒^−𝜇1𝜏)(𝜇+𝑐/(𝑚+ℎ))≜𝛽₀.

Let 𝐸^∗(𝑆^∗,𝐼^∗,𝑦^∗) be any arbitrary equilibrium. Then the characteristic equation about 𝐸^∗ is given by 𝑟−2𝑟𝑆^∗𝐾−𝛽𝑒^−𝜇1𝜏𝐼^∗−𝜆−𝛽𝑒^−𝜇1𝜏𝑆^∗𝑒^−𝜆𝜏0𝛽𝑒^−𝜇1𝜏𝐼^∗𝛽𝑒^−𝜇1𝜏𝑆^∗𝑒^−𝜆𝜏−𝑐𝑚𝑦^∗2𝑚𝑦^∗+𝐼^∗2−𝜇−𝜆−𝑐𝐼^∗2𝑚𝑦^∗+𝐼^∗20𝛿ℎ𝑦^∗2𝐼^∗2𝛿−2𝛿ℎ𝑦^∗𝐼^∗−𝜆=0. (3.1)

For equilibrium 𝐸₁, (3.1) reduces to −𝑟𝑆₁𝐾−𝛽𝑒^−𝜇1𝜏𝐼₁−𝜆−𝛽𝑒^−𝜇1𝜏𝑆₁𝑒^−𝜆𝜏0𝛽𝑒^−𝜇1𝜏𝐼₁𝛽𝑒^−𝜇1𝜏𝑆₁𝑒^−𝜆𝜏−𝜇−𝜆−𝑐00𝛿−𝜆=0. (3.2)

It is easy to see that the equilibrium 𝐸₁ is unstable.

For equilibrium 𝐸₂, (3.1) reduces to 𝑃(𝜆,𝜏)+𝑄(𝜆,𝜏)𝑒^−𝜆𝜏=0 (3.3) with 𝑃(𝜆,𝜏)=𝜆³+𝐴₁(𝜏)𝜆²+𝐴₂(𝜏)𝜆+𝐴₃(𝜏),𝑄(𝜆,𝜏)=𝐵₁(𝜏)𝜆²+𝐵₂(𝜏)𝜆+𝐵₃(𝜏), (3.4) where 𝐴₁(𝜏)=𝛿+𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝑟𝑆₂𝐾,𝐴₂(𝜏)=𝛿𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝑟𝑆₂𝐾𝛿+𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝛿ℎ𝑐𝐼²₂𝑚𝑦₂+𝐼₂²,𝐴₃(𝜏)=𝑟𝑆₂𝐾𝛿𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝛿ℎ𝑐𝐼²₂𝑚𝑦₂+𝐼₂²,𝐵₁(𝜏)=−𝛽𝑒^−𝜇1𝜏𝑆₂,𝐵₂(𝜏)=−𝛿𝛽𝑒^−𝜇1𝜏𝑆₂−𝑟𝛽𝑒^−𝜇1𝜏𝑆²₂𝐾+𝛽𝑒^{−𝜇1𝜏2}𝑆₂𝐼₂,𝐵₃(𝜏)=−𝑟𝛿𝛽𝑒^−𝜇1𝜏𝑆²₂𝐾+𝛿𝛽𝑒^{−𝜇1𝜏2}𝑆₂𝐼₂. (3.5) For 𝜏=0, the transcendental equation (3.3) reduces to (3.4) 𝜆³+𝐴₁(0)+𝐵₁(0)𝜆²+𝐴₂(0)+𝐵₂(0)𝜆+𝐴₃(0)+𝐵₃(0)=0. (3.6) We can easily get 𝐴₁(0)+𝐵₁(0)=𝑟𝑆₂𝐾+𝛿−𝛿^∗=𝑟𝛽₀𝛽+𝛿−𝛿^∗𝐴>0,₂(0)+𝐵₂(0)=𝑟𝑆₂𝐾𝛿−𝛿^∗+𝛽²𝑆₂𝐼₂𝑐=𝑟𝜇+𝑚+ℎ1+𝛿−𝛿^∗−𝐾𝛽𝛽₀𝛽,𝐴₃(0)+𝐵₃(0)=𝛿𝛽²𝑆₂𝐼₂𝑐=𝑟𝛿𝜇+𝑚+ℎ1−𝛽₀𝛽>0 (3.7) if 𝛿>𝛿^∗ and 𝛽>𝛽₀ (obviously, 𝛽₀>𝛽₀), where 𝛿^∗=𝑐ℎ/(𝑚+ℎ)² and 𝛽₀=(1/𝐾)(𝜇+𝑐/(𝑚+ℎ)).

By Routh-Hurwitz criterion, we can get the following theorem.

Theorem 3.1. When 𝜏=0, the positive equilibrium 𝐸₂ is locally asymptotically stable provided that the conditions(𝐻₁)𝑟𝛽₀>𝛿^∗𝛽, (𝐻₂)𝛿>𝛿^∗, (𝐻₃)𝛽>𝛽₀.

In the following, we investigate the existence of purely imaginary roots 𝜆=𝑖𝜔(𝜔>0) to (3.3). Equation (3.3) takes the form of a third-degree exponential polynomial in 𝜆, which all the coefficients of 𝑃 and 𝑄 depending on 𝜏. Beretta and Kuang [29] established a geometrical criterion which gives the existence of purely imaginary of a characteristic equation with delay-dependent coefficients.

In order to apply the criterion due to Beretta and Kuang [29], we need to verify the following properties for all 𝜏∈[0,𝜏_max), where 𝜏_max is the maximum value which 𝐸₂ exists. We can easily obtain that 𝜏_max=(1/𝜇₁)ln(𝐾𝛽(𝑚+ℎ)/𝜇(𝑚+ℎ)+𝑐).(a)𝑃(0,𝜏)+𝑄(0,𝜏)≠0;(b)𝑃(𝑖𝜔,𝜏)+𝑄(𝑖𝜔,𝜏)≠0;(c)limsup{|𝑃(𝜆,𝜏)/𝑄(𝜆,𝜏)|∶|𝜆|→∞,Re𝜆≥0}<1;(d)𝐹(𝜔,𝜏)=|𝑃(𝑖𝜔,𝜏)|²−|𝑄(𝑖𝜔,𝜏)|² has a finite number of zeros;(e)Each positive root 𝜔(𝜏) of 𝐹(𝜔,𝜏)=0 is continuous and differentiable in 𝜏 whenever it exists.

Here, 𝑃(𝜆,𝜏) and 𝑄(𝜆,𝜏) are defined as in (3.4).

Let 𝜏∈[0,𝜏_max). It is easy to see that 𝑃(0,𝜏)+𝑄(0,𝜏)=𝐴₃(𝜏)+𝐵₃(𝜏)>0. (3.8) This implies that (a) is satisfied, and (b) is obviously true because 𝑃(𝑖𝜔,𝜏)+𝑄(𝑖𝜔,𝜏)=−𝑖𝜔³−𝐴₁(𝜏)𝜔²+𝑖𝐴₂(𝜏)𝜔+𝐴₃(𝜏)−𝐵₁(𝜏)𝜔²+𝑖𝐵₂(𝜏)𝜔+𝐵₃(𝜏)=−𝐴₁(𝜏)+𝐵₁(𝜏)𝜔²+𝐴₃(𝜏)+𝐵₃(𝜏)+𝑖−𝜔³+𝐴₂(𝜏)+𝐵₂(𝜏)𝜔≠0. (3.9)

From (3.4) we know that lim_{|𝜆|→+∞}𝑄(𝜆,𝜏)𝑃(𝜆,𝜏)=0. (3.10) Therefore (c) follows.

Let 𝐹 be defined as in (d). From |𝑃(𝑖𝜔,𝜏)|²=−𝜔³+𝐴₂(𝜏)𝜔²+−𝐴₁(𝜏)𝜔²+𝐴₃(𝜏)²=𝜔⁶+𝐴²₁(𝜏)−2𝐴₂(𝜏)𝜔⁴+𝐴²₂(𝜏)−2𝐴₁(𝜏)𝐴₃(𝜏)𝜔²+𝐴²₃(𝜏),|𝑄(𝑖𝜔,𝜏)|²=𝐵²₁(𝜏)𝜔⁴+𝑏²₂(𝜏)−2𝐵₁(𝜏)𝐵₃(𝜏)𝜔²+𝐵²₃(𝜏), (3.11) we have 𝐹(𝜔,𝜏)=𝜔⁶+𝑎₁(𝜏)𝜔⁴+𝑎₂(𝜏)𝜔²+𝑎₃(𝜏), (3.12) where 𝑎₁(𝜏)=𝐴²₁(𝜏)−2𝐴₂(𝜏)−𝐵²₁(𝜏)=𝛿²+𝜇²+𝑐²𝑚²(𝑚+ℎ)²+𝑟𝑆₂𝐾²+2𝜇𝑐𝑚(𝑚+ℎ)²−2𝛿𝑐ℎ(𝑚+ℎ)²−𝛽²𝑒^{−𝜇1𝜏2}𝑆²₂,𝑎₂(𝜏)=𝐴²₂(𝜏)+2𝐵₁(𝜏)𝐵₃(𝜏)−2𝐴₁(𝜏)𝐴₂(𝜏)−𝐵²₂(𝜏)=𝛿𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝑟𝑆₂𝐾𝛿+𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝛿ℎ𝑐𝐼²₂𝑚𝑦₂+𝐼₂²×𝛿𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝑟𝑆₂𝐾𝛿+𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝛿ℎ𝑐𝐼²₂𝑚𝑦₂+𝐼₂²−2𝛿+𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝑟𝑆₂𝐾−𝛿𝛽𝑒^−𝜇1𝜏𝑆₂²−𝑟𝛽𝑒^−𝜇1𝜏𝑆²₂𝐾−𝛽𝑒^{−𝜇1𝜏2}𝑆₂𝐼₂²,𝑎₃(𝜏)=𝐴²₃(𝜏)−𝐵²₃(𝜏)=𝑟𝑆₂𝐾²𝛿𝜇+𝑐𝑚𝑦²₂𝑚𝑦₂+𝐼₂²+𝛿ℎ𝑐𝐼²₂𝑚𝑦₂+𝐼₂²²−𝛿𝛽𝑒^{−𝜇1𝜏2}𝑆₂𝐼₂−𝑟𝛿𝛽𝑒^−𝜇1𝜏𝑆²₂𝐾². (3.13)

It is obvious that property (d) is satisfied, and by Implicit Function Theorem, (e) is also satisfied.

Now let 𝜆=𝑖𝜔(𝜔>0) be a root of (3.3). Substituting it into (3.3) and separating the real and imaginary parts yields −𝐴₃(𝜏)+𝐴₁(𝜏)𝜔²=𝐵₂(𝜏)𝜔sin𝜔𝜏+𝐵₃(𝜏)−𝐵₁(𝜏)𝜔²𝜔𝜔cos𝜔𝜏,₃−𝐴₂(𝜏)𝜔=𝐵₂(𝜏)𝜔cos𝜔𝜏+𝐵₃(𝜏)−𝐵₁(𝜏)𝜔²𝜔sin𝜔𝜏. (3.14)

From (3.14) it follows that𝐵sin𝜔𝜏=¹(𝜏)𝜔⁵+𝐴₁(𝜏)𝐵₂(𝜏)−𝐴₂(𝜏)𝐵₁(𝜏)−𝐵₃(𝜏)𝜔³𝐵²₁(𝜏)𝜔⁴𝐵²₂(𝜏)−2𝐵₁(𝜏)𝐵₃(𝜏)𝜔²+𝐵²₃(𝜏)+𝐴₂(𝜏)𝐵₃(𝜏)−𝐴₃(𝜏)𝐵₃(𝜏)𝜔𝐵²₁(𝜏)𝜔⁴𝐵²₂(𝜏)−2𝐵₁(𝜏)𝐵₃(𝜏)𝜔²+𝐵²₃(𝜏),𝐵cos𝜔𝜏=²(𝜏)−𝐴₁(𝜏)𝐵₁(𝜏)𝜔⁴+𝐴₁(𝜏)𝐵₃(𝜏)+𝐴₃(𝜏)𝐵₁(𝜏)−𝐴₂(𝜏)𝐵₂(𝜏)𝜔²𝐵²₁(𝜏)𝜔⁴𝐵²₂(𝜏)−2𝐵₁(𝜏)𝐵₃(𝜏)𝜔²+𝐵²₃(𝜏)−𝐴₃(𝜏)𝐵₃(𝜏)𝐵²₁(𝜏)𝜔⁴𝐵²₂(𝜏)−2𝐵₁(𝜏)𝐵₃(𝜏)𝜔²+𝐵²₃(𝜏). (3.15a)(3.15b)

By the definitions of 𝑃(𝜆,𝜏), 𝑄(𝜆,𝜏) as in (3.3), and applying the property (a), (3.7) can be written as 𝑃()sin𝜔𝜏=Im𝑖𝜔,𝜏𝑄(),𝑃()𝑖𝜔,𝜏cos𝜔𝜏=−Re𝑖𝜔,𝜏𝑄(),𝑖𝜔,𝜏 (3.16) which yields |𝑃(𝑖𝜔,𝜏)|²=|𝑄(𝑖𝜔,𝜏)|². (3.17)

Assume that 𝒥⊆ℝ₊₀={𝜁∣𝜁>0} is the set where 𝜔(𝜏) is a positive root of 𝐹(𝜔,𝜏)=|𝑃(𝑖𝜔,𝜏)|²−|𝑄(𝑖𝜔,𝜏)|² (3.18) and for 𝜏∉𝒥, 𝜔(𝜏) is not definite. Then for all 𝜏 in 𝒥, 𝜔(𝜏) satisfied 𝐹(𝜔,𝜏)=0. (3.19) The polynomial function 𝐹 can be written as 𝐹(𝜔,𝜏)=𝑓(𝜔²,𝜏), where 𝑓 is a third-degree polynomial, defined by 𝑓(𝑧,𝜏)=𝑧³+𝑎₁𝑧²+𝑎₂𝑧+𝑎₃. (3.20)

Depending on the determinant of equation 𝑓(𝑧,𝜏)=𝑧³+𝑎₁𝑧²+𝑎₂𝑧+𝑎₃=0, (3.21)𝑀=(𝑞/2)²+(𝑝/3)³, where 𝑝=𝑎₂−(1/3)𝑎²₁, 𝑞=(2/27)𝑎³₁−(1/3)𝑎₁𝑎₂+𝑎₃, there are three cases for the solutions of (3.21).(i)If 𝑀>0, (3.21) has a real root and a pair of conjugate complex roots. The real root is positive and is given by 𝜈₁=³−𝑞2+√𝑀+³−𝑞2−√1𝑀−3𝑎₁. (3.22)(ii)If 𝑀=0, (3.21) has three real roots, of which two are equal. In particular, if 𝑎₁>0, there exists only one positive root, 𝜈₁=2³−𝑞/2−𝑎₁/3; If 𝑎₁<0, there exists a positive root 𝜈₁=2³−𝑞/2−𝑎₁/3 for ³−𝑞/2>−𝑎₁/3, and there exist three positive roots for 𝑎₁/6<³−𝑞/2<−𝑎₁/3, 𝜈₁=2³−𝑞/2−𝑎₁/3, 𝜈₂=𝜇₃=−³−𝑞/2−𝑎₁/3.(iii)If 𝑀<0, there are three distinct real roots, 𝜈₁=2(|𝑝|/3)cos(𝜑/3)−𝑎₁/3, 𝜈₂=2(|𝑝|/3)cos(𝜑/3+2𝜋/3)−𝑎₁/3, 𝜈₃=2(|𝑝|/3)cos(𝜑/3+4𝜋/3)−𝑎₁/3, where cos𝜑=−𝑞/2(|𝑝|/3)³. Furthermore, if 𝑎₁>0, there exists only one positive root. Otherwise, if 𝑎₁<0, there may exist either one or three positive real roots. If there is only one positive real root, it is equal to max(𝜈₁,𝜈₂,𝜈₃).

Clearly, the number of positive real roots of (3.21) depends on the sign of 𝑎₁. When 𝑎₁≥0, (3.21) has only one positive real root. Otherwise, there may exist three positive real roots.

It is easy to know that 𝑎₁=𝐴²₁(𝜏)−2𝐴₂(𝜏)−𝐵²₁(𝜏)=𝛿²+𝜇²+𝑐²𝑚²(𝑚+ℎ)²+𝑟𝑆₂𝐾²+2𝜇𝑐𝑚(𝑚+ℎ)²−2𝛿𝑐ℎ(𝑚+ℎ)²−𝛽²𝑒^{−𝜇1𝜏2}𝑆²₂>0 (3.23)

if condition(𝐻₄)∶𝛿²+𝜇²+𝑐²𝑚²(𝑚+ℎ)²+𝑟𝑆₂𝐾²+2𝜇𝑐𝑚(𝑚+ℎ)²>2𝛿𝑐ℎ(𝑚+ℎ)²+𝛽²𝑒^{−𝜇1𝜏2}𝑆²₂ (3.24)

holds.

Hence, (3.21) has only one positive real root if condition (𝐻₄) holds. We denote by 𝑧₊ this positive real root. Thus, (3.21) has only one positive real root 𝜔=√𝑧₊, and the critical values of 𝜏 and 𝜔(𝜏) are impossible to solve explicitly, so we will use the procedure described in Beretta and Kuang [29]. According to this procedure, we define 𝜃(𝜏)∈[0,2𝜋) such that sin𝜃(𝜏) and cos𝜃(𝜏) are given by the right-hand sides of (3.15a) and (3.15b), respectively, with 𝜃(𝜏) given by (3.21). This define 𝜃(𝜏) in a form suitable for numerical evaluation using standard software.

And the relation between the argument 𝜃(𝜏) and 𝜔(𝜏)𝜏 in (3.19) for 𝜏>0 must be 𝜔(𝜏)𝜏=𝜃(𝜏)+2𝑛𝜋,𝑛=0,1,2,…. (3.25) Hence we can define the maps: 𝜏_𝑛∶𝒥→ℝ₊₀ given by 𝜏_𝑛(𝜏)𝜃(𝜏)∶=+2𝑛𝜋𝜔(𝜏),𝜏_𝑛>0,𝑛=0,1,2,…, (3.26) where a positive root 𝜔(𝜏) of (3.19) exists in 𝒥.

Let us introduce the functions 𝑀_𝑛(𝜏)∶𝒥→ℝ, 𝑀_𝑛(𝜏)𝜃(𝜏)=𝜏−+2𝑛𝜋𝜔(𝜏),𝑛=0,1,2,…. (3.27) that are continuous and differentiable in 𝜏. Thus, we give the following theorem which is due to Beretta and Kuang [29].

Theorem 3.2. Suppose 𝛽>𝛽₀ and condition (𝐻₄) hold. Assume that 𝜔(𝜏) is a positive root of (3.4) defined for 𝜏∈𝒥, 𝒥⊆ℝ₊₀, and at some 𝜏^∗∈𝒥, 𝑀_𝑛(𝜏^∗)=0 for some 𝑛∈ℕ₀={1,2,3,…}. Then a pair of simple conjugate pure imaginary roots 𝜆=±𝑖𝜔 exists at 𝜏=𝜏^∗ which crosses the imaginary axis from left to right if 𝛿(𝜏^∗)>0 and crosses the imaginary axis from right to left if 𝛿(𝜏^∗)<0, where 𝛿𝜏^∗=sign𝜕𝐹𝜕𝜔𝜔𝜏^∗,𝜏^∗sign𝑑𝑀_𝑛(𝜏)𝑑𝜏_𝜏=𝜏^∗. (3.28)
Since𝜕𝐹()𝜕𝜔𝜔,𝜏=2𝜔𝜕𝑓𝜔𝜕𝑧²,𝜏, (3.29) condition (3.28) is equivalent to sign𝑑Re𝜆𝑑𝜏_{𝜆=𝑖𝜔(𝜏}^∗₎=sign𝜕𝑓𝜔𝜕𝑧²𝜏^∗,𝜏^∗sign𝑑𝑀_𝑛(𝜏)𝑑𝜏_𝜏=𝜏^∗. (3.30)
We can easily observe that 𝑀_𝑛(0)<0. Moreover, for all 𝜏∈𝒥, 𝑀_𝑛(𝜏)>𝑀_𝑛+1(𝜏), when 𝑛∈ℕ. Therefore, if 𝑀₀ has no root in 𝒥, then the 𝑀_𝑛 functions have no root in 𝒥 and, if the function 𝑀_𝑛(𝜏) has positive roots 𝜏∈𝒥 for some 𝑛∈ℕ, there exists at least one root satisfying𝑑𝑀_𝑛(𝜏)𝑑𝜏>0. (3.31)

Using Theorem 3.1, we can conclude the existence of a Hopf bifurcation as stated in the following theorem.

Theorem 3.3. Assume the following.(1)If the function 𝑀₀(𝜏) has no positive root in 𝒥, then the positive equilibrium 𝐸₂ is locally asymptotically stable for all 𝜏≥0.(2)If the function 𝑀₀(𝜏) has at least one positive root in 𝒥, then there exists 𝜏^∗∈𝒥, such that the positive equilibrium 𝐸₂ is locally asymptotically stable for 0≤𝜏<𝜏^∗ and becomes unstable for 𝜏≥𝜏^∗, with a Hopf bifurcation occurring when 𝜏=𝜏^∗, if and only if sign𝜕𝑓𝜔𝜕𝑧²,𝜏>0. (3.32)

4. Numerical Simulations

To determine whether it is possible that there are parameter values such that a stability switch is possible, we choose a set of parameters as follows: 𝑟=0.1, 𝐾=500, 𝛽=0.001, 𝜇=0.03, 𝑐=8, 𝛿=0.2, 𝑚=150, ℎ=0.5, 𝜇₁=0.04.

Using Maple 10, we compute the coefficients 𝑎₁(𝜏), 𝑎₂(𝜏), 𝑎₃(𝜏) for the above values. We can get 𝑎₁(𝜏)=0.2829795477+0.004706297728𝑒^0.04𝜏>0. Hence the polynomial function 𝑓, defined in (3.20), has only one positive real root.

The function 𝑀₀ is drawn for 𝜏∈𝒥 in Figure 1. One can see that there are two critical values of the delay 𝜏 for which stability switches occur, that is, 𝜏^∗=4.69 and 𝜏^∗∗=13.78. In particulary from Theorem 3.3, a Hopf bifurcation occurs when 𝜏 is approximately equal to 4.69. Thus, periodic solutions appear.

Figure 1 

The graph of the function 𝑀0(𝜏). Two critical values of 𝜏 (i.e., 𝜏∗ and 𝜏∗∗), for which stability switches can occur, appear.

Using DDE23, a Matlab Solver for delay differential equation, we can compute the solutions of (1.4) for the above-mentioned values of the parameters, see Figures 2–4. Before the Hopf bifurcation occurs (0≤𝜏<𝜏^∗), solutions are asymptotically stable and converge to the positive equilibrium. When 𝜏^∗≤𝜏≤𝜏^∗∗, periodic solutions appear (See Figure 3). When 𝜏 increases, longer-periods oscillations exist. For 𝜏=𝜏^∗∗(=13.78), a stability switch occurs, that is, the positive equilibrium becomes asymptotically stable again and solutions converge to the positive equilibrium (Figure 4). In Figure 5, using the time delay 𝜏 as the bifurcation parameter, we can obtain the bifurcation diagram by Matlab package DDE-BIFTOOL.

(a)

(b)

(c)

(d)

Figure 2 

Time evolution of all the population for the model (1.4) with 𝜏=3.3 and initial value (50,80,160). (a), (b) and (c) are 𝑆, 𝐼 and 𝑦 versus 𝑡, respectively, and (d) is the projected 𝑆-𝐼-𝑦 phase plane.

(a)

(b)

(c)

(d)

Figure 3 

Time evolution of all the population for the model (1.4) with 𝜏=5.6 and initial value (50,80,160). (a), (b) and (c) are 𝑆, 𝐼 and 𝑦 versus 𝑡, respectively, and (d) is the projected 𝑆-𝐼-𝑦 phase plane.

(a)

(b)

(c)

(d)

Figure 4 

Time evolution of all the population for the model (1.4) with 𝜏=15.2 and initial value (50,80,160). (a), (b) and (c) are 𝑆, 𝐼 and 𝑦 versus 𝑡, respectively, and (d) is the projected 𝑆-𝐼-𝑦 phase plane.

Figure 5 

Bifurcation diagram for 𝜏.

5. Discussion

In this paper, we consider the ecoepidemiology model with time delay and Leslie-Gower schemes. The existence of a Hopf bifurcation at a positive steady-state is obtained through the study of an exponential polynomial characteristic equation with delay-dependent coefficients. That is to say, the positive steady-state of this model depends on the maturation period 𝜏 through a factor 𝑒^−𝜇1𝜏, where 𝜇₁ is the rate of disease-induced death. We can find stability switch occurs from stability to instability. Our computer simulation can easily confirm this. So the time delay plays a stabilizing role for our model.

In addition, if 𝜇₁=0, the steady state does not depend on time delay 𝜏. From [28], we can see that the stability switch occurs toward instability if it is stable at 𝜏=0, and further increase of time delay does not cause the occurrence of stability switch. This implies that both time delay and disease-induced death rate are necessary for the asymptotical stability of system (1.4).

Finally, although we give the reason why we just think of such situation by citing that the rate of the predation on infected prey may be 31 times higher compared to that on susceptible prey [17], it is not usual in real world. In fact, if we deal with a huge population of the prey, then the susceptible population may be too large to ignore their role as prey of the predator. Thus, we can present the more realistic model as follows: 𝑑𝑆(𝑡)(𝑡)𝑆(𝑡)𝑑𝑡=𝑟𝑆1−𝐾−𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼(𝑡−𝜏)−𝜇₁𝑆(𝑡),𝑑𝐼(𝑡)𝑑𝑡=𝛽𝑒^−𝜇1𝜏𝑆(𝑡)𝐼()−(𝑡)𝐼(𝑡)𝑡−𝜏𝑐𝑦𝑚𝑦(𝑡)+𝐼(𝑡)−𝜇₂𝐼(𝑡),𝑑𝑦(𝑡)(𝑡)(𝑡)𝑑𝑡=𝛿𝑦1−ℎ𝑦𝐼(𝑡)(𝑡),+𝜎𝑆 (5.1) where 𝜇₁<𝜇₂ and 0≤𝜎<1. Study on the dynamic behavior of system (5.1) will be obtained in the future.

Acknowledgments

The authors are grateful for the constructive suggestions of the anonymous referees on our original manuscript. This work is supported by the National Natural Science Foundation of China (no. 10771179), Program for Innovative Research Team (in Science and Technology) in University of Henan Province (no. 2010IRTSTHN006) and Program for Key Laboratory of Simulation and Control for Population Ecology in Xinyang Normal University (no. 201004), Natural Science Foundation of the Education Department of Henan Province (no. 2010A110017) and the Cadreman youth teacher of Xinyang Normal University.