Analysis of an Ecoepidemiological Model with Prey Refuges

Abstract

An ecoepidemiological system with prey refuges and disease in prey is proposed. Bilinear incidence and Holling III functional response are used to model the contact process and the predation process, respectively. We will study the stability behavior of the basic system from a local to a global perspective. Permanence of the considered system is also investigated.

1. Introduction

Ecoepidemiology is the branch of biomathematics that understands the dynamics of disease spread on the predator-prey system. Modeling researches on such ecoepidemiological issues have received much attention recently [1–13]. Anderson and May [1] investigated a prey-predator model with prey infection and observed destabilization due to the spread of infectious diseases within animal and plant communities. Chattopadhyay analyzed predator-prey system with disease in the prey [5] and applied the ecoepidemiological study to the Pelicans at risk in the Salton Sea [6]. Bairagi et al. [13] made a comparative study on the role of prey infection in the stability aspects of a predator-prey system with several functional responses. An ecoepidemiological model with prey harvesting and predator switching was investigated by Bhattacharyya and Mukhopadhyay [14]. Kooi et al. [15] studied stabilization and complex dynamics in a predator-prey system with disease in predator. Most of the above-mentioned studies focused on the role of disease in regulating the dynamical consequences of the interacting populations concerned, such as disease-induced stabilization and destabilization of population states [13].

In fact, the dynamical consequences of the predator-prey model can be determined by much ecological effect, such as the Allee effect and prey refuge. Theoretical research and field observations on population dynamics of prey refuges lead to the conclusion that prey refuges have two influences (stabilizing and destabilizing effect) on predator-prey models and prey extinction can be prevented by the addition of prey refuges [13–34]. Here, stabilization (destabilization) of stability refers to cases where an equilibrium point changes from an attractor (a repeller) to a repeller (an attractor) due to increase in the value of a control parameter [17]. Ruxton [16] proposed a continuous-time predator-prey model under the assumption that the rate of prey moving to refuges is proportional to predator density and the results showed that the hiding behavior of prey has a stabilizing effect. The stabilizing effect was also observed in a simple predator-prey system by González-Olivares and Ramos-Jiliberto [17]. Ma et al. [23] formulated a predator-prey model with a class of functional response incorporating the effect of prey refuges and observed the stabilizing and destabilizing effect due to the increases in the prey refuges.

In the present research, we formulate a mathematical model of prey-predator interaction with prey refuges and disease in prey. We mainly study the positivity and boundedness, the stability behavior of the disease-free equilibrium point, and the permanence of the basic model.

2. Model Formulation

The basic model comprises two population subclasses—(i) prey population with density 𝑁(𝑡) and (ii) predators with density 𝑌(𝑡). To formulate our model, we make the following assumptions.(1)The prey population increases logistically with intrinsic growth rate 𝑟 and environmental carrying capacity 𝐾.(2)The prey population is divided into two subclasses—the susceptible prey (𝑆(𝑡)) and the infected prey (𝐼(𝑡)) due to infectious disease. We also assumed that 𝑁(𝑡)=𝑆(𝑡)+𝐼(𝑡) at any instant of time 𝑡.(3)The susceptible prey is capable of reproducing only and the infected prey is removed by death at a rate 𝑑₁.(4)The disease is spread only among the prey population and the disease is not genetically inherited. The infected prey does not become immune.(5)Susceptible prey becomes infected with the simple mass action law 𝛽𝑆𝐼, where 𝛽 measures the force of infection.(6)The predators vanish due to natural death at a constant rate 𝑑₂. They consume susceptible and infected prey following the Holling III functional response with predation coefficients 𝑐₁ and 𝑐₂, respectively. The consumed prey is converted into predator with efficiency 𝑒.(7)It is assumed that there is a quantity 𝛾𝑋 of (the susceptible and the infected) prey population incorporating refuges. 𝛾  (0<𝛾<1) denotes that a constant proportion of (the susceptible and the infected) prey use refuges. All the above-mentioned parameters are assumed to be positive.

With the previously mentioned assumptions, the generalized predator-prey system with prey refuges and disease in prey can be represented by the following equations: 𝑑𝑆𝑑𝑡=𝑟𝑆1−𝑆+𝐼𝐾𝑐−𝛽𝑆𝐼−¹1−𝛾²𝑆²𝑌𝑚+1−𝛾²𝑆²,𝑑𝐼𝑐𝑑𝑡=𝛽𝑆𝐼−²1−𝛾²𝐼²𝑌𝑚+1−𝛾²𝐼²−𝑑₁𝐼,𝑑𝑌=𝑑𝑡𝑒𝑐₁1−𝛾²𝑆²𝑌𝑚+1−𝛾²𝑆²+𝑒𝑐₂1−𝛾²𝐼²𝑌𝑚+1−𝛾²𝐼²−𝑑₂𝑌. (2.1)

