Bifurcation analysis of Leslie-Gower predator-prey system with harvesting and fear effect

1.   Introduction

Since mathematician Lotka ^[1] and Volterra ^[2] first proposed the first predator-prey system to explain the relationship between two populations, the predator-prey system have been widely used by many biologists and mathematicians to describe their dynamic behaviors in populations over time. Many scholars have identified an increasing number of complex dynamical properties in the predator-prey system, such as global stability ^[3], a variety of bifurcation types: degenerate Hopf bifurcation ^[4], Bogdanov-Takens bifurcation of codimension 3, which includes cusp type ^[5], focus type ^[6], saddle type ^[7] and elliptic type ^[8].

The population growth function and functional response are two important factors in a predator-prey system. The Leslie-Gower system ^[9,10] is an improvement of the Lotka-Volterra predator-prey system, which has the general form:

where $P(x)$ is functional response. Obviously, diverse $P(x)$ have various effects on the prey-predator system ^[11,12,13]. Especially, many scholars have proposed Leslie-type predator-prey system, which included simplified Holling type Ⅳ functional response, i.e., $P(x) = \frac{mx}{a+x^2}$. Li and Xiao ^[14] found that the system (1.1) produces not only Hopf bifurcation but Bogdanov-Takens bifurcation of codimension 2 as well. Shang and Qiao ^[15] revealed that there exists degenerate Hopf bifurcation and a cusp of codimension at least 4 in the system (1.1) with strong Allee effect on prey. However, generalized Holling type Ⅳ functional response i.e., $P(x) = \frac{mx}{a+bx+x^2}$ is more practical, which can be used to describe the phenomenon of group defence in prey population and inhibition in predator population. Shan ^[16] studied the interaction produced by Holling type Ⅳ functional response and strong and weak Allee effects. Meanwhile, the system undergoes degenerate Hopf bifurcation of codimension 3.

Many experiments ^[24,25] have shown that predator do not necessarily affect prey population only by killing. Sometimes the existence of predator may affect the behavior and psychology of the prey, which in turn may change the number of prey. However, many scholars have only considered direct predation from predators, ignoring the fact that the presence of predators themselves can have an effect on the density of prey. Therefore, it is necessary to discuss the effect of fear brought to the prey by the predator. Wang et al.^[26] came up with a predator-prey system with fear effect. As the fear level gets higher, the system stabilizes. However, when the level of fear becomes low, the system generates multiple limit cycles. Since then, many scholars have been influenced by this idea to study predator-prey system with fear effect ^{[17,27,28,29,30,31,32,33,34,35,36]}. Zou ^[37] showed that the presence of the fear effect changes the system from a state of chaos to a state of stability. Wang et al.^[38] showed that the presence of the fear effect favors the occurrence of oscillation in the system. Wang et al.^[19] studied the properties of a modified Leslie-Gower system and showed that fear effect enriches the dynamics of the system.

At the same time, people tend to harvest populations of organisms in the predator-prey system for human needs, such as survival and economic interests. Thus, it is necessary to consider harvesting of the population in the system. Specifically, we need to capture predators for some specific purposes. Therefore, the discussion is more relevant for the harvesting of predators. In ^[18], authors considered the constant-yield harvesting in predator and showed that the system produces various types of bifurcation. In ^[12], authors analyzed a system with the nonlinear Michaelis-Menten type harvesting of predator.

In this article, we investigate the predator-prey system with the constant-effort harvesting on predator and fear effect on prey, i.e.,

where $q$ is the catchability co-efficient of the predator; $m(0 < m < 1)$ is the fraction of the stock available for harvesting; $E$ is the capture effort for harvesting; $r$ and $s$ indicate the intrinsic growth rate of prey and predator populations, respectively; $K$ indicates the environmental capacity of the prey; $h$ indicates the measure of the food quality of the prey for translating into predator births; $a$ represents the level of fear; $\alpha$ stands for maximum predation rate; $b$ and $c$ are the half-saturation constant and the functional response constant, respectively.

For simplicity, taking the following transformations:

and dropping the bars, then

where $k = \frac{arK^2}{\alpha}, d = \frac{c}{K}, e = \frac{b}{K^2}, \delta = \frac{sK}{\alpha h}, \beta = \frac{\alpha h}{rK}, \lambda = \frac{qmE}{r}$ are positive constants.

The article is divided into some parts: In the following section, we show the solutions are bounded and the origin of the system (1.3) is unstable. In Sections 3 and 4, the existence and stability of equilibria are discussed. In Section 5, various forms of bifurcations of the system (1.3) are studied. We study transcritical bifurcation, saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation of codimension 2. In Section 6, the numerical simulation experiments of the system (1.3) are presented in order to better confirm the correctness of the theoretical derivations and to demonstrate these dynamical properties. In the end, a short conclusion is stated in the last section.

2.   Boundedness of the solutions and stability at the origin

In this section, positive and bounded solutions of the system (1.3) are explored. Also, we explore the dynamical properties of the system (1.3) at the origin.

Theorem 2.1. The solutions of the system (1.3) are always positive with the initial conditions $(x(0), y(0))\in\Omega = \left\lbrace {(x, y)\in R^{2}\mid x > 0, y > 0}\right\rbrace$.

Proof. Clearly, the solutions of the system (1.3) can be written in the next form:

So due to the initial condition $(x(0), y(0))\in\Omega = \left\lbrace {(x, y)\in R^{2}\mid x > 0, y > 0}\right\rbrace$, all solutions of the system (1.3) must be positive. □

Theorem 2.2. If $\lambda < \delta \beta$, the solutions of the system (1.3) with positive initial values are bounded and are limited in the positive bounded set

Proof. Observing the first equation of the system (1.3), we get

If $x\ge1$, then there are $\dot{x}\le0$. If $x < 1$, then there are $\dot{x} < x(1-x)$, separating variables, we have $x < 1-\frac{1}{C_{1} e^{t}+1}$, where $C_{1}$ is a constant. Thus, we get $\mathop{lim}\limits_{t \to +\infty} supx(t)\le1$.

Observing the second equation of the system (1.3), we get

If $y\ge\beta-\frac{\lambda}{\delta}$, then there are $\dot{y}\le0$. If $y < \beta-\frac{\lambda}{\delta}$, then there are $\frac{dy}{y(\beta-\frac{\lambda}{\delta}-y)} < \delta ydt$, when $\lambda < \delta \beta$, we have $y < \beta-\frac{\lambda}{\delta}-\frac{\beta-\frac{\lambda}{\delta}}{C_{2} e^{\delta t}+1}$, where $C_{2}$ is a constant. Thus, we get $\mathop{lim}\limits_{t \to +\infty} supy(t)\le\beta-\frac{\lambda}{\delta}$.

Therefore, the system (1.3) is bounded and the bounded region is

□

Theorem 2.3. The origin (0, 0) in the system (1.3) is unstable.

Proof. Since the system (1.3) is not defined at the origin, we first introduce a new time variable $\tau$ by

We get the new system

which is topologically equivalent to the system (1.3).

Since the Jacobi matrixs of the system (2.1) at the origin is a null matrix, which means that the origin is non-hyperbolic, we analyze the stability at the origin by the blow-up method. For the system (2.1), if we let $x = 0$, then its second equation becomes $\frac{dy}{d\tau } = -\delta ey^2(1+ky) < 0$, which implies that the $y$-axis is an invariant line converging to origin. Making the following transformations $x = u, y = uv$ and $\tau = t/u$.

