Stability and Bifurcation Analysis of a Predator-Prey Model with Michaelis-Menten Type Prey Harvesting

Wenchang Chen

School of Mathematics and Physics, Wenzhou University, Wenzhou, China.

DOI: 10.4236/jamp.2022.102038                

Abstract

In this paper, we investigated stability and bifurcation behaviors of a predator-prey model with Michaelis-Menten type prey harvesting. Sufficient conditions for local and global asymptotically stability of the interior equilibrium point were established. Some critical threshold conditions for transcritical bifurcation, saddle-node bifurcation and Hopf bifurcation were explored analytically. Furthermore, It should be stressed that the fear factor could not only reduce the predator density, but also affect the prey growth rate. Finally, these theoretical results revealed that nonlinear Michaelis-Menten type prey harvesting has played an important role in the dynamic relationship, which also in turn proved the validity of theoretical derivation.

Keywords

Michaelis-Menten Type Prey Harvesting, Stability, Bifurcation, Limit Cycle

1. Introduction

The predator-prey model is one of the dominant population models, which has been researched extensively to comprehend interactions between various species in a fluctuant natural environment [1] [2] [3] [4] [5]. However, in reality, predators are unable to catch all prey because prey can decrease the risk of predation by using refuge [6] [7] [8]. So the prey refuge will more accord with authentic circumstance, and many scholars have achieved considerable process in this field [9] [10] [11] [12].

However, many researchers only consider directly killing of prey by predator, but ignore the impact of the presence of predator on prey. Zanette et al. [13] experimentally showed that the predation fears can reduce offspring production by 40%, which is even more intensively than the impact of direct hunting. In fact, all prey can exhibit different kinds of anti-predator behaviors, such as changes of habitat, physiological and foraging [14] [15] when they are confronted with predation risks. Many studies in this direction have been carried out and obtained numerous attractive results [16] [17] [18] [19] [20]. Particularly, Wang et al. [19] firstly proposed a predator-prey model with the cost of fear:

$(\begin{array}{l}\frac{\text{d}u}{\text{d}t}=r_{0}uf\left(k,v\right)-du-a{u}^2-g\left(u\right)v,\\ \frac{\text{d}v}{\text{d}t}=-mv+cg\left(u\right)v,\end{array}$

where k is the level of fear, which causes anti-predator behaviors of prey. Biologically speaking, $f\left(k\mathrm{,}v\right)$ can reasonably obey the following conditions:

$f\left(0,v\right)=1,f\left(k,0\right)=1,\underset{k\to \infty }{\lim }f\left(k,v\right)=0,$

$\underset{v\to \infty }{\lim }f\left(k,v\right)=0,\frac{\partial f\left(k,v\right)}{\partial k}<0,\frac{\partial f\left(k,v\right)}{\partial v}<0.$

They concluded that the cost of fear does not change the dynamical behaviors of this model and the unique equilibrium is globally asymptotically stable when the model with the linear functional response.

Inspired by the insightful work [19], Chakraborty [20] proposed a predator-prey model with fear factor, which is given by:

$(\begin{array}{l}\frac{\text{d}u}{\text{d}t}=ru\left(1-\frac{u}{K}\right)\left(u-\theta \right)\frac{1}{1+fv}-auv,\\ \frac{\text{d}v}{\text{d}t}=\alpha auv-mv,\end{array}$ (1)

where u and v represent prey and predator densities at time t, respectively. r is the intrinsic growth rate, k is the carrying capacity of prey, a is the predation rate, $\alpha$ is the conversion factor, m is the mortality rate of predator, $\theta$ is the Allee

threshold. $\frac{1}{1+fv}$ is the fear factor term. They showed that fear can dramatically

lessen the per-capita growth rate, but cannot affect the equilibrium stability, but can generate richer dynamics such as bi-stability.

From the perspective of human needs and long-term progress, the exploitation of natural resources and the storage of renewable energy, harvesting are always one of the most crucial factors in the dynamics of predator-prey model [21]. We find that the predator-prey model with harvesting can lead to more complex properties than the model without it [22] [23] [24], which inspires us to take harvesting into account. Harvesting regimes can be classified into three main types: 1) $h\left(x\right)=h$, constant rate harvesting, 2) $h\left(x\right)=qEx$, constant-effort harvesting, and 3) $h\left(x\right)=\frac{qEx}{m_{1}E+m_{2}x}$, nonlinear harvesting [25] [26]. Huang et al. [27] systematically studied the dynamics of a Leslie-Gower type predator-prey model with constant-yield predator harvesting. They have shown that the model can have various kinds of bifurcations, such as saddle-node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcations. These results reveal that the harvested model can exhibit more complex dynamics compared to the model without harvesting. Refs. [28] [29] also paid close attention to investigating the impact of constant-yield harvesting in predator-prey models. A ratio-dependent predator-prey model with non-constant predator harvesting rate was analyzed by Lajmiri et al. [30]. They investigated the stability of equilibria and some bifurcation behaviors. Gao et al. [31] considered the temporal, spatial and spatiotemporal dynamics of a ratio-dependent predator-prey diffusive model, in which the predator subject to constant effort harvesting. Lenzini and Rebaza [32] discovered several bifurcations and connecting orbits by studying a ratio-dependent predator-prey model with non-constant harvesting. Singh et al. [33] proposed a Leslie-Gower predator-prey model with Michaelis-Menten type predator harvesting to explore the stability and bifurcation behaviors. Liu et al. [34] presented a super cross-diffusion predator-prey model to study its Turing instability and pattern formation. Chen et al. [35] proposed a predator-prey model with nonlinear harvesting to consider the effect of diffusion on pattern selection. The dynamics of a delayed diffusive predator-prey model with nonlinear predator harvesting has been investigated by Liu and Zhang [36], and they obtained the conditions for Turing and Hopf bifurcation, and found that the delay has remarkable impact on the emergent spatial patterns. Liu and Jiang [37] discussed the stability and Hopf bifurcation in detail for a Gause predator-prey model with gestation delay, Michaelis-Menten prey harvesting and Holling type III functional response. Zhang et al. [38] proposed stochastic non-autonomous predator-prey models with and without impulse; they concluded that the model dynamics can be appreciably influenced by the stronger noises and nonlinear harvesting, which also can cause the extinction of the predator species.

Stimulated by the above review of literature, we will propose a predator-prey model with Michaelis-Menten type prey harvesting and prey refuge based on model (1), which can be expressed by the following equation,

$(\begin{array}{l}\frac{\text{d}X}{\text{d}\tau }=rX\left(1-\frac{X}{K}\right)\frac{1}{1+f_{1}Y}-a_1\left(1-m\right)XY-\frac{qEX}{m_{1}E+m_{2}X},\\ \frac{\text{d}Y}{\text{d}\tau }=e_1{a}_1\left(1-m\right)XY-d_{1}Y,\end{array}$ (2)

where q is the catchability coefficient, E is the effort used on harvesting the prey species, $m_1$ and $m_2$ are proper constants. The rest of the parameters have the same meanings as model (1). For simplicity, we will nondimensionalize the model (2) with the substitutions such as $x=\frac{X}{K}$, $y=\frac{a_{1}Y}{r}$ and $t=r\tau$, hence we have

$(\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(1-x\right)\frac{1}{1+fy}-\left(1-m\right)xy-\frac{hx}{c+x}=F\left(x,y\right),\\ \frac{\text{d}y}{\text{d}t}=e\left(1-m\right)xy-dy=G\left(x,y\right),\end{array}$ (3)

where $h=\frac{qE}{r{m}_{2}K}$, $c=\frac{m_{1}E}{m_{2}K}$, $e=\frac{a_1{e}_{1}K}{r}$, $f=\frac{r{f}_1}{a_1}$ and $d=\frac{d_1}{r}$ are positive constants.

The objective of this paper is to provide a stability and bifurcation analysis of model (3). The rest of the paper is arranged as follows. The existence of equilibria and their stability are presented in Section 2. Section 3 deals with the bifurcation, such as transcritical bifurcation, saddle-node bifurcation and Hopf bifurcation. In section 4, we analyze the impact of fear and prey refuge. Numerical simulations for model (3) are given in section 5 to illustrate the theoretical works. Finally, a brief conclusion is drawn in section 6.

2. Existence and Stability of Equilibria

2.1. Existence of Equilibria

Model (3) always has the trivial equilibrium point $E_0=\left(\mathrm{0,0}\right)$, and the other equilibria are the intersection of the horizontal isocline and vertical isocline of model (3), which are given by

$(\begin{array}{l}F\left(x,y\right)=x\left(1-x\right)\frac{1}{1+fy}-\left(1-m\right)xy-\frac{hx}{c+x}=0,\\ G\left(x,y\right)=e\left(1-m\right)xy-dy=0.\end{array}$ (4)

It is obvious that the possible boundary equilibria are the positive roots of the quadratic equation $F\left(x,0\right)=0$, i.e.,

$1-x-\frac{h}{c+x}=0,$ (5)

which discriminant is $\Delta \left(h\right)={\left(1+c\right)}^2-4h$. Denote $h^{*}=\frac{{\left(1+c\right)}^2}{4}$ when $\Delta \left(h\right)=0$. then the equation of (5) has two roots $x_1=\frac{1-c-\sqrt{\Delta }}{2}$ and $x_2=\frac{1-c+\sqrt{\Delta }}{2}$ when $h<h^{*}$ ; a double root $x_{SN}=\frac{1-c}{2}$ when $h=h^{*}$ ; no real root when $h>h^{*}$. We summarize our finding in the following Lemma.

Lemma 1. Besides $E_0$, the existence of boundary equilibria of model (3) is as follows.

1) If $c\in \left(\mathrm{0,1}\right)$, the existence of boundary equilibria is concluded in Table 1.

2) If $c=1$, there exists an equilibrium point $E_3=\left(x_3,0\right)=\left(\sqrt{1-h},0\right)$ when $0<h<1$.

3) If $c>1$, there exists an equilibrium point $E_4=\left(x_4,0\right)=\left(\frac{1-c+\sqrt{\Delta }}{2},0\right)$ when $0<h<c<h^{*}$.

