Dynamics of prey–predator model with strong and weak Allee effect in the prey with gestation delay

. This study proposes two prey–predator models with strong and weak Allee effects in prey population with Crowley–Martin functional response. Further, gestation delay of the predator population is introduced in both the models. We discussed the boundedness, local stability and Hopf-bifurcation of both nondelayed and delayed systems. The stability and direction of Hopf-bifurcation is also analyzed by using Normal form theory and Center manifold theory. It is shown that species in the model with strong Allee effect become extinct beyond a threshold value of Allee parameter at low density of prey population, whereas species never become extinct in weak Allee effect if they are initially present. It is also shown that gestation delay is unable to avoiding the status of extinction. Lastly, numerical simulation is conducted to verify the theoretical ﬁndings.


Introduction
Prey-predator relationship is a very dominant phenomenon that occurs in nature. It has been an issue of attention among ecologists and biologists since last few decades. First model on prey-predator interaction is formulated and proposed by Lotka and Volterra. It contains a pair of first-order nonlinear differential equations frequently used to describe the dynamics of biological systems in which two species interact. After that several attempts have been made to generalize and extend these equations.
A general two dimensional model of interaction between prey and predator is represented by dx dt = xf (x) − yg(x, y), dy dt = y −d + cg(x, y) , where x and y denote prey and predator densities at time t respectively. f (x) is per capita growth rate of prey. g(x, y) and cg(x, y) are functional and numerical response of predator for prey, where c (0 < c < 1) stands for conversion coefficient denoting the number of newly born predators for each captured prey. d is mortality rate of predator population.
One vital factor of the prey-predator interaction is the intake rate of prey by a predator i.e. functional response. It helps to predict about a prey-predator dynamics with more accuracy. There are many types of functional response: Holling type I-III, ratio dependent, Beddington-DeAngelis, Crowley-Martin, Hassel-Verley. Holling type I-III functional responses are prey dependent whereas Beddington-DeAngelis, Crowley-Martin, Hassel-Verley are prey and predator dependent i.e. functional response is function of both the prey and predator's density.
Crowley-Martin [11] assumed that predation will decrease when the predator density is high due to interference among predators. Some investigations have been conducted on prey-predator model including Crowley-Martin functional response [11,25,38]. This type of response function is written as where α, a and b are positive parameters denoting attack rate, handling time and magnitude of interference among predators, respectively. The effect of intraspecific interference among predators has been investigated in preypredator model with Holling type II functional response in [37,42], with Holling type III functional response in [17], with Beddington-DeAngelis type functional response in [23,29].
Allee effect plays a major role in the structure of population. It creates the possibilities of extinction of species [13,42] and has a huge impact in population dynamics [1]. The Allee effect can be classified into two types on the basis of per capita growth rate at low density. These are known as strong Allee effect and weak Allee effect. Strong Allee effect have negative per capita growth rate at low population level and implies the existence of a threshold level of population so that the species become extinct below this level. Recently, Verma and Misra [40] have studied the impact of a constant prey refuge on the dynamics of a ratio-dependent predator-prey system with strong Allee effect in prey growth. They found that if prey refuge is less than the Allee threshold, the incorporation of prey refuge increases the threshold value of the predation rate and conversion efficiency at which unconditional extinction occurs. They also vindicated that the species can be protected by creating safe zones in accordance with the Allee threshold. On the other hand, in weak Allee effect, the per capita growth rate decreases but remains positive at low population level. Sexual selection [3,27], reduced mating efficiency [12] and alleviated foraging efficiency [2] are some other reasons to give rise to Allee effect. Figure 1 makes us more clear that initially, the per capita growth rate is negative in strong Allee effect (blue) while it remains positive in weak Allee effect (brown).
The Allee effects are observed in many natural species. For example, in plants [15,18], insects [22], marine invertebrates [35], in birds and mammals [10]. Recently, many ecologists have paid attention to the Allee effect [9,13,14,16,28,30,36,41]. Some crucial results have been investigated in [7,42] via a comparative analysis of prey-predator  system with and without Allee effect. Some studies have been conducted on strong Allee effect on prey-predator model [17,29,30]. Weak Allee effect has been well studied by ecologists [16,32,33]. Some researchers have found the natural evidence of weak Allee effect by experimental work on flour beetles of the genus Tribolium [1,21,39]. They have shown that the per capita growth of beetles reaches its maximum at a medium density and the rate is positive at low density. Hopf-bifurcation is an important tool which helps to understand the behavior of system. It gives us such a critical value of a parameter that the stability behavior of system is contrasty in both the sides of the critical value [16,23]. Hopf-bifurcation analysis with Allee effect has been carried out in [5,16,29,31].
Time delay occurs in every biological movement. A delay differential equation shows much more complicated behavior than an ordinary differential equation. Delay is capable to change the stability behavior of any system. Due to time lag in conversion of prey population to predator population (gestation delay), dynamics of system changes. The prey-predator population model with gestation delay [6,8,20,26,38] [24,25] has been studied. Some authors [32,33] have considered an eco-epidemiological model with weak Allee effect in prey-predator population. They have concluded that chaotic dynamics can be controlled by the Allee parameter. Further Biswal et al. [6] have applied gestation delay and observed that the system exhibits chaotic oscillation due to increase of the delay.
Some studies have been conducted with Allee effect including gestation delay [4,5,29,31]. Li et al. [23] investigated the stability and Hopf-bifurcation of a delayed density dependent prey-predator system with Beddington-DeAngelis functional response. Pal and Mandal [29] analyzed a modified delayed Leslie-Gower prey-predator model with strong Allee effect and shown that delay is incapable to decrease the risk of extinction. A prey-predator model with Crowley-Martin functional response including prey refuge has been considered by Maiti et al. [25]. They also examined the effect of gestation delay on the dynamics of the system.
To the best knowledge of the authors, a comparative analysis between strong and weak Allee effect in a prey-predator model with Crowley-Martin type functional response and gestation delay has not been studied. The main motive of this work is to analyze the dynamical complexity of Allee effect in the prey-predator model, studied by Tripathi et al. [38], and further to show the impact of gestation time delay on the dynamics of the system. Keeping all these in mind, we reconstruct the model described by Tripathi et al. [38] by incorporating strong and weak Allee effects. Then we introduce gestation delay in the predator population.

