Dynamics of a predator-prey model with strong Allee effect and nonconstant mortality rate

1.   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) - \alpha(P)$, then we get the following equation:

where $e$ and $\alpha(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., $\alpha(P) = u$ ^[18,19]. However, Cavani and Farkas ^[20] introduced a nonconstant mortality rate of predator $\alpha(P) = \frac{\gamma + \delta p}{1 + P}$, where $\gamma > 0$ and $\delta > 0$ represent mortality rate at the low density and the maximal mortality, respectively($\gamma \leq \delta$). The nonconstant mortality $\alpha(P)$ is a bounded and increasing function of $P$. When $\gamma = \delta$, 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$\acute{a}$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, $\theta$ is the refuge protecting of the prey ($0 < \theta < 1$), $e$ and $b$ respectively represent the conversion efficiency of predator by consuming prey and the predation rate, $\gamma$ and $\delta$ 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.

2.   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

According to model (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 = \frac{a}{Kb(1-\theta)}(K-x)(x-m)$ and $H(x) = 0$. Here,

If $\Delta > 0$, then the root of $\frac{{\rm d}H(x)}{{\rm d}x} = 0$ is

where

By direct calculation, we have

For $x = g(y)$, $\left(\frac{\gamma}{eb(1 - \theta)}, 0 \right)$ 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:

(i) $\gamma < eb(1 - \theta)m$;

(ii) $\gamma = eb(1 - \theta)m$;

(iii) $eb(1 - \theta)m < \gamma < eb(1 - \theta)K$;

(iv) $\gamma \ge eb(1 - \theta)K$.

When $\gamma \ge eb(1 - \theta)K$, the model has no positive equilibrium. That is to say, if there exist positive equilibria, then $\gamma < eb(1 - \theta)K$, i.e., $H(K) < 0$.

In conclusion, when $m < {x_A} < K$, $H({x_A}) = 0$ and $\gamma < eb(1 - \theta)m$, the model has a unique positive equilibrium ${E_*} = ({x_*}, {y_*})$, where $x_{*} = x_{A}$ and $y_{*} = \frac{a(K-x_{*})(x_{*}-m)}{Kb(1-\theta)}$; when $m < {x_A} < K$, $H({x_A}) > 0$ and $\gamma < eb(1 - \theta)m$, ${E_3} = ({x_3}, {y_3})$ and ${E_4} = ({x_4}, {y_4})$ are two positive equilibria of the model, where $H({x_3}) = 0$, $H({x_4}) = 0$ and $y_{3, 4} = \frac{a(K-x_{3, 4})(x_{3, 4}-m)}{Kb(1-\theta)}$; when $m < {x_A} < K$, $\gamma = eb(1 - \theta)m$ or $eb(1 - \theta)m < \gamma < eb(1 - \theta)K$, the model has a unique positive equilibrium ${E_4} = ({x_4}, {y_4})$.

About the stability of all the equilibria, we get the following theorems.

Theorem 2.1. (i) $E_{0} = (0, 0)$ is always a stable node.

(ii) $E_{1}$ is a saddle if $\gamma > eb(1 - \theta)m$; $E_{1}$ is an unstable node if $\gamma < eb(1 - \theta)m$; $E_{1}$ is a repelling saddle-node if $\gamma = eb(1 - \theta)m$.

(iii) $E_{2}$ is a stable node if $\gamma > eb(1 - \theta)K$; $E_{2}$ is a saddle if $\gamma < eb(1 - \theta)K$; $E_{2}$ is an attracting saddle-node if $\gamma = eb(1 - \theta)K$.

Proof. The Jacobian matrix of model (3) is

here, $J_{11} = b(1-\theta)(f(x)-y)+b(1 - \theta)xf'(x)$, $J_{12} = -b(1 - \theta)x$, $J_{21} = eb(1 - \theta)y$, $J_{22} = eb(1-\theta)(x-g(y))- eb(1 - \theta)yg'(y)$.

By direct calculation, we can get

(i) The eigenvalues of $J$ at $E_{0} = (0, 0)$ are ${\lambda_1} = -ma < 0$ and ${\lambda_2} = -\gamma < 0$. Therefore, $E_{0} = (0, 0)$ is a stable node.

(ii) The eigenvalues of $J$ at $E_{1} = (m, 0)$ are ${\lambda_1} = ma(1 - \frac{m}{K}) > 0$ and ${\lambda_2} = eb(1 - \theta)m - \gamma$. If $\gamma > eb(1 - \theta)m$, then ${\lambda_2} < 0$, $E_{1}$ is a saddle; if $\gamma < eb(1 - \theta)m$, then ${\lambda_2} > 0$, $E_{1}$ is an unstable node. For the case $\gamma = eb(1 - \theta)m$, we need to discuss further.

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^\infty}$ functions at least of third order with respect to $(X, Y)$ and

Let $X = x - \frac{a_{01}}{a_{10}}y$, $Y = y$ and $\tau = a_{10}t$, the model becomes

By Theorem 7.1 in ^[41], since the coefficient of ${y^2}$ is $\frac{{{b_{02}}}}{{{a_{10}}}} - \frac{{{a_{01}}{b_{11}}}}{{a_{10}^2}} = \frac{{\gamma b(1 - \theta) - ma(1 - \frac{m}{K})(\delta - \gamma)}}{{{m^2}{a^2}{{(1 - \frac{m}{K})}^2}}} \ne 0$ and ${a_{10}} = ma(1 - \frac{m}{K}) > 0$, the equilibrium $E_{1}$ is a repelling saddle-node.

(iii) At $E_{2} = (K, 0)$, the eigenvalues of $J$ are ${\lambda_1} = -a(K-m) < 0$ and ${\lambda_2} = eb(1 - \theta)K - \gamma$. When $\gamma > eb(1 - \theta)K$, $E_{2}$ is a stable node due to ${\lambda _2} < 0$. When $\gamma < eb(1 - \theta)K$, $E_{2}$ is a saddle due to ${\lambda_2} > 0$.

In the following, we will discuss the case $\gamma = eb(1 - \theta)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^\infty }$ functions at least of third order with respect to $(X, Y)$, and

