Dynamics study of nonlinear discrete predator-prey system with Michaelis-Menten type harvesting

1.   Introduction

Population is a collection of individuals of the same species living together within a certain spatial range. The study of population dynamics mainly describes the dynamic relationship between population communities in predation systems and food web systems. Food web refers to the complex relationship network between predation, competition, cooperation and reciprocity between biological populations in biological systems. The study of food web population dynamics can give humans an important understanding of the basic nature of ecosystems, promoting the evolution of life, protecting the natural environment, and maintaining ecological balance at the macro level, and give certain guiding opinions on the protection of endangered species at the micro level. Therefore, studying the interactions between different biological populations through mathematical modeling is an important area of research for ecosystems.

In 1976, May showed in ^[1] that the first-order discrete equation model, although simple, can produce a series of extremely complex dynamic behaviors, such as from stable points to unstable bifurcation levels, and eventually produce chaos. However, in continuous-time models, achieving such complex dynamic behavior requires differential equations of three dimensions and more than three dimensions to cause chaos ^[2]. It can be seen that discrete systems described by the difference equation are richer in dynamical phenomena. Consequently, the discrete dynamical system model has attracted the attention and research of more scholars ^{[3,4,5,6,7,8,9,10,11]}.

In ^[12], Gan et al. studied the stability in a simple food chain system with Michaelis-Menten functional responses and nonlocal delays, using the Lyapunov functional to derive sufficient conditions for global stability of positive steady state and semi-trivial steady state. In ^[13], Clark et al. described mathematical models of exploited fish stocks, assuming that certain stocks can be obtained through dynamic aggregation processes. The effects of aggregation on yield-effort relationships, abundance indices, and fishery dynamics are discussed, as well as various management approaches for these models. On the other hand, with the increasing demand for food and other resources, the exploitation of biological resources is also increasing. More ecologists are very interested in studying these types of models, such as predator-prey models, and consider the impact of exploiting (harvesting) resources to protect the sustainable use of biological resources ^{[14,15,16,17]}. To do this, they applied optimal capture and control strategies to achieve the recyclability of biological resources.

In ^[20], we studied the behavioral analysis of a class of discrete dynamical system with linear harvest rates, and obtained many complex dynamic phenomena. Considering the finite nature of resources and space, the linear capture rate has no upper bound, so we will further study the Michaelis-Menten^[13] type capture on the basis of the original, and this type of capture is a kind of harvest that gradually rises until the saturation state with the increase of the number of captured objects, which is bounded and more in line with the realistic ecological environment. So, in this paper we consider a discrete-time predator-prey model with Michaelis-Menten type harvesting in preys, which is given by

where the meanings of all parameters of model (1.1) are shown in Table 1.

The paper is organized as follows. In Section 2, we analyze the dynamics of system (1.1), including the existence and local stability of the equilibrium points, and the boundedness of system (1.1). In Section 3, using the bifurcation theory and the center manifold theorem, the flip bifurcation and Neimark-Sacker bifurcation are analyzed, and the conditions for determining the bifurcation direction and the stability of the bifurcation periodic solution are obtained. In Section 4, a feedback control strategy is used to control chaos in the system, and an optimal harvesting strategy is introduced to obtain the optimal value of the harvesting coefficient. In Section 5, we verify our analytical results through numerical simulations. In the last Section, the article is ended with a brief conclusion.

2.   Model dynamics

Lemma 1 Solutions of system (1.1) with nonnegative initial conditions remain nonnegative. If $u_0 = 0$, then $u_n = 0$ for all $n\geq0$. If $v_0 = 0$, then $v_n = 0$ for all $n\geq0$. If $u_0 > 0$ and $v_0\geq0$, then $u_n > 0$ for all $n\geq0$. If $u_0\geq0$ and $v_0 > 0$, then $v_n > 0$ for all $n\geq0$.

Proof. It can be directly demonstrated by system (1.1) structure.□

Lemma 2

(I) System (1.1) always has a trivial equilibrium point $E_{0}(0, 0)$. At this moment, the two species will go extinct.

