Bifurcations of a prey-predator system with fear, refuge and additional food

1.   Introduction

In ecology, prey-predator interaction happens at many higher trophic levels, and due to these interactions, predators always have an impact; direct, indirect or both, on prey population. The direct effect can be seen in nature when the predators catch the prey and kill ^[1,2]. In resent years, more and more habitats and ecological communities have been modified by human activities, and this state will continue. Unlike in the past, the predator-prey interactions increasingly occur in other aspects ^[3]. The theories with respect to the evolution of predator-prey interactions are reviewed by Peter A. Abrams ^[4], which includes the dynamics and stability of both populations and traits, as well as the anti-predator traits to environmental changes. The crucial feature of predator-prey interactions is the predator's functional response, which is used to describe the predation of predator on prey ^[5,6]. They are generally classified into two types: the prey-dependent and predator-dependent, namely, the Holling type and Bedding type respectively ^[7,8,9]. Whereas, the Holling-Ⅱ type functional response is most commonly used in mathematical modeling ^[10,11].

Apart from the direct killing of prey by predator, more and more experimental data show that population of prey reduced due to fear of the predator ^{[12,13,14,15,16,17]}. Due to the fear effect, the prey shows a variety of antipredator responses, including different psychological changes, habitat changes, and foraging. For example, the impact of fear has been investigated on the prey's birth rate, Zanette et al. ^[15] carried out an experiment on song sparrows during the whole breeding season in 2011, the experiment was conducted in a free-living area for birds and there were no predators, the experimenter have used playbacks of the predator's voice to manipulate artificial risk, they found a reduction in reproduction by $40\%$ in the number of offsprings due to the fear of predation alone, the experiment showed that prey could respond to the danger of being preyed. Similarly, another experiment revealed the effect of predator's fear for the death rate of prey, the experiment was performed by on the prey (grasshopper nymphs) and the predator (toothless spider) ^[17], the results of the experiment showed that the presence of predator increased the number of the death of prey by $20\%$ compared with no predators, if the number of killing by predators was taken into account, then the death rate of prey was enhanced only by $9\%$, in fact, the indirect predation rate was 2 times larger than the direct killing. Further evidence has been found about Siberian jay in Eggers et al. ^[18], birds in Hua et al.^[19] and lambs in Krapivsky et al. ^[20]. In addition, the recent research shows that the fear effect can induce bifurcations and even stabilize chaos of predator-prey system dynamics^[21,22,23]. The adult prey also shows weaker anti-predator behaviors if the cost of fear is higher, such as reducing outside activity, foraging and starving to prevent from predation. These anti-predation sensitivity of prey should be considered into the mathematical model to reveal the realistic system dynamics ^[24]. Almost all experiment evidences confirmed that the impact of fear is very important in studying the system dynamics, refer to Blumstein ^[25]. Therefore, the impact of fear and the anti-predation of prey should be added to the ecological model, we can reveal how these factors affect the system stability.

In the anti-predation behaviors of prey, the prey usually move to safer place to avoid attack from predator, so the refuge area is very important for prey during their survival ^{[26,27,28,29,30]}. Most of the researchers have shown that refuges can preserve the prey and stabilize the system, which is very important for the survival of prey and predator. Kar ^[31] proposed a prey-predator model with Holling-Ⅱ type functional response, and he added shelter for prey in the model, Kar established some sufficient conditions assuring the existence and stability of the equilibrium points, he found that a stable limit cycles appeared when the equilibrium point was unstable by bifurcation analysis. Recently, the prey-predator models with refuge areas have attracted many researchers's attention ^[32].

Actually, due to the destruction of ecological environment, apart from the construction of refuge zones to preserve the survival of animals, human want to reduce the prey's negative effect of fear by releasing additional food to predator or prey or both of them, this is a very popular strategy in biological control ^{[33,34,35,36]}. Releasing proper additional food for predator can reduce the consumption of predator and conserve the prey, Moorland Working Group ^[34] made an experiment in the spring and summer of 1998 and 1999 years, they put in additional food to the hen harriers (predator) living on Langholm moor in south Scotland, the predation rate of prey reduced from 3.7 chicks to 0.5 chicks per hour. However, it is worth of pointing out that additional food supplement may increase the reproduction rate of predators causing the reduce or extinction of prey ^[35]. Certainly, the "quality" and "quantity" of additional food plays an key role in the system stability ^[36]. The researchers found that good additional food benefited the survival of predators which has led to arising high predation rate, low quality of food would do harm to the predator and benefit the prey ^[37,38], so the strategy of releasing additional food should be based on some theoretical results and be done scientifically and reasonably.

