Steady states and spatiotemporal dynamics of a diffusive predator-prey system with predator harvesting

1.   Introduction

The study of predator-prey systems is increasingly becoming an important topic in biology and mathematics because it helps us to better understand the connections between populations. Turing pointed out in 1952 that stable homogeneous states in reaction-diffusion systems can destabilize under certain conditions and spontaneously generate a wide variety of ordered and disordered patterns ^[1]. Populations do not remain in a fixed space for a variety of reasons, so it is relevant to introduce diffusion into the predator-prey systems. Many scholars have begun to study the effect of diffusion terms on system patterns (see ^[2,3,4,5,6]).

Among many ecological systems, the Leslie-Gower system has the following form

where $v$ and $u$ denote the densities of predators and prey, respectively; $Q(u)$ denotes the functional response; $K$ stands for the environmental capacity; $h$ measures the translation of prey food quality into predator birth rates; and $s$ and $r$ represent the intrinsic growth rate of predator and prey populations, respectively. This system has been studied by many researchers, such as the simplified Holling IV $Q(u) = \frac{\alpha u}{u^2+b}$ ^[7]. Meanwhile, the generalized Holling IV functional response can describe an ecological phenomenon: When the density of the prey population exceeds a critical value, the group defense capability of the prey population can increase, which not only does not promote the increase of the predator population, but also inhibits its increase ^[8,9]. Thus, we use the generalized Holling IV functional response function $Q(u) = \frac{\alpha u}{u^2+cu+b}$ to describe the interaction between predators and prey, where $\alpha$, $c$, $b$ are biologically meaningful positive numbers.

At the same time, according to experiments ^[10], it is known that the fear effect in prey cannot be ignored. The fear effect can have an important impact on the dynamical behavior of the system (see ^{[11,12,13,14,15]}). Also the harvesting of predators is often of great practical importance (see ^[16,17]). Thus, we can continue to add these factors to the system (1.1), and a new system can be represented as follows

where $a$, $q$, $m(0 < m < 1)$, and $E$ are biologically meaningful positive numbers. Due to the inevitability of the diffusion effect, by adding the diffusion to the system (1.2), we can obtain the system (1.3)

where $\Delta$ indicates the Laplacian operator, and $d_1$ and $d_2$ represent the diffusion rates of the prey and predators, respectively.

For simplicity, by taking the following transformations:

One can attain a new diffusive predator-prey system as follows

where $\Omega \subset \mathbb{R}^{N}(N \geq 1)$ denotes a smooth bounded domain and $\partial\Omega$ is its boundary; n and $\partial \textbf{n}$ denote outward unit normal vector and directional derivative, respectively; and $u_0(x)$ and $v_0(x)$ stand for non-negative smooth initial conditions.

Its ODE system is as follows:

whose dynamic behaviors have been studied by us ^[18]. It is obvious that the positive equilibrium points of the system (1.5) are the intersection of two intersecting curves $\frac{1-u}{1+kv}-\frac{v}{u^2+du+e} = 0$ and the straight line $\delta(\beta-\frac{v}{u}) = \lambda$ in the first quadrant. That is, they are the positive roots of the following cubic equation:

Let $\Gamma = p^{2}-4qr,$ $q = m^2-3n,$ $p = mn+9e$, $r = n^2+3em,$ where $m = (\beta-\frac{\lambda}{\delta})^2k+d-1$, $n = \beta-\frac{\lambda}{\delta}-d+e$.

The conditions for the number of positive equilibrium points have been obtained in detail in ^[18], so they are not given here, but they are the basis for subsequent numerical simulations in this paper.

The Jacobi matrix of the system (1.5) at the internal equilibrium point is not difficult to obtain

It is easy to find that the equilibrium point $E_{5}^{\ast}$ is always a saddle for the reason that $Det(J_{E_{5}^{\ast}}) < 0$. However, $Det(J_{E_{i}^{\ast}}) > 0$ at the other equilibrium point $E_{i}^{\ast}$. Therefore, the equilibrium point $E_{i}^{\ast}$ is locally asymptotically stable if the condition $a_{11}+a_{22} < 0$ is satisfied, where

and $F(u, v) = \frac{u(1-u)}{1+kv}-\frac{ uv}{u^2+du+e}$, $G(u, v) = \delta v(\beta-\frac{v}{u})-\lambda v$.

The framework of the article is organized as follows. In Section 2, the boundedness and existence conditions for solutions of the system (1.4) are given. Then, we analyze the existence and non-existence of non-constant steady states of the elliptic system corresponding to the system (1.4) in Section 3, which facilitates the determination of the existence of Turing patterns. In Section 4, we give the conditions for Turing instability and the amplitude equations at the neighborhood of the threshold of Turing instability. The existence and direction of the Hopf bifurcation are explored in Section 5. Finally, numerical simulations and short conclusions are presented in Sections 6 and 7, respectively.

2.   Boundedness and existence of solutions

In this section, the existence and boundedness conditions for solutions of the system (1.4) will be given.

Theorem 2.1. Suppose that $k, d, e, \delta, \beta, \lambda, u_0(x)\ge0, v_0(x)\ge0$, and $d_1 > 0, d_2 > 0$ in $\Omega \subset \mathbb{R}^{N}(N \geq 1)$, then solutions of the system (1.4) are unique and positive, i.e., $u > 0$ and $v > 0$ for $(u, v)\in \Lambda, t\ge 0$ and $x\in \bar{\Omega }$, where $\Lambda = \left \{ (u, v):k(1-u)(u^2+du+e)+(1+kv)^2\ge0 \right \}$. Also, if $kc^2+c-e < 0$ holds, where $c = \beta -\frac{\lambda}{\delta}$, we have

(i) $\lim\limits_{t \to \infty} sup \underset{x\in \bar{\Omega}}{max} u(\cdot, t)\le 1$, $\lim\limits_{t \to \infty} sup \underset{x\in \bar{\Omega}}{max} v(\cdot, t)\le c$,

(ii) $\lim\limits_{t \to \infty} inf \underset{x\in \bar{\Omega}}{min} u(\cdot, t)\ge 1-\frac{c}{e}(1+kc)$, $\lim\limits_{t \to \infty} inf \underset{x\in \bar{\Omega}}{min} v(\cdot, t)\ge c[1-\frac{c}{e}(1+kc)]$.

Proof. Let $(\bar{u}, \bar{v}) = (\tilde{u}(t), \tilde{v}(t))$ and $(\underline{u}, \underline{v}) = (0, 0)$, where $(\tilde{u}(t), \tilde{v}(t))$ is the unique solution to the following system

Then

and

