Dynamic analysis of a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting

1.   Introduction

Fishing is a method used in the industry to acquire fish products from natural or artificial bodies of water. With the development of fisheries, fishing has become more common. We see that harvesting for economic gain is a relatively regular occurrence in nature and that it significantly affects both the ecological balance and system dynamics. It is crucial to develop biological resources at their maximum sustainable yield while preserving the survival of all interacting populations, both ecologically and economically. However, if a species is overharvested, it can lead to ecological problems, as some people may prioritize profit over protecting the environment. Thus, the authors of ^[1,2] built mathematical models to analyze these problems, whose dynamical behaviors have attracted the interest of many scholars. There are three forms of harvesting: $1)$ constant harvesting, $h(x) = h$; $2)$ linear harvesting, $h(x) = qEx$; $3)$ nonlinear harvesting, $h(x) = \frac{qEx}{m_{1}E+m_{2}x}$, which is also called Michaelis-Menten-type harvesting.

Leslie and Gower studied the predator-prey relationship between two species and developed the famous Leslie-Gower predator-prey model ^[3], which has been widely discussed ^[4,5,6]. For example, the Leslie-Gower predator-prey model with the Allee effect and a generalist predator was studied in ^[7], where the authors found that the system exhibits a multi-stability phenomenon and undergoes various bifurcations. Huang et al. ^[8] studied the Leslie-Gower-type predator-prey model with constant-yield harvesting, and they found that the dynamical behavior of the model is very sensitive to the constant yield harvest of the predator.

Gupta et al. ^[9] considered the following Leslie-Gower predator-prey model with Michaelis-Menten-type prey-harvesting:

where $x$ and $y$ are the prey and predator population densities, respectively. They studied the stability and bifurcation (saddle-node bifurcation and Hopf bifurcation) of system (1.1). Also, the existence of bionomic equilibria and optimal singular control were investigated. Based on system (1.1), Gupta and Chandra ^[10] introduced the Holling type Ⅱ functional response and obtained the bistable situation. The model exhibits several local bifurcations (saddle-node, Hopf, homoclinic and Bogdanov-Takens) which are ecologically important. Considering the group defense and nonlinear harvesting in prey, Kumar and Kharbanda ^[11] obtained that the density of the predator increases as the harvest rate of the predator decreases. Caraballo Garrido et al. ^[12] investigated the predator prey model with nonlinear harvesting with both constant and distributed delay by varying parameters. Some scholars ^{[13,14,15,16,17]} have combined other functional responses and harvesting to obtain more complex dynamical behaviors.

As with direct killing, indirect killing also has a significant impact on the dynamic behaviors of the system. The studies mentioned above, however, solely take into account the predator's direct killing. Predation danger may drive the prey to engage in anti-predation behaviors, such as habitat modifications or foraging, which may lower the prey's birth rate. Hence, Wang et al. ^[18] incorporated the fear effect into the reproduction of prey animals and obtained the following prey-predator model:

where $f(k, y) = \frac{1}{1+k_0y}$ accounts for the cost of anti-predator defense due to fear. They studied a model with a linear functional response or Holling type Ⅱ functional response. It was found that the fear effect has no impact on the dynamic behaviors of model (1.2). However, the dynamic behavior of model (1.2) with Holling type Ⅱ functional responses can be affected by the fear effect. Chen et al. ^[19] considered the influence of the fear effect and Leslie-Gower function on the dynamic behavior of the predator-prey model, and they demonstrated that there are many types of bifurcation phenomena, including transcritical bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation. Zhang et al. ^[20] studied a delayed diffusive predator-prey model with spatial memory and a nonlocal fear effect, and they investigated the stability, Hopf bifurcation and Turing-Hopf bifurcation of the system. Many scholars ^{[21,22,23,24,25]} have studied the prey-predator model with a fear effect.

In this paper, we incorporate the fear effect into system (1.1) and obtain the following model:

where $r_0$ is the birth rate of the prey and $d_0$ is the natural death rate of the prey. In the ecological sense, it is clear that $r_0 > d_0$. $a$ represents the intra-species competition, $c_0$ is the maximum predation rate, $s_0$ is the intrinsic growth rates of the predators, $h$ is a measure of the quality of the prey as food for the predator, $k_0$ is the fear parameter, $q_0$ is the catchability coefficient, $E$ is the effort applied to harvest the prey species and $m_1, m_2$ are suitable constants. For simplicity, letting

for $x, y$ be positive, and dropping the bars, system (1.3) becomes

where $r > d$, and $r$, $k$, $d$, $c$, $q$ and $s$ are positive constants.

The key aim of this study on prey-predator models is to discuss the impacts of prey fear and prey harvesting on system dynamics. The bifurcation phenomenon that distinguishes system (1.4) from system (1.1) deserves further discussion. In addition, by analyzing the observed bifurcation phenomena, we can elucidate the benefits and drawbacks of prey harvesting on both populations.

The structure of the article is as follows. In Section 2, we obtain the boundedness of solutions and analyze the dynamical behaviors of origin. In Section 3, we discuss the existence of boundary equilibria and positive equilibria. In Section 4, we analyze the stability of equilibria. In Section 5, we show that system (1.4) undergoes saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation. In Section 6, we have a summary of the article.

2.   Preliminaries

We show that the positive solutions of system (1.4) are ultimately bounded.

Theorem 2.1. All solutions of system (1.4) are bounded for all $t\geq0$.

Proof. Since system (1.3) is equivalent to system (1.4), we now prove that the solution of system (1.3) is bounded. From the first equation of system (1.3), we have

for $t\geq0$, which immediately implies that $\limsup\limits_{t\to\infty}x(t)\le \frac{r_0-d_0}{a}\triangleq M$. Then from the second equation of system (1.3), it follows that

