Study of dynamical properties in a bi-dimensional model describing the prey–predator dynamics with strong Allee effect in prey

. In this work, we analyze a bi-dimensional differential equation system obtained by considering Holling type II functional response in prey–predator model with strong Allee effect in prey. One of the important consequence of this modiﬁcation is the existence of separatrix curve which divides the behaviour of the trajectories in the phase plane. The results show that the origin is an attractor for any set of parameter values. Axial equilibrium points are stable or unstable according to the different parametric restrictions. The unique positive equilibrium point, if it exists, can be either an attractor or a re-peller surrounded by a limit cycle, whose stability and uniqueness are also established. Therefore long-term coexistence of both populations is possible or they can go to extinction. Conditions on the parameter values are derived to show that the positive equilibrium point can be emerged or annihilated through transcritical bifurcation at axial equilibrium points. The existence of two heteroclinic curves is also established. It is also demonstrated that the origin is a global attractor in the phase plane for some parameter values, which implies that there are satisfying conditions where both populations can go to extinction. Ecological interpretations of all analytical results are provided thoroughly. 58F14, 58F21.


Introduction
In population dynamics, classical mathematical models generally assume that an increase of the population density has a negative effect on the reproduction and survival of each individuals. An usual way to describe such an assumption is to consider models that incorporate per capita birth and mortality rates depending on the population size. Several equations are used to describe this density dependence among which we can cite the classical logistic, Bernoulli and Gompertz equations (see [34] for a review and mathematical properties). However, in some specific situations, the density dependence can operate positively in small populations: while an increase in density is negative for high densities, it may be beneficial for low densities.
The Allee effect, first introduced by W. C. Allee in 1931 [3], is one of the important and ecologically relevant factor in population dynamics, which is exhibited by the growth of natural population and in case of Allee effect the per capita growth rate increases at small species density [9]. Allee effects have been studied extensively in recent years, mainly because of their potential role in extinctions of already endangered, rare or dramatically declining species (see [33] and references therein). Several possible mechanisms of Allee effect are well-known, such as mate limitation, satiation of generalists predators or group defense against a predator (see [11,20]). There are two types of Allee effect: the strong Allee effect and weak Allee effect. Strong Allee effect is subjected to the negative population growth rate when population density falls below a critical value but for weak Allee effect, the growth rate decreases yet remains positive at low population density [31].
The prey-predator models having Allee effect in the growth rate of prey population can exhibit variety of dynamic behaviours compared to the analogous models with logistic growth in prey population [4,8,27,31]. In [8], Conway and Smoller investigated the predator-prey model with Holling type I functional response and having Allee effect in the prey growth function. They observed a variety of dynamical behavior such as Hopf bifurcation, existence of stable limit cycle around non trivial unstable steady state and existence of homoclinic loop. The prey-predator model subjected to the strong Allee effect in prey population and with Holling type II functional response was investigated by Berezovskaya et al. and Morozov et al. in [4,27] and they reported various rich dynamics including existence of heteroclinic loop, two limit cycles under certain parametric restrictions. Recently, González-Olivares et al. reported that the system exhibits very rich dynamics when different choices of Allee effect is considered into the prey growth function [2,13].
On the other hand habitat complexity is the structural complexity of habitats. Habitat complexity can strongly mediate predator-prey interactions, affecting not only total predation rates, but also modifying selectivities for different prey species or size classes [32]. Pennings [29] and Grabowski [15] found that habitat complexity reduces encounter rates of predators with prey. For example, aquatic habitat becomes structurally complex in presence of submerged vegetation or aquatic weeds. It is observed that structural complexity of the habitat stabilizes the predator-prey interaction between piscivorous perch (predator) and juvenile perch and roach (prey) by reducing predator foraging efficiency. Luckinbill prolonged the coexistence of paramecium aurelia (prey) and Didinium nasutum (predator) in laboratory system by increasing of habitat complexity using methyl cellulose in the Cerophyl medium (nutrient) [26]. Therefore it is important to incorporate the effect of habitat complexity when predator-prey interaction is studied by means of models.
The goal of this article is to study the impact of predator population control on the dynamics of a prey population that is subject to Allee effect. Here, we have studied the complete global dynamics of the system under different parametric restrictions. We have performed an extensive bifurcation analysis to determine the global dynamics of the model in a two dimensional plane. It has been observed that the model posses very rich dynamics. We have studied the existence of transcritical bifurcation of equilibrium point to determine the change in the number of equilibrium points. Remarkably, a heteroclinic curve exists for some parameter values, which can be used as the boundary of deterministic extinction of both species. In addition, we establish the existence of a unique limit cycle surrounding the positive equilibrium point and also derived its stability condition.
The paper is organized as follows: in Section 2, we present the model and its nondimensionalized version. Section 3 deals with some basic properties of the model like positivity and boundedness of solution. Stability and bifurcation properties of the model are investigated in Sections 4 and 6 respectively. Existence of heteroclinic curves is shown in Section 5. The global dynamics of the model together with detailed numerical simulations is discussed in Section 7 and the conclusion is given in Section 8.