This shows that $(\bar{u}, \bar{v}) = (\tilde{u}(t), \tilde{v}(t))$ and $(\underline{u}, \underline{v}) = (0, 0)$ are the upper and lower solutions of the system (1.4), respectively. Furthermore, we have $0\le u_0(x)\le\tilde{u}$ and $0\le v_0(x)\le\tilde{v}$. By a simple calculation, we can easily find $F_v(u, v) = -u\frac{k(1-u)(u^2+du+e)+(1+kv)^2}{(1+kv)^2(u^2+du+e)}\le0$, $G_u(u, v) = \frac{\delta v^2}{u^2}$ for $(u, v)\in \Lambda$, where $\Lambda$ is written as $\Lambda = \left\lbrace (u, v):k(1-u)(u^2+du+e)+(1+kv)^2\ge0\right\rbrace$. By the well-known conclusion in ^[19], we find that the system (1.4) is a mixed quasi-monotone system, which owns a globally defined unique solution $(u, v)$ that satisfies $0\le u(x, t)\le\tilde{u}(t)$ and $0\le v(x, t)\le\tilde{v}(t)$. Furthermore, we have $u(x, t) > 0$ and $v(x, t) > 0$ for $t\ge 0$ and $x\in \bar{\Omega }$ by the strong maximum principle.

Next, we begin to explore the boundedness of the solutions $u(x, t)$ and $v(x, t)$. We find that (i) is easily obtained by comparison principle. However, (ii) is the one that needs to be worked out. Observing the first equation of the system (1.4), we have

where $c = \beta-\frac{\lambda}{\delta}$. Using the comparison principle again, there exists $T_1 > 0$ and $\varepsilon_1 > 0$ such that $u(x, t)\ge 1-\frac{c}{e}(1+kc)+\varepsilon_1$ holds for $t > T_1$ and $x\in \bar{\Omega}$. Similarly, we still have

where $c = \beta-\frac{\lambda}{\delta}$. Thus, the comparison principle helps us to obtain that there exists $T_2 > 0$ and $\varepsilon_2 > 0$ such that $v(x, t)\ge c[1-\frac{c}{e}(1+kc)+\varepsilon_1]+\varepsilon_2$ holds for $t > T_2$ and $x\in \bar{\Omega}$. By the arbitrariness of $\varepsilon_1$ and $\varepsilon_2$, we complete the proof.□

3.   Non-existence and existence of non-constant steady states

In this section, we will focus on the elliptic equations of the system (1.4):

Next, the conditions for the existence and non-existence of the steady states of the system (3.1) will be given.

3.1.   A priori estimates

Theorem 3.1. Assume that $k, d, e, \delta, \beta, \lambda, d_1, d_2 > 0$ and $\lambda > -\delta(\frac{-1+\sqrt{1+4ke}}{2k}-\beta)$ holds, we have

$1-\frac{c}{e}(1+kc)\le u(x) \le 1$, $c[1-\frac{c}{e}(1+kc)]\le v(x) \le c$, where $c = \beta-\frac{\lambda}{\delta}$.

Proof. Let $(u(x), v(x))$ be a non-negative solution of the system (3.1), and

By maximum principle ^[20], the system (3.1) follows

and

Thus, we obtain a set of upper bounds for $u, v$

Similarly, using maximum principle again, we can derive

and

The above two inequalities indicate that

Next, we would like to explore whether an unrestricted positive lower bound exists in the priori estimates for positive solutions.

Theorem 3.2. Let $\check{d}$ be a given positive constant; then, there exists a positive constant $\hat{C}$, which depends on $k, d, e, \delta, \beta, \lambda, \check{d}$, such that the solution $(u, v)$ of the system (3.1) for $d_1, d_2\ge\check{d}$ satisfies

Proof. It is clear that $0 < u(x)\le 1, 0 < v(x)\le \beta -\frac{\lambda}{\delta}$ with the help of the maximum principle.

Denote

We have

Then

where $c = \beta-\frac{\lambda}{\delta}$ and $D_1$ is a positive constant.

Define

then, $u$ satisfies the condition

Using Harnack inequality ^[21], it can be easily pointed out that

where $D_{\ast}$ depends on $\left | z \right |_{\infty }$.

Therefore, we can get

which means $u(x_1)\ge \hat{C}$.

Furthermore,

which means $v(y_1)\ge \hat{C}(\beta-\frac{\lambda}{\delta})$.□

Next, we will prove the non-existence of non-constant steady states of the system (3.1). We assume that all eigenvalues of the operator $-\Delta$ with zero-flux boundary conditions in $\Omega$ are $0 = \mu_{0} < \mu_{1}\le\mu_{2}\le \cdots\le \mu_{j} < \infty$ and $\lim\limits_{j \to \infty} \mu_{j} = \infty$.

3.2.   The non-existence

Now, we set $\bar{u} = \frac{1}{\left | \Omega \right | }\int_{\Omega } u(x)dx, \bar{v} = \frac{1}{\left | \Omega \right | }\int_{\Omega } v(x)dx$, where $(u(x), v(x))$ is a solution of the system (3.1), and if $\phi = u-\bar{u}, \psi = v-\bar{v}$, then we have $\int_{\Omega } \phi dx = \int_{\Omega } \psi dx = 0$. Thus, we have Theorems $3.3$ and $3.4$.

Theorem 3.3. Assume $\lambda > -\delta(\frac{-1+\sqrt{1+4ke}}{2k}-\beta)$, for $\phi$ and $\psi$, the following inequalities hold

(i) $\int_{\Omega } \phi^2 dx+\int_{\Omega } \left | \bigtriangledown\phi \right |^2dx \le\frac{(1+\mu_{1})e^2\left | \Omega \right | }{16d_{1}^2(1+kc)^2(e-kc^2)^2\mu_{1}^2},$

(ii) $\int_{\Omega } \psi^2 dx+\int_{\Omega } \left | \bigtriangledown\psi \right |^2dx \le\frac{(1+\mu_{1})\delta^2c^4\left | \Omega \right | }{d_{2}^2\mu_{1}^2},$

where $\mu_{1}$ denotes the first positive eigenvalue of the operator $-\Delta$ and $c = \beta-\frac{\lambda}{\delta}$.

Proof. Using the first equation of the system (3.1) and Cauchy-Schwarz inequality, there holds

and

Through Poincaré's inequality, we have

This means

Again with the help of Poincaré's inequality, we have

$\int_{\Omega } \phi^2 dx+\int_{\Omega } \left | \bigtriangledown\phi \right |^2dx \le\frac{(1+\mu_{1})e^2\left | \Omega \right | }{16d_{1}^2(1+kc)^2(e-kc^2)^2\mu_{1}^2},$

$\int_{\Omega } \psi^2 dx+\int_{\Omega } \left | \bigtriangledown\psi \right |^2dx \le\frac{(1+\mu_{1})\delta^2c^4\left | \Omega \right | }{d_{2}^2\mu_{1}^2}.$□