for $t\geq0$, that is, $\limsup\limits_{t\to\infty}y(t)\le hM$. Hence, $x(t)$ and $y(t)$ are bounded. The proof is completed.

Denote

Next, we show the dynamic behaviors of the origin of system (1.4).

Lemma 2.1. The types of origin in system (1.4) are as follows:

1) if $r-d < q$ or $r-d = q\le q_0$, the origin of system (1.4) is a non-hyperbolic attractor;

2) if $r-d > q$ or $r-d = q > q_0$, the origin of system (1.4) is a non-hyperbolic repeller.

Proof. The Jacobian matrix of system (1.4) at the origin is degenerate; then, we apply the blow-up method to analyze the type of origin. Notice that when $x = 0$, we have that $\dot{x} = 0$ and $\dot{y} = -sy^2 < 0$, which means that system (1.4) has the invariant line $x = 0$. Using the horizontal blow-up

system (1.4) can be rewritten as

(ⅰ) The equilibria of system (2.1) in $u = 0$ are $A(0, 0)$ and $B\left(0, 1-\frac{r-d-q}{s}\right)$ when $r-d < q+s$. The Jacobian matrix at the equilibria $A$ and $B$ are, respectively

and

where

The eigenvalues of matrix $J_{A}$ are $\lambda_{J_{A1}} = r-d-q$ and $\lambda_{J_{A2}} = q+s-(r-d) > 0$, and the eigenvalues of matrix $J_{B}$ are $\lambda_{J_{B1}} = r-d-q$ and $\lambda_{J_{B2}} = r-d-(q+s) < 0$. If $r-d < q$, that is, if $\lambda_{J_{A1}} < 0$ and $\lambda_{J_{B1}} < 0$, $A$ is a saddle and $B$ is a stable node (see Figure 1(a)). If $q < r-d < q+s$, that is, if $\lambda_{J_{A1}} > 0$ and $\lambda_{J_{B1}} > 0$, $A$ is an unstable node and $B$ is a saddle (see Figure 1(c)).

If $r-d = q$, that is, if $\lambda_{J_{A1}} = 0$ and $\lambda_{J_{B1}} = 0$, both $A$ and $B$ are degenerate equilibria. First, we consider the degenerate equilibrium $A$. Taking the time variable

we can obtain the Taylor expansion of system (2.1) as follows:

By Theorem 7.1 in ^[26], the degenerate equilibrium $A$ is a saddle-node if $q\neq 1$. If $q = 1$, we have that the equilibrium $A$ is a degenerate saddle from $- \frac{q}{s} < 0$.

Next, for the degenerate equilibrium $B$, make the following transformation:

System (2.1) becomes

Using $(X, Y) = \left(sX_1, \big[(k-1)q+dk+c+1\big]X_1+Y_1\right)$ and $\mbox{d}\tau = -s\, \mbox{d}t$, system (2.2) becomes

where

If $\alpha_{20} > 0$ (or $\alpha_{20} < 0$), that is, if $q < q_0$ (or $q > q_0$), $B$ a saddle-node with a stable parabolic sector on the right (or left). If $\alpha_{20} = 0$, which implies that $k < 1$, we get that $q = \frac{d k+c+1}{1-k}$. Next, substituting $q = \frac{d k+c+1}{1-k}$ into the coefficients of the $X^3$ term of system (2.3), we have

Explicitly, from Theorem 7.1 in ^[26] we know that the degenerate equilibrium $B$ is a stable node if $q = q_0$.

In summary, when $r-d = q\le1$, A and B are as shown in Figure 1(a). When $1 < r-d = q\le q_0$, A and B are as shown in Figure 1(b). When $r-d = q > q_0$, A and B are as shown in Figure 1(c).

After a blow-down, the origin in system (1.4) is an attractor when $r-d < q$, $r-d = q\le 1$ (see Figure 2(a)) or $1 < r-d = q\le q_0$ (see Figure 2(b)). The origin is a repeller when $r-d = q > q_0$ or $q < r-d < q+s$ (see Figure 2(c)).

(ⅱ) When $r-d = q+s$, system (2.1) has only one equilibrium, $A(0, 0)$ at $u = 0$, whose Jacobian matrix is

Expanding system (2.1) in a Taylor series and taking a time variable $\mbox{d}\tau = s\, \mbox{d}t,$ we have

The coefficient of $v^2$ is $-1 < 0$; from Theorem $7.1$ in ^[26], we know that the equilibrium $A$ is a saddle-node (see Figure 1(d)). After a blow-down, we see that the origin is a repeller for system (1.4) (see Figure 2(d)).

(ⅲ) System (2.1) has a unique equilibrium $A(0, 0)$ when $r-d > q+s$. The Jacobian matrix at the equilibrium $A$ is

Obviously, the two eigenvalues are $r-d-q > 0$ and $d-r+q+s < 0$. It shows that the equilibrium $A$ is a saddle (see Figure 1(d)). Similarly, the origin is a repeller for system (1.4) (see Figure 2(d)). The proof is completed.

3.   Existence of equilibria

In this section, we will discuss the existence of the boundary equilibria and the positive equilibria of system (1.4).

First, we analyze the existence of boundary equilibria. When $y = 0$, the first equation of system (1.4) can be simplified into

We have

and the discriminant of $f(x)$ is as follows:

where

The two roots of $f(x) = 0$ can be expressed as

If $r-d > q$, $f(x) = 0$ has only one positive root $x_{2}$. If $r-d = q$, $f(x) = 0$ has only one positive root $q-1$ if $q > 1$.

