Applied Mathematics and Computation 217 (2010) 3542-3556 ELSEVIER Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc APPLIED MATHEMATICS COMPUTATION The impact of Allee effect on a predator-prey system with Holling type II functional response Jian Zua'*, Masayasu Mimurab 'Graduate School of Science and Technology, Meiji University, Kawasaki, Kanagawa 214-8571, Japan bMeiji Institute for Advanced Study of Mathematical Sciences, Kawasaki, Kanagawa 214-8571, Japan ARTICLE INFO ABSTRACT Keywords: Extinction risk Limit cycle Heteroclinic cycle Stability switches Positive interactions In this paper, the Allee effect is incorporated into a predator-prey model with Holling type II functional response. Compared with the predator-prey model without Allee effect, we find that the Allee effect of prey species increases the extinction risk of both predators and prey. When the handling time of predators is relatively short and the Allee effect of prey species becomes strong, both predators and prey may become extinct. Moreover, it is shown that the model with Allee effect undergoes the Hopf bifurcation and heteroclinic bifurcation. The Allee effect of prey species can lead to unstable periodical oscillation. It is also found that the positive equilibrium of the model could change from stable to unstable, and then to stable when the strength of Allee effect or the handling time of predators increases continuously from zero, that is, the model admits stability switches as a parameter changes. When the Allee effect of prey species becomes strong, longer handling time of predators may stabilize the coexistent steady state. Crown Copyright © 2010 Published by Elsevier Inc. All rights reserved. 1. Introduction Both positive and negative interactions among species are common in communities [1,2]. For example, synchronous breeding is an important mechanism by which colonial guillemots (Una aalge) increase reproductive success. Many carnivore species are better able to capture large prey by cooperative hunting and group living may thus be favored in areas with abundant large prey [3]. Historically, attention has focused on negative interactions, such as competition. However, for the last decade, the importance of positive interactions such as the Allee effect [4] has recently been recognized [1,5]. The Allee effect [4], and more recently as depensation [6] or inverse density dependence [7], refers to a decrease in per capita fertility rate at low population densities. Allee effect may occur under several mechanisms, such as difficulties in finding mates when population density is low [8,9], social dysfunction at small population sizes and increased predation risk due to failing flocking or schooling behavior [10,11]. When such mechanisms operate, the per capita fertility rate of a species increases with density, that is, positive interactions among species occur. The primary consequence of the Allee effect is that it creates a threshold density below which a population cannot survive. For example, this might correspond to the density below which it is so difficult to find a mate that reproduction does not compensate mortality. Each population whose density stochastically goes below this threshold is fated to extinction and species experiencing Allee effect are therefore more extinction prone [11]. The phenomenon has received considerable attention from ecologists [12,13]. The importance of this dynamic process in ecology has been under-appreciated and recent evidence [14] now suggests that it might have an important influence on the population dynamics of many plants and animal species. Therefore, the investigation of Allee effects is important for conservation biology [15]. * Corresponding author. E-mail address: jianzumb@gmail.com (J- Zu). 0096-3003/$ - see front matter Crown Copyright © 2010 Published by Elsevier Inc. All rights reserved, doi: 10.1016/j.amc.2010.09.029 J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 3543 Interactions among species in a food web are multiform. Among these, the interactions between predators and prey are important. Predator-prey interactions are ubiquitous in nature and the dynamical behaviors of predator-prey system are very complex. Many forces may influence the dynamical behaviors of the predator-prey system [11,12]. It has been shown that a predator-prey model can lead to oscillatory behavior because of a nonlinear functional response. Interest in the stability of predator-prey system has continued unabated since the theoretical work of Lotka [16], Volterra [17] and the experimental work of Cause [18]. Theoreticians and experimentalists have proceeded to investigate the processes that affect the stability of predator-prey system [19,20]. Among the many processes that the Lotka-Volterra model ignores, the Allee effect may be the most important [21]. The Allee effect increases the likelihood of local and global extinction. But there have been few papers discussing its stabilizing or destabilizing effects on the predator-prey system (except for [21,22,33]). Kent et al. [22] conclude that the predator-prey system is stabilized by an influx of prey in the form of a rescue effect, and destabilized by an out flux of Allee effect. By combining mathematical analysis with numerical simulation, Zhou et al. [21] have shown that the Allee effect may be a destabilizing force in predator-prey system. However, the Lotka-Volterra model they considered does not consider the density-dependent effect of prey and the functional response is linear. Functional response is a double rate: it is the average number of prey killed by per individual predator per unit of time. In general, the functional response can be classified into two types: prey-dependent and predator-dependent. Prey-dependent means that the functional response is only a function of the prey's density, while predator-dependent means that the functional response is a function of both densities of the prey and predators. Functional response that is strictly prey-dependent, such as the Holling family, is predominant in the literature. For example, since 1959, Holling's prey-dependent type II functional response has served as the basis for many literatures on predator-prey theory. Therefore, the predator-prey model with Holling type II functional response is more realistic. However, the impact of the Allee effect on the stability of a predator-prey system with Holling type II functional response is poorly understood both empirically and theoretically. The purpose of this paper is to show that the Allee effect of prey species has significant effects on the dynamics of predator-prey model with Holling type II functional response. We will investigate how predator and prey species are threatened by extinction when the population density of prey becomes low. Moreover, we will show that the predator-prey model with Allee effect and Holling type II functional response undergoes a sequence of bifurcations including supercritical Hopf bifurcation, subcritical Hopf bifurcation and heteroclinic bifurcation. We will also present a global analysis of the model by means of numerical simulations. The organization of this paper is as follows. In the next section, we present the formulations of mathematical model with Allee effect. In Section 3, we present a qualitative analysis of the model. In Section 4, we use numerical simulations to reveal the global bifurcation structures and the influence of Allee effect on the dynamical behaviors of the model. A brief discussion is given in Section 5. 2. Model formulations In this section, we first introduce the Rosenzweig-MacArthur predator-prey model with Holling type II functional response. Based on this model, we will construct a predator-prey model with Allee effect on prey species. 2.1. Rosenzweig-MacArthur predator-prey model The objective of this subsection is to introduce the predator-prey model with Holling type II functional response and summarize its dynamical behaviors. The classical Lotka-Volterra predator-prey model did not contain saturating effect. More realistically, the functional response of predator population is nonlinear. The Rosenzweig-MacArthur predator-prey model with Holling type II functional response is as follows: where N and P denote the population densities of prey and predators at time t, respectively, b is the per capita maximum fertility rate of prey population, d, (i = 1,2) are the per capita death rates of prey and predators respectively, a denotes the strength of intra-competition of prey population, s denotes the effective search rate, h\ denotes the handling time of predators, and C\ denotes the conversion efficiency of ingested prey into new predators. The product, sN/(l + sh\N), represents the predator's functional response, i.e. the relationship between prey density, N, and the amount ingested by an average predator. All the parameters are positive constants. Further, we assume that the growth rate of prey must exceed its death rate, i.e., b > d^, otherwise, both prey and predators will become extinct. In addition, the maximum growth rate of predator population must exceed its death rate, i.e., C\\h\ > d2. If not, prey population will never be able to sustain predators, which will go to extinction. In model (1), if the carrying capacity of prey (b - d\)ja. is low, i.e., 0) 0 < b-d (2) < s(ci -hid2)' a 3544 J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 then the attractor is the trivial equilibrium ((b - di)/a,0), i.e., predator population goes extinct. By contrast, for intermediate value of the carrying capacity d2 < b-d] < d2 + c, /h. s(c]-h]d2) a ^s(c]-h]d2)' stationary stable coexistence occurs at the strictly positive equilibrium while for high value of the carrying capacity of prey, i.e., b -di d2 +ci/hi s(ci -hid2)' predator and prey populations coexist on a limit cycle and the limit cycle is globally stable. Moreover, the trivial equilibrium (0,0) is always a saddle point in model (1). Therefore, the equilibrium (N\,P\) given by (3) is critically stable and the populations are balanced between stationary and cyclic coexistence. Next, based on the predator-prey model (1), we construct a predator-prey model with Allee effect on prey species and study its impact on the stability of predator-prey system. 2.2. Predator-prey model with Allee effect on prey species Because of difficulties in finding mates when prey population density becomes low, Allee effect may occur in prey species [9]. For example, this might correspond to the density below which it is so difficult to find a mate that reproduction does not compensate for mortality. Let/(N) be the fertility rate of a species that has N adults in an isolated patch, then the fertility rate increases with population density, which is described by /(N)= bN A, +AT where b is the per capita maximum fertility rate of the species, Ai is the Allee effect constant of the species [23]. If Ai > 0, the fertility of the species is therefore zero when N is zero and approaches to b when N becomes very large. How fast f[N) increases with N depends on the parameter/^. The bigger A\ is, the stronger the Allee effect will be. Biologically, A\ is the population density at which a species reaches half its maximum fertility {i.e..,f[A\) = b/2). When/^ = 0, the fertility rate is density independent, with f[N) = b. When prey population is subject to Allee effect, the predator-prey model becomes ■dN ..( bN , M\ sNP dt~ \A,+N 1 ) 1+sM' dP dsNP ( ' dt ~ 1 +sh,N " 2 ' where A\ is a positive parameter denoting Allee effect imposed on prey population. If u a u a b-dt d2 b > di,C] > h]d2, -> 0 < A] < - a s(ci - h]d2)' d2(s(b -di)(ci -h]d2)-ad2) ' s(ci - h]d2)(dis(ci - h]d2) + ad2)' then model (4) yields one non-trivial equilibrium (N-,^) = ( d? V+sh,N'2)(bN*2/(A,+N'2)-d,-aN'2)\ \s(ci - h]d2)' s J It is obvious that now the equilibrium density of prey population is the same as that obtained from Eq. (1), however, the equilibrium density of predator population is smaller than that obtained from Eq. (1). That is to say, a decrease of predator population density at the equilibrium is caused by the Allee effect of prey population. Next, we investigate model (4) by way of mathematical analysis and numerical simulations. First, we perform a global qualitative analysis of system (4). 3. Qualitative analysis In order to reduce parameters, we perform the following transformations: J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 3545 For simplicity, we still use variables N, P, t instead of X, Y, z. Then we obtain where A corresponds to the re-scaled strength of the Allee effect and m\ is a re-scaled death rate of prey population. All the parameters are positive constants. Note that b > d^, so 0 < m.\ = d\\b < 1 and we shall assume that this is the case in this paper. For simplicity of computation, we consider the above system (6), which is equivalent to (4). Henceforth, we perform a qualitative analysis of system (6). We start by studying the local stability of equilibriums of system (6). 3.1. Extinction equilibrium First, we consider the stability of trivial equilibrium £0 = (0,0). Because the stability of £0 = (0,0) is important to understand the qualitative behavior of model (6) influenced by the Allee effect of prey population. The Jacobian matrix of system (6) at £0 is -m-i 0 0 -m2 It can be seen that £0 is always a locally stable node, which implies that both prey and predators will become extinct when their population densities lie in the attraction region of £0. In particular, if the population density of prey becomes low, then both predators and prey will go extinct (see Fig. 1). Compared with the predator-prey model without Allee effect, we can see that the Allee effect of prey population increases the extinction risk of both predators and prey. Further, if model (6) does not have interior equilibrium and other boundary equilibriums, i.e., the Allee effect of prey population is very strong, then £0 is globally asymptotically stable. Any positive orbit converges to £0 as t tends to infinity, i.e., prey and predators cannot coexist even if the initial population density of prey is abundant. A typical phase portrait is shown in Fig. la. 3.2. Boundary equilibriums In order to find positive equilibriums for prey population, set A, = (1 -A-m^f -4/Wni =A2 -2A(m, +l) + (nii -l)2. N N Fig. 1. Phase portraits of model (6). Isoclines are shown as dashed lines, (a) E0 is a global attractor when A = 0.0880. (b) Both E0 and E2 are attractors when A = 0.0160. Other parameter values: h = 11.2500, m, =0.5000, c = 2.5000, m2 = 0.2000. 3546 J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 Then we can see that there exist two boundary equilibriums £i = (Ni,0) and £2 = (N2,0) if Ai > 0 and A < 1 - m.\, where Ni = i (\ - A - m, - \J{\ -A-m.])2 - AAm.] 2 N2 = \{^ - A-mi + y/(l-A-m])2-44m We first determine the local stability of E\. Note that N\I(A + Ni) - m.\ - N\ = 0, then the Jacobian matrix of system (6) at E\ is l-A-m]-2N] Ni M, = 0 1 +hN, ■ m2 1 +hN, The determinant of Jacobian matrix M2 is given by det(M2) = 1 -A-m^y -44m]((c-m2h)N] -m2) (A + N,)(l +hN,) It can be seen that the determinant of the Jacobian matrix M2 may be positive or negative. The trace of Jacobian matrix M2 is , 1 - A - m, - 2Ni tr(M2) =-j—I-i- + V(l-*-mi)2 1 +hN, ■m2 -44m, (c- m2h)Ni - m2 A + N] 1 1 + hN, When det(M2) > 0, the trace of Jacobian matrix M2 is always positive. Therefore, we obtain the following results on the stability of E\. Theorem 1. Assume that 0 < m, < ? and 0 < A < (1 - ^TnT)2. Then (i) £j is a saddle point if ' 1 -A-m^f -44mA (7) h > c/m2 - 2J [ \ -A-mi (ii) £j is unstable if 0 c/m2 ~2 J y - A- m + V0 — -A — "ii) -44m, If model (6) does not have interior equilibrium, i.e., the carrying capacity of prey is low, then its asymptotic behavior is determined by the local stability of £0, £i and £2. Note that A corresponds to a re-scaled strength of the Allee effect and h corresponds to a re-scaled handling time of predator population. From Theorems 1 and 2, it can be seen that if the Allee effect of prey is relatively weak and the handling time of predators is relatively long, then £i is a saddle point and £2 is asymptotically stable. The stable manifold of E\, which is the threshold induced by the Allee effect of prey population, divides the first quadrant into two parts, one is the attraction region of £0, the other is the attraction region of £2 (see Fig. lb). From Fig. lb, it can be seen that if the population density of prey becomes low, then both predators and prey will become extinct. This J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 3547 Fig. 2. Phase portraits of model (6). Isoclines are shown as dashed lines, (a) E0 is a global attractor when A = 0.0720. (b) E0 is a global attractor when A = 0.0400, E3 is unstable, (c) An unstable limit cycle when A = 0.0040. (d) Both E0 and E3 are attractors when A = 0.0020. Other parameter values are the same as in Fig. 1 except h = 0.3125. dependence on the initial population densities is very important. A slight difference in the initial population densities makes the asymptotic state very different [24]. On the contrary, if the Allee effect of prey is relatively strong and the handling time of predators is relatively short, then E\ is unstable and E2 is a saddle point, any positive orbit converges to E0 as t tends to infinity, i.e., E0 is globally asymptotically stable, both predators and prey will go extinct independent of the initial population densities (see Fig. 2 a). Therefore, if model (6) does not have interior equilibrium, then depending on the strength of Allee effect, the handling time of predators and the initial population densities, either predator population becomes extinct or both predators and prey go extinct. 3.3. Coexistent equilibrium To obtain the positive equilibrium E3 = (N*, P*) of system (6), we propose the following assumption: n m2(l + (1 - m,)/;) 0—-—t-^-—, 1 - m, (9] 0 m2h, so the sign of tr(M4) is determined by the sign of F(A), where F(A) is a quadratic function of A, i.e., F(A) =aiA2 +a2A + a3. Because ax < 0, in order to determine the sign of F(A), we further assume that a2 > 0 and a3 < 0, i.e., c> m2(2 + h), h(c - m2h) - (c + m2h) (c2 - m\h2) - 2m2(c+ m2h) h(c-m2h) < mi < 2m2h(c - m2h) ' and set A2 = a\ - 4a]a3 = -(c - m2h)3(Acniim2h(c - m2h) + (c + m2h)(Acm2 + ni\h2 - c2)j. When conditions (12) hold and A2 > 0, then F(A) = 0 admits two positive solutions A* and A**, where a2 - .« a2 + (10) (11) (12) (13) A* = A" = 2di 2di In this case, if 0 A**, then F{A) < 0, i.e., tr(M4) < 0, E3 is locally asymptotically stable. When conditions (12) hold and A2 < 0, then F(A) < 0 for all A > 0, i.e., E3 is always locally asymptotically stable. Set „ _ (c + m2h)(c2 - 4cm2 - m\h2) 1 4cm2h(c - m2h) ' m2((c- m2h)(l - nil)- m2) Ac = (c - m2li)(m2 + nii(c - m2h))' then we have the following results on the local stability of E3 if (9) and (12) hold (see Table 1). J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 3549 Table 1 The local stability of the positive equilibrium E3 of model (6). Range of Range of Allee parameter mi effect parameter A The stability of (AT, P*) 0 < mi < 1715 0 < A < A* or A™ < A < Ac Asymptotically stable mi > 1715 777] = 7775 A' 1 - in] and one of the opposite inequalities of (9) holds; (ii) 0 < A < (1 - ^fnT)2, (7), (8) and one of the opposite inequalities of (9) hold; (iii) 0 0, then there is a family of unstable limit cycles in (6) as A decreases from A*. Proof. Suppose A =A*, then tr(M4) = 0. Set m = ^det(M4), then the eigenvalues of M4 are X\ = on and X2 = -an. For simplicity of computation, we consider the following system which is equivalent to (6): + 2ch)(A*)2 + (-2ch + 2m2h2 - 2m^m2h2 + 6m2h + 2cm^h + 2c)A* + m2h - m2 - m2m^h), (16) -r- = N(N(1 + hN) - (m, + N)(A + N)(l + hN) -P(A + N)), (17) Make a transformation of x = N - N*, y = P - P* to translate (N*, P*) to the origin. Then (17) becomes (18) where f{x,y) (1 = 1,2) represent the higher order terms and a„ = N "(1 + 2hN' - (m, +N *)(A + N *)h - (mi +A + 2N *)(1 + hN *) - P*), au = -N'(A + N'), a2i = (c - m2h)(A + N *)P*, a22 = 0. Setting u = x, anx auy (19) v and using tr(M4) flu +022 =0, at2 = det(M4) = -ai2a2i, we obtain (20) 3552 J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 where fiu,v)=f,(u,——y which represent the higher order terms. Let = ~Yq(Juuu + fuvv Suuv + Svvv) + 15(y (fuvifuu ~\~ fvv) ~ SuviSuu ^ Svv) ~ fuuSuu ^ fvv§vv)> where /u„ denotes (52//(&idv))(0,0), etc. Using the fact that A = A*, with the aid of software MAPLE13, we obtain m2((c - m2h)K + m2)a 8m2(c-m2h)4 Since c > m2h, we can see that the sign of \i is determined by a. The conclusion of this Theorem follows from Theorem 3.4.2 and formula (3.4.11) of Guckenheimer and Holmes [25]. □ As an example, we fix mx = 0.5000, h = 0.3125, c = 2.5000, m2 = 0.2000. Then (9), (12) and 0 < m, < m\ hold. By (16), we have a = 0.2799 >0. Theorem 5 shows that there is an unstable limit cycle when A decreases from A* = 0.0071, which is shown in Fig. 2c when A = 0.0040. Therefore, we can see that the stable coexistent region of predators and prey collapses with the disappearance of the unstable limit cycle. In addition, if the handling time of predators is long, such as h = 7.8000, then we have a = -0.2845 < 0. Thus, there is a stable limit cycle when A increases from A* = 0.0362, which is shown in Fig. 3b when A = 0.0500. This analysis reveals that the unstable or stable periodical oscillation can be caused by the Allee effect of prey species. If model (6) has a unique limit cycle, its asymptotic behavior is determined by the stability of £0, £i, E2, E3 and the limit cycle. Specifically, if E\ and E2 are saddle points and E3 is unstable, then there exists a stable limit cycle. The stable manifold of f-i, which is the threshold induced by the Allee effect of prey, divides the first quadrant into two parts, any positive orbit between the stable limit cycle and the stable manifold of £i converges to the limit cycle as t tends to infinity. Positive orbits outside the stable manifold of E^ converge to £0 as t tends to infinity, i.e., prey and predators cannot coexist if the initial population densities of prey and predators lies outside the stable manifold of E^ (see Fig. 3b). If E^ and £2 are saddle points and £3 is stable, there exists an unstable limit cycle. Any positive orbit outside the unstable limit cycle except the stable manifold of Ex converges to £0 as t tends to infinity (see Fig. 2c). Compared with model (1) which does not include Allee effect, we can say that the Allee effect of prey species may be a destabilizing force in the predator-prey system when the handling time of predators is relatively short. Next, we investigate the global bifurcation of system (4) by means of numerical simulations. 4. Global bifurcation analysis In this section, we provide the global bifurcation analysis of model (4) by means of the software packages PPLANE8 and MATCONT [26]. We will depict and discuss the heteroclinic cycles occurring in this model. 4.1. Local bifurcation analysis Previous calculations show that the handling time of predators h\ and the Allee effect of prey species A\ have an important impact on the dynamical behaviors of model (4), so we concentrate on the effect caused by the changes of A\ and h\. Corresponding to parameter values in Figs. 1-4, we fix b = 0.5000, = 0.2500, d2 = 0.1000, a = 0.0040, s = 0.0500, C\ = 0.2000 in model (4) (we always keep four decimal places for a real number in this paper). When A\ and vary continuously from zero, we obtain two-parameter bifurcation diagram of model (4), which is shown in Fig. 5). It can been seen that there exists a Bogdanov-Takens (or double-zero) point BT and a generalized Hopf bifurcation point GH (also called Bautin point). At the generalized Hopf point GH, an equilibrium with an unstable limit cycle H* turns into an equilibrium with a stable limit cycle H~. When A, > (Vb - v7^)27« = 10.7233, i.e., in region VII, the Allee effect of prey species is very strong, then only the trivial equilibrium £0 = (0,0) is biologically relevant equilibrium, which is globally asymptotically stable. All positive orbits converge to this equilibrium £0 as t tends to infinity, that is to say, both predators and prey cannot exist in this region. A typical phase portrait in this region is qualita- J. Zu, M. Mimura/Applied Mathematics and Computation 217 (2010) 3542-3556 3553 tively similar to Fig. la. When/^ = 10.7233, there is a transcritical bifurcation IQ occurring. This bifurcation is also called 'an extinction threshold for prey'. Further, when _05___0.1 10 7233 0.25/(0.1/(0.01 -0.005/;,)) +0.004 0.01 - 0.005/i, < ]< ' boundary equilibriums E, and E2 are also biologically relevant. If 0 < /i, < 1.2275, that is, the handling time of predator is relatively short, i.e., in region VI, then E0 is the globally stable equilibrium, E, is unstable and E2 is a saddle point. A typical phase portrait in this region is qualitatively similar to Fig. 2a. However, if h\ > 1.2275, that is, the handling time of predator is relatively long, i.e., in region VIII, then both equilibriums E0 and E2 are asymptotically stable, E, is a saddle point. The stable manifold of which is the threshold induced by the Allee effect of prey population, divides the first quadrant into two parts, one is the attraction region of E0, the other is the attraction region of E2. A typical phase portrait in region VIII is qualitatively similar to Fig. lb. Therefore, if the Allee effect of prey is relatively strong and the handling time of predator is relatively short, then prey population will become extinct even if the initial population density of prey is abundant. When 0.5 0.1 0.25/(0.1/(0.01 - 0.005/i,))+0.004 0.01 -0.005/!, ' there is another transcritical bifurcation TC2 occurring. This bifurcation is also called 'a predator invasion boundary'. When 0.5 0.1 0