Complex dynamics induced by harvesting rate and delay in a diffusive Leslie-Gower predator-prey model

1.   Introduction

Lotka and Volterra separately proposed two differential equations to provide a description for the relationship between predators and their prey in 1925 and 1926, respectively ^[1,2]. Because of the importance and practicability of the predator-prey model in most biological problems, researchers have worked to study this dynamical behavior over the past hundred years ^[1,2,3,4,5]. From the perspective of human social relations, one key goal is to understand the economic benefits in fisheries, forestry and wildlife management, which involves the development and utilization of living and biotic resources, such as the scientific management of reproducible resources and an economic harvest of the population ^[6,7]. This produces a strong and intense motivation to further study the predator-prey model. Many scholars have extensively studied the depredator model, and the impacted factors on management and re-usability have been consulted ^{[8,9,10,11,12,13]}. In 1979, May et al. put forward two types about harvesting systems ^[14]: (i) constant food production and continuable, uninterrupted harvesting, described as the harvest biomass, where there is no concern with the plant or animal population size, and (ii) continuous efforts that yield benefits (i.e, cations yield benefits), such that the biomass gathered in crops either increases or decreases.

Continuous crop yield gathering and continuous effort reaping are not very realistic and are worse than nonlinear harvesting from the perspective of biological significance and economic benefits. This mainly has two reasons: on one hand, with a constant yield or constant effort, the harvesting rate is not always constant; on the other hand, some unrealistic characteristics and limitations are reflected in the constant-effort harvesting ^[15,16,17]. Based on the achievements mentioned above, we will rigorously consider the predator-prey model together with ratio-dependent Holling type III functional response and nonlinear prey harvesting:

where $u(t)$ represents the prey density, $v(t)$ represent the predator density, $\alpha$ represents the standard good search effort of $v(t)$ versus $u(t)$, $c$ represents the biomass conversion or consumption rate, $h$ represents the maximum harvested rate of the prey species, $\beta$ represents the number of prey captured, which is the time required to calculate the maximum probability of reaching half time of the maximum probability, and $\gamma$ is either the conversion or the consumption probability rate of prey to predator.

Inspired by the literature ^{[6,7,8,9,10,11,12]}, under the homogeneous Neumann boundary condition and circumstances, we propose to study a sort of predator-prey model with a ratio-dependent Holling type functional response and nonlinear predator-prey harvesting, which has not been performed in the existing literature:

where $u(x, t)$ represents the prey densities and $v(x, t)$ represents the predator densities at the location $x$ and at time $t$, $d_1$ and $d_2$ represent the diffusion coefficients of the prey and the predator population, respectively, and $\Delta$ is the Laplace operator; we assume that the habitat of the predator and prey is a bounded domain $\Omega$.

Based on existing research results, a realistic predator-prey model should include a space and time delay. Therefore, we sought to include a time delay, which will lead to more complex dynamical behaviors of the systems, and continue to keep on studying the dynamics of the following systems

where the delay effects are represented by a nonnegative or positive parameter $\tau$.

In this paper, with the right and proper use of the normal form and the use of the center manifold theory, we will consider a delay-induced Hopf bifurcation for the predator-prey system (1.3). This document can be summarized as follows. In Section 2, we consider the Hopf bifurcation of the system (1.1) and extensively investigate the existence of the delay-induced Hopf bifurcation for the predator-prey model with diffusion. In Section 3, we further discuss the dynamical draw near behavior of the Hopf bifurcation value induced by time delay by carefully calculating the normal and regular forms on the central manifold. In Section 4, we present numerical simulations to illustrate and expand our theoretical outcomes and results.

2.   Stability and bifurcation analysis of the system (1.1)

2.1.   Existence and stability of equilibria of triviality and semi-triviality

To better understand the dynamic behavior near the equilibrium points of system (1.1), the zero growth isoclines of the system are given by and endowed with the following formula:

The equilibrium scores and points of intersection of these zero growth isoclines. The trivial and worthless equilibrium points for the system (1.1) are as follows:

(1) The origin $E_0(0, 0)$;

(2) The equilibrium points without a predator are $E_L(u_L, 0)$ and $E_H(u_H, 0)$, where $u_L$ and $u_H$ are the roots of the following quadratic equation:

