Predator-prey models are studied in both applied mathematics [1–3] and ecology [4–7]. The goal of these studies is to describe and analyse the predation interaction between the predator and the prey and predict how they respond to future interventions [8, 9]. These studies often use mathematical models to describe the species' interactions and the time-series behaviors [10, 11]. The models aim to be representative of real natural phenomena capturing the essentials of the dynamics. However, new technologies used to study biological and physical phenomenon reveal that species' interactions are more complex than previously used in the models [7, 12–14]. The importance of these more complex interactions are becoming increasingly apparent as research findings have shown that ecosystem dynamics depend on the particular nature of the interaction processes, such as the functional response and predation rate [6, 7, 15, 16].

A standard approach for using models to understand ecological systems is to design a framework based on simple principles and compare species' abundance that result from the predicting analysis of those models. However, this approach becomes more difficult when additional nuances to standard models are added, making them more complex and difficult to analyse. For instance, Graham and Lambin [17] showed that field-vole (Microtus agrestis) survival can be affected by reducing least weasel (Mustela nivalis) predation. They also demonstrated that weasels were suppressed in summer and autumn, while the vole (Microtus agrestis) population always declined to low density. However, the authors argued that the underlying model was too difficult to study due to a large number of parameters.

The model used by Graham and Lambin in [17] is a Leslie–Gower predator-prey model [18] and is given by (1)dNdt=rN(1-NK)-qNPN+a,dPdt=sP(1-PhN). Here, N(t) > 0 and P(t) > 0 are used to represent the size of the prey and predator population at time t, respectively. The prey population grows logistically with carrying capacity K and intrinsic growth rate r. The growth of the predator is also logistic, with intrinsic growth rate s, but the carrying capacity is prey dependent and h is a measure of the quality of the prey as food for the predator. The functional response is Holling type II where q is the maximum predation rate per capita and a is half of the saturated response level [7].

In population dynamics, many ecological mechanisms are connected with individual cooperation such as strategies to hunt, collaboration in unfavorable abiotic conditions and reproduction [19]. An example of such an ecological mechanism is the cooperative interaction observed in Kakapo (Strigops habroptilus). This interaction can be influenced by reducing the fecundity levels and the population size [20]. Another factor that affects the population dynamics is the species density. For instance, when the predator population size is low they have more resources and benefits. However, there are species that may suffer from a lack of conspecifics, which may be less likely to reproduce or survive in a small-sized population [21]. In these instances, the size of the population is important, as for a smaller size of biomass adaptability may be diminished [20]. When Allee and collaborators [22] analyzed the data of the false weevil (Tribolium confusum) they observed that the highest growth rates of their populations per capita were at intermediate densities [20]. The fact that they were lower in high densities was not surprising, as intraspecific competition is high. In contrast, when fewer males were present, females produced fewer eggs, which is not an obvious correlation for an insect. In this case, optimal egg production was thus achieved at intermediate densities. This effect is now referred to as the Allee effect [22, 23]. It is modeled by incorporating the factor (N − m) into the logistic equation, with the idea that when the prey population N decreases below the Allee threshold m > 0 the prey rate of growth becomes negative, indicating signs of extinction. Therefore, the prey logistic growth r(1 − N/K) is replaced by r(1 − N/K)(N − m). For 0 < m < K, the per-capita growth rate of the prey population with the Allee effect included is negative, but increasing, for N ∈ [0, m), and this is referred to as the strong Allee effect. When m ≤ 0, the per-capita growth rate is positive but increases at low prey population densities and this is referred to as the weak Allee effect [20, 23]. With the Allee effect included, the Leslie–Gower predator-prey model (1) of [17] becomes (2)dNdt=rN(1-NK)(N-m)-qNPN+a=N·W(N,P),dPdt=sP(1-PhN)=P·R(N,P). The aim of this manuscript is to study the Leslie–Gower predator-prey model with weak Allee effect on prey, that is (2) with m < 0. By considering a weak Allee effect in the prey population we complement the results of Ostfeld and Canhan [24]. The authors studied the stabilization of a file-vole population which depends on the variation in reproductive rate and recruitment of the population, i.e., an Allee effect. The weak Allee effect can be also observed in the File-vole species as the survival rate for adults were delayed [24]. However, these phenomena were not considered by the authors as they used model (1). This manuscript also extends some of the results obtained by Arancibia–Ibarra and González–Olivares [25] and González–Olivares et al. [26] for a modified Leslie–Gower model with m = 0, that is, with a specific type of weak Allee effect, see second row of Table 1. Leslie–Gower models with strong Allee effect and different type of functional responses have been extensively studied in [27–29], see first row of Table 1. In these articles the authors showed that the system, for certain system parameters, supports the extinction of both species and it also supports the stabilization of both populations over the time, i.e., coexistence. Moreover, system (2) with m < 0 complements the results of the Leslie–Gower model studied by Courchamp et al. [20]. The authors in [20] argue that the Allee effect in prey generally destabilizes the dynamics between the prey and the predator. Moreover, the Allee effect can prevent the oscillation of both populations. However, in this manuscript, we find that the Leslie–Gower model with weak Allee effect supports the coexistence and oscillation of both populations.

