Bifurcation analysis of the predator-prey model with the Allee effect in the predator.

The use of predator-prey models in theoretical ecology has a long history, and the model equations have largely evolved since the original Lotka-Volterra system towards more realistic descriptions of the processes of predation, reproduction and mortality. One important aspect is the recognition of the fact that the growth of a population can be subject to an Allee effect, where the per capita growth rate increases with the population density. Including an Allee effect has been shown to fundamentally change predator-prey dynamics and strongly impact species persistence, but previous studies mostly focused on scenarios of an Allee effect in the prey population. Here we explore a predator-prey model with an ecologically important case of the Allee effect in the predator population where it occurs in the numerical response of predator without affecting its functional response. Biologically, this can result from various scenarios such as a lack of mating partners, sperm limitation and cooperative breeding mechanisms, among others. Unlike previous studies, we consider here a generic mathematical formulation of the Allee effect without specifying a concrete parameterisation of the functional form, and analyse the possible local bifurcations in the system. Further, we explore the global bifurcation structure of the model and its possible dynamical regimes for three different concrete parameterisations of the Allee effect. The model possesses a complex bifurcation structure: there can be multiple coexistence states including two stable limit cycles. Inclusion of the Allee effect in the predator generally has a destabilising effect on the coexistence equilibrium. We also show that regardless of the parametrisation of the Allee effect, enrichment of the environment will eventually result in extinction of the predator population.

Introduction

Modelling predator–prey interactions has always been a mainstream area in mathematical biology and theoretical ecology. Our models have evolved tremendously since the famous Lotka–Volterra system, with one realistic modification being the introduction of non-monotonous per capita growth rates to the interacting species, as opposed to the monotonically decreasing per capita growth rate seen in the logistic equation. For instance, it is currently well recognised that the growth of natural populations can be subjected to the so-called Allee effect, where the per capita growth rate increases at low species densities (Courchamp et al. 2008; Fowler and Ruxton 2002). The Allee effect can emerge at the population level due to a variety of mechanisms including enhancement in foraging efficiency, reproductive facilitation, collective defense and the modification of environmental conditions by organisms (Berec et al. 2007; Courchamp et al. 2008; Fowler and Ruxton 2002). There exist two types of Allee effect: weak and strong Allee effects. The weak Allee effect describes situations in which the per capita growth rate is increasing at small densities, but which nonetheless remains positive for low, nonzero population densities, while a strong Allee effect is characterized by a negative population growth at low densities since reproduction cannot compensate mortality rate. It has been demonstrated that including the Allee effect in predator–prey models has a strong impact on dynamics, in particular promoting population collapse and a further species extinction (Boukal et al. 2007; Hilker 2010; Lewis and Kareiva 1993; Morozov et al. 2006; Sen et al. 2012). In previous theoretical works, however, the main focus has been the scenario where there is an Allee effect in the growth rate of the prey rather than that of the predator. The scenarios where the predator growth is subject to the Allee effect are explored in the literature only partially. The aim of this paper is to contribute to bridging the gap.

The existing literature on the Allee effect in predators is scarce, and mainly focused on foraging facilitation among predators which occurs as a result of cooperative hunting (Alves and Hilker 2017; Berec 2010; Cosner et al. 1999; Sen et al. 2019). Mathematically, this implies that the functional response of the predator is an increasing function of the predator density. In particular, it was shown that it might be detrimental for cooperative hunters to be too efficient in catching prey since this may cause resource over-exploitation and eventual extinction of the predator (Alves and Hilker 2017; Sen et al. 2019). On the other hand, the Allee effect can occur in predators due to other mechanisms such as low fertilization efficiency, a lack of mating partners, sperm limitation and cooperative breeding mechanisms (Berec et al. 2007; Courchamp et al. 2008; Dennis 1989). From the modelling point of view, including an Allee effect in this case should affect the numerical response of the predator, since the food conversion efficiency becomes an increasing function of predator density, while the functional response remains unchanged. As such, the model properties and ecological predictions will be different compared to the case of the foraging facilitation scenario. Some studies have considered the Allee effect in predators due to non-foraging mechanisms, but none of them have been studied exhaustively in terms of the bifurcation structure, possible dynamical regimes and the role of parameterisations of the Allee effect in the model equations (Costa and dos Anjos 2018; Zhou et al. 2005). The latter problem may be a general issue in ecological modelling and is related to so-called structural sensitivity, which is briefly described below.

In many ecological models, predator–prey systems in particular, there is often an uncertainty regarding which precise mathematical formulation of the model functions we need to implement in the model equations (Adamson and Morozov 2013). It is often impossible to determine which particular function we need to use in the model equations to describe predation, growth, mortality, competition, etc. Several parameterisations can fit available empirical data well, and different mathematical formulations can have a valid biological rationale (Flora et al. 2011). Furthermore, implementation of close mathematical functions (both in terms of functional forms and their derivatives) in the same predator–prey model may result in different outcomes, in particular in topologically distinct bifurcation structures yielding different dynamical regimes (Adamson and Morozov 2013, 2014a). This property is called the structural sensitivity of biological models (Adamson and Morozov 2013, 2014a, b; Seo and Wolkowicz 2018). Structural sensitivity may cause major problems in terms of generality of results obtained using specific concrete formulations of model functions such as growth rates or functional responses (Adamson and Morozov 2014a; Aldebert et al. 2019). A possible way to address structural sensitivity is to allow for an unspecified formulation of some functions in the model equations with other functions being fixed, an approach is known as partially specified modelling (Wood and Thomas 1999). Implementation of the partially specified models approach is especially relevant for systems with the Allee effect in predators since this phenomenon is often caused by a variety of mechanisms, and is thus hard to describe by a single universal functional relation (Courchamp et al. 2008). Moreover, the Allee effect can depend on the spatial scale of modelling, in which case the use of a single specific mathematical formulation for the dependence of the numerical response on the overall predator density is highly questionable (Courchamp et al. 2008).