Defining 𝐺={(𝑆,𝐼,𝑌)∈𝑅³₊∣𝑆+𝐼≤𝐾,𝑆≥0,𝐼≥0,𝑌≥0}, it is easy to show that the set 𝐺 is the positively invariable set of system (2.1).

3. The Positivity and Boundedness

Theorem 3.1. All solutions of system (2.1) initiating 𝑅³₊ are positive and ultimately bounded.

Proof. Let (𝑆(𝑡),𝐼(𝑡),𝑌(𝑡)) be one of the solutions of system (2.1).
Integrating (2.1) with initial conditions (𝑆₀,𝐼₀,𝑌₀), we have 𝑆(𝑡)=𝑆₀exp^𝑡₀𝑟𝑆(𝜉)(𝜉)1−+𝐼𝐾(𝜉)−𝑐−𝛽𝐼¹1−𝛾²𝑆(𝜉)𝑌(𝜉)𝑚+1−𝛾²𝑆²(𝜉)𝐼(𝑡)𝑑𝜉>0,=𝐼₀exp^𝑡₀(𝜉)−𝑐𝛽𝑆²1−𝛾²𝐼(𝜉)𝑌(𝜉)𝑚+1−𝛾²𝐼²(𝜉)−𝑑₁𝑌(𝑡)𝑑𝜉>0,=𝑌₀exp^𝑡₀𝑒𝑐₁1−𝛾²𝑆²(𝜉)𝑚+1−𝛾²𝑆²(𝜉)+𝑒𝑐₂1−𝛾²𝐼²(𝜉)𝑚1−𝛾²𝐼²(𝜉)−𝑑₂𝑑𝜉>0. (3.1)
Hence all solutions starting in 𝑅³₊ remain in 𝑅³₊ for all 𝑡≥0.

Next, we will prove the boundedness of the solutions.

Because 𝑑𝑆/𝑑𝑡≤𝑟𝑆(1−𝑆/𝐾), then we have lim_𝑡→+∞sup𝑆(𝑡)≤𝐾. (3.2)

Let 𝑊=𝑒𝑆+𝑒𝐼+𝑌, then we obtain that 𝑑𝑊𝑆𝑑𝑡=𝑟𝑒𝑆1−𝐾−𝑟𝑒𝐾𝑆𝐼−𝑒𝑑₁𝐼−𝑑₂𝑌≤𝑟𝑒𝑆−𝑒𝑑₁𝐼−𝑑₂𝑌=(𝑟+1)𝑒𝑆−𝑒𝑆+𝑑₁𝑒𝐼+𝑑₂𝑌≤𝑒(𝑟+1)𝐾−𝛿𝑊, (3.3) where 𝛿=min{1,𝑑₁,𝑑₂}.

Hence, we have 𝑑𝑊()𝑑𝑡+𝛿𝑊≤𝑒𝑟+1𝐾. (3.4) That is ()≤𝑒()𝐾0≤𝑊𝑆,𝐼,𝑌𝑟+1𝛿+𝑊(𝑆(0),𝐼(0),𝑌(0))𝑒^𝛿𝑡,𝑒()𝐾0≤𝑊≤𝑟+1𝛿(𝑡→+∞). (3.5)

Thus, all curves of system (2.1) will enter the following region: )𝑒()𝐾𝐵=(𝑆,𝐼,𝑌∶0≤𝑊≤𝑟+1𝛿+𝜀,∀𝜀>0. (3.6)

4. The Equilibrium Point

All equilibrium points of system (2.1) can be obtained by solving the following equations: 𝑟𝑆1−𝑆+𝐼𝐾𝑐−𝛽𝑆𝐼−¹1−𝛾²𝑆²𝑌𝑚+1−𝛾²𝑆²𝑐=0,𝛽𝑆𝐼−²1−𝛾²𝐼²𝑌𝑚+1−𝛾²𝐼²−𝑑₁𝐼=0,𝑒𝑐₁1−𝛾²𝑆²𝑌𝑚+1−𝛾²𝑆²+𝑒𝑐₂1−𝛾²𝐼²𝑌𝑚+1−𝛾²𝐼²−𝑑₂𝑌=0, (4.1) These points are as follows:(1)the trivial equilibrium point 𝐸₀(0,0,0), (2)the equilibrium point 𝐸_𝐾(𝐾,0,0), (3)the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0), (4)the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝑌^∗), where 𝑑𝑆=¹𝛽=̇𝑆_𝑀,𝑟𝐼=𝛽𝐾−𝑑₁𝛽𝑟+𝛽𝐾=̇𝐼_𝑀,𝑆^∗=11−𝛾𝑑₂𝑚𝑒𝑐₁−𝑑₂,𝑌^∗=𝑒𝑟𝑆^∗𝑑₂𝑆1−^∗𝐾. (4.2)