In section 2, we non-dimensionalise the Leslie–Gower model with weak Allee effect and discuss the number of equilibria the model has in the first quadrant. The main mathematical difference between system (1), (2) with m < 0, and (2) with m ≥ 0, is the fact that (2) with m < 0 has at most three positive equilibria¹ in the first quadrant instead of one for system (1) and two for system (2) with m ≥ 0 [25, 26, 29]. These additional equilibrium gives rise to different type of bifurcations, such as saddle-node bifurcations, Bogdanov–Takens bifurcations, and homoclinic bifurcations. The main properties of the equilibria are studied in section 3. In particular, we study the stability of the equilibria and present the conditions for which the model undergoes different type of bifurcations. Finally, in section 4, we summarize the results and discuss the ecological implications of the model.

The Leslie–Gower model with weak Allee effect is given by (2) with m < 0, and we only consider the model in the domain Ω = {(N, P) ∈ ℝ², N > 0, P ≥ 0} and (r,K,q,a,s,h)∈ℝ+6. The axes in system (2) are invariant since (2) is of Kolmogorov type [30] (i.e., dN/dt = N · W(N, P) and dP/dt = P · R(N, P)). The equilibria of system (2) are (K, 0) and (x*, y*), which are the intersections of the nullclines P=hN and P=rq(1-NK)(N+a)(N-m). We follow the non-dimensionalisation approach of [29] to simplify the analysis and remove the N = 0 singularity for the model. We introduce Ω̄={(u,v)∈ℝ2,u≥0,v≥0} and the dimensionless variables (u, v, τ) given by (3)φ:Ω̄×ℝ→Ω×ℝ where φ(u,v,τ)=(NK,PhK,rKtu(u+aK)). Observe that φ (3) is a diffeomorphism which preserve the orientation of time [31, 32]. Next, set A: = a/K ∈ (0, 1), S: = s/(rK), Q: = hq/(rK) and M: = m/K, then (2) transforms into the non-dimensionalised system (4)dudτ=u2((u+A)(1-u)(u-M)-Qv),dvdτ=S(u+A)(u-v)v. The u-nullclines of system (4) are given by u = 0 and v = (u + A)(1 − u)(u − M)/Q, while the v-nullclines of interest are given by v = 0 and v = u. Hence, the equilibria of system (4) are (0, 0)², (1, 0), and the positive equilibria (u*, v*) with v* = u* and where u* is determined by the solution(s) of (u + A)(1 − u)(u − M)/Q = u, or, equivalently, (5)g(u):=u3-T(A,M)u2-L(A,M,Q)u+AM=0, with T(A, M) = 1 − A + M and L(A, M, Q) = A(M + 1)−Q−M. Since g(0) = AM < 0 and g(1) = Q > 0, there is always at least one positive equilibrium (u*, u*) with 0 < u* < 1. To obtain information about potential other positive equilibria, we divide (5) by (u − u*) and obtain the quadratic equation (6)u2+u(u*-T(A,M))+u*(u*-T(A,M))-L(A,M,Q)=0. Note that as system (4) is topologically equivalent to system (2), the solution u* of (5) is related to a population equilibrium of system (2). Therefore, a positive equilibrium point P = (u*, u*) in system (4) correspond to a coexisting equilibrium point in system (2).