In this paper we explore a predator–prey model with an Allee effect in the predator which affects the numerical response of the predator without affecting its functional response. We consider a partially specified model, where the mathematical formulation of a strong Allee effect has only a few generic constraints to its shape. We explore the bifurcation structure of the model including saddle-node, Hopf, generalised Hopf and Bogdanov–Takens bifurcations of co-dimensions two and three. Then we construct and compare full bifurcation portraits obtained for three possible parameterisations of the Allee effect: the hyperbolic (Monod), exponential (Ivlev) and trigonometric formulations. We demonstrate that the model may exhibit structural sensitivity with respect to parameterisation of the Allee effect function. We find that adding the Allee effect results in emergence of multiple non-trival attractors in the system which can potentially explain some empirically observed alternative states in ecosystems. We argue that the Allee effect in the predator growth has a large destabilising effect on population dynamics, which has been somehow neglected previously.

Model formulation and biological rationale

We consider a Gause type prey-predator ODE model with a specialist predator (Bazykin 1998; Hsu et al. 2001; Kuang and Freedman 1988; Turchin 2003). The model equations read as follows where N and P are the population densities of prey and predator, respectively, at time T. Note that in the above model, N and P can be also interpreted as the density of a resource and its consumer, respectively. In other words, the ‘predator’ in the current model can be a herbivore consuming an autotroph.

The function f(N) is the per capita growth rate of the prey which we consider here to be logistic, i.e., and is the intrinsic death rate of the predator which is assumed to be constant. Functional response of predator (the rate of food consumption per predator) which we consider here to be of Holling type II and we use the following parametrisation of g(N) known as the Holling disk equation (Jost et al. 1999)In this model, we incorporate the Allee effect in the numerical response of the predator by assuming that its food conversion efficiency is a function of predator density. This is different from previous models where the Allee effect was also included in the functional response of the predator (Alves and Hilker 2017; Cosner et al. 1999; Sen et al. 2019). The maximum food conversion coefficient is given by and this value is reached at high P. For simplicity, we neglect direct competitive effects and interference within the predator population. We assume that the reduction of the overall growth rate at high predator densities occurs solely due to over-exploitation of food, i.e., due to a decrease in N. At low predator density, the per capita reproduction rate becomes smaller, which is described by the function .

Biologically, inclusion of the Allee effect in the numerical response of the predator only (and not in its functional response) can mimic a multitude of scenarios/mechanisms. A major scenario is that at a low population size it is hard to find a suitable and receptive mate and this largely reduces the reproduction rate. Mate-finding Allee effects have been found in a large number of species ranging from small insects to large birds in terrestrial ecosystems and from zooplankton to whales in the sea (Courchamp et al. 2008). Due to space limitation, here we can only list a small number of empirical examples from the literature. Namely, Mate-finding Allee effects were found in populations of flour beetles (Allee et al. 1949), muskrats (Errington 1940), whales (Ton 1948), box turtles (Mosimann 1958), condors (Mertz 1971), various zooplankton copepods (Gerritsen 1980; Kiørboe 2006), piping plovers (Strauss 1991), primates (Dobson and Lees 1989), various parasitoids (Hopper and Roush 1993), whooping cranes (Wells et al. 1998), pelagic fish (Liermann and Hilborn 2001), and elk (Larkin et al. 2002).

Another related mechanism impeding reproduction at low numbers is sperm limitation, where fertilisation of eggs requires a sufficient amount of sperm. The biological rationale is that a female needs to find a male of an optimal size, or she needs to have a sufficiently large number of males in the surroundings. The corresponding empirical examples include the blue crab (Hines et al. 2003) and the Caribbean spiny lobster (MacDiarmid and Butler 1999), where sperm limitation occurs due to excessive fishery which targets large-size males and selectively removes them from the population. The Allee effect is also possible due to low fertilization efficiency, which is observed in sessile (e.g. corals) or semi-sessile organisms (e.g. echinoderms, polychaete worms). For such broadcast spawners, at low population density, the probability of meeting of sperm and egg in water column becomes largely reduced (Aronson and Precht 2001; Courchamp et al. 2008). Our model would also mimic the Allee effect due to reproductive facilitation mechanisms. According to this scenario, only the presence of a sufficiently large number of conspecifics in the neighborhood can trigger reproduction instincts of individuals. This is the case, for example, of the flour beetle (Allee et al. 1949), queen conch (Stoner and Ray 1993), snails, and lizards (Crews and Fitzgerald 1980; Vernon 1995). Finally, our model can mimic the Allee effect due to cooperative breeding mechanisms, where at low numbers, breeding groups are less efficient in reproduction and raising young animals. Cooperative breeding is well-known in a large number of bird species (Koenig and Dickinson 2004).