Model formulation
Tripathi et al. [38] have analyzed the dynamics of following density dependent nonlinear mathematical model In this model prey population grows logistically and predator is survived only on the prey population. They follow Crowley-Martin functional response to hunt the prey population. r, K, α, a, b, δ 0 and δ 1 are positive constants representing intrinsic rate of prey, carrying capacity, capture rate, handling time, magnitude of interference among predators, natural death rate of predators and crowding effect, respectively. Now at low and sparse population, prey exhibits strong Allee effect. Let θ s be the Allee parameter and f an auxiliary parameter, which shapes the Allee function. The prey-predator dynamics with strong Allee effect in prey population is governed by the following system: On the other hand, the model with weak Allee effect is based on probability of successful mating of prey population. It is incorporated into the population growth model by multiplying the probability P(x) with birth term of prey population, where P(x) is the probability of successful mating for a female prey during the reproductive period and should follow the bellow criteria: 1. No mating occurs at zero population size, P(0) = 0. 2. P (x) > 0 i.e. if population size increases the probability of successful mating increases. 3. Mating is guaranteed when the population is sufficiently large, that is P(x) → 1 as x → ∞.
We consider the probability function as P(x) = x/(θ w + x), θ w > 0 (rectangular hyperbolic) [12,32]. Thus, model (1) with weak Allee effect can be written as In real situation, the conversion of hunted prey into predators growth is not instantaneous process rather, there occurs a time lag for gestation of predator biomass. Therefore, we assume that the reproduction of predator population after hunting prey is arbitrated by a constant time lag, called gestation delay. In order to get the rich dynamics of the system, we introduce gestation delay τ s and τ w in models (2) and (3), respectively. Then model (2) takes the form subject to the nonnegative condition Similarly, in the presence of gestation delay, model (3) can be written as subject to the nonnegative condition

Dynamics of nondelayed systems
In this section, we will study the dynamics of models (2) and (3).

Positivity and boundedness of model system (2)
It is necessary to prove that the model is biologically well behaved before the detailed study.
From model system (2) we can write ds , which shows that all solutions remain within the first quadrant of the xy-plane starting from an interior point. In the following theorem, we show that all solutions of system (2) are bounded, which refers that the model is biologically well behaved.
We also note that dW/dt < 0 if W > 2rK/δ. Hence all solutions of system (2) point towards Ω. Thus, Ω is a positively invariant set and all the solutions of model (2) are bounded.