Then the system (2.1) becomes the following system

Bringing $u = 0$ into the second equation of the system (2.2), we get $v = 0$ or $v = \beta-\frac{\lambda}{\delta}$. Then we get two boundary equilibria $(0, 0)$ and $(0, \beta-\frac{\lambda}{\delta})$ from the system (2.2). In addition, the Jacobi matrix of the system (2.2) at the two equilibria are

and

Clearly, the eigenvalues of this matrix $J(0, 0)$ are $e > 0$ and $e(\delta \beta -\lambda -1)$, which means (0, 0) is unstable. When $\lambda < 1+\delta \beta$, the matrix $J(0, \beta-\frac{\lambda}{\delta})$ has eigenvalues $e > 0$ and $e(\lambda-\delta \beta -1) < 0$, which means $(0, \beta-\frac{\lambda}{\delta})$ is a saddle. When $\lambda > 1+\delta \beta$, the matrix $J(0, \beta-\frac{\lambda}{\delta})$ has eigenvalues $e > 0$ and $e(\lambda-\delta \beta -1) > 0$, which means $(0, \beta-\frac{\lambda}{\delta})$ is a unstable node. So the origin of the system (1.3) is an unstable point.□

3.   Existence of equilibria

Apparently, the system (1.3) always has a unique boundary equilibrium point $E_{1}(1, 0)$. The positive equilibria of the system (1.3) satisfy the next equations:

From the second equation of (3.1), we get $y = (\beta -\frac{\lambda}{\delta})x$. In order to satisfy the existence of the positive equilibria, $\beta-\frac{\lambda}{\delta} > 0$ is required to hold. Substituting $y = (\beta -\frac{\lambda}{\delta})x$ into the first equation of (3.1), we have

Let

$B = \left[ (\beta-\frac{\lambda}{\delta})^2k+d-1\right](\beta-\frac{\lambda}{\delta}-d+e)+9e$, \qquad $C = (\beta-\frac{\lambda}{\delta}-d+e)^2+3e\left[ (\beta-\frac{\lambda}{\delta})^2k+d-1\right].$

First we note that $f(x) = 0$ must have a positive root because $f(0) = -e < 0$ and $f(1) = (\beta-\frac{\lambda}{\delta})^2k+(\beta-\frac{\lambda}{\delta}) > 0$.

Theorem 3.1. For the number of positive equilibria, we have:

1) when $\Delta > 0$, then the system (1.3) has only one positive equilibrium point $E_{1}^{\ast} = (x_{1}^{\ast}, y_{1}^{\ast})$, where $x_{1}^{\ast}$ is a single positive root.

2) when $\Delta = 0$,

(i) $(\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ and $\beta-\frac{\lambda}{\delta}-d+e > 0$, then the system (1.3) has two positive equilibria $E_{2}^{\ast} = (x_{2}^{\ast}, y_{2}^{\ast})$ and $E_{+} = (x_{+}, y_{+})$, where $x_{2}^{\ast}$ is a single positive root and $x_{+}$ is a dual positive root.

(ii) $(\beta-\frac{\lambda}{\delta})^2k+d-1\ge0$ or $\beta-\frac{\lambda}{\delta}-d+e = 0$ or $(\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ and $\beta-\frac{\lambda}{\delta}-d+e < 0$, then the system (1.3) has only a positive equilibrium point $E_{3}^{\ast} = (x_{3}^{\ast}, y_{3}^{\ast})$, where $x_{3}^{\ast}$ is a single positive root.

(iii) $A = 0$ and $(\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$, then the system (1.3) has a positive equilibrium point $E_{++} = (x_{++}, y_{++})$, where $x_{++}$ is a triple positive root.

3) when $\Delta < 0$,

(i) $(\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ and $\beta-\frac{\lambda}{\delta}-d+e > 0$, then the system (1.3) has three unequal positive equilibria $E_{4}^{\ast} = (x_{4}^{\ast}, y_{4}^{\ast})$, $E_{5}^{\ast} = (x_{5}^{\ast}, y_{5}^{\ast})$ and $E_{6}^{\ast} = (x_{6}^{\ast}, y_{6}^{\ast})$, where $x_{4}^{\ast}$, $x_{5}^{\ast}$, $x_{6}^{\ast}$ are all single positive roots and $x_{4}^{\ast} < x_{5}^{\ast} < x_{6}^{\ast}$.

(ii) $(\beta-\frac{\lambda}{\delta})^2k+d-1\ge0$ or $\beta-\frac{\lambda}{\delta}-d+e = 0$ or $(\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ and $\beta-\frac{\lambda}{\delta}-d+e < 0$, then the system (1.3) has only a positive equilibrium point $E_{7}^{\ast} = (x_{7}^{\ast}, y_{7}^{\ast})$, where $x_{7}^{\ast}$ is a single positive root.

Theorem 3.2. For the case of the equilibria in Theorem 3.1, we have $g(x_{i}^{\ast}) > 0$, where $i = 1, 2, 3, 4, 6, 7$. We also have $g(x_{+}) = 0$, $g(x_{++}) = 0$ and $g(x_{5}^{\ast}) < 0$.

4.   Stability of equilibria

In this section, the stability of the equilibria is discussed.

The Jacobian matrix of the system (1.3) is:

Thus, we obtain the next theorem.

4.1.   Stability of the equilibrium point $E_{1}$

Theorem 4.1. The stability of boundary equilibrium point $E_{1}(1, 0)$ is discussed as follows:

1) If $\lambda < \delta\beta$, then $E_{1}$ is a hyperbolic saddle.

2) If $\lambda > \delta\beta$, then $E_{1}$ is a hyperbolic stable node.

3) If $\lambda = \delta\beta$, then $E_{1}$ is a attracting saddle-node. That is, a adequately small domain of $E_{1}$ is split into two parts by two dividing lines along the top and bottom of $E_{1}$ that tend to $E_{1}$. The left half-plane is a parabolic sector and the right half-plane is two hyperbolic sectors.

Proof. The Jacobian matrix of $E_{1}$ is:

Apparently, $J_{E_{1}}$ has two eigenvalues $\mu_{1} = -1$ and $\mu_{2} = \delta\beta-\lambda$. Hence, if $\mu_{2} = \delta\beta-\lambda > 0$, $i.e.$, $\lambda < \delta\beta$, then $E_{1}$ is a hyperbolic saddle; if $\mu_{2} = \delta\beta-\lambda < 0$, $i.e.$, $\lambda > \delta\beta$, then $E_{1}$ is a hyperbolic stable node; if $\mu_{2} = \delta\beta-\lambda = 0$, $i.e.$, $\lambda = \delta\beta$, then the two eigenvalues are $-1$ and $0$. In order to investigate the stability of $E_{1}$ when $\lambda = \delta\beta$, we convert $E_{1}$ to the origin by $(X, Y) = (x-1, y)$. Meanwhile, we expand the system (1.3) near the origin to the third order, then the system (1.3) becomes

where $a_{11} = k+\frac{1-e}{\left(d +e +1\right)^{2}}, a_{21} = k +\frac{\left(d +3\right) e -1}{\left(d +e +1\right)^{3}}$, $P_{1}(X, Y)$ and $Q_{1}(X, Y)$ are $C^\infty$ functions at least of order fourth in $(X, Y)$.

By the transformation

the system (4.1) becomes a standard form

where $b_{02} = \frac{(1+d+e)(a_{11}+\delta)+1}{(1+d+e)^2}, b_{11} = a_{11}+\frac{2}{1+d+e}, b_{03} = \frac{k^2(1+d+e)-\delta+a_{21}}{(1+d+e)^2}, b_{12} = \frac{k^2(1+d+e)-\delta+2a_{21}}{1+d+e}, b_{21} = a_{21}$, $P_{2}(x, y)$ and $Q_{2}(x, y)$ are $C^\infty$ functions at least of order fourth in $(x, y)$.

Introduce a new variable $\tau$ through $\tau = -t$, we have

where $P_{3}(x, y)$ and $Q_{3}(x, y)$ are $C^\infty$ functions at least of order fourth in $(x, y)$. Since the coefficient of $y^2$ for the second equation of the system $(4.3)$ is $-\delta < 0$, the equilibrium point $E_{1}(1, 0)$ is a attracting saddle-node from Theorem 7.1 in Chapter 2 in ^[20], which means a adequately small region of $E_{1}$ is split into two parts by two dividing lines along the top and bottom of $E_{1}$ that tend to $E_{1}$. The left half-plane is a parabolic sector and the right half-plane is two hyperbolic sectors.□

4.2.   Stability of the internal equilibria

In this section, we dicuss the stability of the internal equilibria by simplifying the Jacobi matrix. Obviously, we can find that the internal equilibria satisfy the Eq (3.1).

Hence, we have

where $x^{\ast}$ and $y^{\ast}$ are the horizontal and vertical coordinates of internal equilibria, respectively.

Then, we obtain the Jacobi matrix at the internal equilibria

So the determinant and trace of this matrix are

4.2.1.   Stability of the equilibria $E_{+}$ and $E_{++}$

Lemma 1(^[23]). The model

is equivalent to the following model

by certain nonsingular transformations in the domain of $(0, 0)$.

Theorem 4.2. If $\Delta = 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ and $\beta-\frac{\lambda}{\delta}-d+e > 0$, the system (1.3) has a positive equilibrium point $E_{+} = (x_{+}, y_{+})$, where $x_{+}$ is a dual positive root, then:

1) If $\lambda\neq\delta\beta+x_{+}\frac{3x_{+}^2+(2d-2)x_{+}+e-d}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}$ and $\lambda\neq\delta \beta -\delta\left[ \frac{(x_{+}-1)(-x_{+}^3+3ex_{+}+de)(1+ky_{+})^2+(x_{+}^{2}+dx_{+}+e)^2}{(x_{+}^{2}-e)(1+ky_{+})^3(x_{+}^{2}+dx_{+}+e)-k(x_{+}^{2}+dx_{+}+e)^2}\right]$, $E_{+}$ is a saddle-node.

