Abstract and Applied Analysis

Volume 2013 (2013), Article ID 152459, 19 pages

http://dx.doi.org/10.1155/2013/152459

Research Article

Stability and Bifurcation Analysis of a Delayed Leslie-Gower Predator-Prey System with Nonmonotonic Functional Response

Jiao Jiang¹ and Yongli Song²

¹Department of Mathematics, Shanghai Maritime University, Shanghai 201306, China

²Department of Mathematics, Tongji University, Shanghai 200092, China

Received 21 December 2012; Accepted 7 April 2013

Academic Editor: Kunquan Lan

Copyright © 2013 Jiao Jiang and Yongli Song. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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).

(a)

(b)

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.

(a)

(b)

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,

(a)

(b)

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

with , . It is well known [35] that the sign of determines the direction of the bifurcation (supercritical if , subcritical if ), and the sign of determines the stability of the nontrivial periodic orbits (stable if , unstable if ). Thus, if the coefficients of system (4) are given, then we will analyze the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits at .

4. Numerical Simulations

In this section, we present some numerical simulations to verify and extend our theoretical analysis proved in Section 3 by using the Simulink for Matlab. To investigate the effect of delay, we are interested in the cases when , , and . We would also like to mention that for the case of delay, the initial values should be given by ,   for . However, in the Simulink for Matlab, we only need to specify the initial values as and . The initial values for the case of delay can be selected automatically by the default function of Simulink.

