Dynamics of a predator-prey model with strong Allee e ﬀ ect and nonconstant mortality rate

: In this paper, dynamics analysis for a predator-prey model with strong Allee e ﬀ ect and nonconstant mortality rate are taken into account. We systematically studied the existence and stability of the equilibria, and detailedly analyzed various bifurcations, including transcritical, saddle-node, Hopf and Bogdanov-Takens bifurcation. In addition, the theoretical results are veriﬁed by numerical simulations. The results indicate that when the mortality is large, the nonconstant death rate can be approximated to a constant value. However, it cannot be considered constant under small mortality rate conditions. Unlike the extinction of species for the constant mortality, the nonconstant mortality may result in the coexistence of prey and predator for the predator-prey model with Allee e ﬀ ect.


Introduction
Interactions between a wide variety of species are a feature of ecosystems. Understanding the relationship between predator and prey, which plays a key status in predicting the dynamics of ecosystems. Mathematical models are increasingly influential in theoretical ecology, contributing not only to the quantitative understanding of ecosystems [1][2][3][4][5] but also the development of mathematical capability [6][7][8][9][10][11]. A classical general predator-prey model is shown below [12]: where the densities of predator and prey are respectively expressed as P and N, the per capita growth rate of predator and prey are respectively denoted as h(g(N, P), P) and f (N), g(N, P) is the functional response.
There are many key factors that affect the dynamics in predator-prey models, especially functional response and nonconstant mortality rate of the predator. Functional response is generally of two kinds: ratio-dependent and prey dependent. Holling [13] introduced the concept of functional response and the study of small-mammal predation of the European Pine Sawfly1 revealed the composition of predation. Subsequently, there have been many discussions on the effects of functional response on predator-prey models [14][15][16]. Abrams and Ginzburg [12] studied the nature of predation and compared the effects of different functional responses on population dynamics. Zhang et al. [17] studied Beddington-DeAngelis functional response, discussed the linear stability and obtained the condition for Turing instability.
Consider h(g(N, P), P) = eg(N, P) − α(P), then we get the following equation: dP dt = (eg(N, P) − α(P))P, where e and α(P) stand for the conversion efficiency and the specific mortality rate of predator without prey, respectively. In most researches, the mortality rate of predator is assumed to be constant, i.e., α(P) = u [18,19]. However, Cavani and Farkas [20] introduced a nonconstant mortality rate of predator α(P) = γ+δp 1+P , where γ > 0 and δ > 0 represent mortality rate at the low density and the maximal mortality, respectively(γ ≤ δ). The nonconstant mortality α(P) is a bounded and increasing function of P. When γ = δ, it is a constant mortality rate [21,22]. About the nonconstant death rate in predator-prey model, many researchers have found that nonconstant mortality has a major impact on the population dynamics [23][24][25][26][27].
Experiments have shown that due to the factors, such as mate finding, predator satiation, cooperative defense and habitat attention, Allee effect exists during the growth of prey species [28,29]. Allee effect is any mechanism guiding a positive correlation between a component of individual fitness and population density. This means that population at low density will increase the risk of extinction because of positive correlation between density and growth rate [30][31][32]. In the past several years, many researchers have shown a keen interest in population models with Allee effect and proposed several interesting mathematical models [33,34]. González-Olivares et al. [35] found Allee effect induces the equilibria to loss their stability, and may be a destabilizing force. The species will extinct due to strong Allee effect when the density of prey population is low, while weak Allee effect does not cause species to die out [32]. Huang et al. [30] illustrated that the Allee effect or the fear effect plays a negligible role in the density of prey, but it influences the predator population greatly.
Besides Allee effect, prey refuge should be considered in predator-prey models. Generally, preys have refuges where they can avoid predators. Therefore, to better describe the interactions between predator and prey, the inclusion of prey refuge is necessary in predator-prey models. As many researchers have made rich achievements in the predator-prey model with prey refuge, it has gradually become the focus of research in this area [36][37][38]. Chen et al. [39] found that prey refuge plays a negligible role in the persistence and extinction of the prey and the predator, but it influences the population density greatly. In [40], the authors concluded that the prey refuge parameter affect the dynamic behavior of the model greatly. As the amount of refuge increases, the prey density rises resulting in an outbreak of prey population.
In this paper, under the assumption that the prey population growth rate obeys logistic type, we describe the predation by using the prey dependent functional response. Additionally, we consider nonconstant mortality, Allee effect and prey refuge, and propose the following predator-prey model: where K and a describe the carrying capacity of prey and the intrinsic growth rate without predators, respectively. m is the Allee effect threshold of the prey species, θ is the refuge protecting of the prey (0 < θ < 1), e and b respectively represent the conversion efficiency of predator by consuming prey and the predation rate, γ and δ are the same as above.
The rest of the article is arranged as follows. In Section 2, we discussed the existence of the equilibria and their stability in detail. Subsequently, the bifurcation analysis of model (2) is given in Section 3. Some numerical simulations are presented to verify the theoretical results in Section 4. Finally, conclusions are presented in Section 5.