2) If $\lambda = \delta\beta+x_{+}\frac{3x_{+}^2+(2d-2)x_{+}+e-d}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}$, $E_{+}$ is a cusp. Further, if $\lambda\neq \lambda_{1}$ and $h(\lambda)\neq 0$, then $E_{+}$ is a cusp with codimension 2; if $\lambda = \lambda_{1}$ or $h(\lambda) = 0$, then $E_{+}$ is a cusp with codimension at least 3, where

Proof. Firstly, we convert the equilibrium point $E_{+}$ into the origin by the translation $(X, Y) = (x-x_{+}, y-y_{+})$, then the system (1.3) becomes

where

and $P_{4}(X, Y)$, $Q_{4}(X, Y)$ are $C^\infty$ functions at least of order third in $(X, Y)$.

Case 1: $\lambda\neq\delta\beta+x_{+}\frac{3x_{+}^2+(2d-2)x_{+}+e-d}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}$.

Under this condition, the Jacobi matrix of the system (4.6) is as follows

so the eigenvalues of $J_{E_{+}}$ are $\mu_{1} = 0$, $\mu_{2} = a_{10}^{\prime}+b_{01}^{\prime}$. Next we reduce the system (4.6) to its standard form by following transformation

The system (4.6) becomes the following form

where

and $c = \beta-\frac{\lambda}{\delta}$, $P_{5}(x, y)$, $Q_{5}(x, y)$ are $C^\infty$ functions at least of order third in $(x, y)$.

We import a new variable $\tau$ to the system (4.7) by $\tau = (a_{10}^{\prime}+\lambda-\delta\beta)t$. For formal simplicity, we still use $t$ instead of $\tau$, we have

where $e_{ij}^{\prime} = \frac{c_{ij}^{\prime}}{a_{10}^{\prime}+\lambda-\delta\beta}$, $f_{ij}^{\prime} = \frac{d_{ij}^{\prime}}{a_{10}^{\prime}+\lambda-\delta\beta}$ $(i+j = 2)$, $P_{6}(x, y)$, $Q_{6}(x, y)$ are $C^\infty$ functions at least of order third in $(x, y)$. Hence, if the coefficient of $y^2$ is $f_{02}^{\prime}\neq0$ for the second equation of the system $(4.8)$, i.e., $d_{02}^{\prime}\neq 0$, $E_{+}$ is a saddle-node from Theorem 7.1 in Chapter 2 in ^[20].

By a simple calculation, we get

That is, if

$E_{+}$ is a saddle-node.

Case 2: $\lambda = \delta\beta+x_{+}\frac{3x_{+}^2+(2d-2)x_{+}+e-d}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}$.

Under this condition, the Jacobi matrix of the system (4.6) is as follows

so the eigenvalues of $J_{E_{+}}$ are $\mu_{1} = 0$, $\mu_{2} = 0$. Next we reduce the system (4.6) to its standard form by the following transformation

The system (4.6) becomes the following form again

where

and $c = \beta-\frac{\lambda}{\delta}$, $P_{7}(x, y)$, $Q_{7}(x, y)$ are $C^\infty$ functions at least of order third in $(x, y)$.

We import a new variable $\tau$ to the system (4.9) by $\tau = \delta t$, For formal simplicity, we still use $t$ instead of $\tau$, we have

where $e_{ij} = \frac{c_{ij}}{\delta}$, $f_{ij} = \frac{d_{ij}}{\delta}$ $(i+j = 3)$, $P_{8}(x, y)$, $Q_{8}(x, y)$ are $C^\infty$ functions at least of order third in $(x, y)$.

Then we use the method of Theorem $7.3$ in Chapter 2 in ^[20] to prove $E_{+}$ is a cusp. First we make the next definitions

If $e_{20}\neq0$, from $y+W(x, y) = 0$, we can obtain

then we have

Hence, $k = 2m, m = 1, a_{2m} = f_{20}, N = 1, B_{N} = 2e_{20}+f_{11}$, $E_{+}$ is a cusp.

From Lemma 1, the system (4.10) is equivalent to the following system

where $Q_{9}(x, y)$ are $C^\infty$ functions at least of order third in $(x, y)$.

By a simple calculation, we get

Furthermore, if $f_{20}\neq0$ and $f_{11}+2e_{20}\neq 0$, i.e., $\lambda\neq\lambda_{1}$ and $h(\lambda)\neq 0$, then $E_{+}$ is a cusp with codimension 2; if $f_{20} = 0$ or $f_{11}+2e_{20} = 0$, i.e., $\lambda = \lambda_{1}$ or $h(\lambda) = 0$, then $E_{+}$ is a cusp with codimension at least 3, where