In this work we have proposed a prey-predator mathematical model with the fear effect and the anti-predation. In view of the importance of system stability, we have mainly studied the local and global stability of the equilibrium states and bifurcations. To observe the impact of fear, refuge and additional food on the dynamics of prey–predator system, we study the nonlinear dynamics of proposed mathematical model using various techniques of stability analysis, permanence, persistence and Hopf bifurcation analysis.

Our main works are as follows:

$1)$ The fear effect and their anti-predation sensitivity to predator have been incorporated in the prey-predator system.

$2)$ Biological control strategies such as constructing refuge area and releasing additional food for predator have been included.

$3)$ Some theorems of the local and global asymptotical stability and bifurcations have been obtained by using characteristic equation and Sotomayor's bifurcation theorems.

$4)$ The interaction of the effects of fear, refuge, additional food and the stability of equilibrium points have been investigated by using theoretical and numerical methods.

The paper has been organized as follows. The mathematical model is formulated in Section 2. The stability and bifurcation analysis are carried out in Section 3. The numerical simulations are executed in Section 4 to establish the bifurcation thresholds of some key parameters. Finally, we end the paper with a discussion of our findings.

2.   Model construction

Let $x(t)$ and $y(t)$ ($x$ and $y$ for simplicity) be the density of prey and predator at time $t$ respectively. In the absence of predator, we assume that the birth rate of prey is $a$, the competition rate between prey is $b$. The functional response of predator to prey is the Holling-Ⅱ type, where the meanings of $c$ and $d$ refer to ^[11,12]. The conversion rate of prey to predator is $\beta$, the dearth rate of predator is $e$, then we have the following prey-predator model:

Taking into account the fear of predator and the anti-predation level of prey to protect them from outside risk, the parameter $f$ is the fear level and the parameter $\kappa$ is the anti-predation sensitivity. Meanwhile, we incorporate parameter $\mu$ to describe the refuge area for prey and $\alpha$, $\beta$ and $F$ to describe the additional food supplement in biological control, where $\alpha F$ and $\beta$ represent the quantity (effectual food level) and quality of additional food respectively (the details refer to ^[33,36]), then we get the following model:

The initial conditions are $x(0) > 0$ and $y(0) > 0$. All parameters in (2.1) are supposed to be positive so as to accord with the biological justifications.

3.   Main results

3.1.   Existence of fixed point

By the theory of functional differential equation ^[40], we can obtain that there exists a unique positive solution for system (2.1) which remains in $R_+^2 = \{(x, y):x > 0, y > 0\}.$ It is easy to obtain that the trivial fixed point of system (2.1) is $Q_0(0, 0),$ the boundary fixed point is $Q_1(\frac{a}{b}, 0).$ Now we discuss the coexistence fixed point $\widetilde{Q}(\widetilde{x }, \widetilde{y })$, by the definition of coexistence fixed point, $\widetilde{x}$ satisfy the following equation,