The model
The most general continuous time ODE model describing the dynamics of prey-predator population with prey dependent functional response is given by [16], subject to the initial condition p(0) > 0, z(0) > 0. Here, p(τ) and z(τ) are the densities of prey and predator population at any instant of time τ respectively, g(p) is the per capita net growth rate of prey in absence of predator, f (p) is the functional response of predator, i.e. the rate of prey consumption by an average predator, θ is the conversion rate and d is a destruction rate of the predator due to a control of their population. The functional response can be classified based on its dependency on prey, predator population as: (a) prey-dependent, when the functional response depends only on prey density; (b) predator-dependent, when the prey and predator density both determine the functional response; (c) multispecies-dependent, when the other species except prey and predator influence the functional response [1]. The most commonly used functional response to describe the prey-predator interaction is the Holling type II functional response [17], defined as where f (p) is the amount of food consumed by predator, p is the amount of food offered, a is the proportionality constant describing the attack rate, h is the handling time per food item. Since the existence of habitat complexity reduces the probability of capturing the prey by reducing the searching efficiency of predator, it affects the attack coefficient [35]. Hence, it is reasonable to consider the attack coefficient as is a dimensionless parameter that measures the degree or strength of habitat complexity. The simplest way of modeling habitat complexity is to replace a by a(1 − c). Here, we assume that the complexity is homogeneous throughout the habitat. Therefore, the total number of prey caught is given by [19] where T s = T − hV. T is the total time, T s is the available searching time. Solving for V the modified Holling type II functional response comes as: Since the predator's functional response is defined as the number of prey caught by a predator per unit time, the functional response in presence of habitat complexity is given by In this paper, we assume the prey growth rate is subject to strong Allee effect. Now we consider the following explicit form of g(p) parametrizing the strong Allee effect [24]: where r is the coefficient determining the magnitude of per capita growth rate; K is the environmental carrying capacity of prey; the parameter x 0 is used to refer the prey survival threshold for strong Allee effect (in this case 0 < x 0 ≤ K). Thus at small population densities (0 < p < x 0 ) the prey growth rate becomes negative since the mortality rate becomes larger than the birth rate.
with the initial condition p(0) > 0, z(0) > 0. In this paper, we study what may be the influence of parameter d on the whole dynamics of the model (2.5)-(2.6).
In order to reduce the number of parameters in the above model we find the nondimensionalized version of the model.

Nondimensionalized version:
Then, we get the following transformed system: Therefore modelling a strong Allee effect implies 0 < β ≤ 1. Equations (2.7)-(2.8) contain only four parameters whereas the original equations contain seven parameters. Moreover, the control dz in model (2.5)-(2.6) is totally characterized by the term δ 1 y in (2.7)-(2.8), so that we aim at studying the impact of parameter δ 1 (death rate of predator) in the dynamics of model (2.7)-(2.8).