(II) If $q < r_1m_1$, then system (1.1) always has a positive semi-trivial equilibrium point $E_{1}(\frac{-\alpha+\sqrt{\Delta }}{2}, 0)$, where $\alpha = \frac{m_{1}E}{m_{2}}-K$, $\beta = \frac{qEK-r_{1}EKm_{1}}{r_{1}m_{1}}$, $\Delta = \alpha^2-4\beta > 0$. At this time, the prey population reaches $\frac{-\alpha+\sqrt{\Delta }}{2}$, and the predator tends to go extinct.

(III) If $a > h_{2}$, then system (1.1) always has a positive semi-trivial equilibrium point $E_{2}(0, \frac{a-h_{2}}{b})$. At this time, the prey population tend to go extinct, and the predator population converge to $\frac{a-h_{2}}{b}$.

(IV) System (1.1) has a positive nontrivial equilibrium point $E^{*}(u^*, \frac{1}{b}(a+\frac{r_2du^*}{u^*+c}-h_2))$, where $a > h_2$, and $u$ is the positive solution to the quartic equation of one variabie

where

When system (1.1) has a positive interior equilibrium point $E^{*}$, two species coexist, i.e. two species do not go extinct.

Proof. Direct computation.□

Lemma 3 ^[18] The equilibrium point $(u, v)$ is called

(I) Sink (locally asymptotically stable) if $|\lambda_{1}| < 1$ and $|\lambda_{2}| < 1$;

(II) Source (locally unstable) if $|\lambda_{1}| > 1$ and $|\lambda_{2}| > 1$;

(III) Saddle if $|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$ (or $|\lambda_{1}| < 1$ and $|\lambda_{2}| > 1$);

(IV) Non-hyperbolic if $|\lambda_{1}| = 1$ or $|\lambda_{2}| = 1$.

Jacobian matrix can be evaluated at $E_{0}(0, 0)$ as

The eigenvalues of $J(E_{0})$ are $\lambda_{1} = e^{r_{1}-\frac{q}{m_{1}}}$ and $\lambda_{2} = e^{a-h_{2}}$. The results regarding dynamical behaviors are listed in Table 2.

From Table 2, we can get the following theorem.

Theorem 1. When $r_{1}m_{1} < q$ and $a < h_{2}$ are satisfied, the trivial equilibrium point $E_{0}(0, 0)$ is locally asymptotically stable.

The Jacobian matrix computed at $E_{1}(\frac{-\alpha+\sqrt{\Delta }}{2}, 0)$ is

The eigenvalues of the Jacobian are $\lambda_{1} = 1-\frac{r_{1}(\Delta^2-\alpha)}{2K}+\frac{2qEm_{2}(\Delta^2-\alpha)}{4m_{1}^2E^2+m_{2}^2(\Delta^2-\alpha)^2+4m_{1}m_{2}E(\Delta^2-\alpha)}$ and $\lambda_{2} = {\rm exp}[a+\frac{r_{2}d(\Delta^2-\alpha)}{(\Delta^2-\alpha)+2c}-h_{2}]$. The properties of semi-trivial equilibrium point $E_{1}(\frac{-\alpha+\sqrt{\Delta }}{2}, 0)$ are summarized in Table 3.

From Table 3, we can have the following theorem.

Theorem 2. When $\frac{r_{1}(\Delta^2-\alpha)}{2K}-\frac{2qEm_{2}(\Delta^2-\alpha)}{4m_{1}^2E^2+m_{2}^2(\Delta^2-\alpha)^2+4m_{1}m_{2}E(\Delta^2-\alpha)} > 2$ and $a+\frac{r_{2}d(\Delta^2-\alpha)}{(\Delta^2-\alpha)+2c} > h_{2}$ are satisfied, the semi-trivial equilibrium point $E_{1}(\frac{-\alpha+\sqrt{\Delta }}{2}, 0)$ is locally unstable.

Jacobian matrix can be evaluated at $E_{2}(0, \frac{a-h_{2}}{b})$ as

The eigenvalues of the Jacobian are $\lambda_{1} = {\rm exp}[r_{1}-\frac{r_{2}(a-h_{2})}{bc}-\frac{q}{m_{1}}]$ and $\lambda_{2} = 1-(a-h_{2})$ at semi-trivial equilibrium point $E_{2}(0, \frac{a-h_{2}}{b})$. The results of dynamical behaviors are listed in Table 4.