(i) Taking and such that , system (6) has a unique positive equilibrium . Without delay, the positive equilibrium is asymptotically stable for any . Let and then by (39), we have the critical value of delay and it follows from Theorem 4 that is asymptotically stable for and unstable for , as shown in Figures 4(a) and 4(b) for . Furthermore, by the procedure of Section 3.3, we get , and , which implies that the Hopf bifurcation associated with this critical value is subcritical. So, there exists a sufficiently small positive real number such that for each , system (6) has an unstable periodic solution near the positive equilibrium . Taking , Figures 4(c) and 4(d) show the existence of the subcritical Hopf bifurcation. In this case, a small-amplitude unstable periodic solution, coming from the Hopf bifurcation, and a large-amplitude stable periodic solution coexist. Taking , the positive equilibrium becomes unstable and numerical simulation also shows the existence of large-amplitude stable periodic solution (see Figures 4(e) and 4(f)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: The trajectory graphs and phase portrait of system (6) for with ,  , and . (a)-(b) When , the equilibrium is stable. (c)-(d) Simulations of subcritical Hopf bifurcation. Starting from the initial value , the trajectory slowly gets away from a small-amplitude unstable periodic solution and finally converges to a large-amplitude stable periodic solution. Starting from the initial value , the trajectory converges to . (e)-(f) When , the equilibrium is unstable with a large-amplitude stable periodic solution.

(ii) Taking and such that , system (6) has a unique positive equilibrium . Without delay, the parameter affects the stability of the positive equilibrium . By Lemma 2, we get . The positive equilibrium is stable for , and unstable for and system (7) undergoes Hopf bifurcation near the positive equilibrium at . By (28), we have , which means that the Hopf bifurcation occurring at is subcritical. So, there exists a sufficiently small positive real number such that for any , system (7) has an unstable periodic orbit near the positive equilibrium . Taking , Figures 5(a) and 5(b) illustrate the instability of the positive equilibrium and the existence of a large-amplitude stable periodic orbit. Taking larger than and closer to , Figures 5(c) and 5(d) illustrate the stability of the positive equilibrium and the existence of subcritical Hopf bifurcation. In this case, a small-amplitude unstable periodic solution, coming from the Hopf bifurcation, and a large-amplitude stable periodic solution coexist. Taking larger than and away from , the positive equilibrium is stable and there are no Hopf bifurcating periodic orbits, as shown in Figures 5(e) and 5(f).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: The trajectory graphs and phase portrait of system (7) without delay for with and . (a)-(b) . The positive equilibrium is unstable with a large-amplitude stable periodic orbit. The initial value is . (c)-(d) . The positive equilibrium is stable, and a small-amplitude unstable periodic solution and a large-amplitude stable periodic solution coexist. The initial values are and . (e)-(f) . The positive equilibrium is stable and there are no Hopf bifurcating periodic orbits. The initial value is .

To investigate the effect of delay, we fix ,  , and . For and , the positive equilibrium is stable and there exists an unstable periodic orbit near the positive equilibrium (see Figures 5(c) and 5(d)). The stability of the positive equilibrium changes with the increase of delay. By (39), the critical value of delay is . By the procedure of Section 3.3, we have , and , which implies that the Hopf bifurcation occurring at is subcritical. So, the positive equilibrium is stable with the coexistence of small-amplitude stable and large-amplitude stable periodic orbits for and unstable for . Figures 6(a) and 6(b) and Figures 6(c) and 6(d) are numerical simulations for and , respectively. For and , the positive equilibrium is stable and there are no periodic orbits near the positive equilibrium (see Figures 5(e) and 5(f)). In this case, the critical value of delay is . The positive equilibrium is stable for and unstable for . Figures 7(a) and 7(b), Figures 7(c) and 7(d), and Figures 7(e) and 7(f) are numerical simulations for , , and , respectively. It follows from the procedure of Section 3.3 that and , which implies that the Hopf bifurcation occurring at is also subcritical. Thus, there exists a positive real number such that for each , there exists an unstable periodic orbit, as shown in Figures 7(c) and 7(d). For , the the positive equilibrium becomes unstable and there exists a large-amplitude stable periodic orbit, as shown in Figures 7(e) and 7(f) for .

(a)

(b)

(c)

(d)

Figure 6: The trajectory graphs and phase portrait of system (6) for with ,  , and . (a)-(b) For , the positive equilibrium is stable with small-amplitude unstable and large-amplitude stable periodic orbits. The initial values are and . (c)-(d) For , the positive equilibrium is unstable with a large-amplitude stable periodic orbit. The initial value is .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: The trajectory graphs and phase portrait of system (6) for with ,  , and . (a)-(b) For , the positive equilibrium is stable. The initial values are . (c)-(d) For , the positive equilibrium is stable with an unstable periodic orbit. The initial values are and . (e)-(f) For , the positive equilibrium is unstable with an stable periodic orbit. The initial values are and .

(iii) Taking and such that , system (6) has three positive equilibria ,  , and . The positive equilibrium is always unstable for any and . But the parameters and affect the stability of and . Without delay, the critical values of for and are and . The positive equilibrium (resp., ) is asymptotically stable for (resp., ) and unstable for (resp., ). Taking , all of these three equilibria are unstable. There is a large-amplitude stable periodic orbit enclosing these three equilibria, as shown in Figures 8(a) and (b). For , the positive equilibria and are still unstable, but the positive equilibrium becomes stable, as shown in Figures 8(c) and 8(d) with .

(a)

(b)

(c)

(d)

Figure 8: The trajectory graphs and phase portrait of system (7) without delay for with and . (a)-(b) . The initial value is . All of the three equilibria are unstable and there exists a stable periodic orbit. (c)-(d) . The initial value is . The positive equilibria and are unstable and the positive equilibrium is stable.

When , the positive equilibria and both become asymptotically stable. Taking , Figures 9(a) and 9(b) illustrate this theoretical result. The critical values of the delay for these two positive equilibria and are and , respectively. Thus the positive equilibrium is asymptotically stable for and unstable for , and the positive equilibrium is asymptotically stable for and unstable for . According to the procedure of Section 3.3, we obtain that , and for and , and , and for and . So, the delay-induced Hopf bifurcations for these two positive equilibria are both subcritical. Taking , the positive equilibrium and are both asymptotically stable. Near the positive equilibria , the subcritical Hopf bifurcation occurs. Figures 9(c) and 9(d) illustrate the existence of this subcritical Hopf bifurcation and show that there is an orbit connecting the unstable periodic orbit to the stable equilibrium . When , the positive equilibria and both becomes unstable. There is a large-amplitude stable periodic orbit enclosing the three equilibria, as shown in Figures 9(e) and 9(f) for .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: The trajectory graphs and phase portrait of system (6) for with ,  , and . (a)-(b) When , the positive equilibria and are both asymptotically stable. The initial values are and . (c)-(d) When , the positive equilibria and are both asymptotically stable and there exists a subcritical Hopf bifurcation near . Here, the initial values is . (e)-(f) When . All of the three equilibria become unstable and there exists a large-amplitude stable periodic orbit enclosing the three equilibria. Here, the initial value is .

5. Conclusions

In this paper, we have considered a predator-prey system with the hunting delay. By analysing the concerned characteristic equations, local asymptotic stability of the various positive equilibria of the system with or without delay is discussed. The dynamics of the system can be described in the plane of the parameters and . For the system without delay, the bistability (existence of two stable equilibria) occurs in some parameter regions. We have shown that the ratio () of the intrinsic growth rates of the predator and prey is the key parameter which has an important influence on the dynamics of the system. In the region where the system has only one positive equilibrium and the ratio can affect the stability of this positive, the ratio can induce large-amplitude periodic oscillations of both species and there are two limit cycles (one is stable and the other is unstable) near the critical value of the ratio via subcritical Hopf bifurcation. In the region where the system has three positive equilibria, the dynamics of the system can evolve from a heteroclinic orbit connecting one positive equilibrium to a large-amplitude limit cycle enclosing three equilibria (see Figures 8(a) and 8(b)), a heteroclinic orbit connecting two equilibria (see Figures 8(c) and 8(d)), two stable equilibria with the increase of the ratio (see Figures 9(a) and 9(b)).

When the hunting delay is introduced by the system, its effect on the dynamics of all positive equilibria is investigated. We first show that when (or equivalent to ), the unique positive equilibrium is global stability for any nonnegative hunting delay. Then the absolute investigate the influence of the delay on the local stability of the positive equilibria. For the case when the system has the unique positive equilibrium, depending on the choice of the parameters and , the delay can induce the emergence of two limit cycles with one stable subcritical Hopf bifraction and another unstable subcritical Hopf bifurcation (see Figures 4(c) and 4(d)) and induce the emerge of semistable limit cycle which is inward unstable and outward stable (see Figures 7(c) and 7(d)). For the case when the system has the three positive equilibria, delay can induce a heteroclinic orbit connecting a limit cycle to a stable equilibrium (see Figures 9(c) and 9(d)). The analytical expressions that determine the properties of bifurcating periodic solutions are also obtained by using the normal form theory and the centre manifold theorem. It is also worthy to mention that for the case without delay, the size of limit cycle increases with the decreasing magnitude of the ratio of the intrinsic growth rates of the predator and prey, but for the case with delay, the size of limit cycle increases with the increasing magnitude of the hunting delay.

Acknowledgments

The authors would like to thank the three anonymous referees very much for their constructive comments, which resulted in much improvement of the paper. The authors would also like to express appreciation for the support from the National Natural Science Foundation of China (nos. 11201294 and 11032009), the Scientific Research Foundation for the Returned Overseas Chinese Scholars, the Fundamental Research Funds for the Central Universities, and the Program for New Century Excellent Talents in University (NCET-11-0385).

References

S. Allesina and S. Tang, “Stability criteria for complex ecosystems,” Nature, vol. 483, no. 7388, pp. 205–208, 2012.

A. J. Lotka, “Analytical note on certain rhythmic relations in organic systems,” Proceedings of the National Academy of Sciences of the United States of America, vol. 6, no. 7, pp. 410–415, 1920.

V. Volterra, “Variations and uctuations of the number of individuals in animal species living together,” in Animal Ecology, R. N. Chapman, Ed., McGraw-Hill, New York, NY, USA, 1931.

C. S. Holling, “The components of predation as revealed by a study of small mammal predation of the European pine sawfly,” The Canadian Entomologist, vol. 91, no. 5, pp. 293–320, 1959.

C. S. Holling, “Some characteristics of simple types of predation and parasitism,” The Canadian Entomologist, vol. 91, no. 6, pp. 385–398, 1959.

C. S. Holling, “The functional response of predators to prey density and its role in mimicry and population regulation,” Memoirs of the Entomological Society of Canada, vol. 46, pp. 1–60, 1965.

P. H. Leslie, “Some further notes on the use of matrices in population mathematics,” Biometrika, vol. 35, pp. 213–245, 1948.

P. H. Leslie, “A stochastic model for studying the properties of certain biological systems by numerical methods,” Biometrika, vol. 45, pp. 16–31, 1958.

P. H. Leslie and J. C. Gower, “The properties of a stochastic model for the predator-prey type of interaction between two species,” Biometrika, vol. 47, pp. 219–234, 1960.

R. M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, NJ, USA, 1973.

D. K. Arrowsmith and C. M. Place, Dynamical Systems, Chapman and Hall Mathematics Series, Chapman & Hall, London, UK, 1992.

S. B. Hsu and T. W. Huang, “Global stability for a class of predator-prey systems,” SIAM Journal on Applied Mathematics, vol. 55, no. 3, pp. 763–783, 1995.

E. Sáez and E. González-Olivares, “Dynamics of a predator-prey model,” SIAM Journal on Applied Mathematics, vol. 59, no. 5, pp. 1867–1878, 1999.

S.-B. Hsu and T.-W. Hwang, “Hopf bifurcation analysis for a predator-prey system of Holling and Leslie type,” Taiwanese Journal of Mathematics, vol. 3, no. 1, pp. 35–53, 1999.

P. A. Braza, “The bifurcation structure of the Holling-Tanner model for predator-prey interactions using two-timing,” SIAM Journal on Applied Mathematics, vol. 63, no. 3, pp. 889–904, 2003.

A. Korobeinikov, “A Lyapunov function for Leslie-Gower predator-prey models,” Applied Mathematics Letters, vol. 14, no. 6, pp. 697–699, 2001.

J. B. Collings, “The effects of the functional response on the bifurcation behavior of a mite predator-prey interaction model,” Journal of Mathematical Biology, vol. 36, no. 2, pp. 149–168, 1997.

C. S. Holling, “Principles of insect predation,” Annual Review of Entomology, vol. 6, pp. 163–182, 1961.

W. Sokol and J. A. Howell, “Kinetics of phenol exidation by washed cells,” Biotechnology and Bioengineering, vol. 23, pp. 2039–2049, 1980.

Y. Li and D. Xiao, “Bifurcations of a predator-prey system of Holling and Leslie types,” Chaos, Solitons & Fractals, vol. 34, no. 2, pp. 606–620, 2007.

Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, vol. 191 of Mathematics in Science and Engineering, Academic Press, Boston, Mass, USA, 1993.

S. Ruan, “On nonlinear dynamics of predator-prey models with discrete delay,” Mathematical Modelling of Natural Phenomena, vol. 4, no. 2, pp. 140–188, 2009.

Z. Lu and X. Liu, “Analysis of a predator-prey model with modified Holling-Tanner functional response and time delay,” Nonlinear Analysis: Real World Applications, vol. 9, no. 2, pp. 641–650, 2008.

A. F. Nindjin, M. A. Aziz-Alaoui, and M. Cadivel, “Analysis of a predator-prey model with modified Leslie-Gower and Holling-type II schemes with time delay,” Nonlinear Analysis: Real World Applications, vol. 7, no. 5, pp. 1104–1118, 2006.

S. Yuan and Y. Song, “Bifurcation and stability analysis for a delayed Leslie-Gower predator-prey system,” IMA Journal of Applied Mathematics, vol. 74, no. 4, pp. 574–603, 2009.

J.-F. Zhang, “Bifurcation analysis of a modified Holling-Tanner predator-prey model with time delay,” Applied Mathematical Modelling, vol. 36, no. 3, pp. 1219–1231, 2012.

F. Lian and Y. Xu, “Hopf bifurcation analysis of a predator-prey system with Holling type IV functional response and time delay,” Applied Mathematics and Computation, vol. 215, no. 4, pp. 1484–1495, 2009.

D. Xiao and S. Ruan, “Multiple bifurcations in a delayed predator-prey system with nonmonotonic functional response,” Journal of Differential Equations, vol. 176, no. 2, pp. 494–510, 2001.

Z. Liu and R. Yuan, “Bifurcations in predator-prey systems with nonmonotonic functional response,” Nonlinear Analysis: Real World Applications, vol. 6, no. 1, pp. 187–205, 2005.

F. Shengjin, “A new extracting formula and a new distinguishing means on the one variable cubic equation,” Natural Science Journal of Hainan Teachers College, vol. 2, no. 2, pp. 91–98, 1989.

J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, vol. 42 of Applied Mathematical Sciences, Springer, New York, NY, USA, 1983.

C. V. Pao, “Dynamics of nonlinear parabolic systems with time delays,” Journal of Mathematical Analysis and Applications, vol. 198, no. 3, pp. 751–779, 1996.

C. V. Pao, “Convergence of solutions of reaction-diffusion systems with time delays,” Nonlinear Analysis: Theory, Methods & Applications, vol. 48, no. 3, pp. 349–362, 2002.

T. Faria and L. T. Magalhães, “Normal forms for retarded functional-differential equations with parameters and applications to Hopf bifurcation,” Journal of Differential Equations, vol. 122, no. 2, pp. 181–200, 1995.

S. N. Chow and J. K. Hale, Methods of Bifurcation Theory, Springer, New York, NY, USA, 1982.