Let 𝑅₁=(𝑒𝑐₁(1−𝛾)²𝐾²/(𝑚+(1−𝛾)²𝐾²))/𝑑₂. It is clear to show that the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝑌^∗) has its ecological meaning when 𝑅₁>1.

Let 𝑅₀=𝛽𝐾/𝑑₁, the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is nonnegative when 𝑅₀>1.

5. The Stability Property

In this section, we will study the local and global stability of the equilibrium points of system (2.1).

Theorem 5.1. Let 𝑅₁=((𝑒𝑐₁(1−𝛾)²𝑆²_𝑀/(𝑚+(1−𝛾)²𝑆²_𝑀))+(𝑒𝑐₂(1−𝛾)²𝐼²_𝑀/(𝑚+(1−𝛾)²𝐼²_𝑀)))/𝑑₂, 𝑅₂=max(𝑅₁,𝑅₁), then one has the following. (1)The trivial equilibrium point 𝐸₀(0,0,0) is always unstable. (2)If 𝑅₀<1, 𝑅₁<1, then the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is globally asymptotically stable. If 𝑅₀>1 or 𝑅₁>1, the equilibrium point 𝐸_𝐾(𝐾,0,0) is unstable. (3)If 𝑅₀>1 and 𝑅₁<1, the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is locally asymptotically stable. If 𝑅₂<1, the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is locally asymptotically stable. If 𝑅₁>1, the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is unstable.

Proof. The Jacobian matrix of system (2.1) at the trivial equilibrium point 𝐸₀(0,0,0) is 𝑟000−𝑑₁000−𝑑₂. (5.1)

Clearly, the trivial equilibrium point 𝐸₀(0,0,0) is unstable.

The Jacobian matrix of system (2.1) at the equilibrium point 𝐸_𝐾(𝐾,0,0) is 𝑐−𝑟−𝑟−𝛽𝐾−¹1−𝛾²𝐾²𝑚+1−𝛾²𝐾²0𝛽𝐾−𝑑₁000𝑒𝑐₁1−𝛾²𝐾²𝑚+1−𝛾²𝐾²−𝑑₂. (5.2)

According to the theorem about the local stability, the local stability of the equilibrium point 𝐸_𝐾(𝐾,0,0) is determined only by the sign of (𝑒𝑐₂(1−𝛾)²𝐾²/(𝑚+(1−𝛾)²𝐾²))−𝑑₂ and 𝛽𝐾−𝑑₁.

Therefore, if 𝑅₀<1 and 𝑅₁<1, the equilibrium point 𝐸_𝐾(𝐾,0,0) is locally asymptotically stable.

Next, we will prove the global stability of the equilibrium point 𝐸_𝐾(𝐾,0,0).

Defining the Lyapunov function 𝑉=𝑒𝐼+𝑌, then we obtain that 𝑑𝑉𝑑𝑡=𝑒𝛽𝑆−𝑑₁𝐼+𝑒𝑐₁1−𝛾²𝑆²𝑚+1−𝛾²𝑆²−𝑑₂𝑌≤𝑒𝛽𝐾−𝑑₁𝐼+𝑒𝑐₁1−𝛾²𝐾²𝑚+1−𝛾²𝐾²−𝑑₂𝑌. (5.3)

Hence, if 𝑅₀<1 and 𝑅₁<1, 𝑑𝑉/𝑑𝑡≤0.

Having 𝐸={(𝑆,𝐼,𝑌)∈𝐺∣𝑑𝑉/𝑑𝑡=0}={𝐼=0,𝑌=0}, the maximum invariable set of system (2.1) is 𝑀=𝐸={𝐼=0,𝑌=0}.

According to the LaSalle invariable set theorem, lim_𝑡→+∞𝐼(𝑡)=0, lim_𝑡→+∞𝑌(𝑡)=0.

Thus, the limit equation of system (2.1) is 𝑑𝑆𝑆𝑑𝑡=𝑟𝑆1−𝐾. (5.4)

Clearly, the equilibrium point 𝑆=𝐾 is globally asymptotically stable.

According to the limit system theorem, if 𝑅₀<1 and 𝑅₁<1, the equilibrium point 𝐸_𝐾(𝐾,0,0) is globally asymptotically stable.