Assume that $r-d < q$. When $r-d\le 1$, obviously, $f(x) = 0$ has no positive roots. When $r-d > 1$, obviously, $q^* > r-d$. If $q > q^*$, then $f(x) = 0$ has no positive roots. If $q = q^*$, then $f(x) = 0$ has only one positive root $\frac{r-d-1}{2}$. If $q < q^*$, then $f(x) = 0$ has two positive roots $x_{1}$ and $x_{2}$.

To summarize, we have the following lemma.

Lemma 3.1. The following claims regarding the existence of the boundary equilibria of system (1.4) are true.

1) If $r-d > q$, system (1.4) has a unique boundary equilibrium $E_{2}(x_{2}, 0)$.

2) If $1 < r-d = q$, system (1.4) has a unique boundary equilibrium $E_{2}(q-1, 0)$.

3) If $1 < r-d < q$, we have the following three cases:

(a) if $1 < r-d < q^* < q$, system (1.4) has no boundary equilibrium;

(b) if $1 < r-d < q = q^*$, system (1.4) has a unique boundary equilibrium $\overline{E}\left(\frac{r-d-1}{2}, 0\right)$;

(c) if $1 < r-d < q < q^*$, system (1.4) has two boundary equilibria $E_{1}(x_{1}, 0)$, $E_{2}(x_{2}, 0)$.

Next, we will discuss the positive equilibria $E(x, y)$ of system (1.4). Letting $\dot{x}\, = \dot{y}\, = \, 0$ in system (1.4), we have

We denote

and

where the discriminant of $F^{\prime}(x)$ is

Define

From $F(x) = 0$, we have

The Jacobian matrix of system (1.4) at $E(x, y)$ is

and

Substituting (3.1) into $DetJ_{E}$ and $F^{\prime}(x)$, we have

When $r-d > q$, it is easy to get that the equation $F(x) = 0$ has a unique positive root $x_4$ (see Figure 3(a)).

When $q_0 < r-d = q$, $\Delta_{2} > 0$. We obtain that $F(x) = 0$ has only one positive root $x_4$ (see Figure 3(b)). When $q_0\ge r-d = q$, we find that $F(x) = 0$ has no positive roots.

When $r-d < q$, obviously, $F(x) = 0$ has no positive roots if $q_0\ge r-d$. If $q_0 < r-d$, we obtain that $F(x) = 0$ has no positive roots when $F(x_{*}) > 0$ (see Figure 3(c)). When $F(x_{*}) = 0$, $F(x) = 0$ has only one positive root $x_{*}$ (see Figure 3(d)). When $F(x_{*}) < 0$, $F(x) = 0$ has two positive roots $x_3$ and $x_4$ (see Figure 3(e)).

To summarize, we have the following lemma.

Lemma 3.2. The following claims regarding the existence of the boundary equilibria of system (1.4) are true.

1) If $r-d > q$, system (1.4) has a unique positive equilibrium $E_{4}(x_{4}, y_{4})$.

2) If $q_0 < r-d = q$, system (1.4) has a unique positive equilibrium $E_{4}(x_{4}, y_{4})$.

3) If $q_0 < r-d < q$, we obtain the following results:

(a) if $F(x_{*}) > 0$, system (1.4) has no positive equilibrium;

(b) if $F(x_{*}) = 0$, system (1.4) has a unique positive equilibrium $E_{*}(x_{*}, y_{*})$;

(c) if $F(x_{*}) < 0$, system (1.4) has two positive equilibria $E_{3}(x_{3}, y_{3})$ and $E_{4}(x_{4}, y_{4})$.

4.   Stability of equilibria

In this section, we will discuss the stability of the equilibria.

4.1.   Stability of the boundary equilibria $E_1$ and $E_2$

Theorem 4.1. 1) If $r-d > q$, the unique boundary equilibrium $E_{2}$ is a saddle.

2) If $1 < r-d = q$, the unique boundary equilibrium $E_{2}$ is a saddle.

3) If $1 < r-d < q = q^*$, the unique boundary equilibrium $\overline{E}$ is a saddle-node.

4) If $1 < r-d < q < q^*$, $E_{1}$ is unstable and $E_{2}$ is a saddle.

Proof. 1) The Jacobian matrix of system (1.4) at $E_2$ is

Obviously, the equilibrium $E_{2}$ is a saddle.

2) The Jacobian matrix of system (1.4) at $E_2(q-1, 0)$ is

which implies that $E_{2}$ is a saddle.

3) The Jacobian matrix of system (1.4) at $\overline{E}$ is

which means that $\overline{E}$ is a degenerate equilibrium. First, we transform $\overline{E}$ to the origin by letting $X = x- \frac{r-d-1}{2}, Y = y.$ Then, system (1.4) can be rewritten as

where

Next, applying the following transformation:

system (4.1) becomes

where

Since $c_{20} < 0$, we get that $\overline{E}$ is a saddle-node from Theorem 7.1 in ^[26].

4) The Jacobian matrices of system (1.4) at $E_{1}\left(x_1, 0\right)$ and $E_{2}\left(x_2, 0\right)$, respectively, are

and

This proves that $E_{1}$ is unstable and $E_{2}$ is a saddle point. The proof is completed.

4.2.   Stability of the interior equilibrium when $r-d\geq q$

If $E_{3}$ exists, from $F'(x_{3}) < 0$ and (3.2), we obtain that $DetJ_{E_{3}} < 0$. That is, $E_{3}$ is a saddle if it exists. Hence, in the following discussion, we only study the stability of $E_{4}$.

Define

Theorem 4.2. When $r-d > q$ or $r-d = q > q_0$, system (1.4) has a positive equilibrium $E_{4}$. In addition, the following statements are true.

1) If $s^{*}\le 0$ or $0 < s^{*} < s$, $E_4$ is locally asymptotically stable.

2) If $s < s^{*}$, $E_4$ is an unstable node or focus.

