The Basins of Attraction in a Modified May-Holling-Tanner Predator-Prey Model with Allee Effect
TThe Basins of Attraction in a ModifiedMay–Holling–Tanner Predator-Prey Model with AlleeEffect
Claudio Arancibia–Ibarra a,b a School of Mathematical Sciences, Queensland University of Technology (QUT), Brisbane,Australia. b Facultad de Educaci´on, Universidad de Las Am´ericas (UDLA), Santiago, Chile.
Abstract
I analyse a modified May–Holling–Tanner predator-prey model considering anAllee effect in the prey and alternative food sources for predator. Addition-ally, the predation functional response or predation consumption rate is lin-ear. The extended model exhibits rich dynamics and we prove the existenceof separatrices in the phase plane separating basins of attraction related to os-cillation, co-existence and extinction of the predator-prey population. We alsoshow the existence of a homoclinic curve that degenerates to form a limit cy-cle and discuss numerous potential bifurcations such as saddle-node, Hopf, andBogadonov–Takens bifurcations. We use simulations to illustrate the behaviourof the model.
Keywords:
Predator-prey model, Allee effect, bifurcations, basin of attraction.
Preprint submitted to arXiv March 25, 2019 a r X i v : . [ m a t h . D S ] M a r . Introduction In the last decade interactions between species are appearing in differentfields of Population Dynamics. In particular, predation models have proposedand studied extensively due to their increasing importance in both Biologyand applied Mathematics [1, 2]. Predation models such as the Holling-Tannermodel [3] (or May–Holling–Tanner [4]) are of particular mathematical interestfor temporal and spatio–temporal domains [5, 6, 7]. The Holling–Tanner modelhas been used extensively to model many real–world predator–prey interactions[8, 9, 10, 11]. For instance, it has been used by Hanski et al. [8] to investigatethe predator–prey interaction between the least weasel (Mustela nivalis) and thefield vole (Microtus agrestis). This study was based under the hypothesis thatgeneralist predators predicts correctly the geographic gradient in rodent oscil-lations in Fennoscandia. Additionally, the authors showed that the amplitudeand cycle period decreasing from north to south.The traditional Holling–Tanner model is describing by the following pair ofequations dxdt = rx (cid:16) − xK (cid:17) − H ( x ) y,dydt = sy (cid:16) − ynx (cid:17) . (1)Here x ( t ) is used to represent the size of the prey population at time t , and y ( t ) isused to represent the size of the predator population at time t . Moreover, r is theintrinsic growth rate for the prey, s is the intrinsic growth rate for the predator, n is a measure of the quality of the prey as food for the predator, K is the preyenvironmental carrying capacity, nx can be interpreted as a prey dependentcarrying capacity for the predator. Moreover, the growth of the predator andprey population is a logistic form and all the parameters are assumed to bepositive. In addition, the behaviour of the Holling–Tanner model depends onthe type of the predation functional response chosen. This response function isused to represent the impact of predation upon the prey species. A Holling TypeI functional response provides a mechanism to explain the survival advantagefor animals to form large groups or herds assuming protection from externalthreats. In many cases, clustering reduces the total area (relative to the totalmass) exposed to chemicals, extreme weather, bacteria or predators [12, 13].In (1) the functional response corresponds to H ( x ) = qx . Where q is themaximum predation rate per capita [12]. The resulting Holling–Tanner modelis an autonomous two–dimensional system of differential equations and is givenby Kolmogorov type [14] given by dxdt = rx (cid:16) − xK (cid:17) − qxy,dydt = sy (cid:16) − ynx (cid:17) . (2)Additional complexity can be incorporated into these models in order tomake them more realistic. In the case of severe prey scarcity, some predator2pecies can switch to another available food, although its population growth maystill be limited by the fact that its preferred food is not available in abundance.For instance, least weasel is a generalist and nomadic species [8]. This abilitycan be modelled by adding a positive constant c to the environmental carryingcapacity for the predator [15]. Therefore, we have a modification to the logisticgrowth term in the predator equation, namely (1 − y/nx ) is replaced by (1 − y/ ( nx + c )) as shown below; dxdt = rx (cid:16) − xK (cid:17) − qxy ,dydt = sy (cid:18) − ynx + c (cid:19) . (3)Models such as (3) are known as modified Holling–Tanner models [15, 16, 17,18] as the predator acts as a generalist since it avoids extinction by utilising analternative source of food. Note that the Holling–Tanner predator-prey modelalso considers the case of a specialist predator, i.e. c = 0 [19]. It is assumedthat a reduction in a predator population has a reciprocal relationship withper capita availability of its favorite food [15]. Nevertheless, when c >
0, themodified Holling–Tanner does not have these abnormalities and it enhances thepredictions about the interactions. This model was proposed in [15], but themodel was only analysed partially. Using a Lyapunov function [19], the globalstability of a unique positive equilibrium point was shown.On the other hand, in this manuscript we consider a density-dependent phe-nomenon in which fitness, or population growth, increases as population densityincreases [20, 21, 22, 23, 24]. This phenomenon is called Allee effect or positivedensity dependence in population dynamics [25]. These mechanisms are con-nected with individual cooperation such as strategies to hunt, collaboration inunfavourable abiotic conditions, and reproduction [26]. When the populationdensity is low species might have more resources and benefits. However, thereare species that may suffer from a lack of conspecifics. This may impact theirreproduction or reduce the probability to survive when the population volume islow [27]. The Allee effect may appear due to a wide range of biological phenom-ena, such as reduced anti-predator vigilance, social thermo-regulation, geneticdrift, mating difficulty, reduced defense against the predator, and deficient feed-ing because of low population densities [28]. With an Allee effect included, theHolling–Tanner Type I model (3) becomes dxdt = rx (cid:16) − xK (cid:17) ( x − m ) − qxy ,dydt = sy (cid:18) − ynx + c (cid:19) . (4)The growth function ( x ) = rx (1 − x/K )( x − m ) has an enhanced growth rateas the population increases above the threshold population value m . If (0) = 0and (cid:48) (0) ≥ m ≤ ( x ) represents a proliferationexhibiting a weak Allee effect, whereas if (0) = 0 and (cid:48) (0) < m > ( x ) represents a proliferation exhibiting a strong Alleeeffect [29].The Holling–Tanner model with Allee effect is discussed further in Section 2and a topological equivalent model is derived. In Section 3, we study the mainproperties of the model. That is, we prove the stability of the equilibrium pointsand give the conditions for saddle-node bifurcations and Bogadonov–Takensbifurcations. In Section 4 we study the impact in the basins of attraction ofthe inclusion of the modification. We conclude the manuscript summarising theresults and discussing the ecological implications.
2. The Model
The Holling–Tanner model with Allee effect and alternative food is given by(4), and for biological reasons we only consider the model in the domain Ω = { ( x, y ) ∈ R , x ≥ , y ≥ } . The equilibrium points of system (4) are ( K, m, , c ), and ( x ∗ , y ∗ ), this last point(s) being defined by the intersection ofthe nullclines y = nx + c and y = r (1 − x/K )( x − m ) /q . In order to simplifythe analysis, we follow [16, 30, 31] and convert (4) to a topologically equivalentnondimensionalised model that has fewer parameters. Following [16, 30, 31]we introduce a change of variable and time rescaling, given by the function ϕ : ˘Ω × R → Ω × R , where ϕ ( u, v, τ ) = ( x, y, t ) is defined by x = Ku , y = Knv , dτ = Krdt and ˘Ω = { ( u, v ) ∈ R , u ≥ , v ≥ } . Additionally, we set C := c/ ( Kn ), S := s/ ( Kr ), Q := nq/r and M := m/K , so ( M, S, Q, C ) ∈ Π =( − , × R . By substitution of these new parameters into (4) we obtain dudτ = u ((1 − u ) ( u − M ) − Qv ) ,dvdτ = Svu + C ( u − v + C ) , (5)Note that system (4) is topologically equivalent to system (5) in Ω andthe function ϕ is a diffeomorphism preserving the orientation of time sincedet( ϕ ( u, v, τ )) = Kn/r > du/dτ = uR ( u, v ) and dv/dτ = vW ( u, v ) with R ( u, v ) = (1 − u )( u − M ) − Qv and W ( u, v ) = S ( u − v + C ) / ( u + C ), system(5) is of Kolmogorov type. That is, the axes u = 0 and v = 0 are invariant.The u -nullclines of system (5) are u = 0 and v = (1 − u )( u − M ) /Q , while the v -nullclines are v = 0 and v = u + C . Hence, the equilibrium points for system(5) are (0 , , M, , C ) and the point(s) ( u ∗ , v ∗ ) with v ∗ = u ∗ + C and where u ∗ is determined by the solution(s) of(1 − u )( u − M ) /Q = u + C , or equivalently ,d ( u ) = u − (1 + M − Q ) u + ( M + CQ ) = 0 . (6)Define the functions g ( u ) = (1 − u )( u − M ) /Q and h ( u ) = ( u + C ) and ob-serve that lim u →±∞ g ( u ) = −∞ and g (0) = − M . So, (6) can have at most two4 igure 1: The intersections of the function g ( u ) (red line) and the straight line h ( u ) (bluelines) by changing Q for the three possible cases of strong and weak Allee effect. positive real root, which are depending on the value of M and Q , see Figure 1.Additionally, the solution of the equation (6) are given by u , = 12 (cid:16) M − Q ± √ ∆ (cid:17) , with ∆ = (1 + M − Q ) − M + CQ ) , (7)such that u ≤ E ≤ u <
1, where E = (1 + M − Q ) / Modifying the parameter M and Q impacts ∆ and hence the number ofpositive equilibrium points. In particular,1. Strong Allee effect ( M > M − Q > <
0, then (5) has no equilibrium points in the first quadrant;(ii) ∆ >
0, then (5) has two equilibrium points P , = ( u , , u , + C )in the first quadrant; and(iii) ∆ = 0, then (5) has one equilibrium point ( E, E + C ) of ordertwo in the first quadrant,(b) If 1 + M − Q ≤
0, then (5) has no equilibrium points in the firstquadrant.2. Weak Allee effect ( M ≤ M − Q > M + CQ > <
0, then (5) has no equilibrium points in the first quadrant;(ii) ∆ >
0, then (5) has two equilibrium points P , = ( u , , u , + C )in the first quadrant; and(iii) ∆ = 0, then (5) has one equilibrium point ( E, E + C ) of ordertwo in the first quadrant,5b) If 1 + M − Q > M + CQ < M − Q < M + CQ < P in the first quadrant, since u < < u .(c) If 1+ M − Q > M + CQ = 0, then ∆ = (1+ M − Q ) . Therefore,(5) has one equilibrium point P = ( u , u + C ) in the first quadrant,since u = 0 and u became u = 1 + M − Q .(d) If 1+ M − Q = 0 and M + CQ <
0, then ∆ = − M + CQ ). Therefore,(5) has one equilibrium point P = ( u , u + C ) in the first quadrant,since u < u became u = (cid:112) − ( M + CQ ).(e) If 1 + M − Q ≤ M + CQ ≥
0, then (5) has no equilibrium pointsin the first quadrant.
Remark 2.1.
Observe that none of these equilibrium points explicitly dependon the system parameter S . Therefore, S is one of the natural candidates to actas bifurcation parameter.
3. Main Results
In this section, we discuss the stability of the equilibrium points of system(5) for strong and weak Allee effect.
Theorem 3.1.
All solutions of (5) which are initiated in the first quadrant arebounded and eventually end up in
Φ = { ( u, v ) , ≤ u ≤ , ≤ v ≤ } .Proof. First, observe that all the equilibrium points lie inside of Φ. Additionally,as the system is of Kolmogorov type, the u -axis and v -axis are invariant setsof (5). Moreover, the set Γ = { ( u, v ) , ≤ u ≤ , v ≥ } is an invariantregion since du/dτ ≤ u = 1 and v ≥
0. That is, trajectories enteringinto Γ remain in Γ. On the other hand, by using the Poincar´e compactification[32, 33] which is given by the transformation U = u/v and V = 1 /v then dU/dτ = (1 /v ) du/dτ − ( u/v ) dv/dτ and dV /dτ = − (1 /v ) dv/dτ . Then, byapplying the blowing-up method used in [16], the result follows. To determine the nature of the equilibrium points we compute the Jacobianmatrix J ( u, v ) of (5) J ( u, v ) = (1 − u )( u − M ) − Qv + u (1 + M − u ) − QuSv ( u + C ) S ( u + C − v )( u + C ) , The stability of the equilibrium points (0 , , M, , C ) is Lemma 3.1.
The equilibrium points (0 , and (1 , are a saddle points and ( M, is unstable if < M < . Moreover, the equilibrium point (0 , C ) is stableif − CQ ≤ M < and a saddle point if M < − CQ . Furthermore, if there are nopositive equilibrium points in the first quadrant, i.e. for ∆ < , then (0 , C ) is globally asymptotically stable in the first quadrant. roof. The determinant and the trace of the Jacobian matrix (3.1) evaluatedat (0 ,
0) are det( J (0 , − M S and tr( J (0 , S − M . Similarly, thedeterminant of the Jacobian matrix (3.1) evaluated at (1 ,
0) is det( J (1 , − (1 − M ) S <
0, since
M <
1. Then, the determinant and the trace of theJacobian matrix (3.1) evaluated at ( M,
0) are det( J ( M, SM (1 − M ) > J ( M, M (1 − M ) + S >
0, since
M <
1. Finally, the determinantand the trace of the Jacobian matrix (3.1) evaluated at (0 , C ) are det( J (0 , C )) = S ( M + QC ) and tr( J (0 , C )) = − ( M + CQ + S ). It follows that (1 ,
0) is alwaysa saddle point and ( M,
0) is unstable if the system (5) is affected by strongAllee effect (0 < M < ,
0) is a saddlepoint if 0 ≤ M <
M <
0. Furthermore,the equilibrium point (0 , C ) is stable if − CQ ≤ M <
M < − CQ . Finally, by Theorem 3.1 we have that solutions starting in thefirst quadrant are bounded and eventually end up in the invariant region Γ.Moreover, the equilibrium point (1 ,
0) is a saddle point and, if ∆ < ω -limit of all the trajectoriesis the equilibrium point (0 , C ).Next, we consider the stability of the two positive equilibrium points P , of system (5) in the interior of Φ. These equilibrium points lie on the curve u = v + C such that g ( u ) = u + C (6). The Jacobian matrix of system (5) atthese equilibrium points becomes J ( u, u + C ) = (cid:18) u (1 + M − u ) − QuS − S (cid:19) . Moreover, the determinant and the trace of J ( u, u + C ) are given bydet( J ( u, u + C )) = Su ( − − M + Q + 2 u ) , tr( J ( u, u + C )) = u (1 + M − u ) − S. Thus, the sign of the determinant depends on the sign of − − M + Q + 2 u ,while the sign of the trace depends on the sign of u (1 + M − u ) − S . This givesthe following results. Lemma 3.2.
Let the system parameters of (5) be such that M − Q > and M + CQ > . Then, ∆ > , as defined in (7) , and therefore the equilibriumpoint P is a saddle point.Proof. Evaluating − − M + Q + 2 u at u gives: − − M + Q + 2 u = − − M + Q + (1 + M − Q − √ ∆) = −√ ∆ < . Hence, det( J ( P )) < P is thus a saddle point.Note that if M + CQ = 0 then the equilibrium point P collapses with (0 , C )and if M + CQ < P crosses to the second or thirdquadrant. 7 emma 3.3. Let the system parameters be such that M − Q > and M + CQ > or M − Q > and M + CQ < or M − Q < and M + CQ < . Then, ∆ > and therefore the equilibrium point P is:(i) a repeller if < S < S := 12 (1 + M − Q + √ ∆)( Q + √ ∆) ; and(ii) an attractor if S > S ,Proof. Evaluating − − M + Q + 2 u at u gives: − − M + Q + 2 u = − − M + Q + (1 + M − Q + √ ∆) = √ ∆ > . Hence, det( J ( P )) >
0. Evaluating u (1 + M − u ) − S at u = u gives u (1 + M − u ) = 12 (1 + M − Q + √ ∆)(1 + M − (1 + M − Q + √ ∆)) = S , Therefore, the sign of the trace, and thus the behaviour of P depends on theparity of u (1 + M − u ) − S , see Figure 3. Remark 3.1.
Note that if M + CQ = 0 then ∆ = (1 + M − Q ) and therefore u (1 + M − u ) − S = − (1 + M − Q )(1 + M − Q ) − S . Hence the sign of thetrace depends on the parity of (1 + M − Q )(1 + M − Q ) − S If ∆ = 0 (7) the equilibrium points P and P collapse such that u = u = E = (1 + M − Q ) /
2. Therefore, system (5) has one equilibrium point of ordertwo in the first quadrant given by (
E, E + C ). Consequently, we have that ∆ = 0and therefore C = ((1 + M − Q ) − M ) / Q . Lemma 3.4.
Let the system parameters be such that ∆ = 0 (7) . Then, theequilibrium point ( E, E + C ) is:(i) a saddle-node attractor if S < S := 12 Q (1 + M − Q ) ; and(ii) a saddle-node repeller if S > S .Proof. Evaluating − − M + Q + 2 u at u = E gives − − M + Q + 2 u = − − M + Q + (1 + M − Q ) = 0 . Hence, det( J ( P )) = 0. Evaluating u (1 + M − u ) − S at u = E gives u (1 + M − u ) = 12 Q (1 + M − Q ) = S , Therefore, the behaviour of the equilibrium point (
E, E + C ) depends on theparity of Q (1 + M − Q ) / − S .See Figure 4 for phase portraits related to both cases of Lemma 3.4.8 emma 3.5. Let the system parameters be such that M − Q > and M + CQ = 0 . Then, the equilibrium point P is:(i) a repeller if < S < S := − (1 + M − Q )(1 + M − Q ) ; and(ii) an attractor if S > S ,Proof. Evaluating − − M + Q + 2 u at u = 1 + M − Q gives: − − M + Q + 2 u = 1 + M − Q > . Hence, det( J ( P )) >
0. Evaluating u (1 + M − u ) − S at u = 1 + M − Q gives u (1 + M − u ) = − (1 + M − Q )(1 + M − Q ) = S , Therefore, the sign of the trace, and thus the behaviour of P depends on theparity of u (1 + M − u ) − S , see Figure 5. Lemma 3.6.
Let the system parameters be such that M − Q = 0 and M + CQ < . Then, the equilibrium point P is:(i) a repeller if < S < S := √− M − CQ (cid:0) Q − √− M − CQ (cid:1) ; and(ii) an attractor if S > S ,Proof. Evaluating − − M + Q + 2 u at u = √− M − CQ gives: − − M + Q + 2 u = (cid:112) − M − CQ > . Hence, det( J ( P )) >
0. Evaluating u (1 + M − u ) − S at u = √− M − CQ gives u (1 + M − u ) = (cid:112) − M − CQ ( Q − (cid:112) − M − CQ ) = S , Therefore, the sign of the trace, and thus the behaviour of P depends on theparity of u (1 + M − u ) − S , see Figure 5. Lemma 3.7.
Let the system parameters be such that M − Q > , M + CQ > and ∆ < or M − Q ≤ and M + CQ ≥ . Then, there are nopositive equilibrium points in the first quadrant. Therefore, (0 , C ) is globallyasymptotically stable in the first quadrant.Proof. Finally, by Theorem 3.1 we have that solutions starting in the first quad-rant are bounded and eventually end up in the invariant region Γ. Moreover,the equilibrium point (1 ,
0) is a saddle point and, if ∆ < ω -limit of all the trajectories is the equilibriumpoint (0 , C ), see Figure 3 (a). 9 .2. Bifurcation Analysis In this section, we discuss some of the possible bifurcation scenarios of system(5). Observe that the stability of (0 , , M,
0) and P do not change thestability. Additionally, none of the equilibrium points P , P , P and ( E, E + C )explicitly depend on the system parameter S . Therefore, S is one of the naturalcandidates to act as bifurcation parameter. Theorem 3.2.
Let the system parameters be such that ∆ = 0 (7) . Then, system (5) experiences a saddle-node bifurcation at the equilibrium point ( E, E + C ) (forchanging Q ).Proof. The proof of this theorem is based on Sotomayor’s Theorem [34]. For∆ = 0, there is only one equilibrium point (
E, E + C ) in the first quadrant, with E = (1 + M − Q ) /
2. From the proof of Lemma 3.4 we know that det( J ( E, E + C )) = 0 if ∆ = 0. Additionally, let U = (1 , T the eigenvector correspondingto the eigenvalue λ = 0 of the Jacobian matrix J ( E, E + C ), and let W = (cid:18) S (3 + 3 M − Q )(1 + M − Q ) , (cid:19) T be the eigenvector corresponding to the eigenvalue λ = 0 of the transposedJacobian matrix J ( E, E + C ) T .If we represent (5) by its vector form F ( u, v ; Q ) = (cid:18) (1 − u )( u − M ) − Qvu − v + C (cid:19) , then differentiating F at ( E, E + C ) with respect to the bifurcation parameter Q gives F Q ( E, E + C ; Q ) = (cid:32) −
12 (1 + M − Q + 2 C )0 (cid:33) . Therefore, W · F Q ( E, E + C ; Q ) = − S (1 + M − Q + 2 C )(3 + 3 M − Q )(1 + M − Q ) (cid:54) = 0 . Next, we analyse the expression W · [ D F ( E, E + C ; Q )( U, U )]. Therefore, wefirst compute the Hessian matrix D F ( u, v ; Q )( V, V ), where V = ( v , v ), D F ( u, v ; Q )( V, V ) = ∂ F ( u, v ; Q ) ∂u v v + ∂ F ( u, v ; Q ) ∂u∂v v v + ∂ F ( u, v ; Q ) ∂v∂u v v + ∂ F ( u, v ; Q ) ∂v v v . At the equilibrium point (
E, E + C ) and V = U , this simplifies to D F ( E, E + C ; Q )( U, U ) = − M − Q ) − C S (1 + C ) . W · [ D F ( E, E + C ; Q )( U, U )] = − S (2 − M + Q )(3 + 3 M − Q )(1 + M − Q ) − SC (1 + C ) (cid:54) = 0 . By Sotomayor’s Theorem [34] it now follows that system (5) has a saddle-nodebifurcation at the equilibrium point (
E, E + C ), see Figure 2 and Figure 4. Theorem 3.3.
Let the system parameters be such that ∆ = 0 (7) and S = Q (1 + M − Q ) / . Then, system (5) experiences a Bogdanov–Takens bifurcationat the equilibrium point ( E, E + C ) (for changing ( Q, S ) ).Proof. If ∆ = 0, or equivalently Q = 1 + M − E , and E (1 + M − E ) = S , thenthe Jacobian matrix of system (5) evaluated at the equilibrium point ( E, E + C )simplifies to J ( E, E + C ) = (cid:18) E (1 + M − E ) − E (1 + M − E ) E (1 + M − E ) − E (1 + M − E ) (cid:19) , = 12 Q ( M − Q + 1) (cid:18) − − (cid:19) . So, det( J ( E, E + C )) = 0 and tr ( J ( E, E + C )) = 0. Next, we find the Jor-dan normal form of J ( E, E + C ). The latter has two zero eigenvalues witheigenvector ψ = (1 , T . This vector will be the first column of the matrixof transformations Υ. For the second column of Υ we choose the generalisedeigenvector ψ = (1 , T . Thus, Υ = (cid:18) (cid:19) andΥ − ( J ( E, E + C ))Υ = 12 Q (1 + M − Q ) (cid:18) (cid:19) . Hence, we have the Bogdanov–Takens bifurcation [34], or bifurcation of codi-mension two, and the equilibrium point (
E, E + C ) is a cusp point for ( Q, S ) =( Q , S ( Q )) such that ∆ = 0 and E (1 + M − E ) = S [35], see Figure 2 andFigure 4.The bifurcation curves obtained from Lemma 3.3, Lemma 3.4, and Theorem3.2 divide the ( M, S )-parameter-space into five parts, see Figure 2. Modifyingthe parameter M – while keeping the other two parameters ( Q, C ) fixed – im-pacts the number of positive equilibrium points of system (5). The modificationof the parameter S changes the stability of the positive equilibrium point P ofsystem (5), while the other equilibrium points (0 , M, ,
0) and P donot change their behaviour. There are no positive equilibrium points in system(5) when the parameters M, S are located in the red area where ∆ < , C ) is a global attractor, see Lemma 3.7 andFigure 3 (a). For M = M ∗ , which is the saddle-node curve SN in Figure 2,the equilibrium points P and P collapse since ∆ = 0, see Lemma 3.4 andFigures 4 (a) and (b). So, system (5) experiences a saddle-node bifurcation and11 igure 2: The bifurcation diagram of system (5) with strong ( M >
0) and weak ( M ≤ Q, C ) fixed and created with the numerical bifurcation package MATCONT[36]. The curve H represents the Hopf curve where P changes stability (Lemma 3.3), HOMrepresents the homoclinic curve of P , SN represents the saddle-node curve from Lemma 3.4where ∆ = 0, and BT represents the Bogdanov–Takens bifurcation from Theorem 3.3 where∆ = 0.
12 Bogdanov–Takens bifurcation (labeled BT in Figure 2) along this line, seeTheorems 3.2 and 3.3, and see also Figure 4 (c). When the parameter M islocated in − C/Q < M < M ∗ , system (5) has two equilibrium points P and P . The equilibrium point P is always a saddle point, see Lemma 3.2, while P can be unstable or stable. For ( M, S ) in the grey area the equilibrium point P is unstable, see Figure 3 (a). For ( M, S ) in the blue area the stable equilibriumpoint P is surrounded by a stable limit cycle, see also Figures 3 (b). For ( M, S )in the green area the equilibrium point P is stable, see Figure 3 (c). Finally,for ( M, S ) in the light green area system (5) has only one equilibrium point inthe first quadrant which is always stable, see Lemmas 3.5 and 3.6. Since P collapse with (0 , C ) or crosses to the second or third quadrant, see Figure 5)
4. Basins of Attraction
For system parameters (
Q, M, A, C ) such that the conditions 1.(a) (strongAllee effect) and 2.(a) (weak Allee effect) presented in Subsection 2.1 are met andfor
S > S , system (5) has two attractors, namely (0 , C ) and P . Furthermore,at the critical value S = S , such that tr ( J ( P )) = 0, P undergoes a Hopfbifurcation [32]. Note that S depends on Q and it can actually be negative. Inthat case P is an attractor for all S > > , C ) and P (for S > S ) in Φ (see Theorem 3.1). The stable manifold of the saddle point P , W s ( P ), often acts as a separatrix curve between these two basins of attraction.Let W u,s (cid:37) ( P ) be the (un)stable manifold of P that goes up to the right(from P ) and let W u,s (cid:46) ( P ) be the (un)stable manifold of P that goes downto the left (from P ). From the phase plane and the nullclines of system (5) itimmediately follows that W s (cid:37) ( P ) is connected with ( M,
0) and W u (cid:46) ( P ) with(0 , C ). Furthermore, everything in between of W s (cid:37) ( P ), W u (cid:46) ( P ) and the u -axisalso asymptotes to the origin.For ∆ > M > − CQ and depending on the value of S , there are differentcases for the boundary of the basins of attraction in the invariant region Φ, seeTheorem 3.1. By continuity of the vector field in S , see (5), we get:(i) For S on the grey region in the bifurcation diagram showed in Figure 2,the equilibrium point P is unstable, see lemma 3.3, and W u (cid:37) ( P ) connectswith (0 , C ). Hence, Φ is the basin of attraction of (0 , C ), see Figure 3 (a)and (d).(ii) For S on the blue region in the bifurcation diagram showed in Figure2. There is a stable limit cycle that surrounds P and W u (cid:37) ( P ) connectswith this limit cycle. This limit cycle is created around P via the Hopfbifurcation [37]. Therefore, W s ( P ) forms a separatrix curve between thebasins of attraction of P and (0 , C ) in this parameter regime, see Figure3 (b) and (e). 13iii) For S on the green region in the bifurcation diagram showed in Figure 2.Then, W s (cid:46) ( P ) intersects the boundary of Φ, and W s ( P ) again forms theseparatrix curve in Φ, see Figure 3 (c) and (f).Note that the system parameters ( Q, C ) are fixed at (0 . , .
07) and M =0 .
04 in the left panel of Figures 3 (a)-(c). Consequently, u , are constant.In particular, u ≈ . u ≈ . Q, C ) fixed at(0 . , .
07) and M = − .
01 in the right panel of Figures 3 (a)-(c). Consequently, u , are also constant. In particular, u ≈ . u ≈ . M > − CQ and depending on the value of S , there are threedifferent cases for the boundary of the basins of attraction in the invariant regionΦ, see Theorem 3.1. By continuity of the vector field in S , see (5), we get:(i) For 0 < S < S , the equilibrium point ( E, E + C ) is a saddle-node attrac-tor, see Lemma 3.4, see Figure 4 (a).(ii) For S < S , the equilibrium point ( E, E + C ) is a saddle-node repeller, seeLemma 3.4, see Figure 4 (b).(iii) For S = S , the equilibrium point ( E, E + C ) is a cusp point, see Theorem3.3, see Figure 4 (c).For system parameters ( Q, M, A, C ) such that the conditions 2.(b), 2.(c) and2.(d) presented in Subsection 2.1 are met (weak Allee effect), system (5) hasone attractor in the first quadrant, namely P , .For 1 + M − Q > M + CQ ≤ M − Q < M + CQ ≤ M − Q = 0 and M + CQ < S , there arethree different cases for the boundary of the basins of attraction in the invariantregion Φ, see Theorem 3.1. By continuity of the vector field in S , see (5), weget:(i) For S on the grey region in the bifurcation diagram showed in Figure 2,the equilibrium point P , is unstable (see lemma 3.5 and 3.6). Hence, Φis the basin of attraction of (0 , C ), see Figure 5 (a).(ii) For S on the blue region in the bifurcation diagram showed in Figure 2,the equilibrium point P , is unstable surrounded by a stable limit cycle(see lemma 3.5 and 3.6). Therefore, Φ is the basin of attraction of thestable limit cycle, see Figure 5 (b).(iii) For S on the green region in the bifurcation diagram showed in Figure 2,the equilibrium point P , is stable (see lemma 3.5 and 3.6). Hence, Φ isthe basin of attraction of the equilibrium point P , , see Figure 5 (c).Note that the system parameters ( Q, C ) are fixed at (0 . , .
1) and M = − .
055 in Figures 5 (a)-(c). Consequently, u is constant. In particular, u ≈ . Q, C ) fixed at (0 . , .
1) and M = − . u are also constant. In particular, u ≈ . igure 3: The grey region represent the basin of attraction of P , the yellow region representthe basin of attraction of a stable limit cycle and the orange region represent the basin ofattraction of (0 , C ). Moreover, the blue (red) curve represents the prey (predator) nullcline.For C = 0 .
07, and Q = 0 .
45, such that ∆ > M = 0 .
04 (left panel) and M = − .
01 (right panel) the equilibrium point (0 , C ) is global attractor. (b) M = 0 .
04 (leftpanel) and M = − .
01 (right panel) the equilibrium point P is surrounded by a stable limitcycle (grey curve) and W s ( P ) forms the separatrix curve in Φ. (c) M = 0 .
04 (left panel)and M = − .
01 (right panel) the equilibrium points (0 , C ) and P are attractors and W s ( P )forms the separatrix curve in Φ. Observe that the same color conventions are used in theupcoming figures. Moreover, Tikz and Matlab were used to do the simulations. igure 4: For C = 0 . Q = 0 .
5, and M = 0 . E, E + C ) = (0 . , . S < S = 0 . E, E ) is a saddle-node repeller. (b) For
S > S , the equilibrium point( E, E + C ) is a saddle-node attractor. (c) For S = S , the equilibrium point ( E, E + C ) is acusp point (See Figure 3 for the color conventions.) igure 5: For M = − . C = 0 .
1, and Q = 0 .
55, such that u = 0 and (a) S = 0 . < S k ,the equilibrium point (0 , C ) is global attractor. (b) For S k < S = 0 . < S ∗ k the equilibriumpoint P is surrounded by a stable limit cycle (black curve). (c) For S = 0 . > S ∗ k theequilibrium point P is stable. Similarly, for M = − . C = 0 .
1, and Q = 0 .
55, such that u < S = 0 .
19 the equilibrium point P is stable (See Figure 3 for the colorconventions). . Conclusions In this manuscript, the Holling–Tanner predator-prey model with strong andweak Allee effect and functional response Holling type I was studied. Using adiffeomorphism we analysed a topologically equivalent system (5). This systemhas four system parameters which determine the number and the stability ofthe equilibrium points. We showed that the equilibrium points (1 ,
0) and P arealways saddle points, ( M,
0) is an unstable point. Moreover, the equilibriumpoint (0 ,
0) can be a saddle or unstable equilibrium point and the equilibriumpoint can be stable or saddle point, see Lemmas 3.1 and 3.2. In contrast, theequilibrium point P can be an attractor or a repeller, depending on the trace ofits Jacobian matrix, see Lemma 3.3. Furthermore, for some sets of parametersvalues the equilibrium point P can collapses with (0 , C ) and then crosses thothe second or third quadrant. Therefore, there exist one positive equilibriumpoint ( P i with i = 3 ,