Abstract

A delayed Leslie-Gower predator-prey model with
nonmonotonic functional response is studied. The existence and local stability of the
positive equilibrium of the system with or without delay are completely determined
in the parameter plane. Using the method of upper and lower solutions and monotone
iterative scheme, a sufficient condition independent of delay for the global stability
of the positive equilibrium is obtained. Hopf bifurcations induced by the ratio of the
intrinsic growth rates of the predator and prey and by delay, respectively, are found. Employing the normal form theory, the direction and stability of Hopf bifurcations
can be explicitly determined by the parameters of the system. Some numerical
simulations are given to support and extend our theoretical results. Two limit cycles
enclosing an equilibrium, one limit cycle enclosing three equilibria and different types
of heteroclinic orbits such as connecting two equilibria and connecting a limit cycle
and an equilibrium are also found by using analytic and numerical methods.

1. Introduction

Predation is one of the most fundamental interactions between different species and is universal in nature. Recent studies have shown that predator-prey relationships (not competitive or mutualistic relationships) provide the necessary stability for almost infinite numbers of species to exist in ecosystems and make the rich biodiversity of complex ecosystems possible [1]. The mathematical model and the related dynamics are very useful for qualitatively understanding the interaction between predators and preys. Thus, the predator-prey dynamics has been extensively studied and continue to be of interest to both applied mathematicians and ecologists.

The pioneering predator-prey model was initially proposed by Lotka [2] and Volterra [3] (known as Lotka-Volterra predator-prey model) and successfully captured the oscillations in populations of a predator and its prey. Holling extended this model in [4–6] where various types of functional responses were proposed to better understand the components of predator-prey interactions. The Lotka-Volterra predator-prey model with Holling functional response has been a cornerstone in the field of mathematical biology, and much literature has been devoted to study variants of these equations.

Another important type of well-known predator-prey models was initially proposed and studied by Leslie and Gower [7–9]. Compared with Lotka-Volterra model, the novel feature of Leslie-Gower model is that both interacting species are assumed to grow according to the logistic law, the environmental carrying capacity of the predator is not a constant but a function of the available prey quantity, that there are upper limits to the rates of increase of both prey and predator, which are not recognized in the Lotka-Volterra model. Later, Holling functional response was introduced to Leslie-Gower model by May [10] and the corresponding model is described by the following ordinary differential equations:
where and denote the prey and predator densities, respectively, represents the Holling type functional response [6], the positive constants and are the intrinsic growth rates of the prey and predators, respectively, the carrying capacity for the prey is a positive constant , is a measure of the food quality that the prey provides for conversion into predator births, and then is considered as the carrying capacity for the predator. The model (1) is also known as the Holling-Tanner predator-prey model [11–13]. When the functional response is Holling types I, II, and III functions respectively, it has been proven that the dynamics of model (1) are quite interesting and the global stability of the positive equilibrium, the existence of limit cycles, and so forth have been investigated by many researchers (see, e.g., [12–16] and references therein). Holling type I, II, and III functional responses are similar in that each per capita rate of predation increases with prey density to a maximum [4] and it is shown that the model (1) with the type I, II, and III functional responses exhibits qualitatively similar bifurcation and stability behavior [17]. Holling type IV functional response described by Holling in [18] incorporates prey interference with predation. The main characteristic of Holling type IV functional response is that the per capita predation rate increases with prey density to a maximum at a critical prey density beyond which it decreases, while for other Holling functional responses, the per capita predation rate always increases with prey density. Collings [17] has also numerically shown that the dynamics of the model (1) with Holling type IV functional response is very different from that of the model (1) with other Holling functional responses at higher levels of prey interference. The qualitative analysis for the model (1) with a simplified Holling type IV functional response, proposed by Sokol and Howell [19] as follows:
has been studied by Li and Xiao [20]. They mainly studied the case when the model has two equilibria: one is an unstable multiple focus with multiplicity one and the other is a cusp of codimension 2 and the dynamics of the model near the small neighborhoods of these two equilibria.