Setting 𝑅₀>1, the Jacobian matrix of system (2.1) at the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is −𝑟𝑆𝐾−𝑟𝑆𝐾−𝛽𝑐𝑆−¹1−𝛾²𝑆²𝑚+1−𝛾²𝑆²𝛽𝑐𝐼0−²1−𝛾²𝐼²𝑚+1−𝛾²𝐼²00𝑒𝑐₁1−𝛾²𝑆²𝑚+1−𝛾²𝑆²+𝑒𝑐₂1−𝛾²𝐼²𝑚+1−𝛾²𝐼²−𝑑₂. (5.5)

The characteristic equation of system (2.1) at the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is (𝜆−𝐴)𝜆²−𝐵𝜆+𝐶=0, (5.6) where 𝑐𝐴=¹1−𝛾²𝑆²_𝑀𝑚+1−𝛾²𝑆²_𝑀+𝑒𝑐₂1−𝛾²𝐼²_𝑀𝑚+1−𝛾²𝐼²_𝑀−𝑑₂,𝐵=−2𝑆_𝑀𝐾<0,𝐶=𝛽𝐼_𝑀𝑆_𝑀𝐾+𝛽𝑆_𝑀>0. (5.7)

According to the Routh-Hurwitz rule, the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is locally asymptotically stable when 𝑅₁<1.

Next, we will prove the global stability of the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0).

Defining the Lyapunov function 𝑉=𝑌, then we have 𝑑𝑉=𝑑𝑡𝑒𝑐₁1−𝛾²𝑆²𝑚+1−𝛾²𝑆²+𝑒𝑐₂1−𝛾²𝐼²𝑚+1−𝛾²𝐼²−𝑑₂𝑌. (5.8)

Let Φ(𝑆,𝐼)=(𝑐₁(1−𝛾)²𝑆²/(𝑚+(1−𝛾)²𝑆²))+(𝑒𝑐₂(1−𝛾)²𝐼²/(𝑚+(1−𝛾)²𝐼²)) in which (𝑆,𝐼)∈𝐷={(𝑆,𝐼)∣𝑆+𝐼≤𝐾,𝑆≥0,𝐼≥0}.

Thus Φ^′(𝑡)=2𝑐₁𝑒𝑚1−𝛾²𝑆𝑚+1−𝛾²𝑆²²𝑑𝑆+𝑑𝑡2𝑐₂𝑒𝑚1−𝛾²𝐼𝑚+1−𝛾²𝐼²²𝑑𝐼.𝑑𝑡 (5.9)

Clearly, (2𝑐₁𝑒𝑚(1−𝛾)²𝑆/(𝑚+(1−𝛾)²𝑆²)²)>0 and (2𝑐₂𝑒𝑚(1−𝛾)²𝐼/(𝑚+(1−𝛾)²𝐼²)²)>0 on the set 𝐷^′=𝐷−{𝑆=0∣𝐼=0}.

Thus, Φ^′(𝑡)=0 on the set 𝐷^′, if and only if 𝑑𝑆/𝑑𝑡=0 and 𝑑𝐼/𝑑𝑡=0, then 𝑆=𝑆_𝑀, 𝐼=𝐼_𝑀.

If 𝑆=𝐼=0, then Φ(𝑆,𝐼)=0. If 𝑆=𝑆_𝑀, then 𝐼=𝐼_𝑀. If 𝑆<𝑆_𝑀, then 𝑑𝐼/𝑑𝑡<0; that is, 𝐼<𝐼_𝑀.

Therefore, the maximum value of the function Φ(𝑡) is obtained at the point (𝐾,0) or (𝑆_𝑀,𝐼_𝑀).

Hence, if 𝑅₂<1, 𝑑𝑉/𝑑𝑡≤0; that is, 𝐸={(𝑆,𝐼,𝑌)∈𝐺∣𝑑𝑉/𝑑𝑡=0}={𝑌=0}.

According to the LaSalle invariable set theorem, lim_𝑡→+∞𝑌(𝑡)=0. Thus, the limit equation of system (2.1) is 𝑑𝑆𝑆𝑑𝑡=𝑟𝑆1−𝐾−𝛽𝑆𝐼,𝑑𝐼𝑑𝑡=𝛽𝑆𝐼−𝑑₁𝐼. (5.10)

According to the results of the appendix section, if 𝑅₀>1, the equilibrium point (𝑆_𝑀,𝐼_𝑀) is globally asymptotically stable.

According to the limit system theorem, if 𝑅₀>1 and 𝑅₂<1, the predator-extinction equilibrium point 𝐸(𝑆,𝐼,0) is globally asymptotically stable.