Well posedeness of the problem
The fundamental theory of ordinary differential equations assures the existence and uniqueness of solution (x(t), y(t)) of the model system (2.7)-(2.8) with the given initial condition (x(0), y(0)).
Assume that δ 1 > 1. From (2.8) and x(t) > 0 for all t > 0, we have which implies that y(t) ≤ y(0)e −(δ 1 −1)t and lim t→∞ y(t) = 0. Therefore, the asymptotic behaviour of the prey population is the same as the asymptotic behaviour of the solution of the equation Since 0, β, 1 are the equilibrium states of (3.1), it is easy to see the following.
In the rest part of this paper, we assume that δ 1 < 1.

Existence and stability of equilibria
In this section, we present exhaustive analysis for the number of equilibrium points of the system (2.7)-(2.8) and their asymptotic stability.
Each of the boundary equilibrium point exists without any restriction of the model parameters.
The interior equilibrium point is the point(s) of intersection(s) of two zero-growth isoclines x = δ 1 δ 2 1−δ 1 and y = g 1 (x). The isocline x = δ 1 δ 2 1−δ 1 is feasible as δ 1 < 1, lying in the positive quadrant. The function y = g 1 (x) is positive in the interval (β, 1). Therefore, the system (2.7)-(2.8) admits a unique positive equilibrium state E * (x * , y * ) such that which is a feasible equilibrium point for β < x * < 1 or equivalently for Therefore, the considered system has at least three and at most four equilibrium points.

Local asymptotic stability of equilibrium points
After finding all the equilibrium solutions of the system (2.7)-(2.8), we investigate the local stability nature of each equilibrium solutions. This enable us to study the local bifurcations of the system , which help us to understand the nonlinear dynamics associated with the system. To study the bifurcation analysis, we concentrate on two parameters of the system given by δ 1 and δ 2 .
To discuss about the stability properties we derive the Jacobian matrix of the system (2.7)-(2.8) by using standard linearization technique and analyse its eigenvalues to identify the nature of associated equilibrium point. The general form of the Jacobian matrix of (2.7)-(2.8) at any point (x, y) is given by Hence, the Jacobian matrix at E 0 is given by which has the eigenvalues −γβ, −δ 1 . Thus, the origin is a stable node. The Jacobian matrix at E 1 is and saddle point if δ 1 < 1 1+δ 2 . Moreover, {(x, y)|x > 0, y = 0} is the stable manifold of the equilibrium point E 1 . The Jacobian matrix calculated at E 2 is given by which has the eigenvalues γβ(1 − β) > 0 and β β+δ 2 − δ 1 . Thus, E 2 is a saddle point if β β+δ 2 < δ 1 and unstable node if β β+δ 2 > δ 1 . Therefore from the feasibility condition of interior equilibrium point it is clear that E 2 is always a saddle point in presence of interior equilibrium point. Now, we summarize all the above local stability results in the following proposition. (i) E 0 is always a stable node; 2 and an unstable node if δ 1 < β β+δ 2 . For the stability of E * , the Jacobian matrix evaluated at this equilibrium is given by The determinant of matrix J E * is given by Then, E * is stable if Tr(J E * ) < 0. By direct computations we obtain The discriminant of the quadratic F(x) is Hence F(x) admits two roots given by On the other hand, we have Then one root of F(x) belongs to (−∞, β) and the other in (β, 1). Since x * * < x * * , then the unique root of F(x) that belongs to (β, 1) is x * * . Therefore, the stability results of interior equilibrium point can be summarized in the following proposition: then Tr(J E * ) > 0 and therefore the positive equilibrium point E * is unstable; , then Tr(J E * ) < 0 and therefore E * is a stable equilibrium point.

Local bifurcations
Bifurcation theory plays an important rule to analyse the ODE model [23,30], where we study the asymptotic behaviour of the system (equilibria, periodic solution or limit cycle etc.) based on parameter variation for qualitative changes. The parameter values where the qualitative changes of asymptotic behaviour occur are referred as bifurcation points.
In this section, we present several types of local bifurcations of codimension one like Hopf bifurcation and transcritical bifurcation,which occur due to change in stability properties of various equilibrium points.