3) If $s = s^{*} > 0$, $E_4$ is a center or weak focus.

Proof. The Jacobian matrix at the equilibrium $E_{4}(x_{4}, y_{4})$ is

Notice that $F^{\prime}(x_{4}) > 0$; so, from (3.2), we know that $DetJ_{E_{4}} > 0$.

Obviously,

Hence, it is easy to see that $E_{4}$ is locally asymptotically stable if $s^{*}\le 0$ or $0 < s^{*} < s$, an unstable node or focus if $s < s^{*}$ and a center or weak focus if $s = s^{*} > 0$. The proof is completed.

Define

Theorem 4.3. If $E_4$ is locally asymptotically stable and $0 < q < q_1$, then $E_{4}$ is globally asymptotically stable.

Proof. Noting that $r-d > q$ or $r-d = q > q_0$, and according to Lemmas 3.1 and 3.2, system (1.4) has a boundary equilibrium $E_2$ and positive equilibrium $E_{4}$. By Lemma 2.1 and Theorem 4.1, both the origin and $E_2$ are unstable. We assume that $E_{4}$ is locally asymptotically stable. Now, we want to show that there is no limit cycle around $E_{4}$. Hence, taking the Dulac function $\Phi(x, y) = \frac{1}{x^{2}y^{2}}$, we have

where $F = x^2\left(\frac{r}{1+ky}-d-x-cy\right)-\frac{qx^2}{1+x}$ and $G = sy\left(x-y\right)$.

Define

Obviously, $H > 0$ for $q\le 1+2s$. In what follows, we only consider that $q > 1+2s$. By calculation, the discriminant of $H$ is

By calculation, we obtain that $q_1 > 1+2s$ and

Thus, we have that $\Delta_{3} > 0$ for $1+2s < q < q_1$, which means that $H = 0$ has no positive roots. Then, $H > 0$.

To sum up, $H > 0$ when $0 < q < q_1$. Then, we have

which implies that $E_4$ is globally asymptotically stable. The proof is completed.

Remark 4.1. If $q < min\{1, r-d\}$, from Theorem 4.2, $E_4$ is locally asymptotically stable. Therefore, by Theorem 4.3, $E_4$ is globally asymptotically stable. That is, when the intrinsic growth rate of the prey is high and the catchability coefficient for prey is low, the prey and predator will reach a steady state. Hence, a small catchability coefficient for prey will not lead to the extinction of prey and predator.

Remark 4.2. Assume that $r-d = q\leq q_0$. From Lemma 2.1, the origin is a attractor. By Lemma 3.1 and Theorem 4.1, $E_2$ is unstable if it exists. From Lemma 3.2, system (1.4) has no positive equilibrium. Note that system (1.4) has no limit cycle. Therefore, the origin is globally asymptotically stable. That is, when the intrinsic growth rate of the prey and the catchability coefficient for prey are low, the prey and predator will become extinct.

4.3.   Stability of the interior equilibrium when $r-d < q$

Lemma 4.1 (^[27]). The system

is equivalent to the system

in some small neighborhood of $(0, 0)$ after changes to the coordinates.

Lemma 4.2 (^[27]). The system given by

is equivalent to the system given by

by some nonsingular transformations in the neighborhood of $(0, 0)$, where $G = a_{31}-a_{30}a_{21}$.

By computation, from $F(x_{*}) = F^{\prime}(x_{*}) = 0$, we can express $k$ and $r$ in terms of $x_{*}$, $c$, $d$, $s$ and $q$, as follows:

where $q > (1+x_{*})^2(1+c)$ because $k$ and $r$ are positive.

Notice that

clearly, $q-(r-d) > 0$ and $(r-d)-q_0 > 0$ when $q > (1+x_{*})^2(1+c)$. Thus, $q_0 < r-d < q$.

By computation, we have

From the above discussions, we can obtain the following theorem.

Theorem 4.4. Assume that (4.2) and $q > (1+x_{*})^2(1+c)$ hold; system (1.4) has a unique positive equilibrium $E_{*}(x_{*}, y_{*})$.

1) If $s\le cx_{*}$ and $q > (1+x_{*})^2(1+c)$, or if $s > cx_{*}$ and $q > \frac{(x_{*}+s)(1+x_{*})^{2}}{x_{*}}$, then $E_{*}$ is a saddle-node with an unstable parabolic sector.

2) If $s > cx_{*}$ and $(1+x_{*})^2(1+c) < q < \frac{(x_{*}+s)(1+x_{*})^{2}}{x_{*}}$, then $E_{*}$ is a saddle-node with a stable parabolic sector.

Proof. Obviously, we have that $DetJ_{E_{*}} = 0$ by $(3.2)$. Then, the type of $E_{*}$ depends on the sign of $TrJ_{E_{*}}$, as follows:

First, moving $E_{*}(x_{*}, y_{*})$ to the origin by the transformation $(x, y) = (X+x_{*}, Y+y_{*})$, it follows that system (1.4) becomes

where

The eigenvalues of the Jacobian matrix at point $E_{*}$ are $\lambda_{1} = 0$ and $\lambda_{2} = \hat{a}_{10}+\hat{b}_{01}$. If $q\neq \frac{(x_{*}+s) (1+x_{*})^{2}}{x_{*}}$, then $\lambda_{2}\neq 0$.

Next, taking the transformation

and, by introducing the new time variable

system (4.3) is rewritten as

where

By a simple calculation, we get

where

Note that $q > (1+x_{*})^2(1+c)$; then, $q > (1+x_{*})^2$. By computation, we can obtain