Mathematically, the function modelling the Allee effect in th predator is considered to possess the following properties: Next we reduce the number of parameters in the model by non-dimensionalisation and introduce the following non-dimensional variables , and , we can transform the equations (1) to

, since at very low densities the population cannot reproduce due to lack of mating opportunities;

since represents the proportion of the maximal possible conversion rate e;

is an increasing function of P, so for all , thus we do not include effects of intraspecific competition;

for all which signifies that the increase in the reproductive ability (population fitness), while the population size P is being increased, is monotonically decelerating;

for large P (we assume the highest efficiency of reproduction at large densities).

with the following positive dimensionless parameters , and .

In the above model, the function is transformed into a dimensionless function h(y) with the same constraints as are imposed on . As we mentioned in the Introduction, we will explore the basic properties of the model for an arbitrary mathematical formulation of h(y) (model equilibria, stability, possible generic bifurcation, etc), i.e., considering the above system as a partially specified model. We will also consider some concrete parameterisations of h(y) such as (Monod parametrisation), (Ivlev parametrisation) and (hyperbolic tangent parameterisation) to construct a full bifurcation portrait and explore the sensitivity of the model dynamics to mathematical formulation of the Allee effect. Finally, we verify how sensitive the model is with respect to small perturbations of h(y) which still preserve assumptions (A1)–(A5).

Model equilibria and their stability

Possible equilibria in the system

We start our investigation by exploring the number and the location of system equilibria for an arbitrary formulation of the Allee effect h(y). It is easy to see that model (2) always has one trivial equilibrium point and one axial (predator-free) equilibrium point .

An interior equilibrium point will be a point of intersection of the following two non-trivial nullclines in the interior of first quadrant For the feasibility of , we must have (see Eq. (3a)) and from this condition we can verify that . We solve (3b) for x to obtain the equation for the non-trivial predator nullclineFor the feasibility of , we must have . Since the Allee effect function h(y) is bounded by 1, from we find with . We differentiate Eq. (4) with respect to y to obtainFor the second derivative of the predator nullcline we haveas and . From the above we derive that the predator nullcline (3b) is strictly decreasing as a function of x and it is always convex. Furthermore, the curve lies in the region where and . Taking into account the above properties, a possible shape of the predator nullcline is shown in Fig. 1; the dashed lines represents the vertical and the horizontal asymptotes.

Fig. 1

Qualitative behviour of the predator nullcline in model (2)

From the geometric properties of the nullclines one can see that there will be at most two points of intersection (between the non-trivial nullclines in the interior of the first quadrant) and so there can be at most two interior equilibria. An example of intersection of the model nullclines for the parameterisation of h(y) given by is shown in Fig. 2. One can see that a gradual increase in (which defines characteristic predator densities at which the Allee effect has a pronounced strength) from small to large values results a saddle-node bifurcation which is described in detail in the next sections. Note that this property is observed for the other two parametrisations of h(y) considered. Clearly, in the absence of the Allee effect, only one non-trivial equilibrium is possible, corresponding to the intersection of the vertical line and the prey nullcline (3a).

Fig. 2

Relative position of nullclines in model (2) constructed for the Allee effect parametrisation given by the Monod function with (no equilibrium), (a single equibrium), and (two equilibria points): one is a saddle point (black dot), the other is a topological focus (red dot). Other parameters are , and

Stability of equilibria

Here we explore the stability for the equilibria of model (2). The following proposition defines the stability of the axial equilibria.

Proposition 1

For any choice of h(y) satisfying assumptions (A1)-(A5)

The trivial equilibrium point is a saddle;

The axial equilibrium point is locally asymptotically stable.

Proof

The Jacobian matrix of model (2) at any point is given by

The eigenvalues of the Jacobian matrix at are 1 and . Therefore it is a saddle point irrespective of the choice of h(y), having a stable manifold along y-axis and an unstable manifold along x-axis.

The axial equilibrium point is locally asymptotically stable (a stable node) as the eigenvalues of the Jacobian matrix are and for any choice of h(y).

An important conclusion is that, in the presence of an Allee effect in the predator, achieving a very low population densities by the predator will result in its eventual extinction, so the Allee effect is strong.

Next we explore the stability of the interior equilibria. As follows from the previous section, model (2) admits at most two interior equilibrium points which we denote by and such that . The Jacobian matrix evaluated at can be expressed aswhere we have used and . Since and are smooth functions, we can differentiate both expressions (3) to obtainwhere and are the tangent lines to the nullclines and , respectively. We substitute the above expressions into the Jacobian matrixTherefore, for the determinant of the Jacobian we obtainNow and as . Substituting the above derivatives we have where is the x-coordinate of the point where the nullcline attains its maximum within the first quadrant (for ).

Depending upon the positions of two points of intersections between the two nontrivial nullclines (cf. (3)) with respect to the point of maximum on the prey nullcline, we can consider following two cases,

Case:1  and Case:2  .

Case:1  Suppose and are angles made by the tangents to and at . Then from Fig. 3 one can see that which impliesTherefore from the above inequality and expression (8) we get . Hence is a saddle point.

Fig. 3