Local stability and Hopf-bifurcation
In this subsection, first we will find out all feasible equilibrium points of system (2) and present all possibilities for interior equilibrium. Then a brief description on their local stability has been done and lastly the analysis of Hopf-bifurcation through local stability of the positive equilibrium has been carried out.

Existence of equilibrium points
System (2) has following equilibrium points: 1. The trivial equilibrium point E 0 (0, 0). 2. The axial equilibrium points E 1 (θ s , 0) and E 2 (K, 0). 3. System (2) has a unique positive equilibrium point E * (x * , y * ) if the following condition holds true: Remark 1. The number of positive equilibrium for system (2) depends on values of parameters, which we have chosen. Several possibilities are depicted in Fig. 2.

Local stability analysis
To analyze the local stability behavior of the equilibria, whenever they exist, we compute the Jacobian matrix for the model system (2), and further this matrix is calculated at each of equilibria. Then using Routh-Hurwitz criteria, we get following results: is saddle point having stable manifold along the x-axis and unstable manifold along the y-axis if cαK > δ 0 (1 + aK).
In order to investigate the stability behavior of positive equilibrium, let M (E * ) denotes the variational matrix evaluated at E * Then the characteristic equation of M (E * ) is given by where Using the Routh-Hurwitz criteria, both the eigenvalues of M (E * ) have negative real part if and only if A 1 > 0, A 2 > 0. Thus we can state the following theorem.
Consequently the interior equilibrium E * is asymptotically stable.
In equation (7), if we assume A 2 < 0, then one eigenvalue is real and positive and other eigenvalue is real and negative. Thus, the following theorem follows: Now, assume that A 1 < 0, and A 2 > 0. Then both the eigenvalues are real and positive, or both the eigenvalues are complex conjugate having positive real parts. Thus, we can state the following theorem. Theorem 4. If A 1 < 0, and A 2 > 0, then the interior equilibrium E * is unstable.
Remark 3. We have seen that E 0 is always stable equilibrium and interior equilibrium E * is stable if A 11 < 0. Hence, if A 11 < 0 holds, then system (2) exhibits bistability. Bistability is a phenomenon where the system can converge to different equilibrium in the same parametric region based on the variation of initial condition.

Hopf-bifurcation analysis
Here auxiliary parameter f is an important parameter which shapes the Allee function. In this subsection, we analyze how the dynamics of system (2) changes with respect to f by using Hopf-bifurcation analysis. Now, to show the existence of Hopf-bifurcation, we assume that and it leads us to the following theorem.
Theorem 5. Assume that f = f * and A 2 > 0. Then system (2) has Hopf-bifurcation near the equilibrium point E * (x * , y * ) if following condition is satisfied: Proof. At f = f * , we have tr(M (E * )) = 0 and det(M (E * )) > 0, which shows the eigenvalues are purely imaginary and conjugate to each other at f = f * . We also have Hence, the transversality condition holds under the condition This shows that E * changes its nature from locally asymptotic stable to unstable as parameter f crosses the critical value f = f * . Therefore by the Hopf-bifurcation theorem, system (2) exhibits Hopf-bifurcation near the interior equilibrium point E * .

Dynamics of weak Allee effect model
The analysis of model (3) with weak Allee effect is similar to that of model (2) with strong Allee effect. Thus, in this section, we omit the detail mathematical analysis and present the main results of this model system briefly.
1. All the solutions of model system (3) with initial conditions that initiate in R 2 + are positive invariant and uniformly bounded. (3) has following equilibrium points: E 0 (0, 0), E 2 (K, 0) and the unique interior equilibrium E * (x * , y * ). E * exists if the following condition holds:

System
, and a saddle point having stable manifold along the x-axis and unstable manifold along the y-axis if cαK > δ 0 (1 + aK).
4. In order to analyze the stability behavior of interior equilibrium, let M (E * ) denotes the variational matrix evaluated at E * . Then Then the characteristic equation of M (E * ) is given by Using the Routh-Hurwitz criteria, we have: • System (3) is locally asymptotic stable around the interior equilibrium E * if and only if A 1 > 0, A 2 > 0. • It can be noted that if A 11 < 0, then the interior equilibrium E * is locally asymptotically stable. The per capita growth rate of prey population is negative.
The per capita growth rate of prey population is slower then without Allee effect but remains positive (see Fig.1). 2 Model (2) has locally asymptotically stable trivial equilibrium point. Consequently possibility of extinction is high at low population density.
Species always coexist if they initially exist under weak Allee effect (model (3)). 3 Model (2) shows bistability when E * is locally asymptotically stable. If E * is unstable then population moves to (0, 0).
Model (3) never shows bistability. When interior equilibrium is unstable then populations fluctuate around it. 4 Model (2) does not show the Hopf-bifurcation with respect to Allee parameter θs.
Model (3) shows the Hopf-bifurcation with respect to Allee parameter θs.
5. Now, to show the existence of Hopf-bifurcation, we assume that Thus, we can state the following.
Theorem 6. Assume that θ w = θ * w , A 2 > 0 and x * > θ * w . Then system (3) has Hopfbifurcation near the equilibrium point E * (x * , y * ). Now, we are in position to compare the system with strong Allee effect and weak Allee effect. In Table 1, comparison between both the cases has been carried out.

Local stability and Hopf-bifurcation of delayed models
In this section, we will investigate local stability of the positive equilibrium E * and exhibition of local Hopf-bifurcation. We know that delay does not affect the equilibrium of the system. Therefore equilibrium points are same as nondelayed model system. We omit the proof of the following theorem as it is similar to nondelayed system. To see the effect of time delay on the dynamics of the system, we can rewrite model system (4) as Let x(t) = x * + x (t), y(t) = y * + y (t). Then linearizing system (4) about the interior equilibrium solution E * (x * , y * ), we have Thus, the variational matrix of system (4) at E * is given by and corresponding characteristic equation is where A = a 1 − a 3 + a 4 , B = −ca 2 , C = ca 2 a 3 and D = a 4 (a 1 − a 3 ).
Remark 4. The characteristic equation is same as the characteristic equation (7) of the nondelayed model system (2) studied earlier.
All the roots of characteristic equation (11)

Stability and direction of Hopf-bifurcation
In the previous section, we obtained the condition under which periodic solution bifurcates from the steady state at the critical value of τ s . In this section, we will study the direction of Hopf-bifurcation and stability of the periodic solution by using normal form theory and center manifold theory introduced in Hassard et al. [19]. We assume that system (4) undergoes Hopf-bifurcation at the steady state E * for τ s = τ s0 and ±iω 0 is corresponding purely imaginary roots of the characteristic equation at E * . Let x 1 (t) = x(t) − x * , y 1 (t) = y(t) − y * and still denote x 1 (t), y 1 (t) by x(t), (t). Let τ s = τ s0 + µ, µ ∈ R, so that Hopf-bifurcation occurs at µ = 0, system (4) is transformed into , Here we omit the detailed analysis and write only the results, which are obtained. One can easily derive them by using the computation process similar to that in Song and Wei [34] and Tripathi et al. [38]. The standard results can be computed as where g 20 , g 11 , g 02 and g 21 are evaluated as follows: 110 , 2 ) T ∈ R 2 are constant vectors computed as These expressions give a description of the bifurcating periodic solution in the center manifold of system (4) at critical values τ s = τ s0 which can be stated as follows: 1. µ 2 determines the direction of Hopf-bifurcation. If µ 2 > 0 (< 0) then the Hopfbifurcation is supercritical (subcritical). 2. β 2 determines the stability of bifurcated periodic solution. If β 2 > 0 (< 0) then the bifurcated periodic solutions are unstable (stable). 3. T 2 determines the period of bifurcating periodic solution. The period increases (decreases) if T 2 > 0 (< 0).

Remark 6.
We can also analyze the properties of bifurcating periodic solution for weak Allee case by adopting the same process.

Numerical simulation
In this section, we will present numerical simulations to validate the analytical findings, obtained in previous sections using MATLAB R2017a.

