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. Supplementary Information The online version supplementary material available at 10.1007/s00285-021-01707-x.


Supplementary material 1
Proof of Proposition 2: Suppose that the interior equilibrium undergoes a generalised Hopf (GH) bifurcation atÊ = (x,ŷ) for the parameter threshold (β GH , m GH ). SinceÊ is a non-trivial equilibrium point then it must satisfies (3). The Jacobian matrix atÊ is given by . Now using the condition (GH2), we find from the above Jacobian β GH = 1 − 2x + αŷh (ŷ). Now to find the first Lyapunov number atÊ, we translateÊ to origin by using the transformation u = x −x, v = y −ŷ and we getu = au + bv + P (u, v), v = cu + dv + Q (u, v), , P (u, v) and Q(u, v) are analytic function is given by where a ij and b ij are obtained from the following relations We calculate the coefficients a ij , b ij and a, b, c, d of the first Lyapunov number atÊ: After some algebraic calculation we get (assuming that h(y) is as smooth as that h (y) exist) .
Due to algebraic complexity we are unable to show that l((Ê; β = β GH , m = m GH )) = 0, however it can be verified numerically. Thus we can conclude that the system (2) admits GH bifurcation for the parameter threshold (β = β GH , m = m GH ) atÊ.

Supplementary material 2
Proof of Proposition 3: The conditions (BT1) and (BT2) are equivalent to the following ones It is difficult to find the explicit expressions of the bifurcation parameters for the BT bifurcation thresholds since the above equations contain unknown function h(y) and the coordinates of the equilibrium points contain the parameters implicitly. However, we can define the bifurcation thresholds as follows: Let us consider the small perturbation around the bifurcation threshold (β BT , m BT ) by (β BT + λ 1 , m BT + λ 2 ), where λ i , i = 1, 2 are sufficiently small, and substituting it into the system we get Now we shift the equilibriumĒ = (x,ȳ) to the origin by the coordinate transform X = x −x, Y = y −ȳ and substitute it in model (2) to obtaiṅ and H 1 , H 2 are power series of X and Y with terms X i Y j , i + j ≥ 3. Now by using affine transformation X 1 = X, X 2 =ãX +bY , the system transferred tȯ where, k 0 =ãF 10 +bF 20 , k 10 =bc −ãd, k 01 =ã +d, andH 1 ,H 2 are power series of X 1 and X 2 with terms X i 1 X j 2 , i + j ≥ 3. Now using the C ∞ change of coordinates in the neighbourhood of origin by the following transformation the system reduces tov and R 1 , R 2 are power series of v 1 and v 2 with terms v i 1 v j 2 , i + j ≥ 3. Now again using the C ∞ change of coordinates in the neighbourhood of origin by the following transformation w 20 = t 20 − t 01 s 20 − t 11 s 10 , w 11 = t 11 + 2s 20 , andR 1 ,R 2 are power series of z 1 and z 2 with terms z i 1 z j 2 , i + j ≥ 3. Finally, using the C ∞ change of coordinates in the neighbourhood of origin by the following transformation n 1 = z 1 , n 2 = z 2 +R 1 (z 1 , z 2 ), the system changes toṅ 1 = n 2 , n 2 = P (n 1 , λ 1 , λ 2 ) + n 2 φ(n 1 , λ 1 , λ 2 ) + n 2 2 ξ(n 1 , n 2 , λ 1 , λ 2 ), with P (n 1 , λ 1 , λ 2 ) = w 0 + w 10 n 1 + w 20 n 2 1 + P 1 (n 1 ), φ(n 1 , λ 1 , λ 2 ) = w 01 + w 11 n 1 + φ 1 (n 1 ), ξ(n 1 , n 2 , λ 1 , λ 2 ) = ξ 1 (n 1 , n 2 ).
Here P 1 , φ 1 , ξ 1 are the power series in terms of n i 1 , i ≥ 3, n i 1 , i ≥ 2 and n i 1 n j 2 , i+j ≥ 1 respectively. Now applying Malgrange preparation theorem to the function P , we get P (n 1 , λ 1 , λ 2 ) = n 2 1 + w 10 w 20 n 1 + w 0 w 20 B(n 1 , λ 1 , λ 2 ), where B is a power series of n 1 and B(0, λ 1 , λ 2 ) = w 20 . Now using the transformation u 1 = n 1 , u 2 = n 2 √ B(n 1 ,λ 1 ,λ 2 ) and T = t 0 B(n 1 (s), λ 1 , λ 2 )ds, the system reduces tȯ where Q 1 is a power series in (u 1 , u 2 ) with power u i 1 u j 2 , i + j ≥ 3 and j ≥ 2. Now applying the final transformation y 1 = u 1 + w 10 2w 20 , y 2 = u 2 , the system transforms tȯ and Q 2 is a power series in y 1 and y 2 of the form y i 1 y j 2 , i + j ≥ 3. To check the non-degeneracy conditions for BT bifurcation atĒ we need to check the signs of k 20 and l 20 + k 11 for λ 1 = 0 and λ 2 = 0 i.e. at β = β BT and m = m BT . Now where a, b, c, d, a ij , b ij , i + j ≥ 1 are same as in the Supplementary Material A except one needs to calculate the coefficients atĒ. We have 1 − β BT − 2x < 0, as α and h(y) is positive. Using this with the condition (A4) we can find k 20 > 0 Now sign(k 11 ) may changes from negative to positive or vice versa. Therefore l 20 + k 11 may or may not be zero. If l 20 + k 11 = 0 then the system near the BT bifurcation point is topologically equivalent to the following systemẋ which is the normal form of the BT bifurcation of co-dimension two. The coefficient in front of x 1 x 2 can be either +1 or −1 according as sign(k 20 (l 20 + k 11 )). If the sign is positive then the predatorprey model undergoes a subcritical BT bifurcation atĒ and if the sign is negative then the model undergoes supercritical BT bifurcation atĒ.
In the case where l 20 + k 11 = 0, we have a co-dimension 3 Bogdanov-Takens bifurcation and the system is topologically equivalent to the following onė Note that this type of co-dimension 3 Bogdanov-Takens bifurcation includes a double-equilibrium point.