On the other hand, the introduction of time delay into the population model is more realistic to model the interaction between the predator and prey populations and the population models with time delay are of current research interest in mathematical biology [21, 22]. There is extensive literature about the effects of delay on the dynamics of predator-prey models. We refer to a recent survey paper by Ruan [22] and the references cited therein for studies on delayed Lotka-Volterra predator-prey models with Holling functional response. The effects of delay on the dynamics of the Leslie-Gower predator-prey model with and its modified version have been recently studied by many authors (see, e.g., [23–26] and references therein). However, most of them considered the case of the monotonic functional response, under which the corresponding system has only one positive equilibrium. When the nonmonotonic functional response (Holling type IV) is introduced to the Leslie-Gower predator-prey model, there is a possibility of the existence of one, two or three positive equilibria and the dynamics is much more complex [17, 20]. In [27], the authors considered the following delayed Leslie-Gower model with nonmonotonic functional response:
where is incorporated because of the negative feedback of the predator density. Under the assumption that the positive equilibrium exists and it is stable for , the local stability and delay-induced Hopf bifurcations have been investigated. The supercritical stable Hopf bifurcations were also found by normal form theory and numerical simulation.

In this paper, following the work of Li and Xiao [20] and Lian and Xu [27], we study how the hunting delay affects the dynamics of the Leslie-Gower predator-prey model nonmonotonic functional response. Introducing the hunting delay into the Leslie-Gower predator-prey model with nonmonotonic functional response, we have the following delay differential equations:
where the biologic meaning of the parameters , , is the same as in model (1), is the maximal predator per capita consumption rate, that is, the maximum number of prey that can be eaten by predator in each time unit, and is the number of prey necessary to achieve one-half of the maximum rate . To the best of our knowledge, although the dynamics of the delayed Lotka-Volterra predator-prey model with nonmonotonic functional response has been studied in [28, 29], there are a few research papers on the dynamics of the delayed Leslie-Gower predator-prey model (4). Our goal is to obtain the explicit condition for the existence of the positive equilibrium, to investigate the local and global stability, and to study how the hunting delay affects the stability of the positive equilibrium and induces the Hopf bifurcations and their properties.

The paper is organized as follows. In Section 2, we study the explicit condition for the existence of the positive equilibria and their stability and the Hopf bifurcation induced by the intrinsic growth rate of the predator for the system without delay. In Section 3, we investigate the effect of the delay on the dynamics of the system. The local stability is studied by analyzing the related characteristic equation and the global stability is studied by the method of upper and lower solutions and monotone iterative scheme, and the direction and stability of the delay-induced Hopf bifurcation are determined by using the normal form theory. Numerical simulations are performed to illustrate and extend the obtained results in Section 4. A brief discussion is given in Section 5 to conclude the paper.

2. Linear Stability and Hopf Bifurcation Analysis for the System without Delay

In order to reduce the number of parameters and simplify the calculus in system (4), introducing new variables
and then dropping the bars, we derive
where , , and are positive parameters. When , system (6) becomes the following ordinary differential equations:
Systems (6) and (7) have the same equilibria, but the stability of system (6) is affected by delay. We first consider the existence of equilibria and their stability for system (7).

2.1. Existence of Equilibria and Their Stability

From the biological viewpoint, we are only interested in the dynamics of system (7) in the closed first quadrant in the plane. Letting denote the equilibrium of system (7), then are the positive real roots of the following equations:

We first have the following results on the existence of the positive equilibrium of system (7).

Lemma 1. (i) If , then system (7) has a unique positive equilibrium. (ii) If , then system (7) has three positive equilibria for , two positive equilibria for or , and a unique positive equilibrium for or , where , , with