□

The investigation of equilibrium point $E_{++}$ stability when $\Delta = 0, A = 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ is similar to the discussion of stability of $E_{+}$. Imitating the proof of Theorem 4.2, we have the next theorem.

Theorem 4.3. If $\Delta = 0, A = 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$, the system (1.3) has a positive equilibrium point $E_{++} = (x_{++}, y_{++})$, where $x_{++}$ is a triple positive root, then:

1) If $\lambda\neq\delta\beta+x_{++}\frac{3x_{++}^2+(2d-2)x_{++}+e-d}{(1+ky_{++})(x_{++}^2+dx_{++}+e)}$

and $\lambda\neq\delta \beta -\delta\left[ \frac{(x_{++}-1)(-x_{++}^3+3ex_{++}+de)(1+ky_{++})^2+(x_{++}^{2}+dx_{++}+e)^2}{(x_{++}^{2}-e)(1+ky_{++})^3(x_{++}^{2}+dx_{++}+e)-k(x_{++}^{2}+dx_{++}+e)^2}\right]$, $E_{++}$ is a saddle-node.

2) If $\lambda = \delta\beta+x_{++}\frac{3x_{++}^2+(2d-2)x_{++}+e-d}{(1+ky_{++})(x_{++}^2+dx_{++}+e)}$, $E_{++}$ is a cusp. Further, if $\lambda\neq \lambda_{2}$ and $w(\lambda)\neq 0$, then $E_{++}$ is a cusp with codimension 2; if $\lambda = \lambda_{2}$ or $w(\lambda) = 0$, then $E_{++}$ is a cusp with codimension at least 3, where

4.2.2.   Stability of the equilibrium point $E_{i}^{\ast}$

Suppose the existence conditions of $E_{i}^{\ast}$ ($i = 1, 2, 3, 4, 5, 6, 7$) are satisfied, then we can easily find that the values of $Det(J_{E_{5}^{\ast}})$ is always negative by Theorem 3.2. Therefore, $E_{5}^{\ast}$ is a saddle, regardless of the value of its trace. We can also easily find that the values of $Det(J_{E_{i}^{\ast}})$ ($i = 1, 2, 3, 4, 6, 7$) is always positive by Theorem 3.2. Hence, if $Tr(J_{E_{i}^{\ast}}) < 0$, then $E_{i}^{\ast}$ ($i = 1, 2, 3, 4, 6, 7$) is a sink; if $Tr(J_{E_{i}^{\ast}}) > 0$, then $E_{i}^{\ast}$ ($i = 1, 2, 3, 4, 6, 7$) is a source; if $Tr(J_{E_{i}^{\ast}}) = 0$, then $E_{i}^{\ast}$ ($i = 1, 2, 3, 4, 6, 7$) is a center.

Summarize the above analysis to the following theorem.

Theorem 4.4. If $\Delta < 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$ and $\beta-\frac{\lambda}{\delta}-d+e > 0$, the system (1.3) has three unequal positive equilibria $E_{4}^{\ast} = (x_{4}^{\ast}, y_{4}^{\ast})$, $E_{5}^{\ast} = (x_{5}^{\ast}, y_{5}^{\ast})$ and $E_{6}^{\ast} = (x_{6}^{\ast}, y_{6}^{\ast})$, where $x_{4}^{\ast}$, $x_{5}^{\ast}$, $x_{6}^{\ast}$ are all single positive roots and $x_{4}^{\ast} < x_{5}^{\ast} < x_{6}^{\ast}$. Moreover, $E_{5}^{\ast}$ is always a saddle.

Theorem 4.5. If the existence conditions of $E_{i}^{\ast}$ ($i = 1, 2, 3, 4, 6, 7$) are satisfied, the system (1.3) has a positive equilibrium point, i.e.$E_{i}^{\ast} = (x_{i}^{\ast}, y_{i}^{\ast})$, where $x_{i}^{\ast}$ is single positive root. Then $E_{i}^{\ast}$ may be a sink or a source or a center.

5.   Bifurcation analysis

In order to explore the influence of the harvesting term on the system (1.3), we choose $\lambda$ as the bifurcation control parameter. We explore the conditions for codimension one and codimension two bifurcations in this section, such as transcritical, saddle-node, Hopf and Bogdanov-Takens bifurcation.

5.1.   Transcritical bifurcation

Theorem 5.1. The system (1.3) undergoes a transcritical bifurcation when $\lambda = \lambda_{TC} = \beta\delta$.

Proof. If $\lambda = \lambda_{TC} = \beta\delta$, the Jacobian matrix at $E_{1}$ can be written in the following form:

Then we check whether the corresponding condition for transcritical bifurcation hold true by Sotomayor's theorem. Suppose the eigenvectors of $J_{E_{1}}$ and $J_{E_{1}}^T$ with respect to the zero eigenvalue are $V$ and $W$, they are respectively as follows

By a simple calculation, we have

     $F_{\lambda}(E_{1};\lambda_{TC}) = \left(\begin{array}{cc} 0 \\ -y\\ \end{array} \right) _{(E_{1};\lambda_{TC})} = \left(\begin{array}{cc} 0 \\ 0\\ \end{array} \right),$

     $DF_{\lambda}(E_{1};\lambda_{TC})V = \left(\begin{matrix} 0 & 0\\ 0 & -1\\ \end{matrix} \right) \left(\begin{matrix} \frac{1}{1+d+e}\\ -1\\ \end{matrix} \right) = \left(\begin{array}{cc} 0 \\ 1\\ \end{array} \right),$

     $D^{2}F(E_{1};\lambda_{TC})(V, V) = \left(\begin{array}{cc} \frac{\partial^{2}F_{1}}{\partial x^{2}}v_{1}v_{1}+2\frac{\partial^{2}F_{1}}{\partial x \partial y}v_{1}v_{2}+\frac{\partial^{2}F_{1}}{\partial y^{2}}v_{2}v_{2}\\ \frac{\partial^{2}F_{2}}{\partial x^{2}}v_{1}v_{1}+2\frac{\partial^{2}F_{2}}{\partial x \partial y}v_{1}v_{2}+\frac{\partial^{2}F_{2}}{\partial y^{2}}v_{2}v_{2} \\ \end{array}\right) _{(E_{1};\lambda_{TC})} = \left(\begin{array}{cc} -2\frac{k(1+d+e)^2+d+2}{(1+d+e)^3} \\ -2\delta\\ \end{array} \right).$

Therefore, we get

    $W^{T}F_{\lambda}(E_{1};\lambda_{TC}) = 0,$

    $W^{T}[DF_{\lambda}(E_{1};\lambda_{TC})V] = 1\neq0,$

     $W^{T}[D^{2}F(E_{1};\lambda_{TC})(V, V)] = -2\delta\neq0.$

Hence, if $\lambda = \lambda_{TC} = \beta\delta$, the system (1.3) will experience a transcritical bifurcation.□

In order to verify the existence of a transcritical bifurcation, the next numerical simulation is given in Figure 1. From Figure 1, we can find that the system (1.3) undergoes a transcritical bifurcation.