Let $X = x - \frac{a'_{01}}{a'_{10}}y$, $Y = y$ and $\tau = a'_{10}t$, the model becomes

The coefficient of ${y^2}$ is $\frac{{{b'_{02}}}}{{{a'_{10}}}} - \frac{{{a'_{01}}{b'_{11}}}}{{a'_{10}}^2} = \frac{{\delta - \gamma }}{{a(K - m)}} + \frac{{Ke{b^2}{{(1 - \theta)}^2}}}{{{a^2}{{(K - m)}^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.

Theorem 2.2. If the condition $m < {x_A} < K, H({x_A}) = 0$ and $\gamma < eb(1-\theta)m$ holds, the unique positive equilibrium ${E_*} = ({x_*}, {y_*})$ exists in model (3).

(i) When ${c_{10}} + {d_{01}} = 0$, ${E_*}$ is a degenerate critical point (cusp). In addition, if ${f_{20}}({f_{11}} + 2{e_{20}}) \ne 0$, then $E_{*}$ is a cusp of codimension 2; else $E_{*}$ is a cusp of codimension at least 3.

(ii) When $a_{10}+b_{01} \neq 0$ and $h_{20}\neq 0$, if $f'_{01} > 0$, then $E_{*}$ is a repelling saddle-node; else $E_{*}$ is an attracting saddle-node.

Proof. According to the Jacobian matrix of model (3), we can get

Obviously, $\det(J(E_{*})) = 0$.

Letting $X = x-x_{*}$, $Y = y-y_{*}$ to transform $E_{*}$ to $(0, 0)$ for studying the stability of $E_{*}$. We can rewrite model (3) as

where ${R_1}(X, Y)$ and ${S_1}(X, Y)$ are ${C^\infty}$ functions at least of third order with respect to $(X, Y)$ and

At $E_{*} = (x_{*}, y_{*})$, the Jacobian matrix of model (3) is

where $\det(J(E_{*})) = c_{10}d_{01} - c_{01}d_{10} = 0$. Then the eigenvalues of $J(E_{*})$ are $\lambda_{1} = 0$, $\lambda_{2} = c_{10} + d_{01}$.

Case 1: $c_{10} + d_{01} = 0.$

Obviously, the eigenvalues of $J(E_{*})$ are $\lambda_{1} = \lambda_{2} = 0$. Through the transformation of $\left({\begin{array}{*{20}{c}} X\\ Y \end{array}} \right) = \left({\begin{array}{*{20}{c}} {{c_{01}}} & 0\\ { - {c_{10}}} & 1 \end{array}} \right)\left({\begin{array}{*{20}{c}} x\\ y \end{array}} \right),$ model (9) becomes

By Lemma in ^[44], we can obtain the equivalent model of model(10) as follows:

If ${f_{20}}({f_{11}} + 2{e_{20}}) \ne 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} \neq 0$.

Obviously, the eigenvalues of $J(E_{*})$ are ${\lambda_1} = 0$ and ${\lambda_2} \ne 0$. Through the transformation of $\left({\begin{array}{*{20}{c}} X\\ Y \end{array}} \right) = \left({\begin{array}{*{20}{c}} {{d_{01}}} & {{c_{10}}}\\ { - {d_{10}}} & {{d_{10}}} \end{array}} \right)\left({\begin{array}{*{20}{c}} x\\ y \end{array}} \right),$ model (9) becomes

where $R{'_2}(x, y)$ and $S{'_2}(x, y)$ are ${C^\infty }$ functions at least of third order with respect to $(x, y)$ and

Redefine $\tau$ by $\tau = f{'_{01}}t$, we can get

where $R{'_3}(x, y)$ and $S{'_3}(x, y)$ are ${C^\infty}$ functions at least of third order with respect to $(x, y)$ and ${h_{ij}} = \frac{{e{'_{ij}}}}{{f{'_{01}}}}, {k_{ij}} = \frac{{f{'_{ij}}}}{{f{'_{01}}}}$. The coefficient of ${x^2}$ is

If ${h_{20}} \ne 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.

Theorem 2.3. When $m < {x_A} < K$, $H({x_A}) > 0$ and $\gamma < eb(1-\theta)m$, model (3) has two positive equilibria ${E_3} = ({x_3}, {y_3})$ and ${E_4} = ({x_4}, {y_4})$.

(i) $E_{3} = (x_3, y_3)$ is always a saddle.

(ii) if ${x_4} \in [\frac{{m + K}}{2}, K)$, then ${E_4} = ({x_4}, {y_4})$ is local asymptotically stable; if ${x_4} \in (m, \frac{{m + K}}{2})$, 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.

(i) If ${x_3} \in (m, \frac{{m + K}}{2})$, then $\det(J(E_3)) = e{b^2}{(1 - \theta)^2}{x_3}{y_3} \left(1 - f'({x_3})g'({y_3}) \right) < 0$. Therefore, $E_{3} = (x_3, y_3)$ is a saddle.

(ii) If ${x_4} \in [\frac{{m + K}}{2}, K)$, then $\det(J(E_4)) > 0$ and $tr(J(E_4)) = b(1 - \theta){x_4}f'({x_4}) - eb(1 - \theta){y_4}g'({y_4}) < 0$ due to $f'({x_4}) \le 0, g'({y_4}) > 0$. Thus, ${E_4} = ({x_4}, {y_4})$ is local asymptotically stable. If ${x_4} \in (m, \frac{{m + K}}{2})$, then $\det (J(E_4)) = e{b^2}{(1 - \theta)^2}{x_4}{y_4}(1 - f'({x_4})g'({y_4})) > 0$ since $\frac{1}{{g'({x_4})}} > f'({x_4})$. 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.

3.   Bifurcation analysis

We are interested in various possible bifurcations of model (3) in this section, including transcritical, saddle-node, Hopf and Bogdanov-Takens bifurcation.

3.1.   Transcritical bifurcation

According to Theorem 2.1, the boundary equilibrium $E_{2}$ is unstable when $\gamma < eb(1-\theta)K$, while it is stable when $\gamma > eb(1-\theta)K$. Particularly, when $\gamma = eb(1-\theta)K$, $E_{2}$ coincides with the positive equilibrium $E_{4}$. Therefore, we consider the existence of transcritical bifurcation around $E_{2}$ by taking $\gamma$ as bifurcation parameter.

Theorem 3.1. A transcritical bifurcation occurs in model (3) when $\gamma = eb(1 - \theta)K \triangleq \gamma_{\rm {TC}}$.

Proof. We have $\det(J(E_2)) = 0$ when $\gamma = eb(1-\theta)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 $\gamma = eb(1-\theta)K$ according to Sotomayor's theorem.

This is the end of the proof.

3.2.   Saddle-node bifurcation

From the previous analysis, we can see clearly that under the condition $m < {x_A} < K$ and $\gamma < eb(1-\theta)$, model (3) has two positive equilibria ${E_3} = ({x_3}, {y_3})$ and ${E_4} = ({x_4}, {y_4})$ if $H({x_A}) > 0$. The two positive equilibria coincide if $H({x_A}) = 0$, while they disappear if $H({x_A}) < 0$. Therefore, a saddle-node bifurcation will occur at $E_{*} = (x_{*}, y_{*})$. Saddle-node bifurcation at $E_{*}$ can be summarized into the following theorem.

Theorem 3.2. A saddle-node bifurcation occurs in model (3) when $\gamma = {\gamma_1} \triangleq {\gamma_{\rm {SN}}}$, where ${\gamma_1}$ satisfies $H({x_A}) = 0$.

Proof. According to Theorem 2.2, there is a zero eigenvalue $\lambda_{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 $\gamma = \gamma_{\rm {SN}}$ according to Sotomayor's theorem.

This is the end of the proof.

3.3.   Hopf bifurcation

The positive equilibrium $E_{4}$ may be a focus from Theorem 2.3. In this subsection, the conditions for Hopf bifurcation at $E_{4}$ will be discussed. Obviously, $\det (J{({E_4})}{|_{\gamma = \gamma _{\rm {H}}}}) > 0$, and there exists $\gamma \triangleq \gamma_{\rm {H}}$ such that $tr(J(E_4){|_{\gamma = \gamma_{\rm {H}}}}) = 0$.

Theorem 3.3. A Hopf bifurcation occurs at $E_{4}$ in model (3) when $\gamma = \gamma_{\rm {H}}$.

Proof. According to $tr(J(E_4))$, we have

which satisfies the transversality condition for Hopf bifurcation.

For the stability of the limit cycle, we need to calculate the first Lyapunov number $l_{1}$ at $E_{4}$ first.

Letting $X = x-x_{4}$, $Y = y-y_{4}$ to transform $E_{4}$ to $(0, 0)$ and we rewrite model (3) as

where

where $P_{4}(X, Y)$ and $Q_{4}(X, Y)$ are the power series in $(X, Y)$ with the term $X^iY^j$ satisfying $i+j \geq 4.$

For the first Lyapunov number $l_{1}$, which is given by the formula

where $\Delta = \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$.

3.4.   Bogdanov-Takens bifurcation

According to Theorem 3.2, positive equilibrium appears via saddle-node bifurcation. The stability of $E_{4}$ may change by Hopf bifurcation in Theorem 3.3. When saddle-node bifurcation crosses Hopf bifurcation, a Bogdanov-Takens (BT) bifurcation will occur under a small parameter perturbation in model (3). By choosing $\gamma$ and $m$ as two bifurcation parameters, $\gamma_{\rm {BT}}$ and $m_{\rm {BT}}$ are the threshold value, i.e., $\det (J(E_*){|_{(\gamma, m) = (\gamma_{\rm {BT}}, m_{\rm {BT}})}}) = 0$ and $tr(J(E_*){|_{(\gamma, m) = (\gamma_{\rm {BT}}, m_{\rm {BT}})}}) = 0$).

Theorem 3.4. A Bogdanov-Takens bifurcation of codimension 2 will occur in a small neighborhood of $E_{*}$ when $(m, \gamma)$ varies near $(m_{\rm {BT}}, \gamma_{\rm {BT}})$ with $a_{10}+b_{01} = 0$ and $g_{20}g_{11}\neq 0$.

For the proof of the theorem, please see Appendix 5.

4.   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$, $\theta = 0.1$, $\delta = 0.9$, $K = 0.9$, $m = 0.2$. The parameter $\gamma$ 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 $\gamma$, there exist equilibria $E_{0}$, $E_{1}$ and $E_{2}$. As $\gamma$ 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 $\gamma$, the positive equilibria $E_{3}$ and $E_{4}$ appear as the taking place of a saddle-node 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 $\gamma$ increases and crosses the line $L_{2}$ ($\gamma_{\rm {H}} = 0.05400367556$), a Hopf bifurcation occurs. The first Lyapunov number $l_{1} = 6.414584574\pi$ 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 ($\gamma < \gamma_{\rm {H}}$, see Figure 3(c)), and becomes a stable focus in region III ($\gamma > \gamma_{\rm {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 $\gamma$ 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 $\gamma$ 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 $\gamma$ and $m$ as the control parameters to numerically study BT bifurcation. Direct numerical computation shows that BT bifurcation may occur when $m = {m_{\rm {BT}}} = 0.1347130155$ and $\gamma = {\gamma_{\rm {BT}}} = 0.003455490134$. The equilibrium $E_{*}$ is a cusp when $m = {m_{\rm {BT}}}$ and $\gamma = {\gamma_{\rm {BT}}}$, which is dispalyed in Figure 5(a). For $m = 0.4$ and $\gamma = 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 $\gamma = 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 $\gamma-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.

For $m = m_{\rm {BT}}$ and $\gamma = \gamma_{\rm {BT}}$, we can obtain $\mathop {\left| {\frac{{\partial (\mathop \zeta \nolimits_1, \mathop \zeta \nolimits_2)}}{{\partial (\mathop \epsilon \nolimits_1, \mathop \epsilon \nolimits_2)}}} \right|}\nolimits_{\mathop \epsilon \nolimits_1 = \mathop \epsilon \nolimits_2 = 0} \ne 0$. Thus, the parametric transformation $\zeta_{1}(\epsilon) = \gamma_{00}(\epsilon)\gamma^4_{11}(\epsilon)$, $\zeta_{2}(\epsilon) = \gamma_{01}(\epsilon)\gamma_{11}(\epsilon)$ is nonsingular. For small $\epsilon_{1}$ and $\epsilon_{2}$, $\theta_{20} > 0$ and $\theta_{11} \neq 0$, which means that BT bifurcation with codimension 2 occurs in model (3). Under the second-order approximation, the bifurcation curves can be locally expressed as:

1) SN: the local expression of the saddle-node bifurcation curve

2) H: the local expression of the Hopf bifurcation curve

3) HL: the local expression of the homoclinic bifurcation curve

A small area at the origin in the $(\epsilon_{1}, \epsilon_{2})$ parameter plane is divided into four regions by the above three local bifurcation curves, as is shown in Figure 7. When $(\epsilon_{1}, \epsilon_{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 $(\epsilon_{1}, \epsilon_{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 $(\epsilon_{1}, \epsilon_{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 $\gamma$ and $\delta$ is shown in Figure 9(a). When $\gamma$ and $\delta$ 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 $\gamma$, as $\delta$ increases from region I to region II, two positive equilibria coexist (see Figure 9(a)). For a large $\delta$, with the increase of $\gamma$ 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 $\gamma = \delta$, when $\gamma$ is small, $(0, 0)$ is a stable node (see the red dashed line in Figure 9(a)). When $\gamma$ 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 $\gamma$ is at the middle region, there exists olny one positive equilibrium in the model (see the red solid line in Figure 9(a)). Unlike $\gamma = \delta$, the dynamics of predator-prey model with $\gamma < \delta$ are complicated when $\gamma$ 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 $\gamma \times \delta = (0, 0.058) \times (0, 1)$, there is always a $m \in (0, 1)$ satisfying the occurrence conditions of saddle-node bifurcation.

5.   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 $\gamma > eb(1 - \theta)K$ and $(m, 0)$ is always unstable whenever exists. As for positive equilibria, when $\gamma < eb(1 - \theta)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 $\gamma = \gamma_{\rm {H}}$, whose stability is determined by the first Lyapunov number $l_1$. Additionally, the positive equilibrium coincides with $(K, 0)$ at $\gamma = eb(1 - \theta)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 $\gamma_{\rm {BT}}$ and $m_{\rm {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 $\gamma_{\rm {BT}}$ and $m_{\rm {BT}}$, there are only boundary equilibria.

For $\gamma = \delta$ in model (3), the predator has a constant mortality rate. There is a $\gamma_{*}$ such that model (3) only has boundary equilibria for $0 < \gamma < \gamma_{*}$, and $(0, 0)$ is a stable node. However, when $\gamma < \delta$, the mortality rate of predator is nonlinear. There is a $\delta_{*}$ such that two positive equilibria exist in model (3) and the model undergoes a saddle-node bifurcation for $0 < \gamma < \gamma_{*}$ and $\delta > \delta_{*}$. Moreover, when $0 < \gamma < \gamma_{*}$, 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 $\gamma$.

From the ecological view point, when the mortality $\gamma$ at low density and the limiting, maximal mortality rate $\delta$ 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 $\gamma$ is small enough and the limiting, maximal death rate $\delta$ 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 $\gamma$ 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 $\gamma$ 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.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grants No. 61901303 and No. 61871293), the National Key Research and Development Program of China (Grant No. 2018YFE0103700), and Wenzhou basic agricultural science and technology project (Grant No. S20180006).

Conflict of interest

The author declares that there are no conflicts of interest.

Appendix

The proof of Theorem 3.4

Proof. To obtain the local critical expressions of the bifurcation curves around the BT point, firstly we need to transform model (3) into a normal form of BT bifurcation.

Introducing two small perturbations to the parameters $\gamma$ and $m$ at the critical BT bifurcation value through $m = {m_{\rm {BT}}} + \epsilon_{1}$ and $\gamma = {\gamma _{\rm {BT}}} + \epsilon_{2}$, where $\epsilon \equiv(\epsilon_{1},\epsilon_{2})$ is extremely close to zero vector. Under this operation, model (3) can be rewritten as

Then we can introduce the transformations $n_{1} = x-x_{*}$ and $n_{2} = y-y_{*}$ to let $(0,0)$ be the bifurcation point, it follows that

where

and $T_{1}(n_{1},n_{2},\epsilon)$ and $R_{1}(n_{1},n_{2},\epsilon)$ are at least of the third order with terms $n^i_{1}n^j_{2}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Next, we introduce the following $C^\infty$ variable substitutions near the origin

Model (A2) becomes

where

and $R_{2}(z_{1},z_{2},\epsilon)$ is at least of the third order with terms $z^i_{1}z^j_{2}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Next we let $\tau$ to replace the time variable such that $\rm{d}t = (1-\eta_{02}(\epsilon)z_{1})\rm{d}\tau$. Model (A4) can be written as model (A5) when we rewritten $t$ to denote $\tau$.

Let $\upsilon_{1} = z_{1},\upsilon_{2} = z_{2}(1-\eta_{02}(\epsilon)z_{1})$, then model (A5) can be expressed as

where

and $R_{3}(v_{1},v_{2},\epsilon)$ is at least of the third order with terms $\upsilon^i_{1} \upsilon^j_{2}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

The express of $\theta_{20}(\epsilon)$ is complex enough, we omit it. And we can not get any information about the sign of $\theta_{20}(\epsilon)$. Therefore, the following two situations need to be considered to ensure the following process is meaningful.

Case 1: $\theta_{20}(\epsilon)>0$ for small $\epsilon_{i}(i = 1,2)$, making the following variable substitutions

Retaining $\tau$ as $t$, model (A1) can be expressed as

where

and $R_{4}(\mu_{1},\mu_{2},\epsilon)$ is at least of the third order with terms $\mu^i_{1} \mu^j_{2}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Model (A7) can be expressed as model (A8) under the changes of $\omega_{1} = \mu_{1}+\frac{s_{10}(\epsilon)}{2}, \omega_{2} = \mu_{2}.$

where

and $R_{5}(\omega_{1},\omega_{2},\epsilon)$ is at least of the third order with terms $\omega_{1}^{i}\omega_{2}^{j}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Assuming $\theta_{11}(\epsilon)\neq0$, then $r_{11}(\epsilon) = s_{11}(\epsilon) = \frac{\theta_{11}(\epsilon)}{\sqrt{\theta_{20}(\epsilon)}}\neq0$. Making the following variable substitutions

then a versal unfolding form of model (A1) can be written as

where

and $R_{6}(x,y,\epsilon)$ is at least of the third order with terms $x^iy^j$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Case 2: $\theta_{20}(\epsilon)<0$ for small $\epsilon_{i}(i = 1,2)$, making the following variable substitutions

Retaining $\tau$ as $t$, model (A1) can be expressed as

where

and $R'_{4}(\mu'_{1},\mu'_{2},\epsilon)$ is at least of the third order with terms $\mu'^i_{1} \mu'^j_{2}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Model (A10) can be expressed as model (A11) under the changes of $\omega'_{1} = \mu'_{1}-\frac{s'_{10}(\epsilon)}{2}, \omega'_{2} = \mu'_{2}.$

where

and $R'_{5}(\omega'_{1},\omega'_{2},\epsilon)$ is at least of the third order with terms $\omega'^i_{1}\omega'^j_{2}$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Assuming $\theta_{11}(\epsilon)\neq0$, then $r'_{11}(\epsilon) = s'_{11}(\epsilon) = \frac{\theta_{11}(\epsilon)}{\sqrt{-\theta_{20}(\epsilon)}}\neq0$. Making the following variable substitutions

Then a versal unfolding form of model (A1) can be written as

where $\zeta'_{1}(\epsilon) = -\gamma'_{00}(\epsilon)\gamma'^4_{11}(\epsilon), \zeta'_{2}(\epsilon) = -\gamma'_{01}(\epsilon)\gamma'_{11}(\epsilon).$ $R'_{6}(x',y',\epsilon)$ is at least of the third order with terms $x'^iy'^j$, whose coefficients are smooth functions of $\epsilon_{1}$ and $\epsilon_{2}$.

Retain $\zeta_{1}(\epsilon)$ and $\zeta_{2}(\epsilon)$ to denote $\zeta'_{1}(\epsilon)$ and $\zeta'_{1}(\epsilon)$ to combine the above two situations. In a small area near the origin, the transformations of $\zeta_{1}(\epsilon)$ and $\zeta_{2}(\epsilon)$ are homeomorphisms when $\mathop {\left| {\frac{{\partial (\mathop \zeta \nolimits_1 ,\mathop \zeta \nolimits_2 )}}{{\partial (\mathop \epsilon \nolimits_1 ,\mathop \epsilon \nolimits_2 )}}} \right|}\nolimits_{\mathop \epsilon \nolimits_1 {\rm{ = }}\mathop \epsilon \nolimits_2 {\rm{ = }}0}$ is a nonsingular matrix, in the mean while $\zeta_{1}(\epsilon)$ and $\zeta_{2}(\epsilon)$ are independent parameters. Then we know there takes place a BT bifurcation when $\epsilon = (\epsilon_{1},\epsilon_{2})$ around the origin. Under the second-order approximation, the bifurcation curves can be locally expressed as (``+" for $\theta_{20}(\epsilon)>0$, ``-" for $\theta_{20}(\epsilon)<0$):

1) SN: the local expression of the saddle-node bifurcation curve

2) H: the local expression of the Hopf bifurcation curve

3) HL: the local expression of the homoclinic bifurcation curve