Hopf bifurcation
The Hopf bifurcation refers to the development of periodic orbit ("self-oscillation") through the change in stability of stable equilibrium point when some model parameter crosses a critical value. The appearance of such kind of periodic orbit is interpreted as a "shift of stability" from the original stationary solution (stable equilibrium point) to the periodic one. From the expression of Tr(J E * ) in the previous section, it is clear that if x * = x * * , then Tr(J E * ) = 0. Taking δ 1 as bifurcation parameter. So, we have Moreover, to ensure the existence of a Hopf bifurcation we have to check the transversality condition. Differentiating the expression for Tr(J E * ) with respect to δ 1 , we get ∂ ∂δ 1 Tr(J E * )

Stability of limit cycle
The above result establishes the existence of small amplitude periodic solution, i.e. existence of a limit cycle in phase plane near the interior equilibrium point E * . In order to discuss the stability of the limit cycle we now compute the first Lyapunov coefficient l 1 at the critical parametric value δ 1 = δ * 1 . We first translate the equilibrium E * (x * , y * ) to the origin by using the transformation x = x * + u and y = y * + v. Substituting this transformation in (2.7)-(2.8) and expanding in Taylor series at the critical parametric condition δ 1 = δ * 1 we get and a ij , b ij are given by Hence the first Lyapunov coefficient l 1 for a planar system (as defined in [30]) is given by . Hence, the first Lyapunov coefficient l 1 for system (2.7)-(2.8) is given by We also have Noting that b 20 b 11 = cb 21 , then The stability of the Hopf bifurcating periodic solution depends upon the sign of σ. The limit cycle is stable if σ < 0 and is unstable for σ > 0. Noting that σ < 0 for 1 − 5δ * 1 ≤ 0. Therefore, the Hopf bifurcation is supercritical if σ < 0 and subcritical if σ > 0.

Uniqueness of limit cycle
The existence, uniqueness and multiplicity of limit cycle of planar system is one of the most important mathematical question related to Hilbert's 16th problem, and it is also practically useful in explaining the mechanical and natural oscillatory phenomena. Various techniques has been developed to establish the uniqueness of limit cycles of ecological systems [6,21,25]. One of the main ideas is to transform an ecological system into a generalized Liénard system, which we also follow here. We introduce two new variables z and τ 1 in place of y and t within the system of equations (2.7)-(2.8), defined by Under the said transformation, the system of equations (2.7)-(2.8) takes the following form in terms of the variables (x, z, τ 1 ) x . All the parameters involved in the above system of equations are positive and satisfy the following parametric condition To establish the uniqueness of limit cycle arising through Hopf bifurcation, we apply the technique of Kooij et al. [18].

Lemma 5.2. If the functions involved with general Liénard system
satisfies the following conditions (i) G(x) = d dx g(x) and h(x) are continuously differentiable within the open interval (r 1 , r 2 ) with r 1 < r 2 , (ii) f (z) is continuously differentiable and increasing function of the variable z in R, then if for all real 'c', g(x) − ch(x) has no multiple zeros and g(x 0 ) = 0, then the system (5.2) has at most one limit cycle within the strip r 1 < x < r 2 .
First we verify the conditions (i), (ii) and (iii) for the said functions. Using the expression of g(x) we get, Here, we take the open interval (r 1 , r 2 ) = (β, x * * ). Obviously, the functions G(x) and h(x) are continuously differentiable functions over (β, x * * ). Also, from the definition of the function f (z) it is clear that f (z) is continuously differentiable for all real z and f (z) = e −z which is always positive in R. Hence, f (z) is increasing in R. There exists a unique positive root Therefore, all the conditions (i)-(iii) hold good. Our final task is to prove the non-existence of multiple zeros of G c (x) = g(x) − ch(x) on (β, x * * ) for all c ∈ R. Now, Hence, x * * ∈ (β, β * ) and H 1 (x) is increasing on (β, x * * ). Case 2. If c < 0 H 1 (x) > 0 and H 2 (x) ≤ 0 for (β, x * ]. Therefore, H c (x) = H 1 (x) + cH 2 (x) > 0, which implies that the roots of H c (x), if they exist, belongs to the interval (x * , x * * ). On the other hand H 1 and H 2 are both increasing on (x * , x * * ). Then In addition, H 1 and H 2 are both concave down. Hence, H 1 (x) and H 2 (x) are both decreasing on (x * , x * * ). So, we have Suppose that . Under one of these condition we have, That is H c (x) = 0. So, in both cases H c (x) does not have any multiple root. Thus the system possesses at most one limit cycle within the strip β < x < x * * .