Determining the type and stability of interior equilibria in model (2). For detail see the text

We can proceed in a similar fashion and use the fact that (this is not shown in Fig.  3 for brevity)and to prove that , i.e., is not a saddle point. The stability of is determined by the sign of . Note that Hence is locally asymptotically stable ifCase:2 In this case we can also prove that is a saddle point proceeding in a similar manner as above. For the stability of we have from (9a), . Hence we get from (8), . Also As , hence is unstable.

To conclude, the interior equilibrium is always a saddle point, whereas is a topological focus which depending on parameters can be either stable or unstable.

Local bifurcations in the model

Here we consider possible local bifurcations in model (2).

Saddle-node bifurcation

Suppose is the point at which two non-trivial nullclines touch each other in the first quadrant when a bifurcation parameter of the model is being varied. The slope of the tangents to the curves at become equal. This signifies thatand hence we haveIn this case and becomes a non-hyperbolic equilibrium point. This situation corresponds to a saddle-node bifurcation in the model. We explore this bifurcation in more detail.

As an example, we consider m to be the bifurcation parameter and denote by the bifurcation point. The eigenvectors of both the matrix and its transpose corresponding to the zero eigenvalue are, respectively given by and . We need to check the transversality conditions for a saddle-node bifurcation (Perko 2013). We denote and we further follow the same notation of (Perko 2013) to obtainas . Hence the transversality conditions are always satisfied and variation of m results in a saddle-node bifurcation. Similar results can be obtained by varying other model parameters.

Hopf bifurcation

In the previous subsection we show that the two interior equilibrium points are generated through a saddle-node bifurcation. The non-saddle interior equilibrium () can be stable or unstable depending on model parameters. It loses its stability when the sign of the trace of the Jacobian matrix has changed through zero (from negative to positive) via a Hopf bifurcation. In this section we show that system (2) undergoes a Hopf bifurcation when a model parameter is varied. Here we choose as a bifurcation parameter. The Jacobian matrix at is given byNow let us assume that . This is an implicit expression for as the components of the equilibrium point containg as well. Now we assume that the following three conditions are satisfied at , Then loses its stability through a Hopf bifurcation at .

,

,

If is the complex eigenvalue of then

Assuming is real part of a complex eigenvalue of , we can write,Now is equal to zero when . Differentiating with respect to we find thatNow as satisfies (3a) we have,Finally using the fact at and above result in (17) we getThe above expression should be checked for the given mathematical formulation of the Allee effect h(y). In particular, we have numerically verified that this quantity is non-zero at the Hopf bifurcation threshold for the parameterisations considered here: the Monod, Ivlev and trigonometric functions.

Generalized Hopf (Bautin) bifurcation

In this subsection we consider a co-dimension two bifurcation called a Bautin or generalized Hopf (GH) bifurcation. This bifurcation occurs when the interior non-saddle equilibrium has purely imaginary eigenvalues and the first Liapunov number becomes zero. We consider and m as bifurcation parameters. Therefore in –m parametric plane, there is a critical point which lies on the Hopf bifurcation curve. In the next proposition we will show that the model undergoes a GH bifurcation by choosing and m as bifurcation parameters.

Proposition 2

Model (2) undergoes a Bautin (generalized Hopf) bifurcation around the interior equilibrium point at the bifurcation threshold whenever the following conditions hold where l is the first Liapunov number.

,

,

,

See supplementary material SM1.

Examples of the above type of bifurcation for several parameterisations of h(y) are provided in Sect. 5.

Bogdanov–Takens bifurcation

Another type of co-dimension two local bifurcation observed in model (2) is a Bogdanov–Takens (BT) bifurcation. In a two dimensional parametric plane, this bifurcation occurs at a point where a Hopf bifurcation curve meets a saddle-node bifurcation curve tangentially. In the previous subsection, we chose as the bifurcation parameter for the Hopf bifurcation and m for the saddle-node bifurcation. Therefore we will consider and m as bifurcation parameters for the BT bifurcation and suppose that model (2) exhibits a BT bifurcation at and the parametric thresholds are denoted by . From the general bifurcation theory (Perko 2013) it is known that satisfies the equations of nullclines (3) and also the Jacobian matrix is similar to at for the parameter threshold . The following proposition provides the conditions for model (2) to undergo a Bogdanov–Takens bifurcation.

Proposition 3

If we choose and m as bifurcation parameters, then system (2) undergoes a Bogdanov–Takens bifurcation around the interior equilibrium point whenever the following conditions hold

,

.

See supplementary material SM 2 for detail.

We found that model (2) may undergo a Bogdanov–Takens bifurcation of either co-dimension 2 or co-dimension 3. The latter requires an extra condition of degeneracy given in the supplementary material. Note that a co-dimension 3 Bogdanov–Takens bifurcation, if it exists, should be of the type involving a double equilibrium point (Dumortier et al. 1987). Indeed, the other type of this bifurcation—known as the cusp—would require a triple equilibrium point which is impossible for this model as shown in Sect. 3. For the same reason, a co-dimension 4 Bogdanov–Takens bifurcation is impossible in this system. Examples of Bogdanov–Takens bifurcation of co-dimension 2 and 3 for particular parameterisations of h(y) are provided in the next section.

Parametric diagrams and phase portraits