Theorem 5.2. Let 𝑅^∗₁=(𝛽/𝑑₁(1−𝛾))(𝑑₂𝑚/(𝑒𝑐₁−𝑑₂)), 𝑅^∗₂=(((2𝑑₂−𝑒𝑐₁)(1−𝛾)𝐾)/2𝑑₂)(𝑒𝑐₁−𝑑₂)/𝑑₂𝑚. If 𝑅₁>1, then the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝐼^∗) is nonnegative, (1)if 𝑅^∗₁<1 and 𝑅^∗₂<1, then the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝐼^∗) is locally asymptotically stable; (2)if 𝑅₀<1 and 𝑅^∗₂<1, then the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝐼^∗) is globally asymptotically stable; (3)if 𝑅^∗₁>1 or 𝑅^∗₂>1, then the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝐼^∗) is unstable.

Proof. Assuming 𝑅₁>1, the Jacobian matrix of system (2.1) at the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝑌^∗) is 𝑎₁₁−𝛽𝑆^∗−𝑟𝑆^∗𝐾−𝑐₁1−𝛾²𝑆^∗2𝑚+1−𝛾²(𝑆^∗)²0𝑎₂₂02𝑒𝑐₁𝑚1−𝛾²𝑆^∗𝑌^∗𝑚+1−𝛾²(𝑆^∗)²²00, (5.11) where 𝑎₁₁=𝑟1−2𝑆^∗𝐾−2𝑐₁𝑚1−𝛾²𝑆^∗𝑌^∗𝑚+1−𝛾²(𝑆^∗)²²,𝑎₂₂=𝛽𝑆^∗−𝑑₁. (5.12)

The characteristic equation of system (2.1) at the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝑌^∗) is 𝜆−𝑎₂₂𝜆²−𝑎₁₁𝜆+2𝑐²₁𝑒𝑚1−𝛾⁴𝑆^∗3𝑌^∗𝑚+1−𝛾²(𝑆^∗)²³=0. (5.13)

Clearly, if 𝑅^∗₁<1, 𝑎₂₂=𝛽𝑆^∗−𝑑₁<0.

Again, we have 𝑎₁₁=𝑟1−2𝑆^∗𝐾−2𝑐₁𝑚1−𝛾²𝑆^∗𝑌^∗𝑚+1−𝛾²(𝑆^∗)²²=𝑟1−2𝑆^∗𝐾𝑑−2𝑟1−²𝑒𝑐₁𝑆1−^∗𝐾=𝑟1−2𝑆^∗𝐾−2+2𝑆^∗𝐾+2𝑑₂𝑒𝑐₁−2𝑑₂𝑒𝑐₁𝑆^∗𝐾=𝑟𝑒𝑐₁2𝑑₂−𝑒𝑐₁−2𝑑₂1−𝛾𝐾𝑑₂𝑚𝑒𝑐₁−𝑑₂. (5.14)

Hence, if 𝑅^∗₂<1, then 𝑎₁₁<0.

According to the Routh-Hurwitz rule, if 𝑅^∗₁<1 and 𝑅^∗₂<1, then the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝑌^∗) is locally asymptotically stable.

Next, we will study the global stability of the disease-free equilibrium point 𝐸^∗(𝑆^∗,0,𝑌^∗).

Defining the Lyapunov function 𝑉=𝐼, then we obtain 𝑑𝑉=𝑐𝑑𝑡𝛽𝑆−²1−𝛾²𝑌𝑚+1−𝛾²𝐼−𝑑₁𝐼≤𝛽𝐾−𝑑₁𝐼. (5.15)

Thus, if 𝑅₀<1, 𝑑𝑉/𝑑𝑡≤0 and 𝑑𝑉/𝑑𝑡=0 if and only if 𝐼=0; that is, 𝐸={(𝑆,𝐼,𝑌)∈𝐺∣𝑑𝑉/𝑑𝑡=0}={𝐼=0}.

According to the LaSalle invariable set theorem, 𝐼(𝑡)→0 when 𝑡→+∞. The limit system of system (2.1) is 𝑑𝑆𝑆𝑑𝑡=𝑟𝑆1−𝐾−𝑐₁1−𝛾²𝑆²𝑌𝑚+1−𝛾²𝑆²,𝑑𝑌=𝑑𝑡𝑒𝑐₁1−𝛾²𝑆²𝑌𝑚+1−𝛾²𝑆²−𝑑₂𝑌. (5.16)

Clearly, the equilibrium points of the system (5.16) are 𝐸^′₀(0,0),𝐸^′_𝐾(𝐾,0),𝐸^′_∗𝑆^∗,𝑌^∗ (5.17) in which 𝑆^∗ and 𝑌^∗ are similar as the equilibria expression of system (2.1).

It is easy to show that the equilibrium point 𝐸^′_∗ has its ecological meaning when 𝑅₁>1.

According to the Routh-Hurwitz rule, the equilibrium point 𝐸^′₀ is unstable. If 𝑅₁>1, then 𝐸^′_𝐾(𝐾,0) is unstable.