Equilibria and their stability
Firstly, model (2) can be rewrited as below: here, By direct calculation, we can get According to the biological significance of the model variables, set of definitions for model (3) is (3), we can draw the following results. About the boundary equilibria of model (3), we have (i) The origin E 0 = (0, 0); (ii) The predator free equilibria E 1 = (m, 0) and E 2 = (K, 0). The positive equilibria should satisfy the following equation Obviously, the above equation is equivalent to y = a Kb(1−θ) (K − x)(x − m) and H(x) = 0. Here, If ∆ > 0, then the root of dH(x) dx = 0 is By direct calculation, we have .
, 0 is where it intersects the X-axis, and for y = f (x), (m, 0) and (K, 0) are where it intersects the X-axis. We need to consider the following situations to find all the positive equilibria: When γ ≥ eb(1 − θ)K, the model has no positive equilibrium. That is to say, if there exist positive equilibria, then γ < eb(1 − θ)K, i.e., H(K) < 0.
Let X = x − m and Y = y, the model can be written as where P 1 (X, Y) and Q 1 (X, Y) are C ∞ functions at least of third order with respect to (X, Y) and Let X = x − a 01 a 10 y, Y = y and τ = a 10 t, the model becomes By Theorem 7.1 in [41], since the coefficient of y 2 is b 02 a 10 − a 01 b 11 0 and a 10 = ma(1 − m K ) > 0, the equilibrium E 1 is a repelling saddle-node. (iii) At E 2 = (K, 0), the eigenvalues of J are λ 1 = −a(K − m) < 0 and λ 2 = eb(1 − θ)K − γ. When γ > eb(1 − θ)K, E 2 is a stable node due to λ 2 < 0. When γ < eb(1 − θ)K, E 2 is a saddle due to λ 2 > 0.
In the following, we will discuss the case γ = eb(1 − θ)K.
According to the translation X = x − K and Y = y to transform E 2 to (0, 0), we have where P 1 (X, Y) and Q 1 (X, Y) are C ∞ functions at least of third order with respect to (X, Y), and y, Y = y and τ = a 10 t, the model becomes xy + Q 2 (x, y).
The coefficient of 2 > 0 and a 10 = −a(K −m) < 0, so the equilibrium E 2 is an attracting saddle-node by Theorem 7.1 in [41].
That is the end of the proof. (ii) When a 10 + b 01 0 and h 20 0, if f 01 > 0, then E * is a repelling saddle-node; else E * is an attracting saddle-node.
Obviously, the eigenvalues of J(E * ) are λ 1 = λ 2 = 0. Through the transformation of By Lemma in [44], we can obtain the equivalent model of model(10) as follows: If f 20 ( f 11 + 2e 20 ) 0, then E * is a cusp of codimension 2; else E * is a cusp of codimension at least 3. Case 2: c 10 + d 01 0.
Obviously, the eigenvalues of J(E * ) are λ 1 = 0 and λ 2 0. Through the transformation of where R 2 (x, y) and S 2 (x, y) are C ∞ functions at least of third order with respect to (x, y) and Redefine τ by τ = f 01 t, we can get where R 3 (x, y) and S 3 (x, y) are C ∞ functions at least of third order with respect to (x, y) and h i j = .
If h 20 0, the equilibrium E * is a repelling saddle-node if f 01 > 0 by Theorem 7.1 in [41]; otherwise it is an attracting saddle-node.
That is the end of the proof.
, then E 4 = (x 4 , y 4 ) may be a focus or a center or a node.
Proof. According to det(J(E)) and tr(J(E)) in Theorem 2.2. ( . The sign of tr(J(E 4 )) is uncertain, so E 4 = (x 4 , y 4 ) may be a focus or a center or a node.
That is the end of the proof.
For the stability of the unique positive equilibrium E 4 in other cases, the analysis is the same as E 4 in Theorem 2.3, and we omit it here.

Bifurcation analysis
We are interested in various possible bifurcations of model (3) in this section, including transcritical, saddle-node, Hopf and Bogdanov-Takens bifurcation.
Proof. We have det(J(E 2 )) = 0 when γ = eb(1 − θ)K. Hence, there is a zero eigenvalue for J(E 2 ). For zero eigenvalue, the eigenvectors corresponding to matrices J(E 2 ) and J T (E 2 ) are denoted as V and W , respectively. After a simple calculation, we have Moreover, For the transversality conditions, we have A transcritical bifurcation occurs at E 2 when γ = eb(1 − θ)K according to Sotomayor's theorem. This is the end of the proof.

Saddle-node bifurcation
From the previous analysis, we can see clearly that under the condition m < x A < K and γ < eb(1 − θ), model (3) has two positive equilibria E 3 = (x 3 , y 3 ) and Therefore, a saddlenode bifurcation will occur at E * = (x * , y * ). Saddle-node bifurcation at E * can be summarized into the following theorem.
Proof. According to Theorem 2.2, there is a zero eigenvalue λ 1 for J(E * ). For zero eigenvalue, the eigenvectors corresponding to matrices J(E * ) and J T (E * ) are denoted as V and W, respectively.
Furthermore, we can get V and W satisfy the transversality conditions A saddle-node bifurcation occurs at E * when γ = γ SN according to Sotomayor's theorem. This is the end of the proof.
For the first Lyapunov number l 1 , which is given by the formula where ∆ = det(J(E 4 )) = l 10 n 01 − l 01 n 10 . The limit cycle is stable via a supercritical Hopf bifurcation if l 1 < 0, and it is unstable via a subcritical Hopf bifurcation if l 1 > 0.  For the proof of the theorem, please see Appendix 5.

Numerical results
In this section, several numerical simulations are provided to verify the reliability of the theoretical results under the parameters conditions as a = 0.8, b = 0.6, e = 0.8, θ = 0.1, δ = 0.9, K = 0.9, m = 0.2. The parameter γ is chosen as the bifurcation parameter.
The possible equilibria of model (3) are discussed in the Section 2 and the qualities of them are also given, through which we plot a bifurcation diagram (Figure 1) and the phase portraits (Figures 2-5). It can be seen from Figure 2(a) that for a small γ, there exist equilibria E 0 , E 1 and E 2 . As γ increases and reaches the line L 1 , E * is the unique interior equilibrium in model (3) (see Figure 2(b)). With the increase of the parameter γ, the positive equilibria E 3 and E 4 appear as the taking place of a saddlenode bifurcation in model (3). From Figure 2(c), it is clearly that both E 3 and E 4 are unstable near the line L 1 , the former is a saddle and the later is a node.
As the parameter γ increases and crosses the line L 2 (γ H = 0.05400367556), a Hopf bifurcation occurs. The first Lyapunov number l 1 = 6.414584574π is greater than zero, which indicates that the direction of Hopf bifurcation is subcritical. Thus, an unstable limit cycle emerges bifurcating from the positive equilibrium E 4 , shown in Figure 3(a),(b). Additionally, the positive equilibrium E 4 recovers its stability owing to the occurrence of Hopf bifurcation. The equilibrium E 4 is an unstable focus in region II (γ < γ H , see Figure 3(c)), and becomes a stable focus in region III (γ > γ H , see Figure 3(d)). In region IV, E 4 is the unique positive equilibrium in model (3), which coincides with the boundary equilibrium E 2 when γ reaches the line L 4 . Specially, on the left of line L 4 (region IV), the boundary equilibrium E 2 is a saddle (see Figure 4(a)), while on the right (region V), E 2 is a stable node (see Figure 4(c)). As γ increases and crosses the line L 4 , E 2 changes from unstable to stable. Therefore, we can conclude that transcritical bifurcation occurs. Moreover, Theorem 3.4 indicates that a BT bifurcation with codimension 2 occurs around the equilibrium E * in model (3). We choose γ and m as the control parameters to numerically study BT bifurcation. Direct numerical computation shows that BT bifurcation may occur when m = m BT = 0.1347130155 and γ = γ BT = 0.003455490134. The equilibrium E * is a cusp when m = m BT and γ = γ BT , which is dispalyed in Figure 5(a). For m = 0.4 and γ = 0.1682189296, the positive equilibrium E * is a repelling saddle-node (see Figure 5(b)). However, it becomes an attracting saddle-node when m = 0.1301 and γ = 0.00005717997942 (see Figure 5(c)).  We present the bifurcation diagram of model (3) in Figure 6 to illustrate all the possible dynamic behaviors in the γ − m parameter plane, where the curve SN, H and T divide the parameter plane into four regions. In region I, model (3) has no positive equilibrium. However, two positive equilibria E 3 and E 4 appear through the saddle-node bifurcation in region II. On the curve H, the model undergoes a Hopf bifurcation which leads to the appearance of a limit cycle. In addition, as the parameters go through the transcritical bifurcation curve (curve T) and enter region IV, the positive equilibrium E 4 disappears. Particularly, when the curve SN meets the curve H at the BT point, model (3) may undergo a BT bifurcation of codimension 2.
A small area at the origin in the ( 1 , 2 ) parameter plane is divided into four regions by the above three local bifurcation curves, as is shown in Figure 7. When ( 1 , 2 ) = (0, 0), Figure 8(a) suggests that E * is a cusp with codimension 2, and it is the unique positive equilibrium in model (3). In region I of Figure 7, there exists three boundary equilibria, but no positive equilibrium (see Figure 8(b)). There exists one positive saddle-node in model (3) when ( 1 , 2 ) is on the saddle-node bifurcation curve SN. Passing through the curve SN from region I to region II, an unstable saddle and an unstable node appear via saddle-node bifurcation (see Figure 8(c)). Furthermore, on the Hopf bifurcation curve H, a limit cycle emerges via Hopf bifurcation. The limit cycle bifurcating from the positive equilibrium is unstable, since there occurs a subcritical Hopf bifurcation (see Figure 8(d)). When ( 1 , 2 ) reaches the homoclinic bifurcation curve HL, the limit cycle becomes an unstable homoclinic one shown in Figure 8(e). Finally, the unstable homoclinic cycle disappears in region IV. From Figure 8(f), we can see clearly that a saddle and a stable focus can coexist near the curve HL. The bifurcation diagram of model (3) according to the parameters γ and δ is shown in Figure 9(a). When γ and δ are small, model (3) only has boundary equilibria (see region I in Figure 9(a)), and the equilibrium (0, 0) is a stable node. For a small γ, as δ increases from region I to region II, two positive equilibria coexist (see Figure 9(a)). For a large δ, with the increase of γ from region II to region III, two positive equilibria coincide and become one equilibrium via saddle-node bifurcation (see Figure 9(a)). Due to transcritical bifurcation, the unique positive equilibrium coincides with the boundary equilibrium. In region IV, there are only three boundary equilibria, and (K, 0) is a stable node. In addition, Hopf bifurcation may occur on the curve H. Particularly, a Bogdanov-Takens bifurcation of codimension 2 takes place at the BT point in model (3), once the curves SN and H meet.
For γ = δ, when γ is small, (0, 0) is a stable node (see the red dashed line in Figure 9(a)). When γ is large, (0, 0) and (K, 0) are the only two boundary equilibria, and both of them are stable node (see the red dotted line in Figure 9(a)). However, when γ is at the middle region, there exists olny one positive equilibrium in the model (see the red solid line in Figure 9(a)). Unlike γ = δ, the dynamics of predator-prey model with γ < δ are complicated when γ is small (see Figure 9(a)). We also trace the critical value of m for saddle-node bifurcation in Figure 9(b). From Figure 9(b), we can see that for γ × δ = (0, 0.058) × (0, 1), there is always a m ∈ (0, 1) satisfying the occurrence conditions of saddle-node bifurcation.

Conclusions
A predator-prey model with strong Allee effect and nonconstant mortality rate is studied in this paper. The aim is to explore what dynamic behaviors will occurs in the predator-prey model with a strong Allee effect caused by the nonconstant mortality. Our results indicate that there always exists three boundary equilibria in model (3), where (0, 0) is always a stable node, (K, 0) is stable only if γ > eb(1 − θ)K and (m, 0) is always unstable whenever exists. As for positive equilibria, when γ < eb(1 − θ)m and m < x A < K, there is a unique positive equilibrium when H(x A ) = 0 in model (3), and the occurrence of saddle-node bifurcation brings two positive equilibria when H(x A ) > 0. A limit cycle appears with the positive equilibrium undergoing a Hopf bifurcation at γ = γ H , whose stability is determined by the first Lyapunov number l 1 . Additionally, the positive equilibrium coincides with (K, 0) at γ = eb(1 − θ)K via transcritical bifurcation.
When the curve H meets the curve SN, a BT bifurcation of codimension 2 occurs. The BT bifurcation point is sensitive to the perturbation intensity. The predator-prey model can generate rich dynamics when we introduce two small perturbation into the BT bifurcation parameters γ BT and m BT . Two positive equilibria exist in model (3), one is a saddle, and the other undergoes a Hopf bifurcation then produces limit cycles. A homoclinic cycle emerges bifurcating from the limit cycle via homoclinic bifurcation. While for a large perturbation to γ BT and m BT , there are only boundary equilibria.
For γ = δ in model (3), the predator has a constant mortality rate. There is a γ * such that model (3) only has boundary equilibria for 0 < γ < γ * , and (0, 0) is a stable node. However, when γ < δ, the mortality rate of predator is nonlinear. There is a δ * such that two positive equilibria exist in model (3) and the model undergoes a saddle-node bifurcation for 0 < γ < γ * and δ > δ * . Moreover, when 0 < γ < γ * , a m satisfying the occurrence conditions of saddle-node bifurcation always exists. These results suggest that the predator dies out (no positive equilibrium) under constant mortality rate, but nonconstant mortality can result in the coexistence of predator and prey (positive equilibrium exists). Obviously, nonconstant mortality rate is the cause of the emergence of two positive equilibria. Therefore, nonconstant mortality plays an important role in the dynamics of predator-prey model, and we cannot simply approximate nonconstant mortality rate to a constant, especially for small γ.
From the ecological view point, when the mortality γ at low density and the limiting, maximal mortality rate δ are small, both the prey and the predator will go extinct due to Allee effect. Interestingly, the prey and the predator will reach a coexistence state when the mortality γ is small enough and the limiting, maximal death rate δ is large enough. The reason is that the increase in the mortality rate of predator relieves the feeding pressure on the prey and reduces the competition among predators for prey. However, when γ is large, (K, 0) is a asymptotically stable boundary equilibrium. This means that the predator dies out and the prey increases quickly to reach its carrying capacity K when the mortality rate of predator is too high. Therefore, we can approximate the nonconstant mortality rate to a constant value when γ is large. Moreover, prey refuge plays a negligible role in the persistence and extinction of the prey and the predator, but it influences the population density greatly. It is worth noting that the diffusion can also affect the dynamics of predator-prey model with strong Allee effect significantly [42,43]. A detailed analysis on the effect of diffusion is beyond the scope of this work and will be presented in the future.
In short, the nonconstant mortality of predator do have a significant impact on the dynamics of predator-prey model. We cannot simply approximate the nonconstant death rate to a constant unless the mortality rate is large. We expect that our results provide a new insight into the predator-prey model with nonconstant mortality rate. and R 2 (z 1 , z 2 , ) is at least of the third order with terms z i 1 z j 2 , whose coefficients are smooth functions of 1 and 2 .
Next we let τ to replace the time variable such that dt = (1 − η 02 ( )z 1 )dτ. Model (A4) can be written as model (A5) when we rewritten t to denote τ.