When $k = 0.2, d = 0.05, e = 0.024, \delta = 0.3, \beta = 2$, we can directly calculate $\lambda_{TC} = 0.6$. When the value of $\lambda$ is smaller than 0.6, the system (1.3) has two equilibria $E_{1}^{\ast}$ and $E_{1}$, the internal equilibrium point $E_{1}^{\ast}$ is locally asymptotically stable and the boundary equilibrium point $E_{1}$ is a unstable saddle (see Figure 1(a)). That is to say, the prey and predator populations will eventually coexist in a periodic oscillation. However, choosing other initial values, they will coexist in $E_{1}^{\ast}$. But when the value of $\lambda$ is equal to 0.6, the internal equilibrium point $E_{1}^{\ast}$ disappears and the boundary equilibrium point $E_{1}$ is a saddle-node (see Figure 1(b)). Further, when the value of $\lambda$ is greater to 0.6, the system (1.3) has only one boundary equilibrium point $E_{1}$, which is a globally asymptotically stable node (see Figure 1(c)). That is to say, the predator population will eventually become extinct.

5.2.   Saddle-node bifurcation

From theorem 3.1, we get a saddle-node bifurcation surface:

The bifurcation surface indicates that the system (1.3) will occur a saddle-node bifurcation and the number of positive equilibria ranges from one to three when the parameter $\lambda$ changes the sign of $\Delta$.

The next theorem will state the correctness of this result.

Theorem 5.2. Under the condition of

1) $\beta-\frac{\lambda}{\delta} > 0$,

2) $(\beta-\frac{\lambda}{\delta})^2k+d-1 < 0$,

3) $\beta-\frac{\lambda}{\delta}-d+e > 0$,

the system (1.3) undergoes a saddle-node bifurcation when the parameter $\lambda$ satisfies the condition $\lambda = \lambda_{SN}$. The saddle-node bifurcation parameter threshold $\lambda_{SN}$ is given by $\Delta = 0$.

Proof. We check whether the corresponding condition for saddle-node bifurcation hold true by Sotomayor's theorem again. First, if $\lambda = \lambda_{SN}$, the Jacobian matrix at $E_{+}$ can be written in the following form:

where $c = \beta-\frac{\lambda_{SN}}{\delta}$, $a_{12} = -x_{+}\frac{1+2ky_{+}}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}$.

Obviously, one of the eigenvalues of the matrix $J_{E_{+}}$ is zero. Suppose the eigenvectors of $J_{E_{+}}$ and $J_{E_{+}}^T$ with respect to the zero eigenvalue are $V$ and $W$, then

Therefore, we have

        $F_{\lambda}(E_{+}; \lambda_{SN}) = \left(\begin{array}{cc} 0 \\ -y\\ \end{array} \right) _{(E_{+}; \lambda_{SN})} = \left(\begin{array}{cc} 0 \\ -y_{+}\\ \end{array} \right),$

        $D^{2}F(E_{+}; \lambda_{SN})(V, V) = \left(\begin{array}{cc} \frac{\partial^{2}F_{1}}{\partial x^{2}}v_{1}v_{1}+2\frac{\partial^{2}F_{1}}{\partial x \partial y}v_{1}v_{2}+\frac{\partial^{2}F_{1}}{\partial y^{2}}v_{2}v_{2}\\ \frac{\partial^{2}F_{2}}{\partial x^{2}}v_{1}v_{1}+2\frac{\partial^{2}F_{2}}{\partial x \partial y}v_{1}v_{2}+\frac{\partial^{2}F_{2}}{\partial y^{2}}v_{2}v_{2} \\ \end{array}\right) _{(E_{+}; \lambda_{SN})},$

         $= \left(\begin{array}{cc} \frac{-2(1+ck)}{(1+ky_{+})^3}+\frac{2c(dx_{+}^3+3ex_{+}^2-e^2)}{(x_{+}^2+dx_{+}+e)^3} \\ 0\\ \end{array} \right).$

Finally, we get

Hence, the system (1.3) experiences a saddle-node bifurcation when the parameter $\lambda$ satisfies the condition $\lambda = \lambda_{SN}$.□

In order to demonstrate the existence of a saddle-node bifurcation, the next numerical simulation is given in Figure 2. From Figure 2, we can find that the system (1.3) undergoes a saddle-node bifurcation.

When $k = 0.2, d = 0.05, e = 0.024, \delta = 0.3, \beta = 2$, we can easily work out that $\lambda_{SN} = 0.510306079$, which means that the number of internal equilibria of the system (1.3) will change from one to three as $\lambda$ varies around $\lambda_{SN}$. In Figure 2(a), the system (1.3) has only one internal equilibrium point $E_{1}^{\ast}$ when $\lambda = 0.52$, $E_{1}^{\ast}$ is a unstable focus. In Figure 2(b), the system (1.3) has two internal equilibria $E_{+}$ and $E_{2}^{\ast}$ when $\lambda = 0.510306079$, $E_{+}$ is a saddle-node and $E_{2}^{\ast}$ is a unstable node. In Figure 2(c), the system (1.3) has three internal equilibria $E_{4}^{\ast}$, $E_{5}^{\ast}$ and $E_{6}^{\ast}$ when $\lambda = 0.511257$, both $E_{4}^{\ast}$ and $E_{6}^{\ast}$ are unstable node, $E_{5}^{\ast}$ is a saddle, which implies the validity of the Theorem 4.4. In all three cases above, the prey and predator populations will all have periodic oscillations. But the period of their oscillations is not same.

5.3.   Hopf bifurcation

In the preceding analysis, $E_{5}^{\ast}$ is always a saddle whenever it exists and $E_{i}^{\ast}$ (i = 1, 2, 3, 4, 6, 7) may be a sink or a source or a center. Then we only discuss equilibrium point $E_{1}^{\ast}$, because other positive equilibria $E_{i}^{\ast}$ (i = 2, 3, 4, 6, 7) are similar to $E_{1}^{\ast}$. Considering $\lambda$ as the Hopf bifurcation parameter, the Hopf bifurcation parameter threshold is given by $Tr(J_{E_{1}^{\ast}}) = 0$, i.e. $\lambda = \lambda_{H}$. At the same time, it satisfies $Det(J_{E_{1}^{\ast}})|_{\lambda = \lambda_{H}} > 0$. As $\lambda$ changes and pass the threshold $\lambda_{H}$, the stability of $E_{1}^{\ast}$ changes accordingly. Next we prove this conclusion.

Theorem 5.3. If $\Delta > 0$, the system (1.3) has only one positive equilibrium point $E_{1}^{\ast}$, where $x_{1}^{\ast}$ is a single positive root, then the system (1.3) will experience a Hopf bifurcation near $E_{1}^{\ast}$ at $\lambda = \lambda_{H}$.

Proof. We have obtained that $Det(J_{E_{1}^{\ast}}) > 0$ and $Tr(J_{E_{1}^{\ast}}) = 0$ when $\lambda = \lambda_{H}$. Therefore, we only require to check the correctness of transversality condition for the Hopf bifurcation. That is,

The transversality condition is satisfied by our numerical simulation, which means $E_{1}^{\ast}$ lose its stability through Hopf bifurcation at $\lambda = \lambda_{H}$.

It is necessary to determine the direction and stability of the limit cycle, then first Lyapunov number is calculated.