If $2\sqrt{h}-1 < \beta < \min\left\{1, h\right\}$, or $h < \beta$, then

If $2\sqrt{h}-1 < \beta < \min\left\{1, h\right\}$, then

Theorem 1. (a) The origin $E_0(0, 0)$ is a saddle point, if $h > \beta$ and unstable if $h < \beta$;

(b) The axial equilibrium point $E_L(u_L, 0)$ is all the way and invariably unstable;

(c) The axial equilibrium point $E_H(u_H, 0)$ is at all time and forever a saddle point.

2.2.   Existence, stability and bifurcation analysis of positive equilibria

The interior and internal equilibria are $E_{1*} = \left(u_{1*}, v_{1*}\right)$ and $E_{2*} = \left(u_{2*}, v_{2*}\right)$, where $u_{1*}$ and $u_{2*}$ are the roots of the following quadratic equation:

and

2.2.1.   Case 1: $h > \beta(1-\alpha).$

For the system (1.1), if $h > \beta(1-\alpha)$, then

(1) The two distinct interior equilibrium points $E_{1*} = (u_{1*}, v_{1*})$ and $E_{2*} = (u_{2*}, v_{2*})$ exist whenever $\alpha+\beta < 1$ and $h < \frac{(1+\beta-\alpha)^2}{4}$;

(2) If $\alpha+\beta < 1$ and $h = \frac{(1+\beta-\alpha)^2}{4}$, then the two interior equilibrium points $E_{1*} = (u_{1*}, v_{1*})$ and $E_{2*} = (u_{2*}, v_{2*})$, collide and conflict with each other, and are denoted by the instantaneous and saddle-node equilibrium $E = \left(\tilde{u}, \tilde{v}\right)$, where $\tilde{u} = \frac{1-\beta-\alpha}{2}$;

(3) If $h > \frac{(1+\beta-\alpha)^2}{4}$, then no interior equilibrium point exist.

Theorem 2. For the system (1.1), if $\frac{(1+\eta-\alpha)^2}{4} > h > \beta(1-\alpha)$ and $\alpha+\beta < 1$, then

(a) The equilibrium point $E_{1*} = \left(u_{1*}, v_{1*}\right)$ is always and forever a saddle point;

(b) The equilibrium point $E_{2*} = \left(u_{2*}, v_{2*}\right)$ is stable and steady when $\gamma > \gamma_0 = 1-2u_{2*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{2*}\right)^2}$, which is unstable when $\gamma < \gamma_0 = 1-2u_{2*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{2*}\right)^2}$;

(c) The system (1.1) undergoes and experiences a Hopf bifurcation with enough respect esteem and value to the bifurcation parameter $\gamma$ around the equilibrium point $E_{2*} = \left(u_{2*}, v_{2*}\right)$ when $\gamma = \gamma_0 = 1-2u_{2*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{2*}\right)^2}$.

In an effort to go deeper and investigate the exceedingly intricate dynamical behaviors for the system (1.2), we consider the dynamics of system (1.1) in detail. The accurate linearization of system (1.1) at the positive equilibrium points $E_{j*} = \left(u_{j*}, v_{j*}\right)$ is

where

The characteristic equation of (2.1) is

where

When the value of $J_0 > 0$, the equilibria $E_{j*}, j = 1, 2$ is unstable. When the value of $J_0 < 0$, the equilibria $E_{j*}, j = 1, 2$ is locally asymptotically stable if $T_0 > 0$, and the equilibria $E_{j*}, j = 1, 2$ is unstable if $T_0 < 0$. Obviously and apparently, $u_{1*} < \tilde{u} < u_{2*}$, hence, the equilibrium point $E_{1*} = \left(u_{1*}, v_{1*}\right)$ is always and at all time a saddle point, and the equilibrium point $E_{2*} = \left(u_{2*}, v_{2*}\right)$ is very steady and stable when $\gamma > \gamma_0 = 1-2u_{2*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{2*}\right)^2}$, which is unstable when $\gamma < \gamma_0 = 1-2u_{2*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{2*}\right)^2}$.