Evidently, the interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}=\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is the intersection of the nullclines. By simple calculation, we obtain the abscissa of interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ is $x^{*}=\frac{d}{e\left(1-m\right)}$, and the ordinate $y^{\mathrm{<mo>\; *\; </mo>}}$ is the positive root of the equa-

Table 1. Boundary equilibria besides $E_0$ of model (3) when $0<c<1$.

tion $G\left(x^{*},y\right)=0$ if $H<0$ hold. That is, $y^{*}=\frac{-\left(a+hf\right)+\sqrt{\Delta }}{2af}$, where $a=\left(1-m\right)\left(c+x^{\mathrm{<mo>\; *\; </mo>}}\right)$, $H=h-\left(1-x^{*}\right)\left(c+x^{*}\right)$ and $\Delta ={\left(a+hf\right)}^2-4afH$. Hence, we can have the following result.

Lemma 2. The positive equilibrium point $E^{*}=\left(x^{*},y^{*}\right)$ exists when $H<0$.

Where

$x^{*}=\frac{d}{e\left(1-m\right)},y^{*}=\frac{-\left(a+hf\right)+\sqrt{{\left(a+hf\right)}^2-4afH}}{2af},$ (6)

with $a=\left(1-m\right)\left(c+x^{\mathrm{<mo>\; *\; </mo>}}\right)$, $H=h-\left(1-x^{*}\right)\left(c+x^{*}\right)$ and $\Delta ={\left(a+hf\right)}^2-4afH$.

2.2. Stability of Equilibria

Now we analyze the local stability of equilibria identified above. The general Jacobian matrix of model (3) takes the form of

$J=\left(\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right),$ (7)

where

$J_{11}=\left(1-x\right)\frac{1}{1+fy}-\left(1-m\right)y-\frac{h}{c+x}+x\left[-\frac{1}{1+fy}+\frac{h}{{\left(c+x\right)}^2}\right],$

$J_{12}=-x\left[\left(1-x\right)\frac{f}{{\left(1+fy\right)}^2}+\left(1-m\right)\right],$

$J_{21}=e\left(1-m\right)y,$

$J_{22}=e\left(1-m\right)x-d.$

Theorem 1. The origin $E_0=\left(0,0\right)$ is a saddle point if $h<c$ and a stable node if $h>c$ or $h=c=1$ ; when $h=c\ne 1$, $E_0$ is a saddle-node.

Proof. The Jacobian matrix of model (3) at the equilibrium point $E_0$ is given by:

$J\left(E_0\right)=\left(\begin{array}{cc}1-\frac{h}{c}& 0\\ 0& -d\end{array}\right).$

Obviously, the eigenvalues of $J\left(E_0\right)$ are ${\lambda }_1=1-\frac{h}{c}$ and ${\lambda }_2=-d$. According to the Routh-Hurwitz criterion, $E_0$ is a saddle point if $h<c$ and a stable node if $h>c$. When $h=c$, we have ${\lambda }_1=0$. To determinate the stability of $E_0$, we rescale t by $\tau =-dt$ and expand the obtained model in power series around $E_0$, under which we get

$(\begin{array}{l}\frac{\text{d}x}{\text{d}\tau }=\frac{1}{d}\left(1-\frac{1}{c}\right)x^2+\frac{1}{d}\left(1+f-m\right)xy-\frac{1}{d{c}^2}{x}^3-\frac{f}{d}{x}^{2}y-\frac{f}{d}x{y}^2+p_1\left(x,y\right),\\ \frac{\text{d}y}{\text{d}\tau }=y-\frac{e}{d}\left(1-m\right)xy,\end{array}$ (8)

where $P_1\left(x\mathrm{,}y\right)$ is a power series in $\left(x\mathrm{,}y\right)$ with terms $x^i{y}^j$ satisfying $i+j\ge 4$. Hence, by Theorem 7.1 in [39], if the coefficient of $x^2$ is $\frac{1}{d}\left(1-\frac{1}{c}\right)\ne 0$, i.e., $h=c\ne 1$, then $E_0$ is a saddle-node; if $\frac{1}{d}\left(1-\frac{1}{c}\right)=0$, i.e., $h=c=1$, then $E_0$ is a unstable node due to $\frac{1}{d{c}^2}>0$. However, the orbits with time going in the opposite direction because of the transformation $\tau =-dt$. Hence, $E_0$ is a stable node. □

Theorem 2. The equilibrium point $E_{SN}=\left(\frac{1-c}{2},0\right)$ is a saddle-node if $\frac{e\left(1-m\right)\left(1-c\right)}{2}-d\ne 0$, and a saddle point if $\frac{e\left(1-m\right)\left(1-c\right)}{2}-d=0$.

Proof. By replacing $\left(x\mathrm{,}y\right)$ in matrix (7) with $E_{SN}$, The Jacobian matrix of model (3) at $E_{SN}$ is given by:

$J\left(E_{SN}\right)=\left(\begin{array}{cc}0& -x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]\\ 0& \frac{e\left(1-m\right)\left(1-c\right)}{2}-d\end{array}\right).$

Clearly, we find that the eigenvalues of $J\left(E_{SN}\right)$ are ${\lambda }_1=0$ and ${\lambda }_2=\frac{e\left(1-m\right)\left(1-c\right)}{2}-d$. We therefore consider the following two cases.

1) ${\lambda }_2\ne 0$. Transforming $E_{SN}$ to the origin by the translation $X=x-x_{SN}$, $Y=y$ and expanding model (3) in a power series around the origin, under which model (3) becomes

$(\begin{array}{l}\frac{\text{d}X}{\text{d}t}=-x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]Y-\frac{h{x}_{SN}}{c+x_{SN}^3}{X}^2\\ +\left[f\left(2x_{SN}-1\right)-\left(1-m\right)\right]XY+f{x}_{SN}\left(1-x_{SN}\right)Y^2+P_2\left(X,Y\right),\\ \frac{\text{d}Y}{\text{d}t}=\left[e\left(1-m\right)x_{SN}-d\right]Y-e\left(1-m\right)XY,\end{array}$ (9)

where $P_2\left(X\mathrm{,}Y\right)$ is a power series in $\left(X\mathrm{,}Y\right)$ with terms $X^i{Y}^j$ satisfying $i+j\ge 3$. Then, under the transformation

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}1& \frac{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{e\left(1-m\right)x_{SN}-d}\\ 0& 1\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right).$

Model (9) becomes