From Table 4, we can get the following theorem.

Theorem 3 The semi-trivial equilibrium point $E_{2}(0, \frac{a-h_{2}}{b})$ is always locally asymptotically stable when $r_{1} < \frac{r_{2}(a-h_{2})}{bc}+\frac{q}{m_{1}}$ and $0 < a-h_{2} < 2$ are satisfied.

$J|_{(u, v)}$ evaluated at the positive equilibrium point $E^{*}(u^{*}, v^{*})$ is

Then characteristic equation of $J(E^{*})$ is given by

where

Lemma 4 ^[18] Suppose that $F(\lambda) = \lambda^{2}-T\lambda+D,$ and $F(1) > 0$, $\lambda_{1}$ and $\lambda_{2}$ are roots of $F(\lambda) = 0$. Then the following results hold true:

(I) $\left|\lambda_{1}\right| < 1$ and $\left|\lambda_{2}\right| < 1$ if and only if $F(-1) > 0$ and $D < 1;$

(II) $\left|\lambda_{1}\right| < 1$ and $\left|\lambda_{2}\right| > 1$ (or $\left|\lambda_{1}\right| > 1$ and $\left|\lambda_{2}\right| < 1$) if and only if $F(-1) < 0;$

(III) $\left|\lambda_{1}\right| > 1$ and $\left|\lambda_{2}\right| > 1$ if and only if $F(-1) > 0$ and $D > 1;$

(IV) $\lambda_{1} = -1$ and $\left|\lambda_{2}\right|\neq1$ if and only if $F(-1) = 0$ and $D\neq0, 2;$

(V) $\lambda_{1}$ and $\lambda_{2}$ are complex and $\left|\lambda_{1}\right| = \left|\lambda_{2}\right| = 1$ if and only if $T^{2}-4D < 0$ and $D = 1.$

Lemma 5 ^[11] Let $E^{*}(u^{*}, v^{*})$ be the unique positive equilibrium point of system (1.1), then the following propositions hold:

$(1.)$ It is a sink if and only if

$(2.)$ It is a source if and only if

$(3.)$ It is a saddle if and only if

$(4.)$ It is non-hyperbolic if and only if

To sum up, we have the following theorem.

Theorem 4 System (1.1) at the the positive equilibrium point $E^*(u^{*}, v^{*})$ is local asymptotically stable when the conditions

and

hold.

Proof. According to Lemma 3 and Lemma 4, $E^*(u^{*}, v^{*})$ is local asymptotically stable if and only if $F(1) > 0, \; F(-1) > 0$ and $D < 0$, the conclusion of Theorem 4 obtained by calculation holds.

Lemma 6 ^[19] Assume that $u_{t}$ satisfies $u_{0} > 0$, and $u_{t+1}\leq u_{t} \exp[B(1-C u_{t})]$ for $t\in[t_{1}, \infty)$, where $C$ is a positive constant. Then $\lim\limits_{t \to \infty} \sup\ u_{t} \leq \frac{1}{B C} \exp (B-1)$.

Theorem 5 Every positive solution $\left\{\left(u_{n}, v_{n}\right)\right\}$ of system (1.1) is uniformly bounded.

Proof. Suppose that $\left\{\left(u_{n}, v_{n}\right)\right\}$ be an arbitrary positive solution corresponding to system (1.1). Then, from first part of system (1.1), one has

for all $n = 0, 1, 2, \cdots$. Suppose that $u_{0} > 0$, then according to Lemma 6, we gain

From the second part of system (1.1), we have

Assume that $v_{0} > 0$, then using Lemma 6, we obtain

That is to say that $\lim\limits_{n \to \infty} \sup\ (u_{n}, v_{n}) \leq M,$ where $M = \max\left\{M_{1}, M_{2}\right\}$. This completes the proof.

3.   Bifurcation analysis

3.1.   Flip bifurcation

The characteristic equation related to system (1.1) at the unique positive interior equilibrium point $E^{*}(u^*, v^*)$ is

where

Assume that $T^{2}(u^{*}, v^{*}) > 4D(u^{*}, v^{*})$, that is,