Proof. (i) From (8), we have . Let , and . Notice that and is decreasing as the -value increases for since for . So, when , the two curves and have only one interaction point for and .When , the curve has one minimum at and one maximum at and decreases in the interval and increases in the interval , where
The tangent point of two curves and is the key point for determining how many intersection points of these two curves. The tangent point of two curves and can be determined by the roots of the following equations
which yields
It is easy to verify that the function has a minimum at for . Noticing that is equivalent to , we can deduce that (12) has no positive real roots for , one positive real root for , and two different positive real roots for . Thus, when , the straight line is not tangent to the curve for any and . When , the straight line is tangent to the curve only at and , as shown in Figure 1(a). Consequently, when , the two curves and have only one interaction point for any and . This completes the proof of conclusion (i). (ii) When , it follows from Shengjin’s theorem [30] that the two positive real roots of (12) are and defined by (9). This, together with (11), implies that two curves and are tangent at points , , as shown in Figure 1(b) for . This completes the proof of conclusion (ii).

Figure 1: Diagram of positive equilibrium. (a) When , there exists only one positive equilibrium for any . Here, , and . (b) When , there are three cases: (i) only one positive equilibrium for or , (ii) two positive equilibria for or , and (iii) three positive equilibria for . Here, , , and .

According to Lemma 1, the distribution of the positive equilibria of system (7) in the plane can be illustrated in Figure 2.

Figure 2: Bifurcation diagram of positive equilibrium in the plane. System (7) has three positive equilibria in the inner region enclosed by the curves and , a unique positive equilibrium outside this region and on the point and two positive equilibria on the two curves on the boundary curves and except the point .

The Jacobian matrix of system (7) at the positive equilibrium point is
where
Then the characteristic equation at is

Letting and denote the slopes of the straight line crossing the maximum and minimum points of , respectively, we have
where and are defined by (10). When , we have , as shown in Figure 2. The curves , , , divide the plane into four regions , , as shown in Figure 2. From Lemma 1, we know that system (7) has only one positive equilibrium in the regions and , two positive equilibria in the curves and with , three positive equilibria in the regions and . Notice that
and then by a direct but tedious analysis, we have the following results on the signs of and , which determine the stability of the positive equilibrium.

In the regions and , we have

In the regions and , denote the three positive equilibria by , with . Notice that , , . Then we have
In addition, notice that if , then for any . However, when , we can choose as a parameter such that for and for . Thus, by (18), (19), and (20), we have the following results on the distribution of roots of (15).

Lemma 2. Assume that , are defined by (9) and .(i)In the region , for the unique equilibrium , all roots of (15) have negative real parts for any .(ii)In the region , for the unique equilibrium , all roots of (15) have negative real parts for , a pair of purely imaginary roots at , and (15) has at least a root with positive real part for .(iii)In the region , for the positive equilibrium , (15) has at least a root with positive real part for any ; for other two positive equilibria and , all roots of (15) have negative real parts for , a pair of purely imaginary roots at , and (15) has at least a root with positive real part for .(iv)In the region , for the positive equilibrium , (15) has at least a root with positive real part for any ; for the positive equilibrium , all roots of (15) have negative real parts for , a pair of purely imaginary roots at , and (15) has at least a root with positive real part for ; for the positive equilibrium , all roots of (15) have negative real parts for any .

Taking as a parameter and letting is a root of (15) such that and . A straight calculation yields the following transversality condition:

From Lemma 2 and the transversality condition (21), the following theorem on the stability and bifurcation follows immediately.

Theorem 3. Assuming that , are defined by (9) and .(i)In the region , the unique equilibrium is asymptotically stable for any .(ii)In the region , the unique equilibrium is asymptotically stable for and unstable for , and system (7) undergoes Hopf bifurcation at .(iii)In the region , the positive equilibrium is unstable for any ; other two positive equilibria and are asymptotically stable for and unstable for , the and system (7) undergoes Hopf bifurcation near neighborhood of each of these two equilibria at .(iv)In the region , the positive equilibrium is unstable for any ; the positive equilibrium is asymptotically stable for and unstable for , and system (7) undergoes Hopf bifurcation near neighborhood of at ; the positive equilibrium is asymptotically stable for any .