Convert $E_{1}^{\ast}$ to the origin by making the transformation $u = x-x_{1}^{\ast}, v = y-y_{1}^{\ast}$. Therefore, the system (1.3) in a region of the origin is as follows

where

and $P_{10}(u, v)$, $Q_{10}(u, v)$ are $C^\infty$ functions at least of order fourth in $(u, v)$.

The first Lyapunov number

From Theorem 1 in Section 4.4 in ^[23], we have the following conclusion. If $l_1 > 0$, the limit cycle around $E_1^\ast$ is unstable and the Hopf bifurcation is subcritical. If $l_1 < 0$, the limit cycle around $E_1^\ast$ is stable and $E_1^\ast$ lose its stability by a supercritical Hopf bifurcation.

Because the formula of the first Lyapunov number $l_1$ is too complicated, we cannot easily determine its sign. Therefore, the rationality of the theorem will be verified in next simulation from Figure 3.□ Unlike a typical prey-predator system, the system (1.3) can easily produce Hopf bifurcation. When $k = 0.2, d = 0.05, e = 0.024, \delta = 0.06931640934, \beta = 2$, then we can easily work out that $\lambda_{H} = 0.06931640934$, which means that the system (1.3) will produce a stable or unstable limit cycle as $\lambda$ varies around $\lambda_{H}$. Also, the first Lyapunov number $l_{1} = 1269.436132\pi > 0$, i.e. subcritical Hopf bifurcation occurs. In Figure 3(a), when the value of $\lambda$ is greater to 0.06931640934, the system (1.3) has a unstable focus $E_{1}^{\ast}$. As the value of $\lambda$ becomes smaller and equals 0.06931640934, the unstable focus becomes a center in Figure 3(c). Further, when the value of $\lambda$ is smaller than 0.06931640934, the centre becomes a stable focus and the system (1.3) produces a unstable limit cycle in Figure 3(b). Meanwhile, in Figure 3(a), (c), the predator and prey end up coexisting in the form of periodic oscillations regardless of their initial values. In Figure 3(b), choosing different initial values, predator and prey either coexist at stable focus $E_{1}^{\ast}$ or in the form of periodic oscillations.

5.4.   Bogdanov-Takens bifurcation

Above we have discussed the codimension one bifurcations at nonhyperbolic equilibria, and next we will demonstrate the possibility of Bogdanov-Takens bifurcation of codimension 2 occurring. From Theorem 4.2, we obtain that $E_{+}$ is a cusp of codimension 2 under the condition $\Delta = 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0, \beta-\frac{\lambda}{\delta}-d+e > 0, \lambda = \delta\beta+x_{+}\frac{3x_{+}^2+(2d-2)x_{+}+e-d}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}, \lambda\neq\lambda_{1}$ and $h(\lambda)\neq0$. Similarly, we also get that $E_{++}$ is a cusp of codimension 2 when some conditions are satisfied from Theorem 4.3. Next, we will study Bogdanov-Takens bifurcation with parameters $\lambda$ and $k$.

Theorem 5.4. If we choose bifurcation parameters $\lambda$ and $k$, then the system (1.3) experiences a Bogdanov-Takens bifurcation around $E_{+}$ with changing parameters $(\lambda, k)$ near $(\lambda_{BT}, k_{BT})$ for $\Delta = 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0, \beta-\frac{\lambda}{\delta}-d+e > 0, \lambda = \delta\beta+x_{+}\frac{3x_{+}^2+(2d-2)x_{+}+e-d}{(1+ky_{+})(x_{+}^2+dx_{+}+e)}, \lambda\neq\lambda_{1}$ and $h(\lambda)\neq0$, where $(\lambda_{BT}, k_{BT})$ denotes the bifurcation threshold value i.e.

Proof. In order to gain the accurate expressions for saddle-node, Hopf and homoclinic bifurcation curve in a small neighborhood near the Bogdanov-Takens point, we convert the system (1.3) into the standard form of Bogdanov-Takens bifurcation.

Next, a parameter vector $(\varepsilon_{1}, \varepsilon_{2})$ near $(0, 0)$ is introduced in order to give a perturbation to $\lambda$ and $k$ near their B-T bifurcation values given by $\lambda = \lambda_{BT}+\varepsilon_{1}$ and $k = k_{BT}+\varepsilon_{2}$. Substituting the perturbations into the system (1.3), we have

Let $(0, 0)$ become the bifurcation point by introducing the transformation $\eta_1 = x-x_{+}$, $\eta_2 = y-y_{+}$, we get

where

and $P_{1}(\eta_1, \eta_2, \varepsilon)$, $Q_{1}(\eta_1, \eta_2, \varepsilon)$ is power series in $(\eta_1, \eta_2)$ with terms $\eta_1^i\eta_2^j$ requiring $i+j\ge3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

In the following, we perform the transformation

the system (5.2) becomes

where

and $P_{2}(X, Y, \varepsilon)$, $Q_{2}(X, Y, \varepsilon)$ is power series in $(X, Y)$ with terms $X^iY^j$ requiring $i+j\ge3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Then introducing the next $C^{\infty}$ change of coordinates in a small domain of $(0, 0)$:

the system (5.3) can be written as

where

and $Q_{3}(z_1, z_2, \varepsilon)$ is power series in $(z_1, z_2)$ with terms $z_1^iz_2^j$ requiring $i+j\ge3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Using a new variable $\tau$ by $dt = (1-g_{02}(\varepsilon)z_{1})d\tau$, then

Let

then the system (5.5) can be rewritten as

where

and $Q_{4}(v_{1}, v_{2}, \varepsilon)$ is a power series in $(v_{1}, v_{2})$ with terms ${v_{1}}^{i}{v_{2}}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

We note that $h_{20}(\varepsilon)$ is a very complex number, so it is difficult to discern the sign of $h_{20}(\varepsilon)$ when $\varepsilon_{1}$ and $\varepsilon_{2}$ are sufficiently small. Therefore, in order for the next transformations to make sense, we must go on to discuss the next two cases.

Case 1: For sufficiently small $\varepsilon_{1}$ and $\varepsilon_{2}$, if $h_{20}(\varepsilon) > 0$, then we use the next transformation:

We get

where

and $Q_{5}(q_{1}, q_{2}, \varepsilon)$ is a power series in $(q_{1}, q_{2})$ with terms ${q_{1}}^{i}{q_{2}}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Let

then

where

and $Q_{6}(s_{1}, s_{2}, \varepsilon)$ is a power series in $(s_{1}, s_{2})$ with terms $s_{1}^{i}s_{2}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Assume that $h_{11}(\varepsilon)\neq0$, then $J_{11}(\varepsilon) = I_{11}(\varepsilon) = \frac{h_{11}(\varepsilon)}{\sqrt{h_{20}(\varepsilon)}}\neq0$, next through the next transformation:

then

where

and $Q_{7}(X, Y, \varepsilon)$ is a power series in $(X, Y)$ with terms ${X}^{i}{Y}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Case 2: For sufficiently small $\varepsilon_{1}$ and $\varepsilon_{2}$, if $h_{20}(\varepsilon) < 0$, then we use the next transformation:

We get

where

and $Q_{5}^{\prime}(q_{1}^{\prime}, q_{2}^{\prime}, \varepsilon)$ is a power series in $(q_{1}^{\prime}, q_{2}^{\prime})$ with terms ${q_{1}^{\prime}}^{i}{q_{2}^{\prime}}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Let