Theorem 3.4. Assume that $k, d, e, \delta, \beta, \lambda, d_1, d_2 > 0$, then the system (3.1) has no non-constant steady states if $d_1 > \bar{d_1}, d_2 > \bar{d_2}$, where

Proof. Multiplying the first equation for the system (3.1) by $\phi$ and integrating it by parts

where $W_1: = \int_{\Omega } \phi\frac{u(1-u)}{1+kv}dx, \quad W_2: = -\int_{\Omega } \phi\frac{uv}{u^2+du+e}dx.$

Through Theorem 3.2, we get

and

Therefore,

Similarly, multiplying the last equation of the system (3.1) by $\psi$ and integrating it by parts

With the help of Poincaré's inequality, we finally have

where $\bar{d_1} = \frac{1}{\mu_{1}}(\frac{\delta c^2}{2\hat{C}^2}+\frac{(k+24)e^2+4e(1+2c)+4d+4+8c}{8e^2})$, $\bar{d_2} = \frac{1}{\mu_{1}}(\delta c +\frac{\delta c^2}{2\hat{C}^2}+\frac{ke^2+4e+4d+4}{8e^2})$ and $c = \beta-\frac{\lambda}{\delta}$.

Clearly, as soon as $d_1 > \bar{d_1}$ and $d_2 > \bar{d_2}$ are satisfied, there is $\bigtriangledown\phi = \bigtriangledown\psi = 0$, which means that all solutions of the system (3.1) are constant steady states.□

3.3.   The existence

To simplify the calculation process, we set $\textbf{z} = (u, v)$ and

where $\tau_i = -u_{i}^{\ast}\frac{3u_{i}^{\ast2}+(2d-2)u_{i}^{\ast}+e-d}{(1+kv_{i}^{\ast})(u_{i}^{\ast2}+du_{i}^{\ast}+e)}, \rho_i = u_{i}^{\ast}\frac{1+2kv_{i}^{\ast}}{(1+kv_{i}^{\ast})(u_{i}^{\ast2}+du_{i}^{\ast}+e)}$ and $c = \beta-\frac{\lambda}{\delta}$.

Let $T(\mu_j)$ be the eigensubspace generated by eigenfunctions corresponding to the eigenvalue $\mu_j$, $j = 0, 1, 2...$. Set

(i) $\textbf{X} : = \left\lbrace \textbf{z} \in [C^1(\bar{\Omega })]\times [C^1(\bar{\Omega })]| \frac{\partial u}{\partial \textbf{n}} = \frac{\partial v}{\partial \textbf{n}} = 0\quad on \quad\partial\Omega\right\rbrace$,

(ii) $\textbf{X}_{js} : = \left\lbrace h\phi_{js}|h\in R^2\right\rbrace$, where $\left\lbrace \phi_{js}:s = 1, ..., n(\mu_{j})\right\rbrace$ is an orthonormal basis of $T(\mu_j)$, and $n(\mu_{j}) = dimT(\mu_j)$.

Then, $\mathbf{X} = \bigoplus_{j = 1}^{\infty} \mathbf{X}_{j} \text { and } \mathbf{X}_{j} = \bigoplus_{s = 1}^{n(\mu_{j})} \mathbf{X}_{j s}$.

In addition, the system (3.1) can be re-represented as

Therefore, finding a positive solution to the above system if and only if $\mathbf{z}$ satisfies

where $(\mathbf{I}-\Delta)^{-1}$ is the inverse of $\mathbf{I}-\Delta$. Obviously, simple calculations show that

Furthermore, it is easy to know that $\zeta$ is an eigenvalue of $D_\mathbf{z}g(d_1, d_2;E_i^{\ast})$ on $\mathbf{X}_j$, which is equivalent to $\zeta(1+\mu_j)$, which is an eigenvalue of the following matrix

If we set

then we can get the following equation

where $n(\mu_{j})$ denotes the multiplicity of $\mu_{j}$.

Therefore, we only need to discuss the sign of $M_i(d_1, d_2;\mu_j)$ to obtain the index of $g(d_1, d_2;\cdot)$ at the internal equilibrium point $E_i^{\ast}$. Suppose that

then, $M_i(d_1, d_2;\mu_{j})$ has two real roots,

Next, we will discuss the existence of non-constant steady states of the system (3.1) based on the number of internal equilibrium points.

First, we consider the case when the number of internal equilibrium points is one.

Theorem 3.5. Suppose that $\Gamma > 0$ and $3u_{1}^{\ast 2}+(2d-2)u_{1}^{\ast}+e-d > 0$, then the unique equilibrium point $E_{1}^{\ast}(u_{1}^{\ast}, v_{1}^{\ast})$ is locally asymptotically stable, which implies that no non-constant steady states exist near the neighborhood of $E_{1}^{\ast}$.

Proof. Let

$\zeta$ being an eigenvalue of $J$ in $\mathbf{X}_{j}$ is equivalent to $\zeta$ being an eigenvalue of the following matrix

Thus, the characteristic equation of the above matrix is given as follows

Since $3u_{1}^{\ast 2}+(2d-2)u_{1}^{\ast}+e-d > 0$, then $\tau_1 < 0$, which means $M_1(d_1, d_2;\mu_{j}) > 0$ and $d_1\mu_{j}+d_2\mu_{j}-\tau_1+\delta c > 0$. By the Routh-Hurwitz criterion, both roots $\zeta_{j1}, \zeta_{j2}$ of the characteristic equation have negative real parts and $Re\zeta_{j1}, Re\zeta_{j2} < -\gamma$, where $\gamma$ is a positive constant. This assertion is valid.□

Clearly, for a sufficiently large $d_2$, we define $\lim\limits_{d_2 \to \infty} \mu_{+}^{(1)}(d_1, d_2) = \frac{\tau_1}{d_1}: = \mu_{\ast}^{1}$ and $\lim\limits_{d_2 \to \infty} \mu_{-}^{(1)}(d_1, d_2) = 0$.

Theorem 3.6. Suppose that $\Gamma > 0$ and $3u_{1}^{\ast 2}+(2d-2)u_{1}^{\ast}+e-d < 0$. If $\mu_{\ast}^{1}\in(\mu_{q}, \mu_{q+1})$ for some $q\ge1$, and $\sigma_{q} = \sum_{j = 1}^{q}n(\mu_{j})$ is odd, then there is a positive constant $\tilde{d}$ such that the system (3.1) has at least one non-constant positive steady state solution when $d_2 > \tilde{d}$.