which implies that $M < 0$ for $q > (1+x_{*})^2(1+c)$. Therefore, the sign of $\hat{c}_{20}$ is determined by $TrJ_{E_{*}}$. Considering the time transformation, and by using Theorem $7.1$ in ^[26], if $s > cx_{*}$ and $(1+x_{*})^2(1+c) < q < \frac{(x_{*}+s)(1+x_{*})^{2}}{x_{*}}$, that is, if $TrJ_{E_{*}} < 0$, then $E_{*}$ is a saddle-node with a stable parabolic sector (see Figure 4(a)). If $s\le cx_{*}$ and $q > (1+x_{*})^2(1+c)$, or if $s > cx_{*}$ and $q > \frac{(x_{*}+s)(1+x_{*})^{2}}{x_{*}}$, that is, if $TrJ_{E_{*}} > 0$, then $E_{*}$ is a saddle-node with an unstable parabolic sector (see Figure 4(b)). The proof is completed.

From $F(x_{*}) = F^{\prime}(x_{*}) = TrJ_{E_{*}} = 0$, we can express $k$, $r$ and $q$ in terms of $x_{*}$, $c$, $d$ and $s$, as follows:

where $s > cx_{*}$.

Theorem 4.5. Assume that (4.4) and $s > cx_{*}$ hold.

1) If one of the following conditions holds: (1.1) $x_*\geq 1$; (1.2) $0 < x_*\leq\frac{c}{c+2}$; (1.3) $\frac{c}{c+2} < x_* < 1, \, \, s\neq\frac{2 x_{*}^{2}}{1-x_{*}}$, $E_{*}$ is a cusp of codimension two.

2) If $\frac{c}{c+2} < x_* < 1$ and $s = \frac{2 x_{*}^{2}}{1-x_{*}}$ hold, $E_{*}$ is a cusp of codimension three.

Proof. 1) Let $X = x-x_{*}$ and $Y = y-y_{*}$; then, system (1.4) can be rewritten as follows:

Applying the transformation $(u, v) = (- \frac{1}{sx_{*}}X, -X+Y)$, system (4.5) becomes

where

By Lemma 4.1, system (4.6) is equivalent to the following:

where

Substituting $s = cx_{*}$ into $T_1$, we get

So, $T_1 < 0$ for $s > cx_{*}$, that is, $D_{20}\neq 0$.

Obviously, if $x_*\geq 1$, we have that $D_{11} > 0$, that is, $E_{*}$ is a cusp of codimension two by the result in ^[28] (see Figure 4(c)). When $x_* < 1$, from $D_{11} = 0$, we have that $s = \frac{2 x_{*}^{2}}{1-x_{*}}$. Noting that $s > cx_{*}$, we have

Therefore, if $0 < x_*\leq\frac{c}{c+2}$ or $\frac{c}{c+2} < x_* < 1, \, \, s\neq\frac{2 x_{*}^{2}}{1-x_{*}}$ holds and $E_{*}$ is a cusp of codimension two.

2) If $\frac{c}{c+2} < x_* < 1$ and $s = \frac{2 x_{*}^{2}}{1-x_{*}}$ hold, that is, $D_{11} = 0$; we will show that $E_{*}$ is a cusp of codimension three. When $\frac{c}{c+2} < x_* < 1, \, \, s = \frac{2 x_{*}^{2}}{1-x_{*}}$ and (4.4) reduces to the following:

Note that $\frac{c}{c+2} < x_* < 1$; then, $k, r, q$ in (4.7) are positive.

Then, system (4.5) becomes

where

Let $x_2 = y_1$ and $y_2 = \dot{y}_1$; then, system (4.8) becomes

where

Let $x_3 = x_2$ and $y_3 = (1-i_{02}x_2)y_2$; then, system (4.9) becomes

where

To delete the $y_3^3$-term, $x_3y_3^3$-term and $y_3^4$-term in system (4.10), we do the following two transformations:

Hence, system (4.10) becomes

where

Using $\frac{c}{c+2} < x_* < 1$, we have that $m_{20} = - \frac{2((c+4)(c+1)x_*-c^2+d+1)x_{*}^{5}}{(cx_{*}^2+(-c+d-3)x_{*}-d-1)(x_{*}-1)} < 0$. Letting

system (4.11) becomes (still denoting $\tau$ by $t$)

where

By Lemma 4.2, system (4.12) is equivalent to the following system:

where

and

here, the coefficients of $P_i, i = 0, \cdots, 9$ are given in Appendix A.

Using $\frac{c}{c+2} < x_* < 1$, the sign of $G$ is determined by $\delta(x_*)$. By computation, we have that $\delta(\frac{c}{c+2}) = \frac{(16 (29 c^2+80 c+56)) (3 c^2+c d+5 c+2 d+2)^4}{(c+2)^9} > 0$ and $\delta(1) = 384c+128d+384 > 0$. Using Lemma 3.1 in ^[29], the number of roots for $\delta(x_*)$ in $\frac{c}{c+2} < x_* < 1$ is equal to that of positive roots for

in $\frac{c}{c+2} < x_* < 1$, and the coefficients of $M_i, i = 0, \cdots, 9$ are given in Appendix A. Obviously, $M_i, i = 0, \cdots, 9$ are positive. Hence, $\mu(x_*) > 0$ in $\frac{c}{c+2} < x_* < 1$, which implies that $\delta(x_*)$ has no positive zeros in $\frac{c}{c+2} < x_* < 1$. Then, $\delta(x_*)\neq 0$, that is $G\neq0$ in $\frac{c}{c+2} < x_* < 1$, which means that $E_{*}$ is a cusp of codimension three (see Figure 4(d)). The proof is completed.

Theorem 4.6. Assume that $q_0 < r-d < q$ and $F(x_{*}) < 0$; system (1.4) has two positive equilibria $E_{3}(x_{3}, y_{3})$ and $E_{4}(x_{4}, y_{4})$, where $E_{3}$ is always a saddle point. Moveover,