Again, the Jacobian matrix of system (5.16) at the equilibrium point 𝐸^′_∗ is 𝑎₁₁−𝑐₁1−𝛾²𝑆^∗2𝑚+1−𝛾²(𝑆^∗)²2𝑒𝑐₁𝑚1−𝛾²𝑆^∗𝑌^∗𝑚+1−𝛾²(𝑆^∗)²²0, (5.18) where 𝑎₁₁=𝑟1−2𝑆^∗𝐾−2𝑐₁𝑚1−𝛾²𝑆^∗𝑌^∗𝑚+1−𝛾²(𝑆^∗)²². (5.19)

The characteristic equation of system (5.16) at the equilibrium point 𝐸^′_∗ is 𝜆²−𝑎₁₁𝜆+2𝑐²₁𝑒𝑚1−𝛾⁴𝑆^∗3𝑌^∗𝑚+1−𝛾²(𝑆^∗)²³=0. (5.20)

According to the above study, if 𝑅^∗₂<1, then 𝑎₁₁<0.

Hence, if 𝑅^∗₂<1, then the equilibrium point 𝐸^′_∗ is locally asymptotically stable in the region 𝐷 by the Routh-Hurwitz rule.

It is easy to note that the globally asymptotically stability of the equilibrium point 𝐸^′_∗ implies that there is no close orbit in the region 𝐷 for the considered system.

Let 𝑆=𝑑₂𝑚𝑒𝑐₁−𝑑₂1𝑆,𝑌=𝑐₁𝑑₂𝑚𝑒𝑐₁−𝑑₂𝑑𝑌,𝑑𝑡=²1−𝛾²𝑆²+𝑒𝑐₁−𝑑₂𝑑₂𝑒𝑐₁−𝑑₂𝑑𝑡,̃ (5.21) and rewrite 𝑆, 𝑌, and 𝑡̃ into 𝑆, 𝑌, and 𝑡, then system (5.16) becomes as follows: 𝑑𝑆𝐴𝑑𝑡=𝑆₀+𝐴₁𝑆+𝐴₂𝑆²+𝐴₃𝑆³−1−𝛾²𝑆²𝑌≡𝑃(𝑆,𝑌),𝑑𝑌=𝑑𝑡1−𝛾²𝑆²𝑌−𝑌≡𝑄(𝑆,𝑌), (5.22) where 𝐴₀=𝑟/𝑑₂>0, 𝐴₁=−𝑟/(𝑑₂𝐾)𝑑₂𝑚/(𝑒𝑐₁−𝑑₂)<0, 𝐴₂=(1−𝛾)²𝑟/(𝑒𝑐₁−𝑑₂)>0, 𝐴₃=−(1−𝛾)²𝑟/(𝑒𝑐₁−𝑑₂)𝑑₂𝑚/(𝑒𝑐₁−𝑑₂)<0.

Thus, the positive equilibrium point 𝐸^′_∗(𝑆^∗,𝑌^∗) of system (5.16) becomes the positive equilibrium point 𝐸_𝑃(𝑆_𝑃,𝑌_𝑃) of system (5.22), where 𝑆_𝑃=1,𝑌1−𝛾_𝑃=11−𝛾𝑒𝑐₁𝑑₂𝑒𝑐₁−𝑑₂𝑟𝑟−1−𝛾𝐾𝑑₂𝑚𝑒𝑐₁−𝑑₂=𝑟𝑒𝑐₁𝑑₂𝑒𝑐₁−𝑑₂11−𝛾1−1−𝛾𝐾𝑑₂𝑚𝑒𝑐₁−𝑑₂. (5.23)

Considering the Dulac function 𝐵(𝑆,𝑌)=𝑆⁻²𝑌^𝑛−1, then we have 𝜕(𝐵𝑃)+𝜕()𝜕𝑆𝐵𝑄𝜕𝑌=𝑆⁻²𝑌^𝑛−12𝐴₃𝑆³+𝐴₂+𝑛1−𝛾²𝑆²−𝐴₀+𝑛≡𝑆⁻²𝑌^𝑛−1Φ(𝑆,𝑛). (5.24)

In order to prove the global stability, we will prove only that there exists a real number 𝑛 such that Φ(𝑆,𝑛)≤0.

Clearly, if 𝑆=0 or 𝑆=−(𝐴₂+𝑛(1−𝛾)²)/3𝐴₃, then 2𝑆[3𝐴₃𝑆+𝐴₂+𝑛(1−𝛾)²]=0 since 𝜕Φ(𝑆,𝑛)/𝜕𝑆=2𝑆[3𝐴₃𝑆+𝐴₂+𝑛(1−𝛾)²].