and $T(u^{*}, v^{*})+D(u^{*}, v^{*}) = -1,$ that is to say

Then eigenvalue of $F(\lambda) = 0$ are $\lambda_{1} = -1$ and $\lambda_{2} = 2+\Phi+\Psi-\frac{r_1u^{*}}{K}$. The condition $\left|\lambda_{2}\right|\neq1$ indicates that

Consider the following set

When the perturbation parameter changes within a small field of $A_{F}$, the system (1.1) will have flip bifurcation at $E^*$. Let parameters $(a, b, c, d, K, r_1, r_2, m_1, m_2, h_2, q, E) \in A_{F}$ and consider the following systems:

where $r^*$ is a small perturbation parameter and $\left|r^*\right|\ll1.$

Let $x = u-u^{*}$ and $y = v-v^{*}.$ Then we gain

where

Construct a nonsingular matrix $D_{1}$ and consider the following translation:

where

Taking $D_{1}^{-1}$ on both sides of Eq (3.5), we obtain

where

The center manifold theorem $W^{c}$ is applied in a small field of $r^* = 0$ at the equilibrium point $E_{0}$. Then there exists $W^{c}(0)$ as follows:

and satisfies

and we have

Therefore, consider the following map on the center manifold $W^{c}(0)$:

where

By flip bifurcation, we define two non-zero real numbers $\delta_1$ and $\delta_2$, where

Based on the above analysis, the following theorem can be obtained.

Theorem 6 System (1.1) undergoes a flip bifurcation at the positive internal equilibrium point $E^{*}(u^{*}, v^{*})$ if $\delta_1\neq0$, $\delta_2\neq0$ are satisfied and when parameter $r^*$ changes within a small field of $r_1$. Moreover, if $\delta_{2} > 0$ (resp., $\delta_{2} < 0$), then the period-two orbits that bifurcate from equilibrium point $E^{*}(u^{*}, v^{*})$ are stable (resp., unstable).

3.2.   Neimark-Sacker bifurcation

Consider the characteristic equation at $E^*$, then $\mathbb{F}(\lambda) = 0$ has two complex conjugate roots with modulus one if the following conditions are satisfied:

and

Let

When the parameter changes in a small field of $A_{NS}$, system (1.1) will have Neimark-Sacker bifurcation at the unique positive equilibrium point $E^{*}(u^{*}, v^{*}).$ Select parameter $(a, b, c, d, K, r_1, r_2, \\ m_1, m_2, h_2, q, E)\in A_{NS}$ and analyze the following system:

where $\overline{r}$ is a small perturbation parameter and $\left|\overline{r}\right|\ll1.$

Let $x = u-u^{*}$ and $y = v-v^{*}.$ Then we gain

where $W_{ij}(i = 1, 2, \; 1\leq j\leq 9)$ are given in (3.6) by substituting $r_{1}$ for $r_{1}+\overline{r}$.

The characteristic equation of system (3.9) at $(x, y) = (0, 0)$ is as follows:

where

Since parameters $(a, b, c, d, K, r_1, r_2, m_1, m_2, h_2, q, E)\in A_{NS}$, the roots of the characteristic equation are

and we have

Suppose that

In addition, it is required that $\overline{r} = 0, \; \lambda^{l}_{1, 2} \neq 1\; (l = 1, 2, 3, 4)$ which is equal to $T(0)\neq$ -2, 0, -1, 2. Because $(a, b, c, d, K, r_1, r_2, m_1, m_2, h, q, E)\in A_{NS}$, thus $T(0)\neq$ -2, 2. We only require $T(0) \neq 0, -1,$ so that

Let $\eta = \frac{T(0)}{2}, \; \omega = \frac{\sqrt{4D(0)-T^{2}(0)}}{2}$, we use the following transformation:

and system (3.10) becomes into

where

System (3.9) undergoes the Neimark-Sacker bifurcation if the following quantity $\Lambda$ is not zero

where

If $\Lambda\neq0$, Neimark-Sacker bifurcation will occur in system (1.1), and the following theorem holds:

Theorem 7 System (1.1) undergoes a Neimark-Sacker bifurcation at the positive equilibrium point $E^{*}(u^{*}, v^{*})$ if conditions (3.10) are satisfied and $\Lambda\neq0$. Moreover, if $\Lambda < 0$(resp., $\Lambda > 0$), an attracting (resp., repelling) invariant closed curve bifurcates from the steady state for $r_1 > \overline{r}$ (resp., $r_1 < \overline{r}$).

4.   Chaos control and optimal harvesting policy

4.1.   Chaos control

In this section, we will adopt the feedback control method^[23,24,25] to stabilize the chaotic orbit at an unstable equilibrium point by adding a feedback control term to the system (1.1). Therefore, system (1.1) makes the following form:

where $h(u_{n}, v_{n}) = q_{1}(u_{n}-u^{*}) + q_{2}(v_{n}-v^{*})$ is feedback controlling force, $q_{1}$ and $q_{2}$ are feedback gains. Furthermore, $\widetilde{f}(u^{*}, v^{*}) = u^{*}$, and $\widetilde{g}(u^{*}, v^{*}) = v^{*}$.

The Jacobian matrix corresponding to system (4.1) at interior equilibrium point $(u^{*}, v^{*})$ is as follows:

Thus, the characteristic equation related to $A(u^{*}, v^{*})$ is:

Let $\lambda_{1}$ and $\lambda_{2}$ be the eigenvalues of characteristic equation (4.2), then

Next, we must solve equations $\lambda_{1} = \pm1$ and $\lambda_{1}\lambda_{2} = 1$ to gain the critical stability line. At the same time, it also ensures that the absolute value $\lambda_{1}$ and $\lambda_{2}$ are less than one.

Suppose that $\lambda_{1}\lambda_{2} = 1$, then we gain

Assume that $\lambda_{1} = 1$, then we have

Assume that $\lambda_{1} = -1$, then we obtain

Thus, the stable eigenvalues lie within the triangular region with the boundaries of the straight lines $L^{'}, \; L^{''}, \; L^{'''}$. In addition, when the control parameters $q_{1}$ and $q_{2}$ take values in the triangular region, system (4.1) will not create chaos phenomena.

4.2.   Optimal harvesting policy

For the sustainable use of biological resources and the protection of the natural environment on which human beings depend. Therefore, the development of renewable resources must be reasonable and proportionate. Under the premise of achieving sustainable development of biological resources, pursue maximum yield or best economic benefits. The biological and economic equilibrium combines to form the bioeconomic equilibrium. Biological equilibrium^{[14,15,16,17]} can be obtained by solving $u_{n+1} = u_{n}, \; v_{n+1} = v_{n}$ and when the economic rent is equal to zero (meaning that the total income equals the total cost), the economic equilibrium can be obtained. If $h_1, \; p$ are the cost of harvest per unit and unit price of the prey population, respectively, then the total cost is $TC = h_1E$ and total income as $TR = \frac{pqE}{m_1E+m_2u_n}$. Then the economic rent at the moment $t$ can be expressed as $\Xi = TR-TC = (\frac{pq}{m_1E+m_2u_n}-h_1)E$. The bioeconomic equilibrium can be obtained by solving the following simultaneous equations:

At present, if $\frac{pq}{m_1E+m_2u_n} < h_1$, then we stop capturing because the cost of harvesting is greater than the revenue which also means losses. Similarly, if $\frac{pq}{m_1E+m_2u_n} > h_1$, then we will continue to capture because of the harvest cost less than revenue which means profit. In order to find the bioeconomic equilibrium point $(u^*, v^*, E^*)$ from system (4.4), we can perform the following steps: first, we solve the value of $u^*$ from the third equation, and then substitute the value of $u^*$ into the second equation to get the value of $v^*$. Finally, substitute the values of $u^*$ and $v^*$ into the first equation to get the value of $E^*$.

Next, we aim to maximize net income while maintaining ecological balance. Define the net income function $J = \sum {\rm exp}(-\delta t)(\frac{pq}{m_1E+m_2u_n}-h_1)E$, where $\delta$ is the discount rate and ${\rm exp}(-\delta t)$ is the discount factor. At the same time, we use the discrete Pontryagin maximum principle^[21] to acquire the optimal capture effort. So the optimal capture problem is