In this section, we construct global parametric diagrams for the considered model for three different mathematical formulations of the Allee effect h(y) given by the Monod, Ivlev and trigonometric tangent functions. Note that all of them satisfy assumptions (A1)–(A5). Examples of all three curves constructed for are shown in Fig.  4. Note that for the plotted functions the initial slopes and their asymptotic values for large y are the same.

Fig. 4

Graphs of parameterisations of the Allee effect in the predator given by the Monod (), the Ivlev () and the trigonometric () functions;

The model contains 4 parameters, so it is convenient to present our results in a 3 dimensional parametric space and then explore the alteration to the portrait by varying a fourth parameter. We construct portraits in the space with a further variation of . For all considered formulations of h(y), the parameter can be interpreted as the intensity of the Allee effect. In particular, in the case where vanishes the system becomes the classical Rosenzweig–MacArthur predator–prey model.

Examples of parametric portraits for the three functional forms of h(y) are given in Fig. 5, in each case is kept fixed. We show the skeletons of the parametric portraits given by local bifurcations: to avoid overloading the diagram, we do not include the non-local bifurcations which are shown in the corresponding cross sections in next figures. The saddle-node bifurcation surface is denoted by the blue curves. These curves show intersections of the saddle-node bifurcation surface with boundaries of the parameteric diagram. The intersection of the Hopf bifurcation surface with the boundaries is denoted by red curves. The saddle-node and Hopf surfaces intersect along the Bogdanov–Takens bifurcation curve which consists of green and magenta coloured parts: the magenta colour corresponds to a Bogdanov–Takens bifurcation of codimension 2 with a positive product of the state variables in the normal form (see SM 2 for detail) and is denoted as ; the green part of the curve gives Bogdanov–Takens bifurcation of codimension 2 with this product having negative sign and is denoted as . The black curve represents the location of Generalised Hopf points on the Hopf bifurcation surface. This curve emerges from the point of Bogdanov–Takens bifurcation of codimension 3 (denoted as ). From comparison of the diagrams in Fig. 5, we conclude that the global bifurcation structure in the parametric space remains the same topologically for all three formulations of h(y).

Fig. 5

Three dimensional diagram of model (2) for the Allee effect parameterised by: a the Monod response ; b the Ivlev response ; c the trigonometric tangent function . In each diagram, . The explanation of the surfaces and curves is in the text

To better understand the parametric structure and feasible phase portraits in the model, we explored two-dimensional cross sectional diagrams for a constant and . In the main text, we present the diagrams for the Monod formulation of h(y). The diagrams for the other functional forms of h(y) are shown in the supplementary material (SM 3). An example of a diagram constructed for above the point is shown in Fig.  6a; the other two diagrams in the same figure are constructed for below the point. Figure 6b describes the situation where the plane does not intersect the GH bifurcation curve, the opposite case is shown in Fig. 6c. The corresponding phase portraits of the model are given in Fig. 7.

Fig. 6

a–c Examples of bifurcation diagrams in –m parametric plane for the Monod parameterisation of the Allee effect in predator plotted for above and below the point shown in Fig. 5a. The parameter values are: a ; b and c ; in all cases . d Two dimensional diagram constructed for and . In each diagram, the dark blue curve is a saddle-node bifurcation curve. The red curve is a Hopf bifurcation curve. The green curve is the curve of a fold bifurcation of limit cycles. The cyan curve describes a homoclinic bifurcation. The meaning of the regions is explained in Fig.7 and in the main text

Fig. 7

Phase portraits for the model with the Monod form of h(y). a–f are constructed for , ; the other parameters are a , , region (); b , , region (); c , , region (); d , , region (); e , , region (); f , , region (). g is plotted for ; ; ; , region (); h is plotted for , , , , region ()

From Fig. 6a–c one can see that for large values of m and (region ) there are no coexistence equilibria in the system: the only (global) attractor is the state , where only the prey survives, whereas the predator goes to extinction. The corresponding phase portrait is shown in Fig. 7a.

Reduction in the strength of the Allee effect (small and high rates of mortality m) results in the emergence of a pair of equilibrium points: a saddle and a node (region ). The non-saddle point is only locally stable: its basin of attraction is limited by that of the axial equilibrium point , which is shown in Fig. 7b. For large values of m (on the right hand side of a BT point), a decrease in will result in a saddle-node bifurcation where the non-saddle point will be unstable (region ). In this case, the global attractor will be the prey only state (the phase portrait is shown in Fig. 7d). The loss of stability of when crossing the Hopf bifurcation curve around the BT point depends on the sign of the BT point. For (Fig. 6a) transition from to occurs via region by crossing the homoclinic loop bifurcation curve. A locally stable interior equilibrium bcomes surrounded by an unstable cycle which forms its basin of attraction (Fig. 7e). All trajectories starting outside this cycle will be attracted to the prey only state . The transition from region to region occurs via a homoclinic loop bifurcation. For (Fig. 6b), the transition from to occurs via region by crossing a supercritical Hopf bifurcation curve. In region , an unstable internal equilibrium is surrounded by a stable limit cycle (Fig. 7c). One can see that in Fig. 6a, for smaller m, a decrease in from region results in a fold bifurcation of limit cycles. In region we have two limit cycles: the inner cycle is stable, the outer cycle is unstable (Fig.7f). The outer cycle forms the boundary of the basin of attraction for the state .