To discuss its fixity, stability, steadiness of the positive equilibrium $E_{2*}$ of system (1.1) more accurately and intuitively, the mathematical relation between $\gamma$ and $h$, which appeared and yielded in the previous equation. The Hopf bifurcation line of that system (1.1) is represented as the following:

Then, the stability region is $D = \left\{(\gamma, h)|\gamma_0 < \gamma\right\}$ of the positive and nonnegative equilibrium $E_{2*} = \left(u_{2*}, v_{2*}\right)$ to the system (1.1), and moreover $T_0(h, \gamma_0(h)) = 0$.

In the following substance, what taken as the bifurcation parameter, the existence of the Hopf bifurcation at the interior equilibrium $E_{2*}$ is the parameter $\gamma$. As a matter of fact, the parameter $\gamma$ can be looked upon as the percent conversion or the consumption rate of prey to predator, is fully represented by the predator, and plays a necessary role in determining the stability of the interior equilibrium, and in deeply impacting and influencing the existence of the Hopf bifurcation.

The equation (2.3) will have a pair of opposite and contrary imaginary eigenvalues, $\omega = \pm \sqrt{J_0}$, if we choose or select to treat the parameter $\gamma$ as a bifurcation parameter. Additionally, the parameter $\gamma$ is $\gamma = \gamma_0$. System (1.1) should be a non-constant periodic solution with a very small amplitude that diverges from the positive equilibrium point $E_{2*}$ when the parameter $\gamma$ crosses through $\gamma_0$ if the cross-sectional condition is met.

Let and allow $\lambda(\gamma) = \alpha(\gamma)+i\omega(\gamma)$ be the root of (1.1), then

Hence, $\alpha(\gamma_0) = 0$ and

This may mean that the system (1.1) will undergo and suffer from the Hopf bifurcation at $E_{2*}$ as $\gamma$ crosses through the $\gamma_0$ as if the transversal condition (2.4) is content and satisfied.

2.2.2.   Case 2: $h\leq \eta(1-\alpha)$.

Under this circumstance, only one internal equilibrium point exists and is represented by $E_* = (u_*, v_*)\equiv(u_{2*}, v_{2*})$.

Theorem 3. For the system (1.1), if $h\leq \beta(1-\alpha)$, then

(a) The equilibrium point $E_* = (u_*, v_*)$ is stable and steady when $\gamma > \gamma_0 = 1-2u_{*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{*}\right)^2}$, which is unstable when $\gamma < \gamma_0 = 1-2u_{*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{*}\right)^2}$;

(b) Concerning the bifurcation parameter $\gamma$, the system (1.1) goes through a Hopf bifurcation around the equilibrium point $E_* = (u_*, v_*)$ when $\gamma = \gamma_0 = 1-2u_{*}-\frac{2\alpha c}{1+c}-\frac{\beta h}{\left(\beta+u_{*}\right)^2}$.

In the next Section, we will continue to consider the harvesting rate-induced Turing-Hopf bifurcation and delay-induce Turing-Hopf bifurcation in order to understand the diffusive predator-prey model.

3.   Spatial-temporal dynamics for the diffusive predator-prey model

Let

The linearization of (1.3) at the positive and nonnegative equilibrium $E_* = (u_*, v_*)$ is

with

where $a_{11}, \; a_{12}, \; a_{21}$ and $a_{22}$ were already abandoned (2.2).

Hence, one can see that the characteristic equation of (3.1) is

where $I_2$ is the $2\times2$ identity matrix and $P_k = -k^2diag\{d_1, d_2\}, k\in N_0$, which can imply that

with

3.1.   The diffusive predator-prey model without delay

When the value of $\tau$ is $\tau = 0$,

in there

Consequently, we are able to receive numerous Hopf bifurcation branching lines $H_k$ as follows

3.2.   The diffusive predator-prey model with delay

We assume and posit that $\lambda = i\omega$, substitute $i\omega$ into (3.3), and separating the real part from the imaginary part, fancy and notional part when the parameter $\tau\neq0$, we can get

which is solved by

Then, we obtain

the roots of (3.7) are

We present the following hypothesis and assumptions