Let −(𝐴₂+𝑛(1−𝛾)²)/3𝐴₃>0  (𝐴₂>0, 𝐴₃<0), that is, 𝐴₂+𝑛(1−𝛾)²>0, then we obtain that if 0<𝑆<−(𝐴₂+𝑛(1−𝛾)²)/3𝐴₃, then 𝜕Φ(𝑆,𝑛)/𝜕𝑆>0. Otherwise, 𝜕Φ(𝑆,𝑛)/𝜕𝑆<0.

Therefore, the function Φ(𝑆,𝑛) has the maximum value at the point 𝑆=−(𝐴₂+𝑛(1−𝛾)²)/3𝐴₃ and Φ(0,𝑛)=−(𝐴₀+𝑛)<0.

Hence, there exists only one real number 𝑛, such that 𝐴₂+𝑛(1−𝛾)²>0 and max_0<𝑆<+∞Φ(𝑆,𝑛)=Φ(𝑆,𝑛)∣_{𝑆=−(𝐴2+𝑛(1−𝛾)}²_)/3𝐴3=𝐴₂+𝑛1−𝛾²³27𝐴²₃−𝐴₂+𝑛1−𝛾²1−𝛾²+𝐴₂1−𝛾²−𝐴₀=0. (5.25)

Thus, we will prove only that there exists 𝐴₂+𝑛(1−𝛾)²>0, such that 𝐴₂+𝑛1−𝛾²³+𝑀𝐴₂+𝑛1−𝛾²+𝑁=0, (5.26) where 𝑀=−27𝐴²₃/(1−𝛾)²<0, 𝑁=27𝐴²₃[𝐴₂/(1−𝛾)²−𝐴₀]>0.

The discriminant of the cubic equation is 𝑁Δ=²4+𝑀²=2727𝐴⁴₃4𝐴₂1−𝛾²−𝐴₀²−4𝐴²₃1−𝛾⁶=𝐴₂𝑆²_𝑝−𝐴₀+2𝐴₃𝑆³_𝑝𝐴₂𝑆²_𝑝−𝐴₀−2𝐴₃𝑆³_𝑝. (5.27)

It is easy to show that if 𝑑₂<𝑒𝑐₁<2𝑑₂, then 𝐴₂𝑆²_𝑝−𝐴₀=𝑟(2𝑑₂−𝑒𝑐₁)/𝑑₂(𝑒𝑐₁−𝑑₂)>0.

Again, if 𝑅₁>1 and 𝑅^∗₂<1, then 𝐴₂𝑆²_𝑝−𝐴₀−2𝐴₃𝑆³_𝑝>0 and 2𝐴₃𝑆³_𝑝+𝐴₂𝑆²_𝑝−𝐴₀≤0. If 2𝐴₃𝑆³_𝑝+𝐴₂𝑆²_𝑝−𝐴₀=0, then Δ=0.

According to Shengjin's distinguishing means, the cubic equation has one negative real root and two positively real roots.

If 2𝐴₃𝑆³_𝑝+𝐴₂𝑆²_𝑝−𝐴₀<0, the cubic equation has three roots which are not equal.

According to the Descartes rule of signs, the cubic equation has two positively real roots and one negatively real root at most.

Therefore, the cubic equation has at least one positively real root. That is to say, there exists a number 𝑛, such that 𝐴₂+𝑛(1−𝛾)²>0.

Furthermore, we obtain that 𝜕(𝐵𝑃)+𝜕()𝜕𝑆𝐵𝑄𝜕𝑌≤0. (5.28)

According to the Bendixson-Dulac theorem, there does not exist the limit cycle for the limit system.

Hence, the equilibrium point 𝐸^′_∗ is globally asymptotically stable.

Therefore, if 𝑅₀<1 and 𝑅^∗₂<1, then the equilibrium point 𝐸^′_∗ is globally asymptotically stable according to the limit system theorem.

6. Permanence

Theorem 6.1. If 𝑅₀>1, 𝑅₁>1, 𝑅₁>1, and 𝑅^∗₁>1, then system (2.1) is permanent.

Proof. Considering the average Lyapunov function 𝑉(𝑆,𝐼,𝑌)=𝑆^𝛼1𝐼^𝛼2𝑌^𝛼3, where 𝛼_𝑖(𝑖=1,2,3) is positive, then in the region 𝑅³₊, we have 𝑉̇𝑉=𝛼₁𝑟1−𝑆+𝐼𝐾𝑐−𝛽𝐼−¹1−𝛾²𝑆𝑌𝑚+1−𝛾²𝑆²+𝛼₂𝑐𝛽𝑆−²1−𝛾²𝐼𝑌𝑚+1−𝛾²𝐼²−𝑑₁+𝛼₃𝑒𝑐₁1−𝛾²𝑆²𝑚+1−𝛾²𝑆²+𝑒𝑐₂1−𝛾²𝐼²𝑚+1−𝛾²𝐼²−𝑑₂≡Ψ(𝑆,𝐼,𝑌). (6.1)