Lemma 2.1. The positive roots of g(u) (5) lie in (0, 1), and

Proof: The cubic equation g(u) (5) always has at least one root in (0, 1), since g(0) = AM < 0 and g(1) = Q > 0, and the first case of the lemma, case 1, thus immediately follows from Descartes' rule of signs ([33], e.g.).

Case (II)(i) also follows directly from Descartes' rule of signs, which for T(A, M) > 0 and L(A, M, Q) < 0 states that g(u) (5) has one or three positive roots (counting multiplicity), and the observation that the quadratic equation 5 has two complex roots as its discriminant is negative.

In the last case, case (II)(ii), g(u) (5) has three real roots: u* and (8)u-*=12(T(A,M)-u*-Δ),u+*=12(T(A,M)-u*+Δ). The latter two being the roots of (6). What remains to show is that u±*∈(0,1). To show this, we look at the derivative of g(u) g′(u)=3u2-2T(A,M)u-L(A,M,Q). The assumptions T(A, M) > 0 and L(A, M, Q) < 0 imply that g′(u) ≥ 0 for u ≤ 0. Furthermore, T(A, M) ∈ (0, 1) and, for u ≥ T, we thus also get g′(u)=3u2-2T(A,M)u-L(A,M,Q)>2u(u-T(A,M))                    +u2>2u(u-T(A,M))≥0. In other words, under the assumptions T(A, M) > 0 and L(A, M, Q) < 0, the cubic function g(u) is strictly increasing outside the interval (0, T(A, M)) (with T(A, M) < 1). Since g(0) = AM < 0 and g(1) = Q > 0, it thus follows that u±*∈(0,1).     □

We observe that none of these equilibria explicitly depend on the system parameter S. Moreover, only L(A, M, Q) in g(u) (5) depends directly on Q. Therefore, S and Q are natural candidates to act as bifurcation parameters. For instance, when two of the three roots of g(u) coincide, say u-* and u+*, we change from having three roots to one root (upon changing Q). So, at this point we have u-*=u+*⇒u*=u*,±               :=13(T(A,M)±2T(A,M)2+3L(A,M,Q)). Implementing this into g(u^*,±) = 0 results in an analytic implicit expression Q = Q^±(A, M) at which the number of roots changes. For example, for (A, M) = (1/10, −1/10) as in Figure 1 we get Q±(1/10,-1/10)=(73±5)/200.³

In case (II)(ii) of Lemma 2.1 we find that g(u) has three roots in (0, 1): u* and u±* (8), with u-*≤u+*. The latter two roots depend on the first root u*, but, a priori, we do not know anything about parity of u*-u±*. However, the roots of g(u) do not depend on which of the three roots we pick a priori as u* in the lemma. Hence, we can, without loss of generality, assume in the remainder of this manuscript that 0<u*≤u-*≤u+*<1. However, see the proof of Corollary 3.3 where we utilize this independence of our initial pick.

In this section, we discuss the main results related to system (4). That is, we discuss the nature of the equilibria (section 3.1) and their bifurcations (section 3.2). First we observe that from [29, Theorem 2] it instantly follows that all solutions of (4) which are initiated in the first quadrant are bounded and end up in (9)Φ={(u,v), 0<u≤1, 0≤v≤1}.

In this manuscript, we studied a Leslie–Gower predator-prey model with weak Allee effect and functional response Holling type II, i.e., system (2) with m < 0. We simplified the analysis by studying a topologically equivalent system (4). The topologically equivalent system (4) has two equilibria on the axis which are always (non-hyperbolic) saddle points, see Lemma 3.1. In addition, system (4) has at most three positive equilibria in the first quadrant, see Lemma 2.1 and Figure 1. As the function φ is a diffeomorphism preserving the orientation of time, the dynamics of system (2) is topologically equivalent to system (4) ([29], Theorem 1). Therefore, we can, for instance, conclude from the results of section 3 that there are conditions on the system parameter for which the predator and prey can coexist without oscillations or with oscillations, see Figure 6 and this behavior depends intrinsically and non-linearly on the system parameters including the predation rate (q) and the intrinsic growth rate of the predator (s).