where $u_t, \; v_t$ are state variables and $E_t$ is the control variable. The Hamiltonian function of the correlation at this moment is

Where $\lambda_{1(t+1)}$ and $\lambda_{2(t+1)}$ are adjoint variables. In addition, the necessary condition for the optimal problem is that $\frac{\partial H_t}{\partial u_t} = 0$, $\frac{\partial H_t}{\partial v_t} = 0$ and $\frac{\partial H_t}{\partial e_t} = 0$ are valid at the same time. Optimal harvest $E^{*}_{t}$ is available at optimal population size level $(u^{*}_{t}, v^{*}_{t})$.

As a consequence, we first solve the value of $\lambda_{1(t+1)} = \frac{{\rm exp}(-\delta t)[pqm_2E^{*}_{t}-h_1(m_1E^{*}_{t}+m_2u^{*}_{t})^2]}{qm_2u^{*2}_{t}}$ from $\frac{\partial H_t}{\partial E^{*}_{t}} = 0$, and then substitute the value of $\lambda_{1(t+1)}$ into the second equation $\frac{\partial H_t}{\partial v^{*}_{t}} = 0$ to get the value of $\lambda_{2(t+1)} = \frac{{\rm exp}(-\delta t)[pqm_2E^{*}_{t}-h_1(m_1E^{*}_{t}+m_2u^{*}_{t})^2]r_2}{qm_2u^{*}_{t}(c+u^{*}_{t})(1-bv^{*}_{t})}$. Finally, substitute the values of $\lambda_{1(t+1)}$ and $\lambda_{2(t+1)}$ into the first equation $\frac{\partial H_t}{\partial u^{*}_{t}} = 0$ to obtain the value of $E^{*}_{t}$.

5.   Numerical simulations

This section will show the bifurcation diagram, phase diagram and maximum Lyapunov exponent diagram of system (1.1) to verify the correctness of theoretical analysis.

Suppose that the parameters $(a, \; b, \; c, \; d, \; E, \; K, \; r_2, \; q, \; m_1, \; m_2, \; h_2) = (0.4, \; 0.2, \; 2, \; 0.1, \; 0.4, \; 5, \; 0.5, \; 0.1, \; 0.1, \; 0, \; 0.01)\in A_F$, $r_1$ as the bifurcation parameter, and the initial value is (3, 1). Meanwhile, when $r_1 < 1.5$, the interior equilibrium point does not exist; when $r_1 > 1.5$, there is a unique interior equilibrium point, and when $r_1 = 1.5$, the system (1.1) has a transcritical bifurcation at the boundary equilibrium point $E_2$. According to Theorem 6, when $r_1 = 3.32$, the system (1.1) will have a flip bifurcation at the interior equilibrium point (3.188, 2.104). The bifurcation diagram and the maximum Lyapunov exponent diagram are shown from Figure 1. Combined with Figures 1 and 2, when $r_1 < 3.32$, the interior equilibrium point is stable. However, when $r_1 > 3.32$, the interior equilibrium loses its stability, and orbits with periods of 2, 4, 8 appear. As the $r_1$ increases, the maximum Lyapunov exponent value is greater than zero, and it can be seen from Figure 1(c) and Figure 2(c), (f) that system (1.1) will generate chaos.

Assume that the parameters $(a, \; b, \; c, \; d, \; E, \; K, \; r_2, \; q, \; m_1, \; m_2, \; h_2) = (1.7, \; 3.5, \; 1.2, \; 0.3, \; 2, \; 3, \; 2.6, \; 0.1, \; 5, \; 2, \; 0.01)$, $r_1$ as the bifurcation parameter, and the initial value is (2, 3). The bifurcation diagram and the maximum Lyapunov exponent diagram are shown from Figure 3. Combined with Figures 4 and 5, it is clear from the figure that for smaller values of $r_1$ the system (1.1) is stable, with the increase of $r_1$ the system (1.1) stability disappears and a flip bifurcation with a period of 2 occurs, which subsequently disappears and tends to stabilize. Furthermore, when $2.43 < r_1 < 2.59$, the equilibrium point is stable. However, when $r_1 > 2.59$, the system loses its stability, and a stable invariant loop appears. At this moment, the system (1.1) produces Neimark-Sacker bifurcation and periodic solution. When $r_1$ increases, system (1.1) produces quasi-periodic solutions and chaotic phenomena. As the $r_1$ continues to increase, the maximum Lyapunov exponent value is greater than zero, the system (1.1) will generate chaos.