The diagram in Fig.6c is more complicated as compared to Fig.6a,b. In particular, a new region emerges, where three limit cycles can coexist: the inner limit cycle is stable, the middle cycle is unstable and the outer one is stable. The corresponding portrait is shown in Fig.7h.

For the Ivlev and trigonometric formulations of h(y), the bifurcation diagrams in the –m plane constructed for above and below the point are topologically equivalent (see supplementary material SM3). However, for fixed and the location of the bifurcation curves as well as the types of bifurcation (e.g., versus type of bifurcation) in the –m plane may be substantially different, especially when comparing the Monod parametrisation with the other two functional forms. This indicates sensitivity of the model to the functional form of h(y). We explore the structural sensitivity in more detail in the next section.

Consider now variation of the fourth model parameter . A decrease in results in a shift of the saddle-node and Hopf bifurcation surfaces in the space closer to the plane. The length of the curve of Bogdanov–Takens points (the intersection between the saddle-node and Hopf bifurcation) will be shortened and it moves upwards on the saddle-node bifurcation surface. The codimension 3 Bogdanov–Takens bifurcation is still observed. The above properties hold true for all three functional forms of h(y) considered.

Finally, we consider the case where gradually decreases and the other parameters are kept fixed. This corresponds to the ecologically important scenario in which the environment undergoes gradual eutrophication: an increase in the carrying capacity K in the original model (1) corresponds to a proportional decrease in in the dimensionless model. An example of a diagram in the plane is constructed for the Monod parametrisation of h(y) (Fig. 6d). The parameter regions in the diagram have the same meanings as in Fig. 6a-c. A gradual decrease of results in destabilisation of the coexistence equilibrium and a further collapse of the population of predator (transition from region to region ). Thus in a eutrophic environment, the only stable equilibrium in the model is being a predator free equilibrium. Note that in this diagram we have a new region denoted as in which a locally stable equilibrium is surrounded by two limit cycles: the inner cycle is unstable whereas the outer one is stable. The corresponding phase portrait is shown in Fig. 7g. Similar bifurcation behavior is observed for the other two parameterisations of h(y).

Structural sensitivity of the model

An important part of our investigation is exploring the dependence of model behaviour on the choice of parametrisation of the Allee effect given by h(y). In the previous section, we show that the skeleton of the bifurcation diagram is topologically robust to the mathematical shape of h(y) (Fig. 5) when we use three different parameterisations given by the Monod, Ivlev and hyperbolic tangent functions. The relative positions of bifurcation surfaces and possible dynamical regimes remain the same. On the other hand, we also find that for a fixed set of parameters the parametric diagrams can differ considerably, even for close functions h(y). This property is known as structural sensitivity of biological models.

Structural sensitivity of models is an important issue in ecological modelling with a large number of insightful examples provided (Adamson and Morozov 2013, 2014a, b; Aldebert et al. 2019; Flora et al. 2011; Wood and Thomas 1999). The biological rationale behind this idea is that any realistic ecological dependence h(y) can be some combination of the three parameterisations considered above, or, more generally, any other mathematical functions. Here we explore the sensitivity of model predictions to a small variation of h(y) by considering the entire set of functions which satisfy assumptions (A1)–(A5), and we use the methodology from (Adamson and Morozov 2013, 2014a). Here we investigate the sensitivity of the stability of the interior non-saddle stationary state to small but finite perturbation of h(y) starting from one of the functional forms considered above. Mathematically, we consider the following situation.

Let the base function h(y) be of Monod type. We consider small deviations from the base function h(y) such that and the second derivative of is negative and bounded, , (note that the base function also satisfies this condition). For simplicity we slightly relax our assumption (A5) and allow , , although this fact is not crucial for the general outcome and conclusions. We conduct the sensitivity analysis in the same way as when investigating the role of the functional response of the predator in stability dynamics in the Rosenzweig–MacArthur model (Adamson and Morozov 2013). To quantify the structural sensitivity to the choice of we use the degree of structural sensitivity introduced in (Adamson and Morozov 2013). represents the probability that for two randomly chosen functions of , the coexistence stationary state will have the same stability properties (Adamson and Morozov 2013), to be stable or unstable. The maximal degree of sensitivity possible in the system is equal to one, signifying maximal uncertainty in the system

Figure 8 shows the degree of sensitivity as the parameter is varied. From the figure one can conclude that a decrease in eventually results in destabilization of the system for any parametrisation of : we gradually move from the region of stability to the region of instability shown in the figure. However, depending on the particular choice of , this destabilisation may occur within a wide range of , thus the system has a plasticity to resist destabilisation caused, for example, by eutrophication (a high value of K in the original system signifies a low in the dimensionless model). Another important observation is that the system exhibits large uncertainty even if the deviation from the base function is small (%). Considering the Ivlev and the hyperbolic tangent as the base function provides similar results.

Fig. 8

The degree of structural sensitivity in model (2) to variation of the parameterisation of the Allee effect in the predator h(y), shown for different . The -neighbourhood of h(y) is defined using the relative difference between the base function and its perturbations. Destabilisation of the equilibrium for the base function (Monod parameterisation) occurs at For further detail see the text and also (Adamson and Morozov 2013, 2014a). The other parameters are , , and