we have

where

and $Q_{6}(s_{1}^{\prime}, s_{2}^{\prime}, \varepsilon)$ is a power series in $(s_{1}^{\prime}, s_{2}^{\prime})$ with terms ${s_{1}^{\prime}}^{i}{s_{2}^{\prime}}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

Suppose that $h_{11}(\varepsilon)\neq0$, then $J_{11}^{\prime}(\varepsilon) = I_{11}^{\prime}(\varepsilon) = \frac{h_{11}(\varepsilon)}{\sqrt{-h_{20}(\varepsilon)}}\neq0$, next we use the next transformation:

we have

where

and $Q_{7}^{\prime}(X^{\prime}, Y^{\prime}, \varepsilon)$ is a power series in $(X^{\prime}, Y^{\prime})$ with terms ${{X}^{\prime}}^{i}{{Y}^{\prime}}^{j}$ requiring $i+j\geq3$, their coefficients depend on $\varepsilon_{1}$ and $\varepsilon_{2}$ smoothly.

To reduce the number of cases to be considered, we retain $\rho_{1}(\varepsilon)$ and $\rho_{2}(\varepsilon)$ to stand for $\rho_{1}^{\prime}(\varepsilon)$ and $\rho_{2}^{\prime}(\varepsilon)$ in $(5.9^{\prime})$. If the matrix $\left|\frac{\partial (\rho_{1}, \rho_{2})}{\partial (\varepsilon_{1}, \varepsilon_{2})}\right|$ is nonsingular, then transformation is a homeomorphism in a adequately small domain of the $(0, 0)$. At the same time, under the above condition, $\rho_{1}$, $\rho_{2}$ are two independent variables. From the conclusions in ^[21,22] and ^[23], it is obtained that B-T bifurcation is produced when $\varepsilon = (\varepsilon_{1}, \varepsilon_{2})$ is located at a fully small domain of the (0, 0). Therefore, the local formulas near the origin of the bifurcation curves can be written as ("+" denotes $h_{20}(\varepsilon) > 0$ and "-" denotes $h_{20}(\varepsilon) < 0$):

1) The saddle-node bifurcation curve can be represented as

     $SN = \left\lbrace (\varepsilon_{1}, \varepsilon_{2}):\rho_{1}(\varepsilon_{1}, \varepsilon_{2}) = 0, \rho_{2}(\varepsilon_{1}, \varepsilon_{2})\neq0\right\rbrace$;

2) The Hopf bifurcation curve can be represented as

     $H = \left\lbrace (\varepsilon_{1}, \varepsilon_{2}):\rho_{2}(\varepsilon_{1}, \varepsilon_{2}) = \pm\sqrt{-\rho_{1}(\varepsilon_{1}, \varepsilon_{2})}, \rho_{1}(\varepsilon_{1}, \varepsilon_{2}) < 0\right\rbrace$;

3) The homoclinic bifurcation curve can be represented as

     $HL = \left\lbrace (\varepsilon_{1}, \varepsilon_{2}):\rho_{2}(\varepsilon_{1}, \varepsilon_{2}) = \pm\frac{5}{7}\sqrt{-\rho_{1}(\varepsilon_{1}, \varepsilon_{2})}, \rho_{1}(\varepsilon_{1}, \varepsilon_{2}) < 0 \right\rbrace$.

In order to study how the parameters $k$ and $\lambda$ synergistically affect the dynamic behavior of the system (1.3), the numerically simulation of B-T bifurcation with $d = 0.05, e = 0.02, \beta = 0.7, \delta = 2.047273228$ will be carried out in Figure 4. By a simple calculation, we obtain $k_{BT} = 0.25, \lambda_{BT} = 0.8362958526$. In Figure 4(a), the system (1.3) has a cusp of codimension 2 and a unstable foucus when $(\varepsilon_{1}, \varepsilon_{2}) = (0, 0)$. Then as $\varepsilon_{1}, \varepsilon_{2}$ varies, the system (1.3) changes from having only one internal equilibrium point to producing three internal equilibria. In Figure 4(b), the system (1.3) has only one internal equilibrium point $E_{2}^{\ast}$ when $(\varepsilon_{1}, \varepsilon_{2}) = (0.01169, 0.0001)$. In Figure 4(c), the system (1.3) produces two internal equilibria near $E_{+}$ when $(\varepsilon_{1}, \varepsilon_{2}) = (0.01167, 0.00077031)$, one of which is the saddle $E_{5}^{\ast}$ and one of which is the unstable focus $E_{6}^{\ast}$. In Figure 4(d), the system (1.3) produces a Hopf bifurcation, which causes the focus $E_{6}^{\ast}$ to become stable and produces an unstable limit cycle. In Figure 4(e), the unstable limit cycle gets bigger and goes through the saddle $E_{5}^{\ast}$ gradually and converts to an unstable homoclinic orbit finally when $(\varepsilon_{1}, \varepsilon_{2}) = (0.01169, 0.0007755)$. In Figure 4(f), the homoclinic orbit disappears and there exists an stable focus $E_{6}^{\ast}$ and a saddle $E_{5}^{\ast}$ when $(\varepsilon_{1}, \varepsilon_{2}) = (0.01169, 0.005)$. Meanwhile, in Figure 4(a)–(c), both the predator and prey end up coexisting in the form of periodic oscillations regardless of their initial values. In Figure 4(d)–(f), choosing different initial values, predator and prey either coexist at stable focus or in the form of periodic oscillations.

In the numerical simulation from Figure 4, we find that the fear effect and harvesting has a crucial role in influencing the dynamical properties of the system. Due to fear effect and harvesting, the system (1.3) undergoes multiple bifurcations and changes in stability near Bogdanov-Takens point.

In particular, we find that when the system (1.3) has three internal equilibriua $E_{4}^{\ast}, E_{5}^{\ast}, E_{6}^{\ast}$ and at least one of $E_{4}^{\ast}$ and $E_{6}^{\ast}$ is locally stable, then the existence of a homoclinic bifurcation can be inferred, since $E_{5}^{\ast}$ is a saddle. Therefore, the system (1.3) must satisfy the following conditions in order to ensure that homoclinic bifurcation occurs: $\lambda < \delta\beta+x_{4}^{\ast}\frac{3x_{4}^{\ast2}+(2d-2)x_{4}^{\ast}+e-d}{(1+ky_{4}^{\ast})(x_{4}^{\ast2}+dx_{4}^{\ast}+e)}$ or $\lambda < \delta\beta+x_{6}^{\ast}\frac{3x_{6}^{\ast2}+(2d-2)x_{6}^{\ast}+e-d}{(1+ky_{6}^{\ast})(x_{6}^{\ast2}+dx_{6}^{\ast}+e)}$. In order to verify the rationality of the conclusion, corresponding numerical simulation is shown in Figure 5. From Figure 5, we can find that the system (1.3) undergoes a homoclinic bifurcation.

6.   Simulation analysis and results

In order to visualize the effect of constant-effort harvesting on predator in the system (1.3), we explore the dynamical phenomenons of the system (1.3) with different harvesting values through numerical simulations.

We fix the following parameters for convenience:

and assume that initial values for prey and predator densities were $(x(0), y(0)) = (0.1, 0.18)$.