Using restrictive condition $e d \kappa < \beta c < \frac{e(1+\delta\alpha F}{\alpha F}$, then

Denote $\frac{c(1-\mu)}{1+ d\kappa(1-\mu) \widetilde{x}+\delta\alpha F} = \Upsilon$ and put $\widetilde{x}$ into the first equation of Eq (2.1), we have

That is

Due to $b\widetilde{x}-a < 0$, it is clear that the Eq (3.1) has a unique positive solution, say $\widetilde{y}$. Consequently, (2.1) has a unique positive coexistence fixed point $\widetilde{Q}(\widetilde{x }, \widetilde{y }).$ Next we will establish some sufficient conditions assuring the stability of these fixed points.

3.2.   LAS and bifurcation analysis

In this part, we focus on the local asymptotic stability (LAS for simplicity) and bifurcation analysis of the fixed points. For (2.1), let $X = (x, y)^T,$ we rewrite (2.1) as the form $X' = G(X), G = (G_1, G_2)^T.$ For the convenience of the reader, we compute the partial derivatives of the operator $G$ as follows:

Next, we begin to discuss the stability of above fixed points one by one.

● Stability of $Q_0$

By the representation of $DG$, we can easily obtain the variational matrix at $Q_0$ is

By using the characteristic equation, we know that there is a positive characteristic root $\lambda = a > 0$, and hence it is unstable. If $\frac{b\alpha F}{1+\delta\alpha F}-e = 0$, then zero is a characteristic root, the eigenvectors of $\lambda = 0$ corresponding to the matrix $M_0$ and $M_0^T$ are $V_1 = (0, 0)^T$ and $V_2 = (0, 1)^T$ respectively. Take $\delta$ as a parameter and denote the solution of $\frac{b\alpha F}{1+\delta\alpha F}-e = 0$ as $\delta = \hat{\delta}.$ We investigate the conditions of the existence of bifurcations ^[39].

(ⅰ) $V_2^T[G_\delta(Q_0, \delta = \hat{\delta})] = 0$,

(ⅱ) $V_2^T[DG_\delta(Q_0, \delta = \hat{\delta})V_1] = (0, 1)\left[\left(\begin{array}{cc} 0 & 0\\ 0 & \frac{\beta c\alpha F}{1+\delta\alpha F} \\ \end{array} \right)\left(\begin{array}{c} 0 \\ 0 \\ \end{array} \right) \right] = 0,$

then there is neither transcritical bifurcation nor saddle-node bifurcation around the trivial fixed point $Q_0$, so we can get the following theorem.

Theorem 3.1. System (2.1) undergoes no bifurcation around the trivial fixed point $Q_0(0, 0)$, which is unstable.

● Stability of $Q_1$

We proceed to study the dynamics of the predator free equilibrium state $Q_1(\frac{a}{b}, 0)$. By same manner as before, we obtain the variational matrix at $Q_1(\frac{a}{b}, 0)$ is,

Denote $m_1 = e- \frac{\beta c[(1-\mu)a +b\alpha F]}{ad\kappa+b(1+\delta\alpha F)}$, then the characteristic equation corresponding to matrix $M_1$ is

The solutions of above equation are $\lambda_1 = -a, \lambda_2 = -m_1.$ If $m_1 > 0$, we obtain that the characteristic roots are all negative, and hence system (2.1) is LAS around $Q_1(\frac{a}{b}, 0)$. If $m_1 = 0$, there is a zero characteristic root, then eigenvectors of zero root w.r.t matrix $M_1$ and $M_1^T$ are $V_1 = (\nu, 1)^T$ and $V_2 = (0, 1)^T$ respectively, where $\nu = - \frac{a f\kappa}{b}-\frac{p(1-\mu)}{ad\kappa+b(1+\delta\alpha F)}.$ Take $\mu$ as a parameter, we verify the following transversality conditions of the existence of transcritical and saddle-node bifurcations respectively. We denote the $\mu$ satisfying $m_1 = 0$ as $\overline{\mu}.$ It is not difficult to compute that

By computation, we have

(ⅰ) $V_2^T[G_\mu(Q_1, \mu = \overline{\mu})]\neq0$.

(ⅱ) $V_2^T[DG_\mu(Q_1, \mu = \overline{\mu})V_1] = (0, 1)\left[\left(\begin{array}{cc} 0 & D_{\mu1}\\ D_{\mu2} & 0 \end{array} \right)\left(\begin{array}{c} \nu \\ 1 \end{array} \right) \right]\neq0,$

where $D_{\mu1} = \frac{c(1+\delta\alpha F)}{(1+ d\kappa(1-\mu) x+\delta\alpha F)^2}, D_{\mu2} = \frac{\beta c[\alpha F-(1+\delta\alpha F)](1-d\kappa(1-\mu) x+\delta\alpha F)}{(1+d\kappa(1-\mu) x+\delta\alpha F)^2}.$

(ⅲ) By the definition of $D^2G$, then

Accordingly,

Therefore, by Sotomayor's bifurcation theorem, there is a transcritical bifurcation around the predator-free fixed point $Q_1$. It is easy to get the representation of the threshold of parameter $\mu$ as,

Similarly, we can obtain the transcritical bifurcations and get their representations of thresholds w.r.t. parameters $\kappa, \delta$ and $\beta.$ We summarize our findings as the following theorem.

Theorem 3.2. System (2.1) is LAS around the predator-free fixed point $Q_1$ if $\frac{\beta c[(1-\mu)a +b\alpha F]}{ad\kappa+b(1+\delta\alpha F)} < e$; and if $\frac{\beta c[(1-\mu)a +b\alpha F]}{ad\kappa+b(1+\delta\alpha F)} = e$, then system (2.1) undergoes transcritical bifurcations w.r.t. parameter $\mu, \kappa, \delta$ and $\beta$ with the following parameter thresholds respectively.

● Stability of $\widetilde{Q}$

We begin with the analysis of the stability of system (2.1) near the coexistence equilibrium state $\widetilde{Q}(\widetilde{x }, \widetilde{y })$. Repeating the above computation leads to the following Jacobian matrix at $\widetilde{Q}(\widetilde{x }, \widetilde{y })$ as

where

The characteristic equation of the matrix $\widetilde{M}$ reads,

Obviously, if $m_{11} > 0, m_{21} > 0,$ then the characteristic roots are both negative and system (2.1) is LAS around the interior fixed point $\widetilde{M}$. If $m_{11} > 0, m_{21} = 0,$ then $\lambda_1 = 0, \lambda_2 = -m_{11}$. We continue to study whether the bifurcation appears around $\widetilde{M}$ or not. By computation, the eigenvectors of $\lambda = 0$ corresponding to $\widetilde{M}$ and $\widetilde{M}^T$ are $\widetilde{V}_1 = (\widetilde{\nu}, 1)^T, \widetilde{V}_2 = (0, 1)^T$ respectively, where $\widetilde{\nu} = -\frac{m_{11}}{m_{12}}.$ Take $\mu$ as a parameter and denote $\widetilde{\mu}$ such that $m_{21} = 0$. By same manner as before, it is not difficult to verify that

(ⅰ) $V_2^T[G_\mu(\widetilde{Q }, \mu = \widetilde{\mu})] = \frac{\beta c\widetilde{x}[\alpha F-(1+\delta\alpha F)]}{(1+d\kappa(1-\mu) \widetilde{x}+\delta\alpha F)^2}\neq0$.

(ⅱ) $V_2^T[D^2G (Q_1, \mu = \widetilde{\mu})(V_1 V_1)]\neq0.$

Therefore Sotomayor's theorem for saddle-node bifurcation ^[39] indicates that there is a saddle-node bifurcation around the coexistence fixed point $\widetilde{Q}$ with the threshold $\widetilde{\mu} = 1-\frac{d\kappa\alpha F}{1+\delta\alpha F}.$ Similarly, we can find that the saddle-node bifurcation around $\widetilde{Q}$ for parameters $\kappa, \delta$ and $\beta.$ We will make a summary in the conclusion part.

Theorem 3.3. System (2.1) is LAS around the coexistence fixed point $\widetilde{M}$ if the following $C_1$ and $C_2$ hold.

$C_1$: $2b\widetilde{x} + \frac{c(1-\mu)(1+\delta\alpha F)\widetilde{y} }{(1+d\kappa(1-\mu)\widetilde{ x}+\delta\alpha F)^2} - \frac{a }{1+f\kappa \widetilde{y} } > 0,$

and

$C_2$: $(1-\mu)(1+\delta\alpha F)-d\kappa\alpha F > 0$.

Further, if $C_1$ as well as $C_3$ hold,

$C_3$: $(1-\mu)(1+\delta\alpha F)-d\kappa\alpha F = 0$,

then system (2.1) undergoes saddle-node bifurcations w.r.t. parameter $\mu, \kappa$ and $\delta$ with the following thresholds respectively.

Next we investigate the existence of Hopf bifurcation. If $m_{11} = 0, m_{21} > 0$, then there are two imaginary roots $\lambda = \pm\sqrt{ m_{12}m_{21}}i$ for Eq (3.2). By the Hopf bifurcation theorem ^[39], we only need to verify the transversality conditions. Take $f$ as a bifurcation parameter and denote $f^H$ such that $m_{11} = 0$. Suppose the root of (3.2) has the general form as $\lambda = \phi_1(f)+\phi_2(f)i$. Substituting it into (3.2) leads to

Differentiating it w.r.t parameter $f$, then

That is

Consequently we have

where $X_1 = 2\phi_1+m_{11}, X_2 = 2\phi_2, X_3 = 2\phi_2-m_{11}'\phi_1-m_{12}'m_{21}-m_{12}m_{21}', X_4 = -\phi_2m_{11}'.$ By the definition of $f^H$ and the necessary condition of bifurcation appearance, at $f = f^H$, we have $\phi_1 = 0, m_{11} = 0$, and

That is, the transversality conditions of Hopf bifurcation hold, and hence we get the following result.

Theorem 3.4. System (2.1) undergoes a Hopf bifurcation around the coexistence equilibrium $\widetilde{M}(\widetilde{x}, \widetilde{y})$ provided the conditions $C_3$ and the following $C_4$ hold,

$C_4$: $2b\widetilde{x} + \frac{c(1-\mu)(1+\delta\alpha F)\widetilde{y} }{(1+d\kappa(1-\mu)\widetilde{ x}+\delta\alpha F)^2} - \frac{a }{1+f\kappa \widetilde{y} } = 0,$

and hence, the threshold is

Similarly, we can find the Hopf bifurcation thresholds about parameters $\mu, \delta$ and $\kappa$, which are to be validated in numerical analysis.

3.3.   GAS analysis

In this section, by using the Lyapunov functional stability theory ^[34], we investigate the global asymptotic stability (GAS) of the equilibrium points of system (2.1).

Theorem 3.5. Assume that system (2.1) is LAS near the predator-free fixed point $Q_1(\frac{a}{b}, 0)$, then it is GAS provided that

Proof. Consider the following function,

where $x_0 = \frac{a}{b}$, then it is positive on the interval $x > 0, y > 0.$ Computing its derivative w.r.t. $t$ along the solution of system (2.1), we have

Therefore, according to the stability theorem of functional differential equation ^[40], system (2.1) is GAS around the fixed point $Q_1$.

Theorem 3.6. Assume that system (2.1) is LAS near the coexistence fixed point $\widetilde{ Q}$, then it is GAS if the following conditions hold.

(ⅰ) $c (1-\mu)\frac{ \widetilde{y}}{\widetilde{x}} < b,$

(ⅱ) $\beta(1+\delta\alpha F-d\kappa\alpha F) < 1+\delta\alpha F+d\kappa(1-\mu)\widetilde{x}.$

Proof. Define

then it is positive on $x > 0, y > 0$. Similarly we have

Using the Lyapunov stability theory of differential equation again, we conclude that system (2.1) is GAS around $\widetilde{ Q}$.

4.   Numerical analysis

In this section, we study the bifurcations and stability of (2.1) from numerical angle intuitively by depicting some pictures.

4.1.   Bifurcation analysis

We begin with our discussion of bifurcation to reveal the impacts on the system stability and explore the dynamics of (2.1) in detail. We classify it into three different scenarios as the bifurcations of fear, refuge and additional food supplement. For system (2.1), we fix a set of parameter values as below

By computation, the predator extinction equilibrium point is $\overline{Q}(20, 0)$ and the interior equilibrium point is $\widetilde{Q}(0.7, 3.8547)$.

● Bifurcation of fear and anti-predation sensitivity

The parameter $f$ represents the fear level of predator on prey, meanwhile, the anti-predation sensitivity parameter $\kappa$ is induced by the fear of predator, which shows the sensitivity of prey that preventing themselves from the predation risk. These two indexes jointly reflect the effect of fear to the system dynamics. Consequently we analyze the bifurcations of parameters $f$ and $\kappa$ numerically. By use of Matlab2014a and continuation software Matcont ^[41], we get the trajectories of equilibrium point of (2.1), see Figure 1. It is clear that the prey is oscillatory if $f < 0.305109$, when $f = 0.305109$, the periodic fluctuation loses and it becomes stable when $f > 0.305109$ (see Figure 1(a)). For the parameter $\kappa < 0.19927$, the equilibrium point is stable, at $\kappa = 0.199270$, it becomes unstable and Hopf bifurcation occurs, while at $\kappa = 0.844747$, the oscillatory disappear and it becomes stable again from unstable state. A bubble formulates between the interval $\kappa\in(0.19927, 0.844747)$, which means that system (4.1) changes its stability from stable to unstable and to stable again, that is, multiple stable statuses appear, see also References [24]. Then as the increase of $\kappa$, the density of prey increases quickly, and at $\kappa = 0.9720$, a transcritical bifurcation occurs, which means the predator is extinct and the prey attains the peak value with $x = 20$ and keeps stable at the equilibrium point $\overline{Q}(20, 0)$.

To clearly see the interactive effect of the fear level $f$ and the anti-predation sensitivity level $\kappa$ on the stability of the equilibrium state, we carry our the two-parameter bifurcation diagrams. See Figure 2(a), where HC curve represents the Homoclinic bifurcation curve produced at the BT point (Bogdanov-Takens point for short) $(f, \kappa) = (-0.016576, 0.959873)$. The Hopf curve and the HC curve separate the whole plane into three regions denoted by $R_1, R_2$ and $R_3$ respectively, where $R_1$ is the region of stable coexistence, $R_2$ is the region of bi-stability between stable coexistence and oscillatory coexistence, $R_3$ is the region of oscillatory coexistence. Figure 2(b) shows the interactive effect of $\kappa$ and the transfer coefficient $\beta$ on the system dynamics. The BT points are $(0.179112 0.740001)$ and $(1.233338 0.740002)$ respectively which produce the same HC curve. Similarly the HC and Hopf curves separate it into three regions $R_i(i = 1, 2, 3)$ with same meanings as before. Without other stated, the meanings of $R_i(i = 1, 2, 3)$ keep fixed later.

Figure 2 shows that the bi-stability phenomena appear, which means system (2.1) will converge to different attractors under different initial data. In order to better visualize the bi-stability, we depict the graph of trajectories with $\kappa = 0.199$ (other parameter values keep unchanged) and different initial conditions, see Figure 3. With two different initial data (see the caption of Figure 3), Figure 3(a), (b) shows (2.1) converges to the stable coexistence and oscillatory coexistence; Figure 3(c), (d) indicates that it converges to two different oscillatory cycles; Figure 3(e), (f) shows that it converges to the same oscillatory cycle.

● Bifurcation of refuge

Due to the fear effect of predator, the prey will seek for some refuges to prevent from the predation of predator. Sometimes some refuge zones are constructed by people to preserve the prey from the risk of being predated. Undoubtedly, the refuge can protect the prey and keep them from extinction, so we continue to study the effect of refuge parameter $\mu$ on the stability of equilibrium point of system (2.1). Numerically we carry out the trajectories of the equilibrium point about parameter $\mu$, see Figure 4. System (2.1) is unstable when $\mu$ is small, but at $\mu = 0.50117$, it becomes stable from periodic fluctuation. The protecting role of refuge on the prey is clearly showed. With the increasing of $\mu$, the prey goes to the maximum value at $\mu = 0.9825$ and a transcritical bifurcation appears. When $\mu > 0.9825$, the predator is extinct and the prey keeps stable in the equilibrium state $\overline{Q}(20, 0)$. Compared with Figure 4(a), Figure 4(b) shows the oscillatory behavior of equilibrium point. From Figure 4, we can be aware the importance for refuge on the protection the prey.

Similarly we give the two-parameter graphs (Figure 3), which explain the interplay effect of $\mu$ and $f$ or $\beta$ on the system stability respectively, where the LPC curve (saddle-node bifurcation of limit cycle for simplicity) of Figure 5(a) originates from the GH point (generalized Hopf point) $(\mu, \beta) = (-0.125707 0.741974)$, and the BT point of Figure 5(b) is $(\mu, \beta) = (0.538642 0.740012)$, which produce the HC curve. The whole plane is divided into three same regions $R_i(i = 1, 2, 3)$ as before, and the bi-stability appears.

● Bifurcation of additional food supplement

For a prey-predator system with fear effect of predator, one effective way to protect the prey is to build refuge zone, another way is to supply additional food for the predator, when there are enough food for predator, they will reduce their attacks on prey so that the prey can survive naturally, the additional food supplement can not only keep the predator's survival, but also benefit to the living of prey, which is applied usually as a biological control strategy. Parameters $\alpha$ and $\delta$ are two crucial indexes representing the quantity and quality of additional food. So we proceed with the study of the effect of $\alpha$ and $\delta$ on the system dynamics. We carry out some numerical simulations about them. See Figure 6. Figure 6(a) shows the bifurcation occurs at $\alpha = 0.6130$ where $y = 3.8764$. Figure 6(b) shows Hopf bifurcation appears at the point $(y, \delta) = (3.860116, 0.402164)$. When $\delta = 0.166667$; A transcritical bifurcation occurs where the prey is extinct and the predator $y = 3.2184$, meanwhile a limit point appears at $(y, \delta) = (3.2184, 0.166667)$. At $\delta = 6.83333$, the transcritical bifurcation appears again where the predator is extinct and prey gets the peak value with $x = 20$. Figure 6 shows that when the account of additional food is relatively low, both the quantity and quality bring positive roles to the survival of predator, while Figure 6(b) shows that if $\delta$ is large enough (for example $\delta > 2$), then the quality of additional food $\delta$ produces negative role, and even leads to the extinction of predator, which reveals that the control strategy should be applied rationally in a suitable region.

To see the interference of $\alpha$ and $\delta$, we depict the two parameter graphs, see Figure 7. The HC curve in $\alpha-\delta$ plane originates from the BT point $(\alpha, \delta) = (0.801951, 0.376521)$, while in $\alpha-\beta$ plane, HC curve originates from the BT point $(\alpha, \beta) = (0.75341, 0.65322)$. The meaning of $R_i(i = 1, 2, 3)$ is the same as before.

Numerically, we find that the fear $f$, anti-predation sensitivity $\kappa$, the refuge $\mu$, and the additional food parameters $\alpha$ and $\delta$ bring large influence to the system dynamics including both positive and negative effect, which result in the system complexity of system (2.1) and should be controlled rationally.

4.2.   Stability analysis

Now we give some numerical examples to validate the criteria of the global asymptotic stability. We start with the verification of Theorem 3.5. For system (2.1), we fix a set of parameter values as below:

Then, the predator extinction equilibrium point is $\overline{Q}(20, 0)$. We compute that the conditions of Theorem 3.5 are satisfied, so it is GAS, see Figure 8.

Next we fix a set of parameter values as below:

By computation, the interior equilibrium point of (2.1) is $\widetilde{Q}(1.75, 7.4754)$. We verify that the conditions of Theorem 3.6 hold, and it is GAS. See Figure 9.

5.   Discussion and conclusions

Considering the joint role of fear effect and direct consumption of predator populations on prey populations, we formulated prey-predator mathematical model with fear and Holling-Ⅱ type functional response. Due to the existence of counter-predation behavior of prey, we introduced an anti-predation sensitivity parameter to describe the level of anti-predation behavior. In addition, in order to prevent the prey from being predated too quickly and protect the prey for biological balance, we add the additional food supplement for predator. By use of Lyapunov stability theory, we investigated the existence of equilibrium points and their stability. See Table 1. Finally, the bifurcation analysis is executed by applying Sotomayor's bifurcation theorem.

Numerically, we explored the effects on the system stability of such parameters as fear level, anti-predation sensitivity, refuge level and the account of addition food. Particularly, we obtained the bifurcation thresholds of these parameters, which are summarized in the following Table 2.

The main findings show that the fear significantly affect the system stability, the construction of refuge area and the releasing of additional food for predator are two kinds of effective ways to reduce the effect induced by the fear of predator and can preserve both predator and prey for their natural survival. While these measures have both positive and negative effects, it can not only keep the continuous survival of species, but also make the species be extinct if they are wrongly used. Therefore, rational use of these measures is very crucial for people to control the ecological balance and continuous development. If we take into the complexity of systems, there are still many factors should be considered, for example, the delay exists in almost every process and the fear delay induced by predator is inevitable. On the other hand, the environmental fluctuation, like earthquake, always appear, and the disaster caused by which can affect the system dynamics seriously. All these are left for us to study in the near future.

Acknowledgments

The paper is supported by NSF of China (Nos. 11861027, 12161062), Science and Technology Program of Inner Mongolia Autonomous Region (Nos. 2022YFHH0017, 2022YFHH0063) and Natural Science Foundation of Inner Mongolia (No. 2021MS01005).

Conflict of interest

The authors declare there is no conflict of interest.