Discussion

The role of the Allee effect in population dynamics has been largely addressed in both empirical and theoretical literature. Surprisingly enough, there has been almost no thorough mathematical investigation into the bifurcation structure of any predator–prey model with an Allee effect in the predator, in particular, this concerns the realistic scenario, where the Allee effect is included in the numerical response of the predator without affecting its functional response (Berec et al. 2007; Courchamp et al. 2008; Dennis 1989). This is in a striking contrast to the situation with single species population models or classical predator–prey models with an Allee effect in the prey growth, which have been discussed in all detail and are now included in standard student textbooks in mathematical ecology (Kot 2001). The current study is intended to partially bridge the existing gap. Importantly, our results are not based on a particular mathematical formulation of the function describing the Allee effect, rather we consider various parameterisations of h(y) which satisfy only few qualitative constraints (A1)–(A5). We have also addressed (for the first time) the issue regarding the sensitivity of the model with respect to parameterisation of the Allee effect (known as the structural sensitivity).

In this section, we will mostly focus on ecological implications of the mathematical results obtained in the previous sections. We should stress, however, that it is crucial to define the way of how we should formally assess the consequences of the Allee effect on the population dynamics. Indeed, this question is far to be a trivial one since distinct paradigms exist in the literature (Alves and Hilker 2017; Berec et al. 2007; Courchamp et al. 2008; Dennis 1989). In fact, evaluation of the role of the Allee effect in population success should largely depend on the choice of the initial density, and this fact is still somehow disregarded in the literature. Indeed, consider two populations, where one possesses a self-accelerating per capita reproduction rate, and the other one which is characterised by a constant per capita growth rate: for simplicity we assume the growth rates of both populations at some low density to be the same. Then an increase in the population density would result in an increase in the per capita growth rate of the species with an Allee effect, and this will clearly indicate the benefits of possessing an Allee effect. For example, hunting cooperation is considered to be beneficial for predators up to certain level of population density (Alves and Hilker 2017). On the other hand, considering higher population densities as the starting point for comparison (e.g. population densities where the Allee effect is not pronounced) will show a different outcome. For the same initial per capita growth rates, the species without an Allee effect will be more advantageous since a sudden drop of the population size would not largely affect its per capita growth whereas the population with an Allee effect may exhibit a significant decline in reproduction rate with a threat of extinction. For example, low fertilization efficiency, a lack of mating partners, and sperm limitation are usually considered as negative and undesirable features for population persistence (Berec et al. 2007; Courchamp et al. 2008; Dennis 1989).

Arguably, for many species, their reproduction rate is often empirically estimated at densities which are away from the extinction threshold, where the Allee effect is not well-pronounced. As such, we suggest that in theoretical models the impact of the Allee effect on dynamics should be assessed via comparison with a scenario without an Allee effect, where for both scenarios per capita reproduction rates are assumed to be the same at some ‘safe’ densities. For the current theoretical study, this mathematically signifies that we need to compare the Allee effect in model (2) with the same model with since for large y the value of h(y) tends to unity, which corresponds to the classical Rosenzweig–MacArthur predator–prey model (Kot 2001). Note that the model with the Allee effect becomes the Rosenzweig–MacArthur model in the case .

Following the above philosophy, our first important conclusion is that introducing the Allee effect in predator’s numerical response generally acts as a destabilising factor of a stable coexistence of the prey and the predator in the case where the predator is specialist. Moreover, the Allee effect can result in extinction of predators regardless of initial density. This can be seen from the bifurcation diagrams, when the parameter , characterising the strength of the Allee effect, increases from small values , corresponding to the Rosenzweig–MacArthur predator–prey model, to some large values. Interestingly, destabilising influence of the Allee effect in the predator is observed even for a linear functional response as well (see supplementary material SM4 for detailed illustration). It is well-known that the coexistence equilibrium of the classical Rosenzweig–MacArthur predator–prey model with a linear functional response and a constant h(y) is globally stable. Introducing an Allee effect into the predator growth results in stability loss, with either generating sustained oscillations (in this case both prey and predator still persist in the system in an oscillatory mode) or leading to extinction of predator (via different scenarios) for any initial population density. Destabilisation of the system is facilitated with a pronounced saturation in the functional response of the predator (small ) and with an increase in the strength of the Allee effect—interpreted in the model as a gradual increase in . Biologically, this signifies that having density-dependent h(y) rather than a constant h(y) efficiency of predator impedes control over the prey population: the response of the predator to variation of the prey density in the system with the Allee effect becomes delayed and this allows the prey to escape from the control.

The predicted by the model destabilising role of the Allee effect demonstrates that cyclic population dynamics should occur in predator–prey systems more frequently as it was suggested earlier. For example, in the case of a Holling type I functional response, which is well-known to be stabilising, empirically observed oscillations of population numbers are usually attributed to factors as environmental/demographic noise, seasonal forcing, complex age structure of the population, complexity of the food webs, etc (see (Barraquand et al. 2017) and the references therein). On the other hand, an Allee effect can be an alternative explanation of oscillatory dynamics of a large number of case studies which occur in non-eutrophic environments with a stabilising functional response of the predator/consumer.