In order to further illustrate Theorems 4.2 and 4.3, we change the values of some parameters in (6.1) to get better simulation results. Different values of $\lambda$ will affect the dynamical properties of $E_{+}$, which may be either a saddle-node or a cusp. In Figure 6(a), the system (1.3) has two internal equilibria $E_{+}$ and $E_{2}^{\ast}$, $E_{+}$ is a saddle-node and $E_{2}^{\ast}$ is a unstable node. That is to say, the prey and predator populations will have periodic oscillation. In Figure 6(c), the form of the trajectory diagram better shows the periodic oscillation. In Figure 6(b), the system (1.3) has two internal equilibria $E_{+}$ and $E_{2}^{\ast}$, $E_{+}$ is a cusp of codimension 2 and $E_{2}^{\ast}$ is a stable focus. That is to say, the prey and predator populations will coexist in $E_{2}^{\ast}$. In Figure 6(d), the form of the trajectory diagram better shows their coexistence.

In Figure 7(a), the system (1.3) has only one internal equilibrium point $E_{++}$, $E_{++}$ is a saddle-node. That is to say, the prey and predator populations will have periodic oscillation. In Figure 7(b), the system (1.3) has only internal equilibrium point $E_{++}$, $E_{++}$ is a cusp of codimension 2. That is to say, the prey and predator populations will have periodic oscillation after a certain time. In Figure 7(c), (d), the form of the trajectory diagram better shows the periodic oscillation. Moreover, we can find that although both Figure 7(c), (d) reach coexistence in the form of periodic oscillation, their periods are different. In Figure 7(d), the increase in the number of codimension makes their periods longer and amplitudes smaller, respectively.

Next we will explore more visually the impact of harvesting on the system (1.3). In Figure 8(a), we can see that two populations eventually converge to the internal equilibrium point $E_{1}^{\ast}(0.012168, 0.024337)$ due to the absence of harvesting i.e. $\lambda = 0$. In Figure 8(b), the two populations still coexist in $E_{1}^{\ast}(0.02511, 0.02511)$ when $\lambda = 0.3$. In Figure 8(c), the two populations still coexist in $E_{1}^{\ast}(0.027533, 0.025287)$ when $\lambda = 0.3244725$. Comparing the above three charts Figure 8(a)–(c), we can find that the densities of predator and prey reaching coexistence respectively increases with the increase of $\lambda$ and the value of coexistence is smaller than the initial value. As the harvesting increases, the system (1.3) will generate multiple stability shifts. As $\lambda$ increases to 0.52, the system (1.3) loses stability and the stable focus becomes an unstable focus. That is to say, two populations coexist in the form of periodic oscillation when $\lambda = 0.52$ (see Figure 8(d)). As $\lambda$ increases to 0.52223, the system (1.3) is stabilized again and the unstable focus becomes an stable focus. In Figure 8(e), when $\lambda = 0.52223$, harvesting has a delayed effect on population densities, which means that the influence of harvesting for the prey and predator populations will not be evident from the outset. That is to say, two populations will eventually coexist, although they will coexist in a periodic oscillation for a long time in the beginning. As $\lambda$ slowly increases, the system (1.3) begins to stabilize and the periodic oscillation of the two populations slowly disappears. Because the stable focus becomes a stable node, the two populations will coexist and periodic oscillation have disappeared when $\lambda = 0.58$ (see Figure 8(f)). In Figure 8(g), when $\lambda = 0.6$, predator and prey populations gradually converge to the boundary equilibrium point $E_{1}(1, 0)$, but they don't exactly overlap. That is to say, the predator is close to the state of extinction. In Figure 8(h), when $\lambda = 0.7$, predator population will become extinct and prey population will reach its peak in a very short period of time. Comparing the above four charts Figure 8(e)–(h), we can find that the densities of coexistence gradually approach $E_{1}(1, 0)$ as $\lambda$ increases. Further, we can also find that the density of predator reaching coexistence decreases with the increase of $\lambda$. However, the change for density of prey is opposite.

7.   Conclusions

Many articles have demonstrated fear effect or harvesting can influence the dynamic behaviors of predator-prey system. However, the joint effects of fear effect and harvesting in the predator-prey system are rarely discussed by researchers. Hence, we investigated the dynamics of a predator-prey system where prey is in possess of fear effect and predator is provided with constant-effort harvesting. Qualitative analysis exposes that predator harvesting term has a crucial role in determining the dynamic behaviors of the system (1.3). Firstly, the system (1.3) at the beginning is reduced to an equivalent system with only six parameters through a suitable transformation. Secondly, we see that the system (1.3) has bounded solutions when the harvesting satisfies certain conditions and the origin is always unstable. Then, the parameter $\lambda$ standing for harvesting can influence the stability and number of equilibria. Also, we get that the system (1.3) always possesses a boundary equilibrium point, which may be a stable node, a saddle or a repelling saddle-node for different $\lambda$. The system (1.3) owns positive equilibria only when the boundary equilibrium is hyperbolic saddle. For positive equilibria, the sign of determinant of their Jacobian matrixs is easy to identify. Therefore, there are a saddle-node or a cusp with codimension 2 choosing different parameters values. Additionally, a cusp with codimension at least 3 also exists. Furthermore, other internal equilibria may be saddle, sink, source or center. Meanwhile, we use Sotomayor's theorem to prove the existence of saddle-node and transcritical bifurcation. In particular, since the determinant of the jacobian matrix is always positive for most of the equilibria, the system (1.3) is prone to produce Hopf bifurcation. Also, we computed the first Lyapunov number to investigate whether the limit cycle generated by the Hopf bifurcation is stable. Next, we chose $k$ and $\lambda$ as bifurcation parameters and proved that the system (1.3) went through Bogdanov-Takens bifurcation with codimension 2 by showing the standard unfolding around the cusp. At the same time, we give the conditions under which the system (1.3) undergoes a homoclinic bifurcation when there are three internal equilibria in the system (1.3).

In the numerical simulation, we fixed a series of parameters and varied the value of $\lambda$ to observe the effect of the harvesting on the two populations. If harvesting exceed a critical threshold, predator will eventually go extinct. When the harvesting is less than the critical threshold and slowly increases, the equilibrium point of the system (1.3) will lose stability and the population will undergo periodic oscillations. Then, this equilibrium point regained stability again. In addition, the system (1.3) is highly susceptible to periodic oscillation.

The harvesting term in the system (1.3) is biologically important, affecting the number and stability of equilibria and the generation of various bifurcations. When the harvesting in predators is too great, i.e. $\lambda > \delta\beta$, the predator will quickly become extinct and the prey will reach its peak. In reality, we need to make measures to control the size of the harvesting to ensure that both species do not become extinct. Also, selecting different harvesting sizes will lead to coexistence or periodic oscillation in the two populations, which will help protect the survival of predator. Therefore, policy makers should develop optimal harvesting programs to achieve ecologically sustainable development. This article will provide some help on how predator harvesting affects the density of two populations.

In the next research work, since neither prey nor predator will remain in a fixed space due to various reasons such as food replenishment, mate choice, and population migration, it will be meaningful to investigate how fear effect for prey and predator harvesting together affect its spatial dynamic distribution ^[39,40,41]. In addition, we can go on to explore the system (1.3) will have a Bogdanov-Takens bifurcation of codimension 3.

Acknowledgments

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

Conflict of interest

The authors declare there is no conflict of interest.