Existence of heteroclinic orbit
In this section, we assume that β β+δ 2 < δ 1 < 1 1+δ 2 , and hence the existence of unique equilibrium point in the interior of first quadrant is guaranteed. Moreover, the equilibrium point E 0 is local attractor and the other two equilibrium points E 1 and E 2 are saddle points. Let W s (β, 0) and W u (1, 0) be the stable and unstable manifolds of E 2 and E 1 respectively.

Lemma 5.3.
There exists a set of parameter values for which W s (β, 0) = W u (1, 0), giving rise to a heteroclinic orbit joining the saddle points E 1 and E 2 .
The proof of the above lemma can be done in a similar fashion like Lemma 3 of [28]. Hence we omit the proof here. Also, it can be proved that there exists a set of parameter values for which W u (β, 0) = W s (1, 0), giving rise to a heteroclinic curve joining the saddle points E 1 and E 2 . We have shown the existence of heteroclinic curve for the given system numerically in Fig. 5.2. The parameter values are given in the caption of the figure. In Fig. 5.2, the red solid circle is origin (attractor) and two black circles are axial equilibrium points which are saddle. Let γ 1β denote the heteroclinic curve for which W s (β, 0) = W u (1, 0) and γ β1 denote the heteroclinic curve for which W u (β, 0) = W s (1, 0). Therefore, we have identified only one heteroclinic orbit γ h , where γ h is given by

Transcritical bifurcation
In this subsection, we show the occurrence of transcritical bifurcation at the equilibrium points E 1 and E 2 . An application of Sotomayer's theorem [30] gives the confirmation of occurring the transcritical bifurcation at E 1 and E 2 when δ 1 = 1 1+δ 2 and δ 1 = β β+δ 2 respectively. Theorem 5.4. The system (2.7)-(2.8) undergoes transcritical bifurcation between E 1 and E * and E 2 and E * when δ 1 crosses the threshold δ 1 = 1 1+δ 2 and δ 1 = β β+δ 2 . Proof. Let Ω(x, y; δ 1 ) represent the vector Differentiating the above function with respect to δ 1 we get The Jacobian matrix J E 1 of the system (2.7)-(2.8) around the axial equilibrium point E 1 = (1, 0) is given by For δ 1 = δ 1 * = 1 1+δ 2 , the above Jacobian matrix becomes Clearly, the matrix P has a simple zero eigenvalue and one negative eigenvalue. Let V = (v 1 , v 2 ) T be the eigenvector of the matrix P corresponding to zero eigenvalue. Then V is given by The eigenvector of P T (transpose of the matrix P) corresponding to zero eigenvalue is given by From equation (5.7) we have, and hence Now we consider (5.14) Hence, and therefore we have The expansion of the term D 2 Ω(x, y; δ 1 )(V, V) is given by Hence, Using (5.10) and (5.11) we get, Using Sotomayer's theorem and the equations (5.13), (5.15), (5.16), we conclude that the preypredator system (2.7)-(2.8) experiences transcritical Bifurcation at E 1 .
In a similar fashion, we can also prove the occurrence of transcritical bifurcation at E 2 for the predator-prey system (2.7)-(2.8).
The unique interior equilibrium point is emerged or annihilated due to the occurrence of two transcritical bifurcations at the axial equilibrium points E 1 and E 2 respectively.

