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

: In this paper, under homogeneous Neumann boundary conditions, the complex dynamical behaviors of a di ff usive Leslie-Gower predator-prey model with a ratio-dependent Holling type III functional response and nonlinear prey harvesting is carefully studied. By scrupulously analyzing and comprehending the distribution of the eigenvalues, the existence and stability (balance) of the extinction and coexistence equilibrium states are determined, and the bifurcations exhibited by the system are investigated by a mathematical analysis. Additionally, based on the theoretical analysis and numerical simulation, (Harvesting rate-induced, Delay-induced), Turing-Hopf bifurcations points are derived. Our results show that delay and nonlinear prey harvesting rates can create spatially inhomogeneous periodic solutions.


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, α 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, β 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 γ 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 ∆ is the Laplace operator; we assume that the habitat of the predator and prey is a bounded domain Ω.
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 τ.
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)

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: 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.

Existence, stability and bifurcation analysis of positive equilibria
The interior and internal equilibria are E 1 * = (u 1 * , v 1 * ) and E 2 * = (u 2 * , v 2 * ), where u 1 * and u 2 * are the roots of the following quadratic equation: and , 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- , then no interior equilibrium point exist.
> h > β(1 − α) and α + β < 1, then (a) The equilibrium point E 1 * = (u 1 * , v 1 * ) is always and forever a saddle point; (b) The equilibrium point E 2 * = (u 2 * , v 2 * ) is stable and steady when γ > γ 0 = 1 − 2u 2 * − 2αc 1+c − βh (β+u 2 * ) 2 , which is unstable when γ < γ 0 = 1 − 2u 2 * − 2αc 1+c − βh (β+u 2 * ) 2 ; (c) The system (1.1) undergoes and experiences a Hopf bifurcation with enough respect esteem and value to the bifurcation parameter γ around the equilibrium point 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 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 * <ũ < u 2 * , hence, the equilibrium point E 1 * = (u 1 * , v 1 * ) is always and at all time a saddle point, and the equilibrium point E 2 * = (u 2 * , v 2 * ) is very steady and stable when γ To discuss its fixity, stability, steadiness of the positive equilibrium E 2 * of system (1.1) more accurately and intuitively, the mathematical relation between γ 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 = {(γ, h)|γ 0 < γ} of the positive and nonnegative equilibrium E 2 * = (u 2 * , v 2 * ) to the system (1.1), and moreover T 0 (h, γ 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 γ. As a matter of fact, the parameter γ 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, ω = ± √ J 0 , if we choose or select to treat the parameter γ as a bifurcation parameter. Additionally, the parameter γ is γ = γ 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 γ crosses through γ 0 if the cross-sectional condition is met.

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 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 × 2 identity matrix and P k = −k 2 diag{d 1 , d 2 }, k ∈ N 0 , which can imply that
Consequently, we are able to receive numerous Hopf bifurcation branching lines H k as follows

The diffusive predator-prey model with delay
We assume and posit that λ = iω, substitute iω into (3.3), and separating the real part from the imaginary part, fancy and notional part when the parameter τ 0, we can get Then, we obtain the roots of (3.7) are We present the following hypothesis and assumptions (H1) (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 ±iω + k at τ j+ k , in the event of (H2) is satisfied, with (H3) is satisfied, the Eq (3.7) has two positive roots, whereupon the Eq (3.3) has a pair of purely virtually imaginary roots ±iω ± k at τ j± k , with (3.10) Permitting λ(τ) = ν(τ) + iδ(τ) be the roots of the Eq (3.3) near from τ = τ j± k which is satisfying ν τ j± k = 0, δ τ j± k = ω ± k . After that, we can get transversality condition as following.
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 ≤ β(1 − α) and γ > γ 0 hold, ω j± k and τ j± k 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 τ * = min k, j τ j± k . (a) The positive equilibrium E * (u * , v * ) of system (1.3) is asymptotically and steadily stable for the parameter τ ∈ (0, τ * ); (b) System (1.3) is at the receiving end the Hopf bifurcations drawing close to the positive equilibrium E * (u * , v * ) at τ j+ k or τ j− k ( j ∈ N 0 ); (c) System (1.3) undergoes a Hopf-Hopf bifurcation approaching the positive equilibrium E * (u * , v * ) at τ j+ k = τ j− k ( j ∈ N 0 ). 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).

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.

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.