Proof. It is easily found that there exists a sufficiently large $\tilde{d}^{'}$ such that $\tau_1d_2-\delta cd_1 > 0$ holds for $d_2 > \tilde{d}^{'}$. Thus, we get $\mu_{+}^{(1)}(d_1, d_2) > 0, \mu_{-}^{(1)}(d_1, d_2) > 0$. By $\mu_{+}^{(1)}(d_1, d_2)\to\mu_{\ast}^{1}, \mu_{-}^{(1)}(d_1, d_2)\to0$, then there exist two positive constants $\tilde{d} > \tilde{d}^{'}$ and $q\ge1$, and we have

when $d_2 > \tilde{d}$. In addition, there exists $d_2^{\star} > \tilde{d}$, which can satisfy

Therefore, it is valid that there exists $d_1^{\star} > d_2^{\star}$ such that $\frac{\tau_1}{d_1^{\star}} < \mu_{1}$; we get

Next, we will assume that the assertion is not valid and introduce the contradiction by means of a homotopy argument.

We define

for $t\in[0, 1]$, and consider the problem

Therefore, $\mathbf{z}$ is a solution to the above problem if and only if $\mathbf{z}$ needs to satisfy the following equation

By Theorem 3.2, we set

Then, we can get $f(\mathbf{z}; t)\ne0$ when $\mathbf{z}\in \partial\Theta$ and $0\le t\le 1$. Homotopy invariance of the Leray-Schauder degree indicates

So we have

This gets the contradiction, which shows that the assertion is correct.□

We now proceed to explore the existence of non-constant steady state solutions for the system (3.1) when the number of internal equilibrium points is three. Therefore, it is always guaranteed that $\Gamma < 0, (\beta-\frac{\lambda}{\delta})^2k+d-1 < 0, \beta-\frac{\lambda}{\delta}-d+e > 0$. Next, we give different conditions to prove it. At the same time, assume that $e < d$ is valid, in order to better distinguish the sign of $\tau_i$.

There are three cases that need to be considered: $(1)u^{\star} < u_{4}^{\ast}, (2)u_{4}^{\ast} < u^{\star} < u_{6}^{\ast}, (3)u^{\star} > u_{6}^{\ast}$ where $u^{\star} = \frac{1-d+\sqrt{(d-1)^2-3(e-d)}}{3}$.

Theorem 3.7. Suppose that $u^{\star} < u_{4}^{\ast}$ holds, if $\mu_{+}^{5}(d_1, d_2)\in (\mu_{m}, \mu_{m+1})$ and $\sigma_{m} = \sum_{j = 0}^{m}n(\mu_{j})$ is even for $m\ge1$, then the system (3.1) has at least one non-constant positive steady state solution.

Proof. Since $u^{\star} < u_{4}^{\ast}$, we get $\tau_{4}, \tau_{5}, \tau_{6} < 0$. In addition, since $\rho_{4}\delta c^2-\tau_{4}\delta c > 0, \rho_{6}\delta c^2-\tau_{6}\delta c > 0, \rho_{5}\delta c^2-\tau_{5}\delta c < 0$, we have

Next, assume that the assertion is not valid. We still use the homotopy argument to introduce a contradiction. For $0\le t\le 1$, we set

and think about the problem

Therefore, $\mathbf{z}$ is a solution to the above problem if and only if $\mathbf{z}$ needs to satisfy the following equation

Similarly to Theorem 3.6, we set

Then, we get $p(\mathbf{z}; t)\ne0$ when $\mathbf{z}\in \partial\Xi$ and $t\in[0, 1]$. Homotopy invariance of the Leray-Schauder degree shows

Therefore,

Moreover, if $t = 0$, $p(\mathbf{z}; 0) = 0$ only owns a positive solution $\mathbf{z}^{\ast}(1, \beta-\frac{\lambda}{\delta})$, by repeating the previous work, we get

This gets the contradiction, which shows that the assertion is correct.□

Theorem 3.8. Suppose that $u^{\star} > u_{6}^{\ast}$ holds, if $\mu_{+}^{5}(d_1, d_2)\in (\mu_{m}, \mu_{m+1}), \mu_{-}^{4}(d_1, d_2)\in (\mu_{p}, \mu_{p+1}), \mu_{+}^{4}(d_1, d_2)$$\in (\mu_{q}, \mu_{q+1}), \mu_{-}^{6}(d_1, d_2)\in (\mu_{s}, \mu_{s+1}), \mu_{+}^{6}(d_1, d_2)\in (\mu_{g}, \mu_{g+1})$ and $\sigma_{m} = \sum_{j = 0}^{m}n(\mu_{j}), \sigma_{q} = \sum_{j = p+1}^{q}n(\mu_{j})$, $\sigma_{g} = \sum_{j = s+1}^{g}n(\mu_{j})$ as long as there are at least two odd numbers or all even numbers, then there exists a positive constant $d^{\ast}$ such that the system (3.1) has at least one non-constant positive steady state solution for $d_2 > d^{\ast}$.

Proof. Since $u^{\star} > u_{6}^{\ast}$, we get $\tau_{4}, \tau_{5}, \tau_{6} > 0$. In addition, since $\rho_{4}\delta c^2-\tau_{4}\delta c > 0, \rho_{6}\delta c^2-\tau_{6}\delta c > 0, \rho_{5}\delta c^2-\tau_{5}\delta c < 0$, we have that there exists a positive constant $d^{\ast}$ such that there is $\mu_{-}^{i}(d_1, d_2) > 0, \mu_{+}^{i}(d_1, d_2) > 0\quad (i = 4, 6)$ when $d_2 > d^{\ast}$.

In addition, we have

Repeating the proof of Theorem 3.7, we have

Obviously,

Our assertion is ultimately proven correct.□

Theorem 3.9. Suppose that $u_{4}^{\ast} < u^{\star} < u_{6}^{\ast}$ holds, if $\mu_{+}^{5}(d_1, d_2)\in (\mu_{m}, \mu_{m+1}), \mu_{-}^{4}(d_1, d_2)\in(\mu_{p}, \mu_{p+1}),$ $\mu_{+}^{4}(d_1, d_2) \in (\mu_{q}, \mu_{q+1})$ and $\sigma_{m} = \sum_{j = 0}^{m}n(\mu_{j}), \sigma_{q} = \sum_{j = p+1}^{q}n(\mu_{j})$, then when $\sigma_{m}+\sigma_{q}$ is even, there exists a positive constant $d^{\ast}$ such that the system (3.1) has at least one non-constant positive steady state solution for $d_2 > d^{\ast}$.

Proof. Since $u_{4}^{\ast} < u^{\star} < u_{6}^{\ast}$, we get $\tau_{4} > 0, \tau_{6} < 0$. Regardless of whether the sign of $\tau_{5}$ is positive or negative, we have

Moreover, there exists a positive constant $d^{\ast}$ such that there is $\mu_{-}^{4}(d_1, d_2) > 0, \mu_{+}^{4}(d_1, d_2) > 0$ when $d_2 > d^{\ast}$. Then, we get

Repeating the proof of Theorem 3.7, we have

Obviously,

The assertion is valid. □

4.   Turing instability and weakly nonlinear analysis

In this section, we derive the conditions for Turing instability. At the same time, the amplitude equation derived from the weak linear analysis ^{[22,23,24,25,26]} facilitates the differentiation of different patterns. For simplicity, we first set the spatial region $\Omega$ as a one-dimensional interval $(0, \pi)$, then the system (1.4) becomes the following system

Now we can consider the following eigenvalue problem

which has eigenvalue $\mu_{j} = j^2$ and eigenfunction $\omega_{j}(x) = cos(jx)$ for $j\in\left\lbrace 0, 1, 2\dots\right\rbrace$.

Suppose that $\begin{pmatrix} \sigma \\ \eta \end{pmatrix} = \sum_{j = 0}^{\infty } \begin{pmatrix} A_{j}\\ B_{j} \end{pmatrix}cos(jx)$ is the eigenfunction corresponding to the eigenvalue $\xi_{j}$ of the matrix $J_{j}$, thus we have

where $I$ is the identity matrix and $J_{j} = \begin{pmatrix} \tau_{i}-j^2d_{1} & -\rho _{i}\\ \delta c^2 & -\delta c-j^2d_{2} \end{pmatrix}$, $\tau_i = -u_{i}^{\ast}\frac{3u_{i}^{\ast2}+(2d-2)u_{i}^{\ast}+e-d}{(1+kv_{i}^{\ast})(u_{i}^{\ast2}+du_{i}^{\ast}+e)}, \rho_i = u_{i}^{\ast}\frac{1+2kv_{i}^{\ast}}{(1+kv_{i}^{\ast})(u_{i}^{\ast2}+du_{i}^{\ast}+e)}, c = \beta-\frac{\lambda}{\delta}$.

Clearly, the characteristic equation has the form

where

The roots of the above equation are

Based on the characteristic equation, we quickly have the following theorem.

Theorem 4.1. If $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d > 0$, then the positive equilibrium point $E_{i}^{\ast}(u_{i}^{\ast}, v_{i}^{\ast})$ is locally asymptotically stable, expect for $E_{5}^{\ast}$.

Proof. Since $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d > 0$, then $\tau_{i} < 0$ holds. In addition, we have $\rho_{i}c-\tau_{i} > 0$, expect for $i = 5$, thus it is easy to find that $T_{j} < 0$ and $D_{j} > 0$ hold, which means $E_{i}^{\ast}$ is locally asymptotically stable.□

4.1.   Existence of Turing instability

Theorem 4.2. Suppose that $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d < 0$, $\tau_{i}-\delta c < 0$ and $d_{2} > \sigma^{1} d_{1}$ are valid, where $\sigma^{1} = \delta c\frac{-\tau_{i}+2\rho_{i} c+2\sqrt{\rho_{i} c(\rho_{i} c-\tau_{i})}}{\tau_{i}^{2}}$, then $E_{i}^{\ast}$ is unstable for the reaction-diffusion system apart from $E_{5}^{\ast}$, which implies that the system suffers Turing instability at $\lambda = \lambda_{T}$ with the wave number $j^2 = j_{T}^{2} = \sqrt{\frac{\delta c(\rho_{i}c-\tau_{i})}{d_{1}d_{2}}}$, where $\lambda_{T}$ is a root that can satisfy the equation $\left(-\delta^{2} d_{1}^{2}+4 d_{1} d_{2} \delta \rho_{i}\right) (\beta-\frac{\lambda}{\delta})^{2}-2 d_{1} d_{2} \tau_{i}(\beta\delta-\lambda)-d_{2}^{2} \tau_{i}^{2} = 0$.

Proof. Since $\rho_{i}\delta c^2-\tau_{i}\delta c > 0$, expect for $i = 5$, we can get that the internal equilibrium point $E_{i}^{\ast}$ of the ODE system is locally asymptotically stable if $\tau_{i}-\delta c < 0$, i.e., $a_{11}+a_{22} < 0$. Next, we will explore the conditions for Turing instability under the assumption that $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d < 0$ and $\tau_{i}-\delta c < 0$, i.e., $a_{11} > 0$ and $a_{11}+a_{22} < 0$. It is well known that we only need to ensure that there exists $j\ge1$ such that $D_{j} < 0$, which causes the equation $\chi(\xi)$ to have a positive real root and a negative real root.

Obviously, if

then $D_{j}$ will reach a minimum value

where $j^{\ast2} = \frac{\tau_{i}d_{2}-\delta cd_{1}}{2d_{1}d_{2}} > 0$.

Let $\sigma = \frac{d_2}{d_1}$ and

we have the following equivalent condition

It is not difficult to find that

because $\rho_{i} c-\tau_{i} > 0$ and $\rho_{i} > 0$.

Therefore, $G(\sigma)$ has two positive real roots, which are

In addition, we have $G(\frac{\delta c}{\tau_{i}}) < 0$, which means $\sigma^{2} < \frac{\delta c}{\tau_{i}} < \sigma^{1}$. So, from the above analysis, $d_2 > \sigma^{1} d_{1}$ can be obtained.□

Remark 4.1. To verify the validity of the Theorem $4.2$, the relation between $Re(\xi_{j})$ and wave number $j$ is described in Figure 1. It is easy to find that there exists a threshold $\lambda = \lambda_{T} = 0.4004$ for Turing instability when $k = 0.63, d = 9, e = 0.01, \delta = 0.1, \beta = 9, d_{1} = 0.118, d_{2} = 0.6$. This means that we should control $\lambda < \lambda_{T}$ to induce the Turing instability.

Remark 4.2. We fix the parameters $k = 0.8, d = 0.8, e = 0.2, \delta = 0.4, \beta = 6, \lambda = 1.6, d_{1} = 0.15$ to explore the effect of the diffusion coefficient $d_{2}$ on Turing instability; then, we can obtain $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d = -0.6066 < 0$ and $\tau_{i}-\delta c = -0.6078 < 0$, which means that the positive equilibrium point $E_{i}^{\ast}$ is locally asymptotically stable in ODE system (1.5). Furthermore, if $d_{2} = 7$ is valid, then $\sigma^{1} = 66.3106$ and $d_{2}-\sigma^{1}d_{1} = -2.9466 < 0$, so the positive equilibrium point $E_{i}^{\ast}$ remains locallyasymptotically stable in PDE system (see Figure 2). If $d_{2}$ increases to $20$, then $\sigma^{1} = 66.3106$ and $d_{2}-\sigma^{1}d_{1} = 10.0534 > 0$, so the positive equilibrium point $E_{i}^{\ast}$ will become unstable in PDE system (see Figure 3).

4.2.   Weakly nonlinear analysis for Turing pattern

In this subsection, in order to distinguish the different patterns, multiple-scale analysis will be used to obtain the amplitude equation near $\lambda = \lambda_{T}$.

Since the weak linear analysis method will be used next, we need to rewrite the system at the positive equilibrium point $E_{i}^{\ast}(u_{i}^{\ast}, v_{i}^{\ast})$ and still use $(u, v)^{T}$ to represent the perturbation solution $(u-u_{i}^{\ast}, v-v_{i}^{\ast})^{T}$ of the system (1.4)

where

Let $\mathbf{U} = \begin{pmatrix} u\\ v \end{pmatrix}$, then the above system can be equivalently written as

where $P$ and $Q$ are linear and nonlinear operators, respectively,

with $n_{11} = -\frac{d\tau_{i}}{d\lambda }\bigg|_{\lambda = \lambda _{T}}$ and $n_{12} = \frac{d\rho_{i}}{d\lambda }\bigg|_{\lambda = \lambda _{T}}$.

Then, the variables are expanded using the small parameter $\varepsilon$:

where

Substituting the variables from the above expansion into the equation and combining the terms about $\varepsilon$, for the order $\varepsilon$, $\varepsilon^2$ and $\varepsilon^3$, we have

By solving the equation, we obtain that: $\begin{pmatrix} u_{1}\\ v_{1} \end{pmatrix} = \begin{pmatrix} \phi \\ 1 \end{pmatrix}(\sum_{j = 1}^{3}W_{j}exp(i\mathbf{k_{j}}\cdot \mathbf{r}) +c.c)$ where $\phi = \frac{\rho_{i}(\lambda_{T})}{d_{1}j_{T}^{2}-\tau_{i}(\lambda_{T})}$, $W_{j}$ denotes the amplitude of $exp(i\mathbf{k_{j}}\cdot \mathbf{r})$ under first-order perturbation, and its form is controlled by the higher-order perturbation term.

According to the Fredholm solvability condition, the vector function at the right end of the equation needs to be orthogonal to the zero eigenvectors of the operator $L_{T}^{+}$, which can guarantee the existence of nontrivial solutions. A simple calculation gives the zero eigenvector of $L_{T}^{+}$ as

where $\varphi = \frac{d_{1}j_{T}^{2}-\tau_{i}(\lambda_{T})}{\delta c^2(\lambda_{T})}$. Applying again the Fredholm solvability condition, one has

with $s_{1} = \frac{F_{uu}}{2}\phi^2+F_{uv}\phi+\frac{F_{vv}}{2}, s_{2} = \frac{G_{uu}}{2}\phi^2+G_{uv}\phi+\frac{G_{vv}}{2},$ and

where

Substituting $\begin{pmatrix} u_{2}\\ v_{2} \end{pmatrix}$ into the equation, the Freudian solvability condition shows

where

and the other two equations can be obtained by transforming the subscripts of $V$ and $W$.

Since the amplitude can be expressed as follows

we obtain the amplitude equation by combining the variables

with

Next, we can decompose the amplitude into $D_{j} = \varrho_{j}exp(i\theta _{j})(j = 1, 2, 3)$, where $\varrho_{j}$ and $\theta _{j}$ denote the mode length and phase angle, respectively. Then, we will substitute it into the amplitude equation, and separate the real and imaginary parts of the equation; thus, we can yield the following result

with $\theta = \theta _{1}+\theta _{2}+\theta _{3}$. Meanwhile, we are only interested in the stable solutions of equations, then we get the following equations

On the basis of the theory of ^[23], the above amplitude equations have four solutions, which will imply four different patterns.

(i) Spotted pattern:

It always exists and is stable if $\xi < \xi_{2} = 0$ and unstable if $\xi > \xi_{2} = 0$.

(ii) Stripe pattern:

It exists when $\xi > 0$ and is stable if $\xi > \xi_{3} = \frac{s^2I_{1}}{(I_{2}-I_{1})^2}$ and unstable if $\xi < \xi_{3} = \frac{s^2I_{1}}{(I_{2}-I_{1})^2}$.

(iii) Hexagonal pattern $(H_{0}, H_{\pi})$:

It exists when $\xi > \xi_{1} = \frac{-s^2}{4(I_{1}+2I_{2})}$. Furthermore, solution $\varrho^{+} = \frac{\left | s \right | + \sqrt{s^2+4\xi (I_{1}+2I_{2})} }{2(I_{1}+2I_{2})}$ is stable if $\xi < \xi_{4} = \frac{s^2(2I_{1}+I_{2})}{(I_{2}-I_{1})^2}$ and solution $\varrho^{-} = \frac{\left | s \right | - \sqrt{s^2+4\xi (I_{1}+2I_{2})} }{2(I_{1}+2I_{2})}$ is always unstable.

(iv) Mixed pattern:

It exists when $I_{2} > I_{1}, \xi > \xi_{3} = \frac{s^2I_{1}}{(I_{2}-I_{1})^2}$ and is always unstable.

5.   Existence, direction, and stability of the Hopf bifurcation

In this section, we will explore the existence of spatial Hopf bifurcation of the system (1.4) with $\Omega \in(0, \pi)$ and further derive the stability and direction of the Hopf bifurcation.

5.1.   The existence

In the previous discussion of Turing instability, we have obtained the characteristic equation of the PDE system as

where

In order to find the Hopf bifurcation value $\lambda_{H}$ and verify the transversality conditions, we need to explore whether the PDE system satisfies the following conditions ^[27], i.e., there exists $j\ge 0$ such that:

and $\eta^{'}(\lambda_{H})\ne 0$ for complex eigenvalues $\eta(\lambda)\pm \mathbf{i}\gamma (\lambda)$.

For the existence of $T_{j}(\lambda _{H}) = 0$, it is only necessary to satisfy $\tau_{i}(\lambda _{H})-\delta c(\lambda _{H})-j^2(d_{1}+d_{2}) = 0$, i.e., $\lambda _{H} = \lambda _{H}^{j} = \frac{\tau_{i}(\lambda _{H})-j^2(d_{1}+d_{2})}{\delta } > 0$. Obviously, we must make it valid so that $\tau_{i} > 0$, i.e., $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d < 0$.

Therefore, there exists a positive integer $j^{\ast}\ge 1$ such that

and it is clear that there is $T_{l}(\lambda _{H})\ne 0$. Furthermore, $\eta ^{'}(\lambda _{H}) = \frac{\tau_{i}^{'}(\lambda _{H})+1 }{2}\ne 0$ can be verified by numerical simulations. With the expression for $D_{j}(\lambda _{H})$, we can easily find out that

where $E_{5}^{\ast}$ is not considered.

Thus, if the following inequality holds

then $D_{j}(\lambda _{H}) > 0$. Similarly, we have $D_{l}(\lambda _{H}) > 0$.

Therefore, assuming $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d < 0$ and $d_{1}d_{2}-(\tau_{i}d_{2}-\delta cd_{1}) > 0$, we know that all roots of the characteristic equation for $\lambda = \lambda _{H}^{0} = \lambda_{0}$ have negative real parts except for the imaginary roots $\pm \mathbf{i} \sqrt{D_{0}(\lambda _{H}^{0})}$. However, at least one of the roots of the equation for $\lambda = \lambda _{H}^{l}\quad (l = 1, \dots, j^{\ast}-1)$ has a positive real part.

Theorem 5.1. Assuming that $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d < 0$ and $d_{1}d_{2}-(\tau_{i}d_{2}-\delta cd_{1}) > 0$ are valid, then the system (1.4) undergoes Hopf bifurcation at $E_{i}^{\ast}$ except $E_{5}^{\ast}$ with $\lambda = \lambda _{H}^{j}\quad (j = 0, \dots j^{\ast}-1)$. Furthermore, for $\lambda = \lambda _{H}^{j}\quad (j = 1, \dots j^{\ast}-1)$, the bifurcating periodic solutions are non-homogeneous, and for $\lambda = \lambda _{H}^{0} = \lambda_{0}$, the bifurcating periodic solutions are homogeneous, which means that it can coincide with the periodic solution of the ODE system.

5.2.   The direction and stability

First, make the following definition: $\mathbf{U_{t}} = R\mathbf{U}$, where $\mathbf{U} = (u, v)^{T}, R = D\Delta +J(E_{i}^{\ast}), J(E_{i}^{\ast}) = \begin{pmatrix} \tau_{i} & -\rho_{i} \\ \delta c^2 & -\delta c \end{pmatrix}, D = diag(d_{1}, d_{2}).$ Meanwhile, we set $R^{\ast }$ as the conjugate operator of $R$, which is defined as

We let $\varpi(\lambda)$ be the imaginary part of the roots of the characteristic equation

which has the following form

Meanwhile, we set

where $\varpi_{0} = \varpi(\lambda _{0}), \quad \zeta = \delta ^2c_{0}^{2}+\varpi_{0}^{2}, \quad c_{0} = \beta -\frac{\lambda_0 }{\delta }.$

It is easy to get that $\left \langle R^{\ast} \nu, \mu \right \rangle = \left \langle\nu, R \mu \right \rangle, R(\lambda _0)q = \mathbf{i}\varpi _{0}q, R^{\ast}(\lambda _0)q^{\ast} = -\mathbf{i}\varpi _{0}q^{\ast}, \left \langle q^{\ast}, \bar{q} \right \rangle = 0, \left \langle q^{\ast}, q \right \rangle = 1$, where $\left \langle \nu, \mu \right \rangle = \int_{0}^{\pi}\bar{\nu }^{T}\mu dx$ indicates the inner product. From ^[28], the complex space is decomposed into $X = X^{c}\oplus X^{s}$, where $X^{c} = \left \{ zq+\bar{z}\bar{q}|z\in C \right \}$ and $X^{s} = \left \{ w\in X|\left \langle q^{\ast}, w \right \rangle = 0 \right \}$.

For any $U = (u_{0}, v_{0})^{T}$, we have that there exist $z\in C$ and $w\in (w_{0}, w_{1})$ such that

Apparently,

Thus, the system (1.4) is represented as

where

From ^[28], $\tilde{g}$ can be written as

where $H, P$ have a complex symmetrical form, and direct calculations show that

Thus,

where

Meanwhile, we get $G_{20} = G_{11} = G_{02} = (0, 0)^{T}$ and $G(z, \bar{z}, w) = O(\left | z \right |^3, \left | z \right |\cdot \left | w \right |, \left | w \right | ^{2})$. From ^[28], it is clear that the system has a center manifold, and we can write it as

Then, we have

Thus, the system that is confined to the center manifold can be represented as

where

and

A straightforward calculation demonstrates

Therefore,

Based on the above analysis, we have the following conclusion.

Theorem 5.2. Assuming that $3u_{i}^{\ast 2}+(2d-2)u_{i}^{\ast}+e-d < 0$ and $d_{1}d_{2}-(\tau_{i}d_{2}-\delta cd_{1}) > 0$ are valid, then the system (1.4) undergoes Hopf bifurcation at $\lambda = \lambda _{H}^{0} = \lambda_{0}$ for $E_{i}^{\ast}$ except $E_{5}^{\ast}$.

(i) The direction of the Hopf bifurcation is supercritical (resp. subcritical) if

(ii) The bifurcating periodic solutions are unstable (resp. stable) if $Re(c_{1}(\lambda _{H}^{0})) > 0$ $(resp. < 0)$.

6.   Numerical simulations

In this section, in order to verify the correctness of the previous theoretical derivations and explore the impact characteristics of harvest on the ecological relationship of predator populations, we will perform numerical experiments.

6.1.   Turing pattern

In this subsection, we set the bounded region $\Omega = [0, 200]\times [0, 200]$, while the time step is limited to $\bigtriangleup t = 0.01$ and the spatial step is determined to be $\bigtriangleup h = 0.8$. The initial condition is set as a random perturbation at the positive equilibrium point $E_{i}^{\ast}$. Moreover, we only give different spatiotemporal pattern formations for prey $u$, since the spatiotemporal patterns of $v$ are similar to $u$. Now we will fix the following parameters

then, we have

Now it is not very difficult to get $\lambda = 0.39 < \lambda_{T}$ and $\xi = 0.02588\in(\xi_{2}, \xi_{3})$. Thus, it is easy to find from Figure 4 that the pattern formation of prey $u$ is spot patterns and stripe patterns coexisting with each other when time is short. As time continues to increase, the spot patterns dominate until the stripe patterns disappear; the spot patterns are the final form and no other structures appear. This result suggests that the prey populations ultimately form a high-density interconnected spatial distribution trend, while the predator populations have its own capture space and are not interconnected, which means that they have a separate capture space.

If we choose the parameter $\lambda = 0.35$ and $\xi = 0.12579\in(\xi_{3}, \xi_{4})$, it is obvious from Figure 5 that the whole region appears as an irregular patterns, in which the spot and stripe patterns are in competition with each other. Thereafter, as time grows, the spot and stripe patterns have a stable distribution until both coexist; finally, mixed patterns are presented in Figure 5.

Furthermore, we similarly give the formation of spot patterns of prey $u$ with the parameter $\lambda = 0.21$ and $\xi = 0.4755 > \xi_{4}$ in Figure 6. As the number of iterations increases, we can clearly observe that eventually, only spot patterns exist, which is not consistent with the theoretical analysis. This phenomenon occurs for the following reasons. In the relationship of $\xi > \xi_{4} > \xi_{3} > \xi_{2} > \xi_{1}$, the value of $\lambda$ is far away from the threshold $\lambda_{T}$, which means that some active modes will dominate compared with the primary slave modes. Consequently, they are very difficult to be adiabatically eliminated in the derivation of the amplitude equation. In addition, during the transition from uniform state modes to active modes, the amplitude equation of $D_1$ has an extra third-order term $D_0 \overline{D_2}\cdot \overline{D_3}$. Similarly, the amplitude equations of $D_2$ and $D_3$ have extra terms $D_0 \overline{D_1}\cdot \overline{D_3}$ and $D_0 \overline{D_1}\cdot \overline{D_2}$. The inclusion of these terms leads to the re-stabilization of the spot patterns, which is why the numerical simulations do not match the theoretical analysis. Similar numerical simulation results can be found in ^[29,30]. More details of the theory can be found in ^[23].

By comparing and analyzing the results of Figures 4–6, it can be concluded that the spatial distribution of prey and predator populations undergoes essential changes as the predator harvesting parameter values decrease, and the spatial distribution density of the prey populations gradually decreases. Furthermore, the final spatial pattern transitions from a spot pattern to stripe and spot mixed patterns, which means that, as the spatial distribution density of the prey populations decreases, the predator populations must spread to the predation domain in order to survive. Therefore, it is worth emphasizing that the size of the predator harvesting not only affects the predation dynamics between predatory populations, but also affects the density spatial distribution characteristics of the populations.

6.2.   Hopf bifurcation

In this subsection, we will fix the following parameters

By a simple calculation, we can get $\lambda_{H}^{0} = 0.1346$, and the parameters can satisfy the conditions in Theorem 5.2. In addition, we can get $Re(c_{1}(\lambda_{H}^{0}))\approx -1.3526 < 0$ and $\eta ^{'}(\lambda _{H}^{0}) > 0$, thus it is worth pointing out that the PDE system undergoes a supercritical Hopf bifurcation at $\lambda = \lambda_{H}^{0}$ and produces stable bifurcated periodic solutions (see Figure 7). This result means that appropriate predator harvesting behavior can promote the formation of a stable periodic growth coexistence mode between prey and predator populations, which is beneficial for their sustainable survival.

7.   Conclusions

This paper mainly proposed a predator-prey system with harvesting and diffusion to explore how harvesting affects predatory ecological relationships. Within the framework of theoretical analysis, we first give the existence of solutions of the system (1.4) by using the methods in ^[19] and boundedness by using the comparison principle, and prove that the solutions of the elliptic system (3.1) are upper and lower bounded. Second, with the help of Poincaré's inequality, the non-existence conditions for the non-constant steady states of the elliptic system (3.1) are investigated. At the same time, the existence of the non-constant steady states is analyzed by homotopy invariance of the Leray-Schauder degree. Moreover, we obtain the condition for Turing instability, and derive the amplitude equation at the threshold of Turing instability by weak linear analysis, which gives different patterns such as spot patterns, mixed patterns, and so on. Finally, the existence, direction, and stability of Hopf bifurcation are analyzed through theories like central manifolds. Under the framework of numerical simulation, we first validate the effectiveness and feasibility of the theoretical analysis results and dynamically display the trend of spatial distribution changes in population density. Second, through comparative analysis of numerical simulation results, the impact characteristics of harvesting behavior on predatory ecological relationships and spatial changes in population density are pointed out. Finally, based on the numerical simulation results, it is clear that appropriate harvesting behavior can promote the formation of a stable periodic growth coexistence mode between the predator and prey populations. Based on the above results, it can be clearly emphasized that harvesting has a significant impact on predatory ecological relationships.

One innovation of this paper is the introduction of the generalized Holling IV functional response to describe the interaction between predator and prey, which can not only enrich the dynamic properties of the system (1.4) but also make it more suitable for exploring the spatial distribution trends of prey and predator in natural ecosystems. Another innovation of this paper is to reveal the spatial coexistence mode of prey and predator during the gradual enhancement of group defense in the prey populations from the perspective of the dynamic evolution process of Turing patterns. Furthermore, it is also worth emphasizing that prey and predator populations have a stable periodic oscillation growth coexistence mode, which can indicate that appropriate harvesting behavior can not only effectively control the growth of prey populations but also maintain sustainable survival between prey and predator. This research result can be applied to the control of Microcystis aeruginosa bloom outbreaks by monitoring the its population density and of filter-feeding fish. If the dynamic change law of Microcystis aeruginosa population density is a stable periodic oscillation mode under low density, and the dynamic change law of filter-feeding fish population density is a stable periodic oscillation mode, this can indicate that filter-feeding fish can effectively control the outbreak of Microcystis aeruginosa blooms.

To emphasize the feasibility of the theoretical and numerical results in this paper, we conducted comparative analysis on the published papers separately. The paper ^[18] has thoroughly explored the bifurcation dynamics behaviors of the system (1.5), and we have compared and analyzed the research results of this paper with those of ^[18]. It is worth pointing out that the systems (1.4) and (1.5) have stable periodic solutions, which means that the predator and prey populations can eventually form a stable periodic oscillation growth coexistence mode with time. This indicates that the system (1.4) continues some of the dynamic characteristics of the system (1.5) in the time state. Furthermore, these papers ^[2,31,32,33] have obtained some excellent research results on steady states and spatiotemporal patterns of predator-prey system with generalized Holling IV functional response, Holling type I functional, and Beddington-DeAngelis functional response. Under the same theoretical analysis and numerical simulation framework, we have investigated all possible stationary distributions of prey and predator in two dimension habitats (for example, spots and mixture of spots and stripes). These research results are similar to those presented in the papers ^[2,31,32,33]. Based on the above, it is necessary to demonstrate that the theoretical and numerical results of this paper are reliable.

In summary, although this paper obtained some theoretical and numerical results in steady states and spatiotemporal patterns, there is still much to be explored in subsequent work, such as using laboratory or field monitoring data to calibrate the values of system parameters, investigating the impact of population migration behavior on the dynamic relationship between prey and predator, etc. Finally, we hope that the research results of this paper can provide some theoretical support for the control of Microcystis aeruginosa blooms.

Author contributions

Rongjie Yu, Hengguo Yu and Min Zhao: Conceptualization, Methodology, Validation, Software, Writing-original draft, Writing-review & editing. All authors contributed equally to the manuscript. All authors have read and approved the final version of the manuscript for publication.

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 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), by the Zhejiang Province College Student Science and Technology Innovation Activity Plan (New Talent Plan) (Grant No: 2024R429A010).

Conflict of interest

The authors declare there is no conflict of interest.