2.2. Direction and Stability of Hopf Bifurcation Associated with

It follows from Theorem 3 that when , , Hopf bifurcation occurs at . When Hopf bifurcation occurs, the corresponding characteristic equation has a pair of purely imaginary roots at and thus the linearization is not sufficient to determine the stability of the corresponding equilibrium, say, . In this subsection, by calculating the focal value we determine the stability of the equilibrium at , which is directly associated with the direction and stability of Hopf bifurcation associated with .

Assume that are a pair of complex conjugates of (15) such that , . Then the eigenvector corresponding to is
Under the following coordinate transformation:
where , , system (7) becomes
where
with
where
Find and at (31). By the formula for the third focal value of a multiple focus in [31], we have
where

It is well known from the result in [31] that if , then the equilibrium of system (7) is a stable multiple focus of multiplicity 1 and the corresponding Hopf bifurcation is supercritical, and if , then the equilibrium of system (7) is an unstable multiple focus of multiplicity 1 and the corresponding Hopf bifurcation is subcritical.

3. Stability and Hopf Bifurcation Analysis for the System with Delay

In this section, we focus on the investigation of the local stability and Hopf bifurcation criteria of the positive equilibrium for the system (6). Letting and , we can rewrite the system (6) by Taylor series expression about as the following system:
where

3.1. Linear Stability and Delay-Induced Hopf Bifurcation

The characteristic equation at takes the form
When the time delay , the equilibrium is asymptotically stable if
Now for , assuming that is a root of (32), then we have
Separating the real and imaginary parts, we obtain
which leads to the following fourth degree polynomial equation:
It is easy to see that if
then (36) has only one positive root
Substituting the value of in (35) and solve for , we get
then when , the characteristic equation (32) has a pair of purely imaginary roots .

Let be a root of (32) such that the following two conditions hold:
Substituting into (32) and taking the derivative of resulting expression with respect to , we get
which, together with (35) and (38), leads to

Therefore, the transversality condition holds, and hence Hopf-bifurcation occurs at . Because the quartic equation (36) has only one positive root, system (6) has no switching stability. Delay-induced small-amplitude period solution bifurcates from interior equilibrium point when passes through its critical magnitude . Here, is the smallest positive value of given in (39).

To study the effect of the delay on the stability of the interior equilibrium in regions , , we just need decide the sign of by (18), (20), and (37). Notice that . So we first determine at which point the slope of curve is negative of the slope of line ; that is, we need to find the solution of the following equations:
which yields
Let . When , , which implies that the function is increasing for any . Notice that and , and we can deduce that (44) has only one positive real root for . When , which implies that the function has a maximum value at and a minimal value at . Furthermore, we can deduce that and . So (44) has also only one positive real root for . By direct computation, the positive real root of (44) is
for any . Then divides the regions and into two parts, respectively; see Figure 3. When , , and when , . Therefore,

Figure 3: Bifurcation diagram of positive equilibrium in the plane. divides and into , and , , respectively.

The above results are summarized in the following theorems.

Theorem 4. (i) In the region , for any , the unique equilibrium is asymptotically stable for and unstable for , and system (6) undergoes Hopf bifurcation near for , . (ii) In the region , for any , the unique equilibrium is asymptotically stable for any .

Theorem 5. In the region ,(i)when , the unique equilibrium is asymptotically stable for and unstable for , and system (6) undergoes Hopf bifurcation near for , ; (ii)when , the unique equilibrium is unstable for any , and system (6) undergoes Hopf bifurcation near for , .

Theorem 6. In the region ,(i)for any , the positive equilibrium is unstable for any ; (ii)for other two positive equilibria and , (a) when , they are asymptotically stable for and unstable for , and system (6) undergoes Hopf bifurcation near each of these two equilibria for , ; (b) when , they are unstable for any , and system (6) undergoes Hopf bifurcation near each of these two equilibria for , .