In more detail, when the system parameters are such that system (4) has only one positive equilibrium, i.e., Q < Q⁻(A, M) or Q > Q⁺(A, M), see Regions I, II, VII and VIII in Figure 7, then it is a repeller surrounded by a stable limit cycle or an attractor, but not a saddle point, see Corollary 3.2 and Figure 3. Note that when the equilibrium is an attractor it is not necessary a global attractor since there is a small region in parameter space near the Hopf bifurcation where the equilibrium is surrounded by two limit cycles, see Figure 8A. The observed behavior for system (4), and hence system (2) with m < 0, in this case of one positive equilibrium – global attractor, attractor with two limit cycles, or repeller with one limit cycle – is similar to the behavior of the original Leslie–Gower predator-prey (1) with logistic growth for the prey [40]. In other words, for q smaller than some q⁻ or larger than some q⁺ both systems (1) and (2) with m < 0 present qualitatively similar behaviors, see Figure 3B.

In contrast, for Q⁻(A, M) < Q < Q⁺(A, M), system (4) has three equilibria in the first quadrant, which is not possible in the original Leslie–Gower predator-prey (1) nor in the model with strong Allee effect, i.e., system (2) with m > 0. In this case, the middle one is always a saddle point, while the outer two are repellers or attracters, see Corollary 3.3 and Figures 4, 6. If the two outer equilibria are attractors, then either the stable manifold of the saddle point determines a separatrix curve which divides the basins of attraction of the two attracters, or there exists an unstable limit cycle dividing the basins of attraction, see Figures 6A–D. When both outer equilibria are repellers the system necessarily possesses a stable limit cycle, see Figure 6F. When the outer equilibria are an attractor and a repeller, the attractor can be a global attractor or there are two limit cycles surrounding the equilibria, see Figures 6E, 8B. The latter one again happens near the Hopf bifurcation. In other words, there are regions in parameter space for which system (4) has two attractors or an attractor and a stable limit cycle. As such, a modification of one, or both, of the species could have the result that you end up in a different basin of attraction and you will thus have significantly different dynamics as you will approach the other attractor, see Figures 6A–D, 8B. We also performed an initial (numerical) bifurcation analysis which revealed the existence of saddle-node and Bogdanov–Takens bifurcations, see Theorems 3.1 and 3.2 and Figure 7.

In summary, we showed that the weak Allee effect in the Leslie–Gower model (2) presents—due to the presence of a region in parameter space where we have three positive equilibria—richer dynamics than the original Leslie–Gower predator-prey model (2) studied, for example, by Saez and Gonzalez-Olivares [40]. The original Leslie–Gower predator-prey model (2) does not possess the bifurcations studied in section 3.2. That is, a saddle-node bifurcation and a Bogdanov–Takens bifurcation. Moreover, the original model only has one positive equilibrium point, which can be stable surrounded by two limit cycles. Therefore, the populations could switch between a stable state or an oscillatory behavior. However, system (2) with a weak Allee effect can have two stable equilibria and thus depending on the initial condition the population could end up in two stable states. From an ecological point of view, we can also conclude that the species in system (2) could coexist or oscillate but never go extinct, see Figure 6. This behavior was also observed by Ostfeld and Canhan [24]. The authors found that the stabilization of vole (Microtus agrestis) populations in Southeastern New York depend on the variation in reproductive rate and thus this impact the reproductive activity of the species. This is again similar to the original Leslie–Gower predator-prey model (2), but different from the Leslie–Gower model (2) with strong Allee effect, i.e., m > 0. In the latter case, the model supports coexistence and extinction of the species [29].

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

CA-I implemented the methodology, performed the analytics and numerical analysis, composed all the figures, drafted the manuscript, wrote the cover letter and acted as the corresponding author. JF provided the technical assistance with the interpretation of the results and critically reviewed the manuscript. PH directed the research, provided the technical assistance with the interpretation of the results and critically reviewed the manuscript. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.