$(\begin{array}{l}\frac{\text{d}X}{\text{d}t}=a_{20}{x}^2+a_{11}xy+a_{02}{y}^2+P_3\left(x\mathrm{,}y\right)\mathrm{,}\\ \frac{\text{d}Y}{\text{d}t}=b_{01}y+b_{11}xy+b_{02}{y}^2\mathrm{,}\end{array}$ (10)

where $P_3\left(x\mathrm{,}y\right)$ is a power series in $\left(x\mathrm{,}y\right)$ with terms $x^i{y}^j$ satisfying $i+j\ge 3$ and

$\begin{array}{c}{a}_{11}=-\frac{2h{x}_{SN}^2\left[\left(1-x_{SN}\right)f+1-m\right]}{{\left(c+x_{SN}\right)}^3\left[e\left(1-m\right)x_{SN}-d\right]}+f\left(2x_{SN}-1\right)-\left(1-m\right)\\ +\frac{e\left(1-m\right)x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{e\left(1-m\right)x_{SN}-d},\end{array}$

$\begin{array}{c}{a}_{02}=-\frac{h{x}_{SN}^3{\left[\left(1-x_{SN}\right)f+1-m\right]}^2}{{\left(c+x_{SN}\right)}^3\left[e\left(1-m\right)x_{SN}-d\right]}+f{x}_{SN}\left(1-x_{SN}\right)\\ +\frac{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{e\left(1-m\right)x_{SN}-d}\left[f\left(2x_{SN}-m-2\right)\right]\\ +\frac{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{e\left(1-m\right)x_{SN}-d}\left[\frac{e\left(1-m\right)x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{e\left(1-m\right)x_{SN}-d}\right],\end{array}$

$a_{20}=-\frac{h_{SN}}{c+x_{SN}^3},b_{01}=e\left(1-m\right)x_{SN}-d,b_{11}=e\left(1-m\right),$

$b_{20}=-\frac{e\left(1-m\right)x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{e\left(1-m\right)x_{SN}-d}.$

Introducing a new time variable $\tau =\left[e\left(1-m\right)x_{SN}-d\right]t$, since the coefficient of $x^2$ is $-\frac{h{x}_{SN}}{{\left(c+x_{SN}\right)}^3\left[e\left(1-m\right)x_{SN}-d\right]}\ne 0$, by Theorem 7.1 in [39] again, $E_{SN}$ is a saddle-node.

2) ${\lambda }_2=0$. Let $\tau =-x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]t$, then model (9) can be rewritten as

$\{\begin{array}{l}\frac{\text{d}X}{\text{d}t}=Y+\frac{h}{{\left(c+x_{SN}\right)}^3\left[\left(1-x_{SN}\right)f+1-m\right]}{X}^2-\frac{f\left(2x_{SN}-1\right)-\left(1-m\right)}{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}XY\\ -\frac{f{x}_{SN}\left(1-x_{SN}\right)}{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}+o\left({\left|X,Y\right|}^2\right)=Y+P_4\left(X,Y\right),\\ \frac{\text{d}Y}{\text{d}t}=-\frac{e\left(1-m\right)x_{SN}-d}{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}XY=Q_4\left(X,Y\right).\end{array}$ (11)

From $Y+P_4\left(X,Y\right)=0$, we obtain an implicit function

$\begin{array}{c}\phi \left(X\right)=-\frac{h}{{\left(c+x_{SN}\right)}^3\left[\left(1-x_{SN}\right)f+1-m\right]}{X}^2\\ +\frac{h\left[f\left(2x_{SN}-1\right)-\left(1-m\right)\right]}{x_{SN}\left[\left(1-x_{SN}\right)f+1-m\right]}{X}^3+\cdots ,\end{array}$

and we have

$\begin{array}{c}\psi \left(X\right)=-\frac{he\left(1-m\right)}{x_{SN}{\left(c+x_{SN}\right)}^3{\left[\left(1-x_{SN}\right)f+1-m\right]}^2}{X}^3\\ +\frac{he\left(1-m\right)\left[f\left(2x_{SN}-1\right)-\left(1-m\right)\right]}{x_{SN}^2{\left(c+x_{SN}\right)}^3{\left[\left(1-x_{SN}\right)f+1-m\right]}^4}{X}^4+\cdots ,\end{array}$

$\delta \left(X\right)=-\frac{2h_{SN}+e\left(1-m\right){\left(c+x_{SN}\right)}^3}{x_{SN}{\left(c+x_{SN}\right)}^3\left[\left(1-x_{SN}\right)f+1-m\right]}X+{\left[X\right]}_2.$

Using the notations of Theorems 7.2 in [39], we have $k=2m+1$, $m=1$ ;

$a_k=\frac{he\left(1-m\right)}{x_{SN}{\left(c+x_{SN}\right)}^3\left[\left(1-x_{SN}\right)f+1-m\right]}>0$. Hence, $E_{SN}$ is a saddle point. □

Theorem 3. The equilibrium point $E_1$ is always unstable. Especially, $E_1$ is a unstable node if $e\left(1-m\right)x_1-d>0$ and a saddle if $e\left(1-m\right)x_1-d<0$.

Proof. The Jacobian matrix of model (3) at $E_1$ is given by:

$J\left(E_1\right)=\left(\begin{array}{cc}\left(1-x_1-\frac{h}{c+x_1}\right)+x_1\left(-1+\frac{h}{{\left(c+x_1\right)}^2}\right)& -x_1\left[\left(1-x_1\right)f+1-m\right]\\ 0& e\left(1-m\right)x_1-d\end{array}\right).$

We note that one of the eigenvalues of $J\left(E_1\right)$ is ${\lambda }_1=\left(1-x_1-\frac{h}{c+x_1}\right)+x_1\left(-1+\frac{h}{{\left(c+x_1\right)}^2}\right)$, which follows that

${\lambda }_1>\sqrt{\Delta \left(h\right)}+c-\frac{4hc}{{\left(1+c\right)}^2}>\sqrt{\Delta \left(h\right)}+c\left(1-\frac{h}{h^{*}}\right)>0.$

So the stability of $E_1$ is determined by the eigenvalue ${\lambda }_2=e\left(1-m\right)x_1-d$ and the conclusion is tenable. □

As discussed in the previous section, we know that $E_2$ exist when $0<c<1$ and $h^{\ast }>h>c>0$ ; $E_3$ exist when $c=1$ and $0<h<1$ ; $E_4$ exist when $c>1$ and $h^{\ast }>c>h>0$. We now discuss the stability of the boundary equilibrium $E_i\left(i=2,3,4\right)$.

Theorem 4. The boundary equilibrium point $E_i\left(i=2,3,4\right)$ is a stable node if $e\left(1-m\right)x_i-d<0$ and a saddle point if $e\left(1-m\right)x_i-d>0$.

Proof. Similar to the proof of Theorem 3, the Jacobian matrix of model (3) at $E_i$ is

$J\left(E_i\right)=\left(\begin{array}{cc}\left(1-x_i-\frac{h}{c+x_i}\right)+x_i\left(-1+\frac{h}{{\left(c+x_i\right)}^2}\right)& -x_1\left[\left(1-x_i\right)f+1-m\right]\\ 0& e\left(1-m\right)x_i-d\end{array}\right).$

Then the eigenvalues of $J\left(E_i\right)$ are ${\lambda }_1=\left(1-x_i-\frac{h}{c+x_i}\right)+x_i\left(-1+\frac{h}{{\left(c+x_i\right)}^2}\right)$

and ${\lambda }_2=e\left(1-m\right)x_i-d$. One can easily prove that the following two expressions

${\left[{\lambda }_1\right]|}_{x_i=\frac{1-c+\sqrt{\Delta }}{2}}=x_i\left[-1+\frac{4h}{{\left(1+c+\sqrt{\Delta }\right)}^2}\right]<0$, for $i=2$ and 4,

${\left[{\lambda }_1\right]|}_{x_3=\sqrt{1-h}}=x_2\left[-1+\frac{4h}{{\left(1+\sqrt{1-h}\right)}^2}\right]<0$,

hold on the basis of existence conditions of $E_i$ from Lemma 1, which means that the stability of $E_i\left(i=2,3,4\right)$ is determined by the sign of $e\left(1-m\right)x_i-d$. Therefore, this theorem follows. □

Next, we study the local stability of the positive equilibrium $E^{*}=\left(x^{*},y^{*}\right)$, and give sufficient conditions on its global stability.

Theorem 5. The interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ is locally asymptotically stable when $0<h<h_{Hp}$ with

$h_{Hp}=-\frac{\left[2\sqrt{e\left\{\left[\left(c-1\right)f+\frac{1}{4}\left(m-1\right)\right]\left(m-1\right)-2df\right\}}+e\left(m-1\right)\right]{\left[c\left(m-1\right)e-d\right]}^2}{2e^{2}f\left(m-1\right)\left[\left(m-1\right)\left(c-1\right)e-2d\right]}$ (12)

and unstable when $h>h_{Hp}$.

Proof. The Jacobian matrix at $E^{\mathrm{<mo>\; *\; </mo>}}$ is:

$J\left(E^{*}\right)=\left(\begin{array}{cc}{x}^{*}\left[-\frac{1}{1+f{y}^{*}}+\frac{h}{{\left(c+x^{*}\right)}^2}\right]& -x^{*}\left[\left(1-x^{*}\right)\frac{f}{{\left(1+f{y}^{*}\right)}^2}+1-m\right]\\ e\left(1-m\right)y^{*}& 0\end{array}\right).$

Then we can easily gain the determinant of the matrix $J\left(E^{*}\right)$ is given by

$DetJ\left(E^{*}\right)=e\left(1-m\right)x^{*}{y}^{*}\left[\left(1-x^{*}\right)\frac{f}{{\left(1+f{y}^{*}\right)}^2}+\left(1-m\right)\right]>0,$

which implies that $E^{\mathrm{<mo>\; *\; </mo>}}$ is an antisaddle. It is easy to see the trace is

$TrJ\left(E^{*}\right)=x^{*}\left[-\frac{1}{1+f{y}^{*}}+\frac{h}{{\left(c+x^{*}\right)}^2}\right].$

By simple calculation, we find that $TrJ\left(E^{*}\right)<0$ is equivalent to $h<h_{Hp}$. Therefore, $J\left(E^{*}\right)$ has two negative eigenvalues when $h<h_{Hp}$, which indicates that $E^{\mathrm{<mo>\; *\; </mo>}}$ is locally asymptotically stable. On the contrary, both characteristic roots of $J\left(E^{*}\right)$ are positive and $E^{\mathrm{<mo>\; *\; </mo>}}$ loses its stability when $h>h_{Hp}$ hold. □

In Theorem 5, we have proved that the interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ of model (3) is locally asymptotically stable when $0<h<h_{Hp}$. Now we are going to give adequate conditions on its global stability.

Theorem 6. Suppose $0<h<h_{Hp}$, then the interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ of model (3) is globally asymptotically stable in the positive quadrant if one of the following conditions holds.

(H₁) $c>1$ and $h^{*}>c>h>0$ ;

(H₂) $c=1$ and $0<h<1$ ;

(H₃) $c<1$ and $h^{*}>c>h>0$.

Proof. Under the above conditions, we find that model (3) has the other two boundary equilibria $E_0$ and $E_2$, which are unstable. Meanwhile, the interior equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ is locally asymptotically stable. It is easy to prove that the interior of $R_{+}^2$ is the invariant set of model (3). Moreover, $F\left(x\mathrm{,}y\right)$ and $G\left(x\mathrm{,}y\right)$, respectively, represent the right hand side function of model (3). Choose the

Dulac function as $B\left(x\mathrm{,}y\right)=\frac{1}{xy}$. After calculations, we have

$D=\frac{\partial \left(BF\right)}{\partial x}+\frac{\partial \left(BG\right)}{\partial y}=\frac{1}{y}\left[\frac{h}{{\left(c+x\right)}^2}-\frac{1}{1+fy}\right]<0$

in the interior of $R_{+}^2$. By the Dulac theorem [40], there exists no periodic orbit in the first quadrant. Thus, $E^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable. □

Theorem 7. Suppose $h>h_{Hp}$ and $e\left(1-m\right)x_2-d>0$, then model (3) exists a limit cycle if one of the conditions of H_i (i = 1, 2, 3) holds.

Proof. By previous calculations, here, $E^{\mathrm{<mo>\; *\; </mo>}}$ is unstable. In addition, $E_0$ and $E_2$ are saddle point. Recall $x^{*}<x_2$ and Let $L_1:x-x_2=0$, then

${\frac{\text{d}{L}_1}{\text{d}t}|}_{\left(3\right)}={\frac{\text{d}x}{\text{d}t}|}_{x=x_2}=x_2\left[\left(1-x_2\right)\frac{1}{1+fy}-\left(1-m\right)y-\frac{h}{c+x_2}\right]<0,\text{for}y\in \left(0,\infty \right).$

Next, let $L_2:y-\lambda =0$ with $\lambda >0$ to be specified later, then

${\frac{\text{d}{L}_2}{\text{d}t}|}_{\left(3\right)}={\frac{\text{d}y}{\text{d}t}|}_{y=\lambda }=\lambda \left[e\left(1-m\right)x-d\right]<0,\text{for}x\in \left(0,x^{*}\right).$

Moreover, let $L_3:2\left(x_2-x^{*}\right)\left(y-\lambda \right)+\lambda \left(x-x^{*}\right)=0$.

Since $0<x<x_2$, we get $x-x^{*}<x_2-x^{*}$ and $y=\lambda -\frac{\lambda \left(x-x^{*}\right)}{2\left(x^2-x^{*}\right)}>\frac{1}{2}\lambda$. For $0<x<x^{*}$, $y>0$. Simply calculation can give

$\begin{array}{c}{\frac{\text{d}{L}_3}{\text{d}t}|}_{\left(3\right)}=2\left(x_2-x^{*}\right)\frac{\text{d}y}{\text{d}t}+\lambda \frac{\text{d}x}{\text{d}t}\\ =2\left(x_2-x^{*}\right)\left[e\left(1-m\right)x-d\right]y+\lambda \left[\left(1-x\right)\frac{1}{1+fy}-\left(1-m\right)y-\frac{h}{c+x}\right]\\ <\frac{1}{\left(1-fy\right)\left(c+x\right)}f\left(\lambda \right).\end{array}$

where

$\begin{array}{c}f\left(\lambda \right)=-\frac{fx\left(1-m\right)\left(c+x\right)}{4}{\lambda }^3\\ +\left\{2f\left(x_2-x^{*}\right)\left(c+x\right)\left[e\left(1-m\right)x-d\right]-\frac{x}{2}\left[hf+\left(1-m\right)\left(c+x\right)\right]\right\}{\lambda }^2\\ +\left\{2\left(x_2-x^{*}\right)\left(c+x\right)\left[e\left(1-m\right)x-d\right]+x\left[hf+\left(1-m\right)\left(c+x\right)\right]\right\}\lambda .\end{array}$

Due to $fx\left(1-m\right)\left(c+x\right)>0$, it follows that ${\frac{\text{d}{L}_3}{\text{d}t}|}_{\left(3\right)}<0$ for sufficiently large $\lambda >0$. Thus, when $h>h_{Hp}$ and $e\left(1-m\right)x_2-d>0$, by Poincare-Bendixson Theorem [40], model (3) exists a limit cycle if one of the conditions of (H_i) (i= 1, 2, 3) holds. □

3. Local Bifurcation Analysis

In this section, we will discuss possible bifurcations of model (3), and derive the conditions for various bifurcations. The following lemma is given for proving saddle-node bifurcation and transcritical bifurcation.

Lemma 3. (Sotomayor’s Theorem [41] ). Consider the system $\stackrel{\dot{}}{x}=F\left(x,\mu \right)$ and assume that $F\left(x_0,{\mu }_0\right)=0$ at equilibrium point $x_0$ holds. The eigenvectors of zero eigenvalues of $n\times n$ Jacobian matrix $J\equiv DF\left(x_0\mathrm{,}{\mu }_0\right)$ and $J^{\text{T}}$ are V and W respectively. Then

1) Suppose

$W^{\text{T}}{F}_{\mu }\left(x_0\mathrm{,}{\mu }_0\right)\ne \mathrm{0,}$

$W^{\text{T}}\left[D^2{F}_{\mu }\left(x_0\mathrm{,}{\mu }_0\right)\left(V\mathrm{,}V\right)\right]\ne 0.$

Hence, when the bifurcation parameter $\mu$ through the thresholds value, that is, $\mu ={\mu }_0$, the system undergoes a saddle-node bifurcation at $x_0$.

2) Suppose

$W^{\text{T}}{F}_{\mu }\left(x_0\mathrm{,}{\mu }_0\right)=\mathrm{0,}$

$W^{\text{T}}\left[D{F}_{\mu }\left(x_0\mathrm{,}{\mu }_0\right)V\right]\ne 0.$

$W^{\text{T}}\left[D^2{F}_{\mu }\left(x_0\mathrm{,}{\mu }_0\right)\left(V\mathrm{,}V\right)\right]\ne 0.$

Hence, when the bifurcation parameter $\mu$ through the thresholds value, that is, $\mu ={\mu }_0$, the system undergoes a transcritical bifurcation at $x_0$.

3.1. Transcritical Bifurcation

From Lemma 1 and Theorem 1, we find that $E_0$ is a saddle when $h<c$, and stable node when $h>c$, which indicates that $E_0$ changes its stability when parameter h over the threshold value $h_{TC}=c$. Particularly, the boundary equilibrium $E_1$ bifurcates from $E_0$ when $h_{TC}=c$. The above analysis implies that the existence of transcritical bifurcation at the trivial equilibrium $E_0$ when the value of the parameter h exceeds the threshold value $h_{TC}=c$.

Theorem 8. Model (3) undergoes a transcritical bifurcation around $E_0$ at $h=h_{TC}=c\left(0<c<1\right)$.

Proof. Clearly, the matrix $J\left(E_0\right)$ has a zero eigenvalue when $h=h_{TC}=c\left(0<c<1\right)$. Let v and w be the eigenvectors associate with the zero eigenvalue for $J\left(E_0\right)$ and $J^{\text{T}}\left(E_0\right)$, respectively. We get

$v=\left(\begin{array}{c}1\\ 0\end{array}\right)\text{and}w=\left(\begin{array}{c}1\\ 0\end{array}\right).$

Defining

$F\left(x,y\right)=\left(\begin{array}{c}{F}_1\left(x,y\right)\\ F_2\left(x,y\right)\end{array}\right)=\left(\begin{array}{c}x\left[\left(1-x\right)\frac{1}{1+fy}-\left(1-m\right)y-\frac{h}{c+x}\right]\\ \left[e\left(1-m\right)x-d\right]y\end{array}\right),$

and calculating, we have

$F_h\left(E_0\mathrm{,}{h}_{TC}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)\mathrm{,}$

$D{F}_h\left(E_0,h_{TC}\right)v=\left(\begin{array}{cc}-\frac{1}{c}& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)=\left(\begin{array}{c}-\frac{1}{c}\\ 0\end{array}\right),$

$D^2{F}_h\left(E_0,h_{TC}\right)\left(v,v\right)=\left(\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{V}_1^2+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{V}_1{V}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{V}_2^2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{V}_1^2+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{V}_1{V}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{V}_2^2\end{array}\right)=\left(\begin{array}{c}2\left(-1+\frac{1}{c}\right)\\ 0\end{array}\right).$

Using the expressions for v and w we get,

$W^{\text{T}}{F}_h\left(E_0\mathrm{,}{h}_{TC}\right)=\mathrm{0,}$

$W^{\text{T}}\left[D{F}_h\left(E_0,h_{TC}\right)v\right]=-\frac{1}{c}\ne 0,$

$W^{\text{T}}\left[D^2{F}_h\left(E_0,h_{TC}\right)\left(v,v\right)\right]=2\left(-1+\frac{1}{c}\right)\ne 0.$

Thus from Lemma 3, we can conclude that model (3) undergoes a transcritical bifurcation from $E_0$ at $h_{TC}=c$. □

3.2. Saddle-Node Bifurcation

It follows from Lemma 1, controlling the parameter h, the coexisting equilibria $E_1$ and $E_2$ collide with each other as h crosses the critical value $h=h^{*}=\frac{{\left(1+c\right)}^2}{2}$ when $\Delta \left(h\right)=0$, and become an unique equilibrium $E_{SN}=\left(\frac{1-c}{2},0\right)$. Then as the value of parameter h changes, the unique equilibrium $E_{SN}$ disappears when $\Delta \left(h\right)<0$. Change in number of equilibria is owing to the occurence of saddle node bifurcation at $h_{SN}=h^{*}=\frac{{\left(1+c\right)}^2}{4}$. Thus, we state the following theorem.

Theorem 9. When $c<1$ and $e\left(1-m\right)\left(1-c\right)-2d\ne 0$, model (3) undergoes a saddle-node bifurcation from $E_{SN}=\left(\frac{1-c}{2}\mathrm{,0}\right)$ at $h_{SN}=\frac{{\left(1+c\right)}^2}{4}$.

Proof. We utilize lemma 3 to verify the transversality condition for the occurence of saddle-node bifurcation at $h=h_{SN}$. It can easily verify that $Det\left[J\left(E_{SN}\right)\right]=0$, which implies that the matrix $J\left(E_{SN}\right)$ has a zero eigenvalue. Let v and w be the eigenvectors corresponding to the zero eigenvalue for $J\left(E_{SN}\right)$ and $J^{\text{T}}\left(E_{SN}\right)$, respectively. Then we obtain

$v=\left(\begin{array}{c}1\\ 0\end{array}\right)\text{and}w=\left(\begin{array}{c}1\\ \frac{\left(1-c\right)\left[\left(1+c\right)f+2\left(1-m\right)\right]}{2\left[e\left(1-m\right)\left(1-c\right)-2d\right]}\end{array}\right).$

Furthermore, using the expression for F in Theorem 8, we can get

$F_h\left(E_{SN}\mathrm{,}{h}_{SN}\right)=\left(\begin{array}{c}-\frac{1-c}{1+c}\\ 0\end{array}\right)\mathrm{,}$

$D^2{F}_h\left(E_{SN},h_{SN}\right)\left(v,v\right)=\left(\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{V}_1^2+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{V}_1{V}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{V}_2^2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{V}_1^2+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{V}_1{V}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{V}_2^2\end{array}\right)=\left(\begin{array}{c}\frac{2\left(c-1\right)}{c+1}\\ 0\end{array}\right).$

It follows that

$W^{\text{T}}{F}_h\left(E_{SN}\mathrm{,}{h}_{SN}\right)=-\frac{1-c}{1+c}\ne \mathrm{0,}$

$W^{\text{T}}\left[D^{2}F\left(E_{SN}\mathrm{,}{h}_{SN}\right)\left(v\mathrm{,}v\right)\right]=-\frac{2\left(c-1\right)}{c+1}\ne 0.$

Therefore, we can conclude that model (3) undergoes a saddle-node bifurcation at $h_{SN}=\frac{{\left(1+c\right)}^2}{2}$. □

3.3. Hopf Bifurcation

From Theorem 5, the steady state of $E^{\mathrm{<mo>\; *\; </mo>}}$ changes as the parameter h crosses the threshold value $h=h_{Hp}$, which implies that the interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ may lose its stability through Hopf bifurcation under certain parametric restrictions. Let us adopt h as the bifurcation parameter, the Hopf bifurcation threshold is a positive root of ${\left[TrJ\left(E^{*}\right)\right]|}_{h=h_{Hp}}=0$, and can satisfy ${\left[DetJ\left(E^{*}\right)\right]|}_{h=h_{Hp}}>0$. Therefore we conclude the result in the following theorem.

Theorem 10. Suppose the conditions of existence of positive equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ stated in Lemma 2 is satisfied, then there is a Hopf bifurcation around $E^{\mathrm{<mo>\; *\; </mo>}}$ when $h=h_{Hp}$, where $h_{Hp}$ is defined in (12).

Proof. The characteristic equation of matrix $J\left(E^{*}\right)$ is ${\lambda }^2-tr\left(J_{E^{*}}\right)\lambda +det\left(J_{E^{*}}\right)=0$, and the conditions for the Hopf bifurcation occurs are stated below

1) ${\left[TrJ\left(E^{*}\right)\right]|}_{h=h_{Hp}}=0$,

2) ${\left[DetJ\left(E^{*}\right)\right]|}_{h=h_{Hp}}>0$,

3) ${\frac{\text{d}}{\text{d}h}\left[TrJ\left(E^{*}\right)\right]|}_{h=h_{Hp}}\ne 0$.

