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

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.


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 N t , 2 the predator, whose population density is denoted by y t .
We make the following assumptions.
A 1 In the absence of infection and predation, the prey population density grows logistically with carrying capacity K K > 0 and an intrinsic birth rate constant r r > 0 13, 14 , In the presence of disease, the total prey population N t are divided into two distinct classes, namely, susceptible populations, S t , and infected populations, I t 13, 14 .Therefore, at any time t, the total density of prey population is We assume that only susceptible prey S t are capable of reproducing with logistic law 1.1 ; that is, the infected prey I t 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 I t still contributes with S t to population growth toward the carrying capacity 13, 14 .
A 4 We assume that the force of infection at time t is given by βe −μ 1 τ S t I t − τ , where β is the average number of contacts per infective per day, μ 1 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 SI model of the infected prey is

1.3
A 5 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 I t .
From the above assumptions we have the following model:

1.4
All the parameters of system 1.4 are positive.The initial conditions for system 1.4 defined in the Banach space where R 3 { S, I, y ∈ R 3 | S ≥ 0, I ≥ 0, y ≥ 0}.We assume that ϕ i 0 > 0 i 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 S t , I t , y t 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.

Positivity and Boundedness of Solutions
It is important to show positivity and boundedness for the system 1.4 because the variables S, I and y 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 where where x 1 S t , x 2 I t and x 3 y t .
Due to lemma in 27 any solutions of 2.2 with 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 M ≥ 0, and we can get that all solutions of system 1.4 with initial condition enter the region B { S t , I t , y t | 0 ≤ S t , I t , y t ≤ M}.

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

exists and remains positive, and E 2 exists and remains positive if
For equilibrium E 1 , 3.1 reduces to It is easy to see that the equilibrium E 1 is unstable.
For equilibrium E 2 , 3.1 reduces to where

3.5
For τ 0, the transcendental equation 3.3 reduces to 3.4 We can easily get if δ > δ * and β > β 0 obviously, β 0 > β 0 , where δ * ch/ m h 2 and β 0 1/K μ c/ m h .By Routh-Hurwitz criterion, we can get the following theorem.Theorem 3.1.When τ 0, the positive equilibrium E 2 is locally asymptotically stable provided that the conditions In the following, we investigate the existence of purely imaginary roots λ iω ω > 0 to 3.3 .Equation 3.3 takes the form of a third-degree exponential polynomial in λ, which all the coefficients of P and Q 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 E 2 exists.We can easily obtain that τ max 1/μ 1 ln Kβ m h /μ m h c .
e Each positive root ω τ of F ω, τ 0 is continuous and differentiable in τ whenever it exists.
Here, P λ, τ and Q λ, τ are defined as in 3.4 .Let τ ∈ 0, τ max .It is easy to see that This implies that a is satisfied, and b is obviously true because Therefore c follows.
Let F be defined as in d .From we have where

3.13
It is obvious that property d is satisfied, and by Implicit Function Theorem, e is also satisfied.Now let λ iω ω > 0 be a root of 3.3 .Substituting it into 3.3 and separating the real and imaginary parts yields

3.14
From 3.14 it follows that

3.19
The polynomial function F can be written as F ω, τ f ω 2 , τ , where f is a third-degree polynomial, defined by f z, τ z 3 a 1 z 2 a 2 z a 3 .
Clearly, the number of positive real roots of 3.21 depends on the sign of a 1 .When a 1 ≥ 0, 3.21 has only one positive real root.Otherwise, there may exist three positive real roots.
It is easy to know 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 ω τ τ θ τ 2nπ, n 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 .

3.30
We can easily observe that M n 0 < 0.Moreover, for all τ ∈ J, M n τ > M n 1 τ , when n ∈ N. Therefore, if M 0 has no root in J, then the M n functions have no root in J and, if the function M n τ has positive roots τ ∈ J for some n ∈ N, there exists at least one root satisfying dM n dτ τ > 0.

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: r 0.1, K 500, β 0.001, μ 0.03, c 8, δ 0.2, m 150, h 0.5, μ 1 0.04.
Using Maple 10, we compute the coefficients a 1 τ , a 2 τ , a 3 τ for the above values.We can get a 1 τ 0.2829795477 0.004706297728e 0.04τ > 0. Hence the polynomial function f, defined in 3.20 , has only one positive real root.
The function M 0 is drawn for τ ∈ J 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.
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.

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 e −μ 1 τ , where μ 1 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 μ 1 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: where μ 1 < μ 2 and 0 ≤ σ < 1. Study on the dynamic behavior of system 5.1 will be obtained in the future.

Figure 1 :
Figure 1: The graph of the function M 0 τ .Two critical values of τ i.e., τ * and τ * * , for which stability switches can occur, appear.

Figure 2 :
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 S, I and y versus t, respectively, and d is the projected S-I-y phase plane.

Figure 3 :Theorem 3 . 3 .12Figure 4 :
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 S, I and y versus t, respectively, and d is the projected S-I-y phase plane.