Strong Allee effect
For the set of values of parameters (15), condition (6) is satisfied. Thus, there exists a positive equilibrium E * (6.4758, 5.1441). It is also noted that A 11 < 0 holds for the set of parameters chosen in (15). So the interior equilibrium E * is asymptotically stable, depicted in Fig. 3. This figure shows that the density of prey and predator species both are increasing initially, then some fluctuations occur and eventually settle down to their respective steady states. The trivial equilibrium point E 0 is always asymptotically stable for system (2). Therefore, under the condition of stability of positive equilibrium E * , system shows bistability. A separatrix lies in the xy-plane,which divides the plane into two regions in such a way  that trajectories starting from different regions approach to different steady states. This phenomenon of system (2) is depicted in Fig. 3(b). From the figure, it is clear that solution curves, which are initiated from left of the separatrix, approach to E 0 (0, 0) and solution curves, which are initiated from right of the separatrix, approach to interior equilibrium point E * . The effect of parameter δ 1 on both the prey and predator species is shown in Fig. 4. In Figs. 4(a) and 4(b), time series analysis is shown for four different values of δ 1 (δ 1 = 0.15, 0.7, 1.5, 2). From Fig. 4(a) it is noted that prey population increases with the parameter δ 1 . The predator population initially grows up with δ 1 but after a threshold value of δ 1 =δ * 1 = 0.825 it starts to decrease and settles down at its equilibrium level (see Fig. 4(b)). Now, we observe the dynamical behavior of the system for the variation of the Allee parameter θ s . It is noted that as we increase the value of parameter θ s , the time of fluctuations for both prey and predator increase. Time series analysis has been shown in Figure 5 for different values of θ s , which refers that the system is stable for θ s = 0.1, 0.5 and 1. At θ s θ * s = 1.4985, interior equilibrium E * becomes unstable and beyond this threshold value of θ s , system converges to the stable trivial equilibrium point E 0 . This shows that as we increase Allee parameter θ s , the life expectancy of both biological species decreases and after a critical value of θ s (= θ * s ), they move to extinction. In the model system (2), the auxiliary parameter f is also a crucial parameter because it shapes the Allee function. Therefore we will analyze how the dynamics of system changes with respect to f by using Hopf-bifurcation analysis. The condition of Hopfbifurcation, which is derived in Theorem 5, is satisfied. The critical value of parameter f , where bifurcation occurs is calculated from equation (8) and it is f = f * = 3.22. In Fig. 3, we depicted time series (Fig. 3(a)) and phase portrait (Fig. 3(b)) for f = 0.002 < f * = 3.22, which refers that the system (2) is stable. Figure 6 shows that the system is unstable for f = 3.7 > f * = 3.22 and periodic solution occurs.

Weak Allee effect
On the other hand for system (3), with weak Allee effect in prey population and with same values of parameters as in (15), condition for the existence of unique equilibrium E * (x * , y * ) is satisfied and it is given by E * (6.4782, 5.1464). The conditions of Theorem 6 are also satisfied. Therefore, Hopf-bifurcation with respect to θ w occurs near the interior equilibrium E * . The threshold value of θ w is evaluated as θ * w = 3.36. Thus the equilibrium point E * is asymptotically stable for θ w = 2.5 < θ * w which is shown in Figure 7 and unstable for θ w = 4.3 > θ * w (Fig. 8) and a periodic solution exists around E * . The bifurcation diagram has been shown in Fig. 9 by taking θ w as a bifurcation parameter. This figure depicts the dynamics of the system as the Allee parameter increases. From the figure, it is clear that for θ w < θ * w , the system (3) is stable but as θ w crosses its critical value, the system loses its stability and undergoes Hopf-bifurcation. We have drawn phase portrait for model (2) and (3) Figure 11. System (4) is locally asymptotically stable when τs = 0.16 < τs 0 , other parameters are same as in (15). parameters as that in (15). We know that introduction of delay does not affect equilibrium of the system. Thus, the interior equilibrium E * (6.4758, 5.1441) remains as it is.
For τ s > 0, we note that conditions (H1) and (H3) are satisfied. So, equation (13) has a unique positive root. Taking i = 0 in equation (14), our computer simulation yields the following: ω 0 = 1.4472, τ s0 = 0.3914, and transversality condition (H5) is satisfied. Here all three conditions of Theorem 8 hold true. Therefore system undergoes a Hopf-bifurcation at τ s = τ s0 = 0.3914. By the algorithm obtained in Section 5, we computed the following: This shows that the Hopf-bifurcation with respect to τ s is supercritical, bifurcating periodic solution is stable and the period increases. Thus, the system is stable for τ s = 0.16 < τ s0 = 0.3914 which is shown in Fig. 11. As τ s passes through its critical value τ s0 , the system loses its stability and a Hopf-bifurcation occurs into the system. In Fig. 12(a), we have shown time series analysis for τ s = 0.45 > τ s0 = 0.3914. Figure 12(b) shows that a periodic solution exists and any solution trajectory initiating from inside and outside the closed trajectory, approaches towards the closed trajectory. This shows the existence of a stable limit cycle. Bifurcation diagram has also been carried out in Fig. 13 by taking τ s as a bifurcation parameter. Figure makes us clear that τ s changes the stable behavior of the system into instable behavior.