In order to prove the permanence of system (2.1), we only indicate the following results: the function Ψ(𝑆,𝐼,𝑌)>0 for all boundary equilibrium points.

Let 𝛼₁>(𝑑₁𝛼₂+𝑑₂𝛼₃)/𝑟, then Ψ𝐸₀=𝛼₁𝑟−𝛼₂𝑑₁−𝛼₃𝑑₂Ψ𝐸>0,_𝐾=𝛼₂𝛽𝐾−𝑑₁+𝛼₃𝑒𝑐₁1−𝛾²𝐾²𝑚+1−𝛾²𝐾²−𝑑₂,Ψ𝐸=𝛼₃𝑒𝑐₁1−𝛾²𝑆²_𝑀𝑚+1−𝛾²𝑆²_𝑀+𝑒𝑐₂1−𝛾²𝐼²_𝑀𝑚+1−𝛾²𝐼²_𝑀−𝑑₂,Ψ𝐸^∗=𝛼₂𝛽𝑆^∗−𝑑₁. (6.2)

Hence, if 𝑅₀>1 and 𝑅₁>1, then Ψ(𝐸_𝐾)>0. If 𝑅₂>1, then Ψ(𝐸)>0. If 𝑅^∗₁>1, then Ψ(𝐸^∗)>0.

Therefore, system (2.1) is permanent when 𝑅₀>1, 𝑅₁>1, 𝑅₁>1, and 𝑅^∗₁>1 by the average Lyapunov function theorem [34].

By simple computation, 𝑅₁>1(<1) is equivalent to (1−𝛾)𝐾(𝑒𝑐₁−𝑑₂)/𝑑₂𝑚>1(<1).

Hence, if 𝑅₁>1, then 𝑅^∗₁>1 implies 𝑅₀>1.

Appendix

Considering the following 𝑆𝐼 model 𝑑𝑆𝑑𝑡=𝑟𝑆1−𝑆+𝐼𝐾−𝛽𝑆𝐼,𝑑𝐼𝑑𝑡=𝛽𝑆𝐼−𝑑₁𝐼, (A.1) where 𝑆(𝑡) and 𝐼(𝑡) are density of susceptible and infected prey population at time 𝑡. The parameters 𝑟, 𝐾, 𝛽, and 𝑑₁ are all positive.

The positive equilibrium point of this system is 𝐸_𝑀(𝑆_𝑀,𝐼_𝑀) in which 𝑆_𝑀=𝑑₁𝛽,𝐼_𝑀=𝑟𝛽𝐾−𝑑₁𝛽𝑟+𝛽𝐾. (A.2)

Let 𝑅₀=𝛽𝐾/𝑑₁.

Clearly, if 𝑅₀=𝛽𝐾/𝑑₁>1, then the equilibrium point 𝐸_𝑀(𝑆_𝑀,𝐼_𝑀) is positive.

Theorem A.1. Let 𝑅₀=𝛽𝐾/𝑑₁. If 𝑅₀>1, then the positive equilibrium point 𝐸_𝑀(𝑆_𝑀,𝐼_𝑀) is globally asymptotically stable.

Proof. The Jacobian matrix of system (A.1) at the positive equilibrium point 𝐸_𝑀(𝑆_𝑀,𝐼_𝑀) is −𝑟𝑆_𝑀𝐾−𝑟𝑆_𝑀𝐾−𝛽𝑆_𝑀0𝛽𝑆_𝑀−𝑑₁. (A.3)

Clearly, if 𝑅₀>1, then the positive equilibrium point 𝐸_𝑀(𝑆_𝑀,𝐼_𝑀) is locally asymptotically stable.

Next, defining the Dulac function 𝐵(𝑆,𝐼)=1/𝑆𝐼, then we have 𝜕(𝐵𝑃)+𝜕()𝜕𝑆𝐵𝑄𝑟𝜕𝐼=−𝐾𝐼<0, (A.4) where 𝑃(𝑆,𝐼)=𝑟𝑆(1−((𝑆+𝐼)/𝐾))−𝛽𝑆𝐼, 𝑄(𝑆,𝐼)=𝛽𝑆𝐼−𝑑₁𝐼.

Therefore, there does not exist closed curve in the region 𝐷 according to the Dulac theorem.

This implies that the positive equilibrium point 𝐸_𝑀(𝑆_𝑀,𝐼_𝑀) is globally asymptotically stable.

Acknowledgments

This work was supported by the Fundamental Research Funds for the Central Universities (no. lzujbky-2011-48) and the National Natural Science Foundation of China (no. 11126183, 30970478, 30970491, 31100306).