Our second important conclusion is that destabilisation of the coexistence equilibrium in the system with an Allee effect in the predator can occur not only via a supercritical Hopf bifurcation scenario (appearance of a small-amplitude stable limit cycle) but via a subcritical one. In the latter case, destabilisation of the equilibrium will leads to eventual extinction of the predator since the only possible attractor of the model is the state with only prey population being present (regime in the model, Fig. 7). On the other hand, eutrophication of the environment—which in the dimensionless model corresponds to a decrease in the parameter —would eventually result in extinction of the predator regardless of the scenario of stability loss of the equilibrium (supercritical and subcritical). Indeed, in the case of a supercritical Hopf bifurcation, the resultant predator–prey cycle grows in size and enters the basin of attraction of the prey-only equilibrium: mathematically this occurs via a homoclinical bifurcation.

Our results may have important implications for the biological control of the pests by predators and parasitoids. It has been reported that in a large number of cases, biological control agents have failed to get established even if under laboratory conditions they could survive by consuming target pest species (Bellows 2001; Roderick and Navajas 2003; Orr 2009; Bompard et al. 2013). A possible explanation is the presence of an Allee effect in the biological control agents which becomes more pronounced in the environment as in a lab. For example, the initial density of the predator can quickly fall below the critical threshold because of dispersal and diffusion. However, there can be a more complicated scenario where the initial density of the predator can be very high, but a pronounced Allee effect will still not allow a long-term persistence of species. Our model mathematically describes this as a globally unstable co-existence state (regime in the model, Fig. 7). As a conclusion, the choice of the appropriate species for an efficient biological control should be made carefully. For example, for parasitoids it is preferred to use haplodiploid species (i.e., where the males are haploid and females diploid) to alleviate the negative demographic consequences of mate-finding Allee effects, which is well-pronounced in diploid species (Hopper and Roush 1993; Bompard et al. 2013).

Thirdly, we found that inclusion of the Allee effect in the predator growth largely increases the complexity of the system as compared to the scenario without an Allee effect, i.e., the original Rosenzweig–MacArthur predator–prey model. One of the most interesting observation is the possibility of non-trivial multiple attractors, which can be interpreted as alternative ecosystem states. The model predicts two such regimes denoted by and , see Fig. 7. They are (i) a stable equilibrium coexisting with a stable limit cycle and (ii) two stable limit cycles, respectively. Note that in the literature, there is an ongoing debate on possible origins of alternative states generated through various ecological mechanisms (Schröder et al. 2005). In the simplest case, alternative stationary states are two contrasting equilibria (Scheffer et al. 2003; Scheffer and Carpenter 2003). However, more interesting patterns with non-equilibrium coexisting attractors have been reported in the literature as well. In particular, empirical observations demonstrate the possibility of alternative attractors where depending on initial condition the population dynamics can show either cyclic or stable coexistence scenario (Zamamiri et al. 2001). Other empirical studies demonstrate the possibility of coexisting cyclic oscillations with contrasting amplitude and periodicity (McCauley et al. 1999; Henson et al. 2002). Interestingly, in the Daphnia-algae predator–prey system reported in (McCauley et al. 1999), the observed coexisting cycles were considered to be a consequence of variability of the available food for zooplankton (e.g. a more complicated prey-dependent functional response of the predator or food-dependent efficiency of reproduction), whereas the density dependence of Daphnia’s vital rates was somewhat intentionally disregarded. On the other hand, the presence of an Allee effect (e.g. an emergent Allee effect is known to be present in Daphnia (de Roos et al. 2003)) can be arguably an alternative explanation of the co-existing cyclic behaviour.

Finally, we should stress that unlike the scenario with the Allee effect in prey—which is currently considered to be straightforward in the literature—including the Allee effect in predators can be somewhat tricky since this can be done either by modifying the functional response or the numeric response of the predator. Biologically, this signifies that distinct mechanisms of emerging the resultant demographic Allee effect—a decrease in reproduction at low population numbers—should be modelled in a different way. In particular, the mechanisms such as low fertilization efficiency, a lack of mating partners, sperm limitation, cooperative breeding and similar mechanisms (see empirical examples in Sect. 2) should be included in the numerical response of predator only, whereas collective exploitation of a resource such as cooperating hunting should be included in both functional and the numerical responses. It is rather surprising that this fact has not been largely explored in the literature yet.

As a first step to fill the existing gap, we compared the stability of two similar predator–prey systems: in one model the Allee effect was due to collective exploitation of resources (as in the study (Alves and Hilker 2017)) and in the other one the Allee effect was due to the lack of mating partners (for details see SM4). Note that unlike the mentioned work (Alves and Hilker 2017) we included saturation in the Allee effect in the functional response of predator (see the model equations in SM4). For simplicity, we considered the case where the functional response of the predator is of Holling type I. We found that a pronounced Allee effect (large saturation in h(y)) has somewhat different consequences for the two systems. Using bifurcation diagrams, we compared the stability regions for both models. We found that an Allee effect due to collective exploitation of resources considered in (Alves and Hilker 2017) facilitates persistence of the predator as compared to the scenario where the Allee effect is included only in numerical response (e.g. due to mate-finding) since the stability region for letter model is smaller in size. However, a more detailed comparison of the two mentioned approaches to modelling the Allee effect in predator should be an interesting separate study.

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 572 KB)