Considering the parameter values $(a, \; b, \; c, \; d, \; E, \; K, \; r_2, \; q, \; m_1, \; m_2, \; h_2) = (1.7, \; 3, \; 1.2, \; 0.3, \; 1.5, \; 3, \; 1.9, \; 0.1, \; 1, \; 2, \; 0.01)\in A_{NS}$ with the initial value is (2, 3), and $r_1$ as the bifurcation parameter. According to Theorem 7, when $r_1 = 2.518$, the system (1.1) has Neimark-Sacker bifurcation at the interior equilibrium point (2.52, 0.69). Figure 6 is the bifurcation and MLE graph corresponding to $u$ and $v$ with $r_1 \in [2.2, \; 3.2]$, and Figure 7 is the phase graph related to Figure 6(a). It can be seen from Figures 6 and 7 that when $r_1 < 2.518$, the equilibrium point is stable; when $r_1 > 2.518$, the equilibrium point loses its stability, and a stable invariant loop appears. At this moment, system (1.1) arises a periodic solution. When $r_1$ increases, system (1.1) generates quasi-periodic solutions and chaotic phenomena. Furthermore, as the $r_1$ continues to increase, the maximum Lyapunov exponent value is greater than zero, the system (1.1) will generate chaos.

To verify the chaotic control theory, we analyze Figure 6 and its numerical simulation parameters. In Figure 6(c) when the bifurcation parameter $r_1 = 3.1$, the maximum Lyapunov exponent value is greater than zero, system (1.1) will produce chaos. When the $q_{1}$ and $q_{2}$ are controlled in the triangular region surrounded by three straight lines $L^{'}$, $L^{''}$, and $L^{'''}$ (see Figure 8), the chaos generated by system (4.1) will be controlled near the equilibrium point and become asymptotically stable state.

Considering the parameter values $(a, \; b, \; c, \; d, \; r_1, \; K, \; r_2, \; q, \; m_1, \; m_2, \; h_2) = (0.95, \; 1.5, \; 1.2, \; 0.3, \; 0.7, \\\; 5, \; 0.5, \; 0.15, \; 1, \; 1, \; 0.1)$ with the initial value is (3, 2), and $E$ as the bifurcation parameter. At this time, the bifurcation phenomenon of system (1.1) will not occur. As the degree of capture effort $E$ increases, the population density of prey and predator will continue to decrease and will not tend to 0 (see Figure 9).

6.   Conclusions

In this paper, we study the stability and bifurcation of equilibrium points in a discrete predator-prey model with Michaelis-Menten type harvesting. The stability analysis indicates that the model has a trivial equilibrium point, two positive boundary equilibrium points and the boundary equilibrium point is always unstable. The bifurcation analysis shows that when $r_1 = 1.5$, the boundary equilibrium point $E_2$ will have a transcritical bifurcation, and when the coexistence equilibrium $E^*$ exists and loses stability, system (1.1) will have a flip bifurcation (see Figure 1). System (1.1) has, in addition, Neimark-Sacker bifurcation occur at the interior equilibrium point $E^*$ when bifurcation parameter $r_1$ changes in $A_{NS}$ small ranges (see Figure 6). Numerical simulation reveals that when the internal growth rate of prey $r_1$ gradually increases, system (1.1) will produce periodic, quasi-periodic windows and chaos.

Finally, we analyze chaos control theory and the existence of bioeconomic equilibrium points. In order to maximize profit in a finite time, we built an optimal control problem with the harvest effort as the control parameter, and theoretically obtained the optimal value of the control variable (harvest effort). As a result, we detected that harvesting efforts for prey and predator populations had a specific value that maximizes net income.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by NNSF of China (No.12161079), the Northwest Normal University Graduate Research Grant Project (No.2022KYZZ-S114), and the Gansu Province Innovation Star Project (No.2023CXZX-325).

Conflict of interest

The authors declare there is no conflict of interest regarding the publication of this paper.