Theorem 7. (i) In the region , for any , the positive equilibrium is unstable for any ; for any , the positive equilibrium is asymptotically stable for and unstable for , and system (6) undergoes Hopf bifurcation near for , ; for the positive equilibrium , (a) when , is asymptotically stable for and unstable for , and system (6) undergoes Hopf bifurcation near for , ; (b) when , is unstable for any , and system (6) undergoes Hopf bifurcation near for , . (ii) In the region , for any , the positive equilibrium is unstable for any , and system (6) undergoes Hopf bifurcation near for , ; the positive equilibrium is asymptotically stable for and unstable for for any ; the positive equilibrium is asymptotically stable for any and for any .

3.2. Global Stability of the Positive Equilibrium

In this section, we give a result on the global stability of the positive equilibrium of system (6).

Theorem 8. If , then the unique positive equilibrium of system (6) is global stability. That is, for any initial values , , , with , the corresponding solution of system (6) satisfies

Proof. By the first equation of system (6), we have , which, together with the comparison principle, implies that . Thus, for any sufficiently small positive number , there exists such that for any ,
This, together with the second equation of system (6) yields that for ,
Again by the comparison principle, we have
which implies that for mentioned above, there exists such that for any ,
It follows from (51) and the second equation of system (6) that for ,
which leads to
So, there exists such that for any ,
Using this lower bound of , we obtain
which, together with the comparison principle, yields
Therefore, there exists such that for any ,
Clearly, , . It follows from (54) and (57) that when , we can choose
such that , . By (48)–(57), it is easy to verify that
which means that and are coupled upper and lower solutions of system (6) (see Definition 2.2 in [32]).Define two sequences of constant vectors and , , as follows:
where is the Lipschitz constant for the vector field of system (6) when .Notice that the above two sequences of constant vectors are recursive sequences starting with and . Then, by the principle of induction, the following monotone property of these sequences holds (see Lemma 2.1 in [33] for more details):
This monotone property implies that there exist positive constants and such that
By (60) and (63), we have
From (64), it is easy to verify that
Subtraction of these two equations gives
Then we can prove by contradiction. In fact, if , then (66) becomes
Equation (65) can be written as follows
Solving (68) gives
Substraction of the two equations of (69) and simplification lead to
Addition of the two equations of (69) and using (70) yield
It follows from (67) and (71) that
which is a contradiction under the assumption . So, the statement is proved.By (64), we also have , . Therefore, since . Letting and , then is the positive solution of (8). In addition, notice that when , (8) has the unique positive solution . Thus, . Then from the result due to Pao (see Theorem 2.2 in [33]), the proof of Theorem is complete.

3.3. Direction and Stability of Hopf Bifurcation Associated with the Critical Values of Delay

In Section 3.1, we have shown that the system (6) undergoes the Hopf bifurcation at critical values . In this subsection, we will derive the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits from the interior equilibrium of system (6) at . We will employ the algorithm of Faria and Magalhães [34] to compute explicitly the normal forms of system (6) on the center manifold.

We rescale the time by to normalize the delay and then introduce the new parameter so that is the critical value of Hopf singularity. System (30) can be rewritten as
in the phase space , where for , we have
By the Riesz representation theorem, there exists a function of bounded variation such that

Setting , then and . Let be the generalized eigenspace associated with the eigenvalues in and let be its dual space. For , the 2-dimensional space of row vectors, define and the adjoint bilinear form on as follows:
Denote the dual bases of and by and , respectively. Consider with
where is a nonzero solution in to the following equation:
and choose a basis for the adjoint space such that , where is a identical matrix. So with
where is the solution to the equation
and satisfies . By direct computation, we can choose
where .

Following the normal form theory of functional differential equations due to Faria and Magalhães [34] and using a very similar procedure as in [29] but a tedious calculation, we can obtain the following normal form at the critical values :
where . The coefficients and can be determined as follows:
where
with
with
The normal form (82) can be written in real coordinates through the change of variables , . Transforming to polar coordinates , , the normal form becomes