$(H1)\; A_k^2-2B_k-C^2 < 0, B_k^2-D_k^2 > 0$, or $\left(A_k^2-2B_k-C^2\right)^2-4\left(B_k^2-D_k^2\right) < 0$;

$(H2)\; B_k^2-D_k^2 < 0$;

$(H3)\; A_k^2-2B_k-C^2 < 0, B_k^2-D_k^2 > 0$, and $\left(A_k^2-2B_k-C^2\right)^2-4\left(B_k^2-D_k^2\right) > 0$.

Lemma 1. (a) If $(H1)$ is satisfied, the Eq (3.7) has no positive root, then the Eq (3.3) has no purely virtually imaginary root;

(b) the Eq (3.7) has one positive root, after that the Eq (3.3) has a couple of purely virtual imaginary roots $\pm i\omega_k^+$ at $\tau_k^{j+}$, in the event of $(H2)$ is satisfied, with

$(c)$ If $(H3)$ is satisfied, the Eq (3.7) has two positive roots, whereupon the Eq (3.3) has a pair of purely virtually imaginary roots $\pm i\omega_k^\pm$ at $\tau_k^{j\pm}$, with

Permitting $\lambda(\tau) = \nu(\tau)+i\delta(\tau)$ be the roots of the Eq (3.3) near from $\tau = \tau_k^{j\pm}$ which is satisfying $\nu\left(\tau_k^{j\pm}\right) = 0, \delta\left(\tau_k^{j\pm}\right) = \omega_k^\pm$. After that, we can get transversality condition as following.

Lemma 2. $\frac{d Re(\lambda)}{d \tau}|_{\tau = \tau_k^{j+}} > 0$, and $\frac{d Re(\lambda)}{d \tau}|_{\tau = \tau_k^{j-}} < 0$.

Proof. It can be proved that after distinguishing the two sides of (3.3), we can chalk up

Thus, by (3.6) and (3.8), we pose and have

□

Theorem 4. Assume that the conditions $h\leq \beta(1-\alpha)$ and $\gamma > \gamma_0$ hold, $\omega_k^{j\pm}$ and $\tau_k^{j\pm}$ is defined by (3.8) and (3.9), distinctively and respectively, and denote the minimum worth of the critical worth to be delayed and postponed by $\tau_* = \min_{k, j}\left\{\tau_k^{j\pm}\right\}$.

(a) The positive equilibrium $E_*\left(u_*, v_*\right)$ of system (1.3) is asymptotically and steadily stable for the parameter $\tau\in(0, \tau_*)$;

(b) System (1.3) is at the receiving end the Hopf bifurcations drawing close to the positive equilibrium $E_*\left(u_*, v_*\right)$ at $\tau_k^{j+}$ or $\tau_k^{j-}$ ($j\in N_0$);

(c) System (1.3) undergoes a Hopf-Hopf bifurcation approaching the positive equilibrium $E_*\left(u_*, v_*\right)$ at $\tau_k^{j+} = \tau_k^{j-}$ ($j\in N_0$).

Consider the system (1.2), fixing $d_1 = 0.05, d_2 = 0.5, \alpha = 0.3, \beta = 0.5, c = 0.2$, we can get hold of the Turing-Hopf bifurcation point $TH(\gamma^*, h^*) = TH(0.0766, 0.3296)$ in the $\gamma-h$ plane, just by a ordinary and casual calculation. (see Figure 1(A)).

Consider the system (1.3), fixing $d_1 = 0.05, d_2 = 0.5, \alpha = 0.3, \beta = 0.5, c = 0.2, h = 0.335, \gamma^* = 0.1314$, through a simple calculation, we can obtain that $\tau_* = 2.4136$, and the point which the resulting Turing-Hopf bifurcation is delayed and induced at $DTH(\gamma^*, \tau_*) = DTH(0.1314, 2.4136)$ is presented in the $\gamma-\tau$ plane. (see Figure 1(B)).

In the following Section, we are going to put out some accurate and precise numerical simulations together with dynamical analysis why it is that harvesting rate-induce Turing-Hopf bifurcation and delay-induced Turing-Hopf bifurcation of these systems (1.3).

4.   Numerical simulations