Supplementary material 3
Here we resent several α cross sections of the (α, δ, m) bifurcation diagrams shown in Fig. 5 of the main text. We consider the case where parametrization of h(y) is given either the Ivlev or the trigonometric response. Parameter values are (a) α = 7.1, β = 0.8; (b) α = 3, β = 0.8. The dark blue curve is saddle-node bifurcation curve, the red curve is a Hopf bifurcation curve, the green curve is the curve of a saddle-node bifurcation of limit cycles. The cyan curve describes a homoclinic bifurcation.
Firstly, we present the bifurcation diagrams for the model with the Ivlev response constructed for α = 7.1 and α = 3 which is shown in Fig. 1S. The bifurcation diagrams contain saddle-node bifurcation curve (blue), Hopf-bifurcation curve (red) and two global bifurcation curves such as the homoclinic bifurcation curve (cyan) and the saddle-node bifurcation of limit cycles (green). Two local bifurcation curves intersect at the Bogdanov-Takens(BT) bifurcation point and generalised Hopf bifurcation (GH) points are located on the Hopf bifurcation curve. The homoclinic bifurcation curve and the saddle-node bifurcation curve of limit cycles emerge from the BT and GH points, respectively. Note that for α = 3, we find two GH points which are connected by a global bifurcation curve. These four curves divide the parametric domain into six and five regions R j respectively. Regions are described in detail in the main text.
Secondly, we present the bifurcation diagrams constructed for the hyperbolic tangent parametrisation of h(y) which are shown in Figs. 2S. We consider the same values of the parameter α. The meaning of the curves and domains are the same as in the previous figure as well as in the main text. One can see that although the mutual position and the shape of bifurcation curves are similar to the Ivlev parameterisations, their shape is slightly different as those in Figs. 1S.

Supplementary Material 4
Here we will discuss the dynamics of two prey-predator models consisting of Allee effect in the predator: one with the Allee function in both functional and numerical responses (System (13)) and another with the Allee function in numerical response only (System (14)). First one describes the Allee effect due to collective exploitation of resources and another describes the same due to the lack of mating partner. We consider RM model with linear functional response and choose Monod parametrisation (h(y) = y y+δ ) for the Allee function. Without loss of generality, we consider dimensionless model and the evolving equations are in the following: To illustrate the dynamics of the above models, we construct bifurcation diagrams for both models in δ-m parametric plane for fixed α = 8. Fig. 3S present the bifurcation portraits of the model (13) consisting three local bifurcation curves: one saddle-node bifurcation curve (blue), two Hopf-bifurcation curves (red) and two non-local bifurcation curves, known as homoclinic bifurcation (cyan). Also the bifurcation diagram consists two BT points, corresponds to Bogdanov-Takens bifurcation of co-dimension 2, where Hopf-bifurcation curve meets tangentially to the saddle-node bifurcation curve and a global bifurcation curve emerges from it. These curves divide the parametric plane into 5 regions (R 1 − R 5 ) with one common region (R 2 ). We label this regions as the same way we did for other bifurcation portraits. Secondly we construct the bifurcation diagram of the model (14) in Fig. 4S when Allee function occurs in the numerical response only. This diagram consists same number of local and global bifurcation curves as previous bifurcation diagram (Fig. 3S), rather these diagrams looks topologically equivalent to each other.    (13) and (14) in δ-m parametric plane at fixed α = 8, (b) same bifurcation diagram with x-axis in logarithmic scale. Solid curves corresponds to system (13) and dashed curve corresponds to system (14) In Fig. 5S we have plotted all the local and global bifurcation curves for the two systems compositely to compare their dynamics. We use solid lines for the system (13) with Allee function in both functional and numerical responses and dashed line for the other. It is clearly seen that diagrams look like same except for the area of some regions (R i ) being increased or decreased. The way of appearance of the local and global bifurcation curves are almost same in both cases which leads the diagrams look like similar. But one noticeable fact is that predator in system (13) can persists in larger domain (in R 1 ) compared to the system (14) with Allee function in numerical response only. Thus the mechanism through which Allee effect appears into the predator growth, either due to exploitation of resources or lack of mating partner or something else, has a strong influence on predator persistence and can be considered as separate interesting study.