1) if $q < \frac{(x_4+s) (1+x_4)^{2}}{x_4}$, $E_4$ is a stable node or focus;

2) if $q > \frac{(x_4+s) (1+x_4)^{2}}{x_4}$, $E_4$ is an unstable node or focus;

3) if $q = \frac{(x_4+s) (1+x_4)^{2}}{x_4}$, $E_4$ is a center or weak focus.

Proof. It is clear that $F^{\prime}(x_{3}) < 0$ and $F^{\prime}(x_{4}) > 0$ (see Figure 3(e)). Combining (3.2) again, we have that $DetJ_{E_{3}} < 0$ and $DetJ_{E_{4}} > 0$. Thus, $E_{3}$ is a saddle point.

In what follows, we consider that

When $TrJ_{E_{4}} < 0$, $E_{4}$ is a stable node or focus; when $TrJ_{E_{4}} > 0$, $E_{4}$ is an unstable node or focus; when $TrJ_{E_{4}} = 0$, $E_{4}$ is a center or weak focus. The proof is completed.

5.   Bifurcation

We will analyze the bifurcations of system (1.4) in this section, including saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation.

5.1.   Saddle-node bifurcation

From Lemma 3.1, when $1 < r-d < q < q^*$, system (1.4) has two boundary equilibria $E_{1}(x_{1}, 0)$ and $E_{2}(x_{2}, 0)$. However, when $q = q_{SN} = q^*$, only the boundary equilibrium $\overline{E}$ exists. Therefore, according to Sotomayor's theorem ^[28], system (1.4) will produce a saddle-node bifurcation at $\overline{E}$.

Theorem 5.1. Assume that $1 < r-d < q$, with $q = q_{SN}$ being the bifurcation parameter; then, system (1.4) will undergo saddle-node bifurcation at $\overline{E}$.

Proof. The eigenvalues of $J_{\overline{E}}$ are $\lambda_{1} = 0$ and $\lambda_{2} = \frac{s (r-d-1)}{2}$. Denote the eigenvectors of $J_{\overline{E}}$ and $J_{\overline{E}}^{T}$ as

and

respectively.

Denote

Then,

It is easy to know that

This proves the transversality conditions, which means that system (1.4) will undergo saddle-node bifurcation at $\overline{E}$. The proof is completed.

Similarly, it follows from Lemma 3.2 that a saddle-node bifurcation occurs at the positive equilibrium $E_{*}$.

Theorem 5.2. Assume that $q_0 < r-d < q$ and $F(x_{*}) = 0$; system (1.4) will undergo saddle-node bifurcation at $E_{*}$.

5.2.   Hopf bifurcation

It follows from Theorem 4.6 that, if $q = \frac{(x_4+s) (1+x_4)^{2}}{x_4}$, then $TrJ_{E_4} = 0$. Noting that $DetJ_{E_4} > 0$, the Jacobian matrix of $E_4$ has a pair of purely imaginary eigenvalues. Thus, system (1.4) may undergo Hopf bifurcation at $E_4$.

For simplicity, similar to the analyses of Dai et al. ^[30] and Lu et al. ^[31], we prove the Hopf bifurcation. Letting

and by dropping the tilde, system (1.4) becomes

where $r > d$ and the other parameters are positive.

Clearly, $\tilde{E}_4(1, 1)$ is an equilibrium of system (5.1), which implies that

Define

Assume that system (5.1) has another positive equilibrium $\tilde{E}_3(\tilde{x}_3, \tilde{y}_3)$. By computation, $\tilde{x}_3$ satisfies the following equation:

where

Note that $\tilde{x}_3 < 1$ is a unique positive root of $\Phi(x)$, which implies that $\alpha k < 1$ and $\tilde{q_0} < q$. Also, substituting $x = 1$ into $\Phi(x)$, we have

that is, $\alpha k < 1$ and $q < \tilde{q}_1$.

The Jacobian matrix of system (5.1) at $\tilde{E_4}$ is

and

where

We have the following results.

Theorem 5.3. Assuming that $\alpha k < 1$ and $\tilde{q_0} < q < \tilde{q_1}$, system (5.1) has the equilibrium $\tilde{E}_4(1, 1)$. Moveover,

1) $\tilde{E_4}(1, 1)$ is a stable hyperbolic node or focus if $q < \tilde{q}$;

2) $\tilde{E_4}(1, 1)$ is an unstable hyperbolic node or focus if $q > \tilde{q}$;

3) $\tilde{E_4}(1, 1)$ is a fine focus or center if $q = \tilde{q}$.

Now, we will study the Hopf bifurcation around $\tilde{E_4}$ in system (5.1). Obviously, the transversality condition

holds. Then, we can determine the stability of the limit cycle around $\tilde{E_4}$ by calculating the first Lyapunov number. First, using the transformation $(\tilde{x}, \tilde{y}) = (x-1, y-1)$, the Taylor expansion of system (5.1) at the origin takes the following form:

where

Next, using the transformation $(\tilde{u}, \tilde{v}) = \left(-\tilde{x}, \frac{\tilde{a}_{10}\tilde{x}+\tilde{a}_{01}\tilde{y}}{\sqrt{D}}\right)$, where $D = \tilde{a}_{10}\tilde{b}_{01}-\tilde{a}_{01}\tilde{b}_{10} = DetJ_{\tilde{E_4}} > 0$, system (5.2) becomes

where the coefficients are given in Appendix B.

According to the results of ^[26], the first-order Lyapunov number can be written as

where

Thus, we can obtain the following theorem about the Hopf bifurcation.

Theorem 5.4. If $\alpha k < 1$, $\tilde{q}_0 < q < \tilde{q}_1$ and $q = \tilde q$, then the following statements hold.