Global dynamics
In the previous section, we have studied some local bifurcations associated with the system (2.7)-(2.8) under several parametric restrictions. We observed that all the parameters except γ, plays an important role in the study of bifurcation theory for considered model. Here we draw a schematic bifurcation diagram in δ 1 δ 2 -parametric space to understand how the bifurcation curves divide the positive quadrant of δ 1 δ 2 plane into different subregions where the given model shows qualitatively different dynamic behaviours.
The schematic bifurcation diagram of model (2.7)-(2.8) for β = 0.2 is presented in Fig. 5.1. This diagram contains Hopf bifurcation curve (green colour curve), two transcritical bifurcation curves δ 1 = 1 1+δ 2 (blue colour curve) and δ 1 = β β+δ 2 (red colour curve). These three curves divide the δ 1 δ 2 -parametric space into four different subregions, named by R 1 , R 2 , R 3 , R 4 and we get four qualitatively different phase portraits for chosen set of parameter values from each region. All boundary equilibrium points exist in each region as their existence do not depend on parametric restrictions. In the subregion R 1 , we have the parametric restriction β β+δ 2 and in the subregion R 4 , δ 1 < min 1 1+δ 2 , β β+δ 2 . Therefore, the interior equilibrium point does not exist in R 1 and R 4 . At δ 1 = 1 1+δ 2 , the interior equilibrium point E * emerges through transcritical bifurcation as δ 1 enters into the domain R 2 . Now as we decrease the value of δ 1 the system will be continued to have the interior equilibrium point for β β+δ 2 < δ 1 < 1 1+δ 2 . Therefore, the system has also the unique interior equilibrium point in R 3 . The interior equilibrium point disappears through transcritical bifurcation at δ 1 = β β+δ 2 . Hence, no interior equilibrium point exists in region R 4 . Now we discuss the phase portraits for a particular choice of parameter set given by β, γ and varying δ 1 and δ 2 to understand the stability and bifurcation analysis very well. We choose β = 0.2, γ = 0.2. These parameter values are chosen from [31]. Values of δ 1 and δ 2 will be chosen in such a way that the point (δ 1 , δ 2 ) belongs to each of the four different domains presented in the schematic bifurcation diagram (see Fig. 5.1). Phase portraits corresponding to each domain are presented in Fig. 5.3. Parameter values for δ 1 and δ 2 for different domains are given in the caption of the Fig. 5.3. In all the phase portraits, stable equilibrium points and stable limit cycle are marked with red colour, unstable equilibrium points are marked with small black circle, manifolds of equilibrium points are marked with green colour but the separatrix between two domains of attraction are marked with magenta colour. In region R 1 , E 1 and E 0 are stable and E 2 is a saddle point. The trajectories with initial conditions above the stable manifold of E 2 have the point origin as ω-limit, whereas the orbits starting below the stable manifold have the point E 1 as their ω-limit. In region R 2 , the unique interior equilibrium point E * and E 0 are stable but E 1 and E 2 are saddle in both R 2 and R 3 . The interior equilibrium point looses it's stability through Hopf bifurcation in region R 3 and hence a limit cycle will be created in this domain around the unstable interior equilibrium point E * . One can easily check that for the given set of parameter values in domain R 3 , the value of σ is −.3964348492 < 0. Hence, the above mentioned limit cycle in R 3 is stable and the Hopf bifurcation through which it occurs is supercritical. The trajectories starting below the stable manifold of E 2 have the point E * and stable limit cycle as their ω-limit in region R 2 and R 3 respectively but for the trajectories initiating above the stable manifold have the point (0, 0) as their ω-limit in both regions R 2 and R 3 . In region R 4 , E 1 is a saddle point and E 2 is unstable and the trivial equilibrium point E 0 is globally asymptotically stable under the parametric restriction δ 1 < min 1 1+δ 2 , β β+δ 2 .