Weak Allee effect
The model system (5), with weak Allee effect in prey population has one interior equilibrium E * (6.4782, 5.1464), with set of parameters (15). We see that the conditions (H1 )  and (H3 ) are satisfied and we obtain and transversality condition (H5 ) is also satisfied. Therefore, system (5) undergoes a Hopfbifurcation around interior equilibrium at τ w = τ w0 = 0.3921 (Theorem 9). Using algorithm derived in previous section, it is obtained This shows that the Hopf-bifurcation with respect to τ w is also supercritical, bifurcating periodic solution is stable and its period increases. Thus, the positive equilibrium E * is asymptotically stable for τ w = 0.22 < τ w0 = 0.3921 which is shown in Figure 14(a) and  unstable for τ w = 0.48 > τ w0 = 0.3921 ( Fig. 14(b)). When τ w = τ w0 , system undergoes a Hopf-bifurcation at the positive equilibrium E * . The phase portrait has been shown in Fig. 14(b), which shows the existence of a stable limit cycle.

Conclusion
In this work, we made an attempt to discuss the impact of Allee effect (strong and weak both) with gestation delay in the model proposed by Tripathi et al. [38]. They analyzed a density dependent nonlinear mathematical model (1). In that model, prey grows logistically and predator fully depends on prey for food that follows Crowley-Martin functional response.
Allee effect plays an important role in the structure of population. The Allee effect increases the possibilities of extinction. Thus, we include Allee effect into model (1). Since there are two types of Allee effect; strong and weak, so we studied both the models separately. In the study, we discussed positivity, boundedness of the solutions, existence of equilibrium points and their stability analysis of both the models. Positivity and boundedness of the solutions refer that the system is well behaved.
We have shown that system (2) and (3) may have more than one interior equilibrium point. Under sufficient condition (6) (for system (2)) and (9) (for system (3)) they have unique interior equilibrium point. We also derived a sufficient condition for asymptotic stability of interior equilibrium point. Then we found that model system (2) is bistable in the presence of positive equilibrium. The existence of periodic solution via Hopfbifurcation with respect to auxiliary parameter f in model (2) and Allee parameter θ w in model (3) have also been shown. In Table 1, we presented that how the dynamics of model (2) differs from model (3). We also observed that the time of fluctuations for both prey and predator increases with increase in Allee parameter for system (2). But after a critical value, interior equilibrium E * becomes unstable. In this situation system converges to stable trivial equilibrium E 0 , which shows extinction of species after a critical value of Allee parameter.
Delay exhibits much more realistic behavior. The reproduction of predator after hunting prey is not instantaneous i.e. there is some time lag for gestation. Therefore to make the model biologically more realistic, we consider gestation delay for predator into both the models (2) and (3). We have analyzed Hopf-bifurcation through local stability considering delay as a bifurcation parameter. When the time delay is small then trajectory of system oscillates around the positive equilibrium for finite time and eventually settle down to equilibrium level. As the time delay increases, time of oscillations also increases and beyond a critical value of gestation delay, then stability of system switches and we obtain periodic solutions. This proves that the time delay can cause a stable equilibrium to become unstable. The stability and direction of Hopf-bifurcation also have been investigated using Normal form theory and Center manifold theory.
The numerical simulation is based on some biologically feasible data to support our theoretical results. We found that Hopf-bifurcation is supercritical and stable with increasing period. Bifurcation diagram with respect to θ w and τ s help us to understand about the stability behavior of the system. This study has some new and significant results that we hope very helpful to understanding the dynamics of prey-predator system with Allee effect and gestation delay.