1) If $l_1 > 0$, then system (5.1) undergoes subcritical Hopf bifurcation and an unstable limit cycle comes out around $\tilde{E_4}$.

2) If $l_1 < 0$, then system (5.1) undergoes supercritical Hopf bifurcation and a stable limit cycle appears around $\tilde{E_4}$.

3) If $l_1 = 0$, then system (5.1) undergoes a degenerate Hopf bifurcation and multiple limit cycles may appear around $\tilde{E_4}$.

By numerical simulation, we show the existence of limit cycles. Letting $k = 0.1, \alpha = 1, d = 1, c = 1, s = 1, q = 8$ and $r = 7.7$, we have that $l_{1} = 0.001984126984$. We perturb $q$ to $q = 8-0.005$; then, there exists an unstable limit cycle around $\tilde{E_4}$ (see Figure 5(a), (b)). On the other hand, letting $k = 0.1, \alpha = 1, d = 1, c = 1, s = 0.7, q = 6.8$ and $r = 7.04$, we obtain that $l_{1} = - 0.06095323795$. We perturb $q$ to $q = 6.8+0.03$; then, there exists a stable limit cycle around $\tilde{E_4}$ (see Figure 5(c), (d)).

Now, we give an example to illustrate the existence of two limit cycles. The parameters are given as follows:

where $l_1 = 0$. We perturb $k$ and $q$ to $k = \frac{7}{116}+\frac{3\sqrt{57}}{116}+0.03$ and $q = \frac{9}{2}+0.01$. Hence, system (5.1) undergoes a degenerate Hopf bifurcation and has two limit cycles (the inner one is stable and the outer is unstable) around $\tilde{E_4}$ (Figure 5(e), (f)).

Remark 5.1. In Figure 5(a), (b), the origin is a stable node and the boundary equilibria are unstable. In addition, system (5.1) has two positive equilibria, where $\tilde{E}_3$ is a saddle point and $\tilde{E}_4$ is a stable point, and an unstable limit cycle appears around $\tilde{E}_4$. The orbits of the phase portraits reveal that the prey and predator tend to coexistent in steady states only when the initial values of system (5.1) lie inside the unstable limit cycle; otherwise, the prey and predator become extinct.

In Figure 5(c), (d), in addition to the origin being stable, the other two boundary equilibria and the two positive equilibria are unstable. System (5.1) has a stable limit cycle that appears around $\tilde{E}_4$. When the initial values lie to the right of the two stable invariant manifolds of the saddle, the prey and predator tend to coexist in periodic orbits. In addition, when the initial values lie to the left of the two stable invariant manifolds of the saddle, the prey and predator tend to go extinct.

Figure 5(e), (f) show that system (5.1) undergoes a degenerate Hopf bifurcation and has two limit cycles (the inner one is stable and the outer is unstable) around $\tilde{E}_4$. Prey and predator will oscillate and coexist if the initial values lie inside of the unstable limit cycle, while the prey and predator will become extinct if the initial values lie outside of the unstable limit cycle.

5.3.   Bogdanov-Takens bifurcation

From Theorem 4.5(1), the unique positive equilibrium $E_{*}$ of system (1.4) is a cusp of codimension two, which means that a Bogdanov-Takens bifurcation of codimension two may occur. Hence, using $q$ and $s$ as the bifurcation parameters, system (1.4) becomes

where $\lambda = (\lambda_1, \lambda_2)$ is a parameter vector in a small neighborhood of the origin.

Theorem 5.5. Assuming that the conditions of Theorem 4.5 (1) hold, system (1.4) undergoes a Bogdanov-Takens bifurcation of codimension two around $E_{*}$.

Proof. First, by initiating the transformation $x_1 = x-x_{*}$ and $y_1 = y-y_{*}$ to move the positive equilibrium $E_{*}$ to the origin, system (5.4) becomes

where

Second, letting

system (5.5) can be written as follows:

where

Taking a new time variable $\tau$ with $\mbox{d}t = (1-j_{02}x_2)\mbox{d}\tau$ and $x_3 = x_2, y_3 = (1-j_{02}x_2)y_2$, system (5.6) becomes

where

From the proof of Theorem 4.4, we have

where $T_1$ is defined in Theorem 4.4. Letting

system (5.7) becomes

where

Next, letting $x_5 = x_4- \frac{m_{10}}{2}$ and $y_5 = y_4$, system (5.8) is equivalent to the following system:

where

From the proof of Theorem 4.4, we obtain

Finally, letting

we obtain the universal unfolding of system (5.4) as follows:

where

Using Maple software, we have

By the results in ^[28], system (1.4) undergoes a Bogdanov-Takens bifurcation of codimension two. The proof is completed.

The local expression of the bifurcation curves are given in ^[28] as follows:

(ⅰ) The saddle-node bifurcation curve

(ⅱ) The Hopf bifurcation curve

(ⅲ) The homoclinic curve

In what follows, we present the phase diagrams of system (5.4), as obtained by some numerical simulations. Choosing $c = 2, d = 1, s = 4, q = 20, k = \frac{1}{6}$ and $r = \frac{49}{3}$, and from Theorem 4.5(1), $E_{*}(1, 1)$ is a cusp of codimension two. Figure 6 shows that system (1.4) undergoes a Bogdanov-Takens bifurcation of codimension two.

6.   Conclusions

In this paper, we consider a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting. Fear of predator and nonlinear harvesting are the main factors affecting the dynamic behavior of system (1.4). Via numerical simulations, we show the influences of the fear effect and nonlinear harvesting on the dynamic behavior of system (1.4).