Conclusions
In this paper, we have studied the dynamics of a prey-predator model with Holling type II functional response, where the prey population is subjected to the most common form of Allee effect (strong Allee effect). We observed that different dynamical behaviours can exist for this type of model such as: stable coexistence of both species, persistent oscillations, predator extinction or extinction of both species. For clear understanding of this paper, we list all properties which we have obtained for our considered model. (ii) The most significant result of our analysis is that the unique interior equilibrium point is either locally asymptotically stable or lead to a stable limit cycle.
(iii) There exist parameter restrictions for which positive equilibrium point is locally stable.
(iv) There are parameter conditions for which one limit cycle exists. This corresponds to the occurrence of Hopf bifurcation when change of stability occurs at the unique positive equilibrium point.
(v) A range of parameters also exists wherein multiple stable states occur, i.e. bistability scenario is observed.
(vi) There exists a separatrix curve determined by the stable manifold of the saddle equilibrium point (β, 0) which divides the behaviour of the trajectories.
(vii) The trajectories which are close to the separatrix curve are highly sensitive to the initial conditions. Those which starts from the below of separatrix curve will converge to the stable equilibrium point or stable limit cycle. Those which starts from the above of separatrix curve will converge to the point (0, 0). The sensitivity of the bifurcation structure on the formulation of the Allee effect is considered in [31], where the ratio dependent model is analyzed.
Some ecological interpretations of these properties are as follows.
(a) From property (i), it is clear that both species can go to extinction at low population densities, irrespective of any parameter values.
(b) However, extinction could occur even for high initial population densities of prey (near carrying capacity K).
(c) Also there exists a region where both populations coexist and this region is determined by the x axis and the stable manifold of the saddle point (β, 0). There, the populations coexist and tend to the unique stable interior equilibrium point (from property (iii)) or they have oscillatory behaviour (from property (iv)).
(d) From property (vii), the behavior of the trajectories originated in the neighbourhood of separatrix curve depend on the initial sizes of both populations,implying that small disturbances due to environmental changes caused by pollution or other reasons, could take the attempt of extinction of both species.
Also, it is observed that the region of extinction of both species gradually increases from R 1 to R 4 and for R 4 it is whole xy plane. It is worthwhile to mention here that the existence of limit cycle does not imply the instability of interior equilibrium point in general. In [14], authors have demonstrated the existence of two limit cycles in a Gauss type predator-prey model with sigmoid functional response and Allee effect on prey. There are many real ecological models which show bi-stability behaviours and investigation of these models that exhibit this behaviour lead to insights about threshold and breakpoint behaviour. This means that the system has capability to model the population explosions (outbreaks) and to crash the stable predator-prey relationship [7].
Also, we would like to explore here briefly the difference between the dynamical properties of models for c = 0 and c = 0. If we take c = 0, the value of δ 2 in system (2.7)-(2.8) will be decreased. As a result the level of interior equilibrium point is going to be changed and the region of existence of interior equilibrium point will be enlarged than c = 0. Also, the value of parameter δ 1 at Hopf bifurcation threshold will be increased for c = 0 compare to c = 0. As the existence of heteroclinic orbit depends on the parameter value δ 2 , hence for c = 0 we can not be ensure about it's existence with the decreasing value of δ 2 . For transcritical bifurcation, the threshold value will be higher for c = 0 compare to c = 0.
As a future direction, we can first mention here that the analytical proof of the absence of stable limit cycle in model (2.7)-(2.8) and this can be approached based on divergence criterion as in [22]. It has been observed that in some other prey-predator models like Leslie-Gower type models can result different number of limit cycles by changing the mathematical formulation of Allee effect [12]. In our considered model we can also expect to have the sensitivity of the bifurcation structure on the mathematical form of Allee effect. The final interesting question is concerned with the explicit outcome of the model (2.7)-(2.8) if we add the spatial dimension in the model. There exists several prey-predator models which show extinction scenario in non-spatial form but they can exhibit persistence properties for a large range of parameter values when they are in spatial form [27]. Thus it is very much important to examine how the parametric range of persistence in model (2.7)-(2.8) will be extended when we allow species to disperse.