The conditions (1) and (2) have already been certified from Theorem 5, we merely need to verify the transversality condition (3). Obviously,

${\frac{\text{d}}{\text{d}h}\left[TrJ\left(E^{*}\right)\right]|}_{h=h_{Hp}}=x^{\mathrm{<mo>\; *\; </mo>}}{\left[\frac{1}{{\left(c+x^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}+\frac{f}{{\left(1+f{y}^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}\frac{\text{d}{y}^{\mathrm{<mo>\; *\; </mo>}}}{\text{d}h}\right]|}_{h=h_{Hp}}\ne 0.$

Consequently, condition (3) is satisfied, which guarantees the occurrence of Hopf bifurcation around $E^{\mathrm{<mo>\; *\; </mo>}}$ at $h=h_{Hp}$.

Furthermore, in order to investigate the stability (direction) of the limit cycle, we are going to calculate the first Lyapunov number $\sigma$ at the equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ of model (3).

We firstly translate $E^{\mathrm{<mo>\; *\; </mo>}}$ to the origin by using the transformation $X=x-x^{*}$, $Y=y-y^{*}$, then model (3) can reduce to

$(\begin{array}{l}\frac{\text{d}X}{\text{d}t}=a_{10}X+a_{01}Y+a_{20}{X}^2+a_{11}XY+a_{02}{Y}^2+a_{30}{X}^3\\ +a_{21}{X}^{2}Y+a_{12}X{Y}^2+a_{03}{Y}^3+P_4\left(X,Y\right),\\ \frac{\text{d}Y}{\text{d}t}=b_{10}X+b_{11}XY,\end{array}$ (13)

where

$a_{01}=-x^{*}\left[\left(1-x^{*}\right)\frac{f}{{\left(1+f{y}^{*}\right)}^2}+1-m\right],$

$a_{20}=\frac{h}{{\left(c+x^{*}\right)}^2}-\frac{1}{1+f{y}^{*}}-\frac{h{x}^{*}}{{\left(c+x^{*}\right)}^3},$

$a_{11}=\frac{f{x}^{*}}{{\left(1+f{y}^{*}\right)}^2}-\frac{f\left(1-x^{*}\right)}{{\left(1+f{y}^{*}\right)}^2}+m-1,a_{02}=\frac{f{x}^{*}\left(1-x^{*}\right)}{{\left(1+f{y}^{*}\right)}^3},$

$a_{30}=-\frac{hc}{{\left(c+x^{*}\right)}^4},a_{21}=\frac{2f}{{\left(1+f{y}^{*}\right)}^2},$

$a_{12}=\frac{2f\left(1-2x^{*}\right)}{{\left(1+f{y}^{*}\right)}^3},a_{03}=\frac{f{x}^{*}\left(1-x^{*}\right)}{{\left(1+f{y}^{*}\right)}^4},$

$b_{10}=e\left(1-m\right)y^{*},b_{11}=e\left(1-m\right),$

$a_{10}=b_{01}=b_{20}=b_{02}=b_{30}=b_{21}=b_{12}=b_{03}=0,$

and $P_4\left(X\mathrm{,}Y\right)$ is power series in $\left(X\mathrm{,}Y\right)$ with term $X^i{Y}^j$ satisfying $i+y\ge 4$. Let ${\Delta }_1=a_{10}{b}_{01}-a_{01}{b}_{10}$. Then the expression of the first Lyapunov number $\sigma$ can be expressed as follows:

$\begin{array}{c}\sigma =\frac{-3\Pi }{2a_{01}{\Delta }^{\frac{3}{2}}}\{[a_{10}{b}_{10}\left(a_{11}^2+a_{11}{b}_{02}+a_{02}{b}_{11}\right)+a_{10}{a}_{01}\left(b_{11}^2+a_{20}{b}_{11}+a_{11}{b}_{02}\right)\\ +b_{10}^2\left(a_{11}{a}_{02}+2a_{02}{b}_{02}\right)-2a_{10}{b}_{10}\left(b_{02}^2-a_{20}{a}_{02}\right)-2a_{10}{a}_{01}\left(a_{20}^2-b_{20}{b}_{02}\right)\\ -a_{01}^2\left(2a_{20}{b}_{20}+b_{11}{b}_{20}\right)+\left(a_{01}{b}_{10}-2a_{10}^2\right)\left(b_{11}{b}_{02}-a_{11}{a}_{20}\right)]\\ -\left(a_{10}^2+a_{01}{b}_{10}\right)\left[3\left(b_{10}{b}_{03}-a_{01}{a}_{30}\right)+2a_{10}\left(a_{21}+b_{12}\right)+\left(b_{10}{a}_{12}-a_{01}{b}_{21}\right)\right]\}\end{array}$

In fact, if $\sigma >0$, the equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ loses its stability through a Hopf bifurcation; if $\sigma <0$, $E^{\mathrm{<mo>\; *\; </mo>}}$ will obtain stability. Since the expression of $\sigma$ is quiet complicated, we cannot differentiate the sign of $\sigma$. Thus, we have presented numerical example in Section 5. □

4. The Effects of the Fear Factor and the Prey Refuge

In this section, we will explore the influence of fear factor (measured by f) and prey refuge (measured by m) on population density of model (3).

4.1. The Impact of the Fear on Population Density

Using the expression of $x^{\mathrm{<mo>\; *\; </mo>}}$ in Equation (6), we note that the density of prey at coexistence equilibrium is independent of f, so the fear effect cannot affect the prey density. However, we find that the prey growth rate is greatly influenced by the fear effect. On the other hand, differentiation of $y^{\mathrm{<mo>\; *\; </mo>}}$ gives

$\frac{\text{d}{y}^{*}}{\text{d}f}=\frac{\left(a+hf\right)h-2a\left(H+\Delta \right)}{4a^{2}f\sqrt{\Delta }}+\frac{1}{2f^2}<0,$

and $\underset{f\to \infty }{\lim }\frac{\text{d}{y}^{\mathrm{<mo>\; *\; </mo>}}}{\text{d}f}=0$, which indicates that $y^{\mathrm{<mo>\; *\; </mo>}}$ is strictly decreasing function with respect to f, i.e., increasing the level of f can decrease the final size of predator.

4.2. The Impact of the Prey Refuge on Population Density

We firstly denote that

$m^{*}=1-\frac{d\left(1-c+\sqrt{\Delta }\right)}{2\left(h-c\right)e}<1$

is the unique solution of the equation $y^{*}=0$ in $m\in \left(\mathrm{0,1}\right)$.

Using the expression in Equation (6) again, we can compute the derivative of $x^{\mathrm{<mo>\; *\; </mo>}}$ with respect to m is

$\frac{\text{d}{x}^{*}}{\text{d}m}=\frac{d}{{\left[e\left(1-m\right)\right]}^2}>0,$

which implies that $x^{\mathrm{<mo>\; *\; </mo>}}$ is a strictly increasing function of $m\in \left(\mathrm{0,}{m}^{\mathrm{<mo>\; *\; </mo>}}\right)$. That is, increasing the value of the prey refuge can rise the prey density.

Similarly, by simple computation, we can obtain

$\frac{\text{d}{y}^{\mathrm{<mo>\; *\; </mo>}}}{\text{d}m}=\frac{g\left(m\right)}{2a^{2}f}\mathrm{,}$

where

$g\left(m\right)=a\left[-{a^{\prime }}_m+\frac{\left(a+hf\right){a^{\prime }}_m-2f\left(H{a^{\prime }}_m+a{H^{\prime }}_m\right)}{\sqrt{\Delta }}\right]-{a^{\prime }}_m\left[-\left(a+hf\right)+\sqrt{\Delta }\right].$

If the equation $g\left(m\right)=0$ has a suitable solution $m_{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,1}\right)$, then it is the unique extreme point of $y^{\mathrm{<mo>\; *\; </mo>}}$ in $\left(\mathrm{0,1}\right)$. Thus, $\frac{\text{d}{y}^{*}}{\text{d}m}>0$ for all $m\in \left(\mathrm{0,}{m}_{\mathrm{<mo>\; *\; </mo>}}\right)$, i.e., strictly increasing; $\frac{\text{d}{y}^{*}}{\text{d}m}<0$ for all $m\in \left(m_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{m}^{\mathrm{<mo>\; *\; </mo>}}\right)$, i.e., strictly decreasing, and $y^{\mathrm{<mo>\; *\; </mo>}}$ reach its maximum at $m=m_{*}$.

5. Simulation Analysis

As we all know, the relationship between prey and predator in the natural environment is mutual restriction and influence. To better comprehend the dynamic relationship between prey and predator, we will perform numerical simulations to exhibit some complex dynamic behaviors of model (3). For convenience, we fix the parameters value as follows:

$f=0.2,c=0.3,e=0.6,d=0.1,m=0.2$ (14)

In order to demonstrate how the nonlinear Michaelis-Menten type prey harvesting affect the bifurcation dynamics of model (3), we have constructed a bifurcation diagram of model (3) in h-x plane by directly continuing the threshold values of $h_{Hp}=0.2100328877$, $h_{TC}=0.3$ and $h_{SN}=0.4225$ with parameters fixed as in (14) (seeing Figure 1(a)), which corresponds to Hopf bifurcation, transcritical bifurcation and saddle-node bifurcation, respectively. Meanwhile, Figure 1(b) is local amplification of (a) for $h\in \left[\mathrm{0.2,0.22}\right]$. Evidently, we note that model (3) can reveal abundant bifurcation dynamics with the value of key parameter h changing in Figure 1(a). Firstly, the interior equilibrium point $E^{*}=\left(x^{*},y^{*}\right)$ is stable when $0.2=h<h_{Hp}=0.2100328877$, and the interior equilibrium point $E^{*}=\left(x^{*},y^{*}\right)$ can change from stable to unstable when the value of h increases and passes through the critical value $h_{Hp}$ with the first Lyapunov number $\sigma =-35.44023532<0$, which implies that supercritical Hopf bifurcation has taken place in model (3). The simulation result has been proved by the phase diagram in Figure 2. Furthermore, it is easy to find that when $h=0.2$ is smaller than $h_{Hp}=0.2100328877$, the interior equilibrium point $E^{*}=\left(0.20833,0.37889\right)$ is locally asymptotically stable, which suggests that the prey and the predator can coexist and form the stable equilibrium state. However, We note that when $h>h_{Hp}$, the interior equilibrium point $E^{*}=\left(0.20833,0.37889\right)$ will lose its stability, a stable limit cycle will be generated, and the periodic orbit is surrounded by a trajectory passer through the

Figure 1. (a) Three vertical dot-dashed lines from left to right represent the critical value of transcritical, Hopf and saddle-node bifurcation, respectively. Two horizontal lines stand for the origin $E_0$ (purple) and the interior equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$ (blue). Two boundary equilibria $E_1$ and $E_2$ are presented by black and red curve, respectively. The dotted line indicates unstable and the solid line implies stable. (b) Local amplification of (a) for $h\in \left[\mathrm{0.2,0.22}\right]$.

Figure 2. Hopf bifurcation of model (3) around $E^{*}=\left(0.20834,0.28788\right)$. (a) Stable periodic orbits bifurcate from Hopf bifurcation around $E^{\mathrm{<mo>\; *\; </mo>}}$ when $h=h_{Hp}=0.2100328877$. (b) Local amplification of (a) for $\left(x\mathrm{,}y\right)\in \left[\mathrm{0.16,0.26}\right]\times \left[\mathrm{0.25,0.32}\right]$. (c) $E^{*}=\left(x^{*},y^{*}\right)$ is local asymptotically stable with $0.2=h<h_{Hp}$. (d) $E^{*}=\left(x^{*},y^{*}\right)$ is a spiral source with $0.22=h>h_{Hp}$.

saddle $E_2=\left(x_2,0\right)$, which suggests that the prey and the predator can coexist and produce periodic oscillation. Secondly, we can see that when the value of $h=0.25$ is less than $h_{TC}=0.3$, the origin $E_0=\left(0,0\right)$ is a saddle. In the meantime, the interior equilibrium point $E^{*}=\left(x^{*},y^{*}\right)$ is a spiral source and the boundary equilibrium point $E_2$ is unstable. However, when the value of h acrosses the threshold $h=h_{TC}=0.3$, the trivial equilibrium point $E_0$ becomes a nodal sink and model (3) will generate the boundary equilibrium point $E_1=\left(x_1,0\right)$, which is a saddle point. At the same time, it is worth noting that $E_1$ is unstable whenever it exists, and the boundary equilibrium point $E_2$ is also a saddle point. This implies that model (3) has occurred a transcritical bifurcation. The simulation result can be exhibited in Figure 3. Finally, if it continues to rise the critical thresholds h, we firstly see that the interior equilibrium point $E^{\mathrm{<mo>\; *\; </mo>}}$ will vanish. Then the two boundary equilibria $E_1$ and $E_2$ will gradually be close to each other, and becomes an unique boundary equilibrium point $E_{SN}=\left(x_{SN},0\right)$ until the annihilation, which is owing to the occurrence of saddle-node bifurcation of model (3) at $h=h_{SN}=0.4225$. The simulation conclusion can be demonstrated by Figure 4. Moreover, it can find from Figure 4 that when $h=0.4$ is less than $h_{SN}=0.4225$, model (3) has two non-trivial boundary equilibria $E_1$ and $E_2$, in which $E_1$ is a nodal source and $E_2$ is a saddle. However, with the

Figure 3. (a) The trivial equilibrium $E_0$ and a boundary equilibrium $E_2$ when $h<h_{TC}$ ; (b) The vector field graph for $h=h_{TC}=0.3$ with the origin $E_0$ collide with $E_1$ ; (c) The boundary equilibrium $E_1$ separate from the origin $E_0$ when $h>h_{TC}$, which is a saddle.

Figure 4. (a) Two boundary equilibria exist, respectively, and $E_1$ is a nodal source and $E_2$ is a saddle, and a trivial equilibrium $E_0$ which is a nodal sink when $h<h_{SN}$ ; (b) Due to the two boundary equilibria $E_1$ and $E_2$ coincide, we get a unique boundary equilibrium $E_{SN}$, which is a saddle-node when $h=h_{SN}$ ; (c) The two boundary equilibria disappear and only the origin $E_0$ when $h>h_{SN}$.

Figure 5. (a) Relationship between the fear and prey growth rate; (b) Relationship between the fear and the final size of the predator.

Figure 6. Relationship between population density with prey refuge: (a) prey and (b) predator.

increase of the thresholds h, there only the trivial equilibrium point $E_0$ when $h=h_{SN}=0.45$, which is a nodal sink.

In addition, in order to deeply investigate their impact mechanism of fear factor (measured by f), we have given the diagram of model (3) in Figure 5. In contrast to the model without fear factor, the increment of f will decrease the prey growth rate (seeing Figure 5(a)) and lessen the maximum of final size of predator density (seeing Figure 5(b)). Besides, the diagram of model (3) in Figure 6 manifests that the population density will change due to the variation of prey refuge (measure by m). Moreover, Figure 6(a) shows that the prey density will rise as the value of m increases, which implies that the increase of prey refuge can protect more prey from predation, Figure 6(b) exhibits that with the increase of m value, the predator density increases first, reaches the maximum value, and then decreases.

6. Conclusion

In this paper, we have studied the dynamics of a predator-prey model with nonlinear Michaelis-Menten type prey harvesting. The main focus of this article is to explore the impact of a Michaelis-Menten type prey harvesting mathematically and numerically. Using an appropriate conversion, the original Michaelis-Menten type prey harvesting term becomes a nonlinear harvesting term with only two parameters. Based on relevant mathematical theory, we reveal that nonlinear prey harvesting term plays an important role in influencing the dynamics of model (3). Note that, the parameters h and c in nonlinear harvesting term can affect the number and stability of boundary equilibria, which can become a saddle point, stable node, unstable node or saddle node for different parameter values. This implies that some possible bifurcation dynamic behaviors will occur, such as saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation. At the same time, comprehensive numerical simulation works have been carried out, which also, in turn, demonstrate the validity of these theoretical results. Moreover, we also discuss the influence of fear factor and prey refuge on the population density, it is easy to obtain that the fear factor can not only reduce the predator density but also affect the prey growth rate, while the prey refuge can affect both prey and predator population density. We hope this type of investigations will be of great help in comprehending the dynamic complexity of ecological system or physical systems in the future when the nonlinear Michaelis-Menten type prey harvesting can interact with prey populations.

Acknowledgements

This work is supported by the National Key Research and Development Program of China (Grant No. 2018YFE103700), the National Natural Science Foundation of China (Grants No. 61901303 and No. 61871293).

Conflicts of Interest

The author declares that there is no conflict of interests regarding the publication of this paper.

References

[1]	Wang, S.T., Yu, H.G., Dai, C.J. and Zhao, M. (2020) Global Asymptotic Stability and Hopf Bifurcation in a Homogeneous Diffusive Predator-Prey System with Holling Type II Functional Response. Applied Mathematics, 11, 389-406. https://doi.org/10.4236/am.2020.115028
[2]	Liu, H.Y., Yu, H.G., Dai, C.J., Ma, Z.L., Wang, Q. and Zhao, M. (2021) Dynamical Analysis of an Aquatic Amensalism Model with Non-Selective Harvesting and Allee Effect. Mathematical Biosciences and Engineering, 18, 8857-8882. https://doi.org/10.3934/mbe.2021437
[3]	Du, Y.H. and Hsu, S.B. (2004) A Diffusive Predator-Prey Model in Heterogeneous Environment. Journal of Differential Equations, 203, 331-364. https://doi.org/10.1016/j.jde.2004.05.010
[4]	Abbakar, K., Yang, Y., Mohamed, A., Xia, S., Abkar, M. and Hassan, O. (2021) Stability Analysis and Hopf Bifurcation for ODE System of Predator-Prey Model with Mutual Interference. Applied Mathematics, 12, 793-802. https://doi.org/10.4236/am.2021.129053
[5]	Liu, H., Dai, C.J., Yu, H.G., Guo, Q., Li, J.B., Hao, A.M., Kikuchi, J. and Zhao, M. (2021) Dynamics Induced by Environmental Stochasticity in a Phytoplankton-Zooplankton System with Toxic Phytoplankton. Mathematical Biosciences and Engineering, 18, 4101-4126. https://doi.org/10.3934/mbe.2021206
[6]	Li, X.X., Yu, H.G., Dai, C.J., Ma, Z.L., Wang, Q. and Zhao, M. (2021) Bifurcation Analysis of a New Aquatic Ecological Model with Aggregation Effect. Mathematics and Computers in Simulation, 190, 75-96. https://doi.org/10.1016/j.matcom.2021.05.015
[7]	Huang, Y.J., Chen, F.D. and Zhong, L. (2006) Stability Analysis of a Prey-Predator Model with Holling Type III Response Function Incorporating a Prey Refuge. Applied Mathematics and Computation, 182, 672-683. https://doi.org/10.1016/j.amc.2006.04.030
[8]	Balram, D., Kumar, A. and Maiti, A.P. (2019) Global Stability and Hopf-Bifurcation of Prey-Predator System with Two Discrete Delays Including Habitat Complexity and Prey Refuge. Communications in Nonlinear Science and Numerical Simulation, 67, 528-554. https://doi.org/10.1016/j.cnsns.2018.07.019
[9]	Kar, T.K. (2005) Stability Analysis of a Prey-Predator Model Incorporating a Prey Refuge. Communications in Nonlinear Science and Numerical Simulation, 10, 681-691. https://doi.org/10.1016/j.cnsns.2003.08.006
[10]	Chen, L.J., Chen, F.D. and Chen, L.J. (2010) Qualitative Analysis of a Predator-Prey Model with Holling Type II Functional Response Incorporating a Constant Prey Refuge. Nonlinear Analysis Real World Applications, 11, 246-252. https://doi.org/10.1016/j.nonrwa.2008.10.056
[11]	Eduardo, G.O. and Rodrigo, R.J. (2003) Dynamic Consequences of Prey Refuges in a Simple Model System: More Prey, Fewer Predators and Enhanced Stability. Ecological Modelling, 166, 135-146. https://doi.org/10.1016/S0304-3800(03)00131-5
[12]	Bairagi, N. and Chakraborty, B. (2018) Complexity in a Prey-Predator Model with Prey Refuge and Diffusion. Ecological Complexity, 37, 11-23. https://doi.org/10.1016/j.ecocom.2018.10.004
[13]	Zanette, L.Y., White, A.F., Allen, M.C. and Clinchy, M. (2011) Perceived Predation Risk Reduces the Number of Offspring Songbirds Produce per Year. Science, 334, 1398-1401. https://doi.org/10.1126/science.1210908
[14]	Cresswell, W. (2011) Predation in Bird Populations. Journal of Ornithology, 152, 251-263. https://doi.org/10.1007/s10336-010-0638-1
[15]	Pettorelli, N., Coulson, T., Durant, S.M. and Gaillard, J.M. (2011) Predation, Individual Variability and Vertebrate Population Dynamics. Oecologia, 167, 305-314. https://doi.org/10.1007/s00442-011-2069-y
[16]	Dai, B.X. and Sun, G.X. (2020) Turing-Hopf Bifurcation of a Delayed Diffusive Predator-Prey System with Chemotaxis and Fear Effect. Applied Mathematics Letters, 111, Article ID: 106644. https://doi.org/10.1016/j.aml.2020.106644
[17]	Tiwari, V., Tripathi, J.P., Mishra, S. and Upadhyay, R.K. (2020) Modeling the Fear Effect and Stability of Non-Equilibrium Patterns in Mutually Interfering Predator-Prey Systems. Applied Mathematics and Computation, 371, Article No. 124948. https://doi.org/10.1016/j.amc.2019.124948
[18]	Pal, S., Pal, N., Samanta, S. and Chattopadhy, J. (2019) Effect of Hunting Cooperation and Fear in a Predator-Prey Model. Ecological Complexity, 39, 1-18. https://doi.org/10.1016/j.ecocom.2019.100770
[19]	Wang, X.Y., Zanette, L. and Zou, X.F. (2016) Modelling the Fear Effect in Predator-Prey Interactions. Journal of Mathematical Biology, 73, 1179-1204. https://doi.org/10.1007/s00285-016-0989-1
[20]	Sasmal, S.K. (2018) Population Dynamics with Multiple Allee Effects Induced by Fear Factors—A Mathematical Study on Prey-Predator Interactions. Applied Mathematical Modelling, 64, 1-14. https://doi.org/10.1016/j.apm.2018.07.021
[21]	Chen, J., Huang, J.C., Ruan, S.G. and Wang, J.H. (2013) Bifurcations of Invariant Tori in Predator-Prey Models with Seasonal Prey Harvesting. SIAM Journal on Applied Mathematics, 73, 1876-1905. https://doi.org/10.1137/120895858
[22]	Chakraborty, S., Pal, S. and Bairagi, N. (2012) Predator-Prey Interaction with Harvesting: Mathematical Study with Biological Ramifications. Applied Mathematical Modelling, 36, 4044-4059. https://doi.org/10.1016/j.apm.2011.11.029
[23]	Zhang, X.B. and Zhang, H.Y. (2020) Global Stability of a Diffusive Predator-Prey Model with Discontinuous Harvesting Policy. Applied Mathematics Letters, 109, Article ID: 106539. https://doi.org/10.1016/j.aml.2020.106539
[24]	Sen, M., Srinivasu, P.D.N. and Banerjee, M. (2015) Global Dynamics of an Additional Food Provided Predator-Prey System with Constant Harvest in Predators. Applied Mathematics and Computation, 250, 193-211. https://doi.org/10.1016/j.amc.2014.10.085
[25]	May, R.M., Beddington, J.R., Clark, C.W. and Holt, S.J. (1979) Management of Multispecies Fisheries. Science, 205, 267-277. https://doi.org/10.1126/science.205.4403.267
[26]	Gupta, R.P., Banerjee, M. and Chandra, P. (2012) Bifurcation Analysis and Control of Leslie-Gower Predator-Prey Model with Michaelis-Menten Type Prey-Harvesting. Differential Equations and Dynamical Systems, 20, 339-366. https://doi.org/10.1007/s12591-012-0142-6
[27]	Huang, J.C., Gong, Y.J. and Ruan, S.G. (2013) Bifurcation Analysis and Control of Leslie-Gower Predator-Prey Model with Michaelis-Menten Type Prey-Harvesting. Discrete and Continuous Dynamical Systems: Series B, 18, 2101-2121. https://doi.org/10.3934/dcdsb.2013.18.2101
[28]	Huang, J.C., Liu, S.H., Ruan, S.G. and Zhang, X.N. (2016) Bogdanov-Takens Bifurcation of Codimension 3 in a Predator-Prey Model with Constant-Yield Predator Harvesting. American Institute of Mathematical Sciences, 15, 1041-1055. https://doi.org/10.3934/cpaa.2016.15.1041
[29]	Yao, Y. (2020) Bifurcations of a Leslie-Gower Prey-Predator System with Ratio-Dependent Holling IV Functional Response and Prey Harvesting. Mathematical Methods in the Applied Sciences, 43, 2137-2170. https://doi.org/10.1002/mma.5944
[30]	Lajmiri, Z., Ghaziani, R.K. and Orak, I. (2018) Bifurcation and Stability Analysis of a Ratio-Dependent Predator-Prey Model with Predator Harvesting Rate. Chaos Solitons and Fractals, 106, 193-200. https://doi.org/10.1016/j.chaos.2017.10.023
[31]	Gao, X.Y., Ishag, S., Fu, S.M., Li, W.J. and Wang, W.M. (2020) Bifurcation and Turing Pattern Formation in a Diffusive Ratio-Dependent Predator-Prey Model with Predator Harvesting. Nonlinear Analysis: Real World Applications, 51, Article ID: 102962. https://doi.org/10.1016/j.nonrwa.2019.102962
[32]	Lenzini, P. and Rebaza, J. (2010) Non-Constant Predator Harvesting on Ratio-Dependent Predator-Prey Models. Applied Mathematical Sciences, 4, 791-803.
[33]	Singh, M.K., Bhadauria, B.S. and Singh, B.K. (2016) Qualitative Analysis of a Leslie-Gower Predator-Prey System with Nonlinear Harvesting in Predator. International Journal of Engineering Mathematics, 2016, Article ID: 2741891. https://doi.org/10.1155/2016/2741891
[34]	Liu, B., Wu, R.C. and Chen, L.P. (2018) Patterns Induced by Super Cross-Diffusion in a Predator-Prey System with Michaelis-Menten Type Harvesting. Mathematical Biosciences, 298, 71-79. https://doi.org/10.1016/j.mbs.2018.02.002
[35]	Chen, M.X., Wu, R.C., Liu, B. and Chen, L.P. (2018) Pattern Selection in a Predator-Prey Model with Michaelis-Menten Type Nonlinear Predator Harvesting. Ecological Complexity, 36, 239-249. https://doi.org/10.1016/j.ecocom.2018.09.004
[36]	Liu, J. and Zhang, L. (2016) Bifurcation Analysis in a Prey-Predator Model with Nonlinear Predator Harvesting. Journal of the Franklin Institute, 353, 4701-4714. https://doi.org/10.1016/j.jfranklin.2016.09.005
[37]	Liu, W. and Jiang, Y.L. (2017) Bifurcation of a Delayed Gause Predator-Prey Model with Michaelis-Menten Type Harvesting. Journal of Theoretical Biology, 438, 116-132. https://doi.org/10.1016/j.jtbi.2017.11.007
[38]	Zhang, Y., Chen, S.H. and Wei, X. (2017) Stochastic Periodic Solution for a Perturbed Non-Autonomous Predator-Prey Model with Generalized Nonlinear Harvesting and Impulses. Physica A: Statistical Mechanics and Its Applications, 486, 347-366. https://doi.org/10.1016/j.physa.2017.05.058
[39]	Zhang, Z.F. (1992) Qualitative Theory of Differential Equation. Science Press, Beijing.
[40]	Meiss, J.D. (2007) Differential Dynamical Systems. Society for Industrial and Applied Mathematics, Philadelphia. https://doi.org/10.1137/1.9780898718232
[41]	Perko, L. (1996) Differential Equations and Dynamical Systems. Springer, New York. https://doi.org/10.1007/978-1-4684-0249-0