First, let $c = 0.1, d = 0.1, s = 0.04, q = 2$ and $r = 2$. When $k = 0.06$, system (1.4) has no positive equilibrium (see Figure 7(a)). In a biological sense, the prey and predator will become extinct when the fear effect is large. When $k = 0.053$, system (1.4) has two positive equilibria, where $E_3$ is a saddle and $E_4$ is an unstable node (see Figure 7(b)). Hence, in this case, the prey and predator are still extinct. When $k = 0.05$, $E_4$ becomes a stable node, and an unstable limit cycle appears around $E_4$ (see Figure 7(c)). Then, for system (1.4), a bistable phenomenon occurs, in which the prey and predator tend to steady states (or extinction), depending on the initial values lying inside (or outside) of the unstable limit cycle. When $k = 0.03$, $E_4$ is still a stable node and an unstable limit cycle disappears (see Figure 7(d)). Then, the prey and predator will survive or become extinct depending on the two stable manifolds of the saddle that act as a separatrix curve. When $k = 0$, that is, without the fear effect, the dynamic behavior of system (1.4) is similar to that shown in Figure 7(d) (see Figure 7(e)). Figure 7 shows that the prey and predator may survive or become extinct when the fear effect is small. With the increase of the fear effect, the survival area of species decreases, until finally, the prey and predator will become extinct if the fear effect is strong enough. Hence, a strong fear effect is not conducive to the survival of the species.

Second, we consider the impact of nonlinear harvesting on system (1.4). Let $c = 0.1, d = 0.1, s = 0.1, k = 0.1$ and $r = 3$. When $q = 3.5$, system (1.4) has no positive equilibrium and the origin is globally asymptotically stable (see Figure 8(a)). When $q = 3.3$, system (1.4) has two positive equilibria, where $E_3$ is a saddle and $E_4$ is an unstable node (see Figure 8(b)). In this case, the prey and predator are still extinct. When $q = 3.287$, system (1.4) has an unstable limit cycle and there is a bistable phenomenon (see Figure 8(c)). That is, an unstable limit cycle acts as a separatrix curve, where the prey and predator will become extinct or survive. When $q = 3.283$, there exists an unstable homoclinic loop in system (1.4) (see Figure 8(d)). When $q = 3.26$, the unstable limit cycle and homoclinic loop disappear. Hence, the prey and predator will tend to steady states (or extinction) if the initial values lies to the right (or left) of the two stable manifolds of the saddle (see Figure 8(e)). When $q = 0$, that is, without nonlinear harvesting, system (1.4) has only one positive equilibrium, which is globally asymptotically stable (see Figure 8(f)). This shows that overfishing can lead to the extinction of the predator and prey, so, maintaining proper harvesting can help the survival of the prey and predator.

By conducting numerical simulations, we were able to clearly observe that, when $k\le0.05$ and $q\le3.281$, the prey and predator tend to coexist around the stable positive equilibrium $E_4$. In other words, by effectively controlling the harvesting, we can ensure that the prey's fear of being caught remains within a smaller range, which benefits the survival of both populations. That is, weaker fear effects and less capture are beneficial to the survival of both predator and prey. We have conducted a theoretical analysis of system (1.4) and obtained some conclusions. However, when it comes to solving practical problems, there are many external factors. The actual application of the model may be difficult to achieve in the short term.

When $r-d > q$, the origin is a repeller, the only boundary equilibrium is a saddle point and the unique positive equilibrium may be stable or unstable. From Remark 4.1, the unique positive equilibrium is globally asymptotically stable if $q < 1$. Then, the prey and predator will tend to a positive coexistent steady state if the birth rate of the prey is high and the catchability coefficient is small. When $r-d = q\leq q_0$, from Remark 4.2, the origin is globally asymptotically stable, which implies that the prey and predator will become extinct. When $r-d = q > q_0$ or $r-d < q$, system (1.4) may have zero, one or two positive equilibria, and these equilibria may be stable or unstable. We show that the unique equilibrium $E_{*}$ is a saddle-node or a cusp of codimension two (or three). Moveover, system (1.4) undergoes saddle-node bifurcation and Bogdanov-Takens bifurcation around $E_{*}$. Also, system (1.4) undergoes a degenerate Hopf bifurcation and multiple limit cycles may appear around $\tilde{E_4}$. In Figure 5, we show that system (1.4) has two limit cycles (the inner one is stable and the outer is unstable) around $\tilde{E}_4$, which implies the bistable phenomenon. That is a large amount of fear and prey harvesting are detrimental to the survival of the prey and predator. Additionally, the prey and predator will reach a steady state if the intrinsic growth rate of the prey is high and the catchability coefficient for the prey is low. However, the prey and predator will become extinct if the intrinsic growth rate for the prey and the catchability coefficient for the prey are small.

In ^[9], the authors studied the stability of the equilibria and demonstrated that there exists a limit cycle in system (1.1). Considering the Holling type Ⅱ functional response, ^[10] showed that a unique equilibrium is a cusp of codimension two and a limit cycle appears. Unlike ^[9,10], we show that a unique equilibrium is a cusp of codimension three and confirm the occurrence of Bogdanov-Takens bifurcation. We have found that system (1.4) has two limit cycles (the inner one is stable and the outer is unstable), which exhibit the bistable phenomenon. Also, we have proven that the origin and equilibrium are globally asymptotically stable under some conditions. The strong fear effect and nonlinear harvesting are not conducive to the survival of the species. These indicate that the dynamic behavior of system (1.4) is more complex than that of the systems in ^[9,10].

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

This work was supported by the Natural Science Foundation of Fujian Province (2021J01613, 2021J011032) and the Scientific Research Foundation of Minjiang University (MJY22027).

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A.   Coefficients in the proof of Theorem 4.5

Appendix B.   Coefficients of system (5.3)