In the following section, for the sake of supporting and developing our previous analysis outcomes, we use the Matlab mathematical software to perform some meaningful numerical simulations.

4.1.   Harvesting rate-induced Turing-Hopf bifurcation

With regard to system (1.2), choosing and selecting the parameters $d_1 = 0.05, d_2 = 0.5, \alpha = 0.3, \beta = 0.5, c = 0.2$, after some simple calculations, we can easily gain the critical value $h = 0.3295, \gamma = 0.0766$ for the Turing-hopf bifurcation. The family in its stable spatial nonhomogeneous periodic solutions is divided into the positive equilibrium point $E_*\left(u_*, v_*\right)$, which are distinctly expressed as a Figure 2.

4.2.   Delay-induced Turing-Hopf bifurcation

In regard to the system (1.3), we all determined to build up the merit $d_1 = 0.05, d_2 = 0.5, \alpha = 0.3, \beta = 0.5, h = 0.335, c = 0.2$. Subsequently, a long list for accurate calculations shows and expresses that $E_*\left(u_*, v_*\right) = (0.4094, 0.3517)$, and the values $\tau_* = 0.9363, \kappa_{01} = 0.0625, \kappa_{02} = -0.1026$. As a consequence, $E_*\left(u_*, v_*\right) = (0.4094, 0.3517)$ is locally stable at the moment of the parameter $\tau\in[0, \tau_*)$. At the time the parameter $\tau$ passed through the critical and primary value $\tau_*$, $E_*\left(u_*, v_*\right) = (0.4094, 0.3517)$, loses its stability and steadiness, the Turing-Hopf bifurcation progress, and a kind of stable spatial inhomogeneous periodic solution will be branched from the positive equilibrium point $E_*\left(u_*, v_*\right) = (0.4094, 0.3517)$, as shown in Figure 3.

5.   Conclusions

Though many researchers have carefully studied the very complex dynamical behavior for a predator-prey model, there was much to discover regarding time and nonlinear harvesting, and given a series of related results, we still need to further study its high codimension bifurcation in this connection. In this subfraction, with ratio dependence and nonlinear predator-prey harvesting, one must study and discuss the spatiotemporal dynamics in the differential Holling-type functional response and the diffusion Leslie-Gower predator-prey model. Concerning this spatial model, we study the characteristics of the roots for the characteristic equation, which is also distributed over an area be equation of the linearized model in the steady-state solution; additionally, we discuss the steadiness of the linear system with the positive and negative roots. Our research shows that under certain conditions, The Turing-Hopf bifurcation is able to emerge in the studied system. We further studied the important dynamic behavior of stable spatial inhomogeneous, where it may be necessary to use the central and major manifold theorem and normal naturally form theory. It showed us that this steadiness and stability or oscillate periodically in this system crossing from the equilibrium between theoretical and numerical results would be controlled and changed by controlling the threshold effect of the nonlinear prey harvesting rate and time delay; therefore, we can easily observe the rich dynamic behavior of the system near the equilibrium point. Some numerical simulation results demonstrated that, a change of the nonlinear prey rate can induce the system to produce spatiotemporal resonance, and the reaction-diffusion system (1.2) will have stable steady spatial inhomogeneous periodic solution (Figure 2). In this reaction-diffusion model with the time delay and postpone equation (1.3), the change about time delay can also trigger the change of system stability, and the system will occur a stable spatial inhomogeneous periodic solution (Figure 3). In future work, we will study the high codimension bifurcation in the reaction-diffusion predator-prey system with time delay and with a nonlinear harvesting rate using the Hopf-Hopf bifurcation and the Turing-Turing bifurcation.

Use of AI tools declaration

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

Acknowledgments

We would like to thank the reviewers for their valuable comments and suggestions, which can significantly improve the quality of our paper indeed. This work is supported by the National Natural Science Foundation of China (No.11701208), the key Projects of Natural Science Research in Colleges and Universities in the Anhui Province (No.2022AH051948), the general Project of Natural Science Research in Colleges and Universities in Anhui Province (No.KJHS2021B14).

Conflict of interest

The authors declare that he has no known competing financial interest or personal relationship that could have appeared to influence the work reported in this paper.