Multi-scale dynamics of predator-prey systems with Holling-IV functional response

1.   Introduction

Predator-prey interactions have a strong influence on population dynamics. In order to describe this effect more accurately, Lotka ^[1] and Volterra ^[2] proposed the predator-prey model, and this model was used to study the dynamics of interacting populations. Further, many more realistic mathematical models have been proposed and studied extensively. Until now, the study of predator-prey systems has been at the center of fields such as ecology and mathematical biology, see ^{[3,4,5,6,7,8,9]} and references therein. To study in depth the evolution of predator-prey interactions and population dynamics, we usually add to the nonlinear system the Holling-type functional response ^{[10,11,12,13]}, the functional form commonly used to describe the response of predators to prey. It is modeled based on the relationship between predator-prey rate and prey density, which acts as a constraint on the growth of prey density. When prey density is high, the prey rate of the predator increases, which reduces the survival and reproductive success of the prey and limits the growth of the prey population. This limiting effect helps to maintain a dynamic balance between predator and prey and prevents excessive growth in prey populations. There are many functional response functions to describe the different relationships between predator and prey, and form of the usual functional response functions are given in Table 1.

Compared with the previous three Holling-type functional response functions, the fourth Holling-type functional response function is more representative and more realistic. Now, some scholars have studied the predator model with a Holling-IV functional response function, and some basic research results have been obtained ^[14,15,16]. However, the simplified Holling-IV functional response ^[17] function parameters in the system are still relatively idealized. The populations that can satisfy this system in the reality of nature are also relatively few. Given this, this paper will study the Holling-IV type functional response function without simplification.

In recent years, predator-prey systems, with respect to functional response functions, have been extensively studied, and most of these studies have focused on the existence of positive equilibrium points, the bifurcation problem, and the existence and number of limit cycles at a single scale. Ma ^[18] et al. studied the Holling-II predator-prey system containing cross-diffusion. Steady state bifurcations and Hopf bifurcations were considered and it was found that when the diffusivity is large enough, the non-constant state will not exist. Ma ^[19] et al. applied the Lyapunov-Schmidt reduction method to study the diffusion Lotka-Volterra model with the fear effect. The effects of population growth rate and fear effect on the stability of the system were considered. The complex dynamical behavior of Holling-IV predator-prey systems, such as Bogdanov-Takens bifurcation of limit cycles, Hopf bifurcation, homoclinic bifurcation, and saddle-node bifurcation, were studied by Huang and Xiao ^[20]. Chuan ^[21] et al. studied a Holling-IV stochastic predator-prey system with antipredator behavior and noise and calculated the stability of the distribution of the system to give sufficient conditions for the continued survival and extinction of the population. Mian ^[22] et al. studied the discrete Holling-IV predator-prey model and investigated the complex dynamical properties of the system using the semi-discretization method, the central prevalence theorem, and bifurcation theory. Bin ^[23] et al. studied the Holling-III predator-prey system containing prey shelters and analyzed the Hopf bifurcation of the system and the limit cycle, as well as the fear and shelter effects.

In contrast, there are fewer multi-scale related studies on predator-prey systems. Multi-scale problems are widely used in modern science to describe various phenomena in mechanics, biology, and chemistry ^{[24,25,26,27,28,29]}. Rinzel ^[30], who separated a system coupled at multiple time scales into a fast subsystem controlled by a fast variable and a slow subsystem controlled by a slow variable, then considers the slow variable as a slow-variable parameter of the fast subsystem and investigated the bursting phenomena of the system by analyzing the bifurcation behavior of the fast subsystem. This method is also the most commonly used method to study the mechanism of bursting phenomena on different time scales. Different slow and fast systems contain significant predictions of complex oscillations ^[31].

The study of dynamic bifurcation problems, such as the variation of bifurcation parameters with time, is an important topic of study in nonlinear dynamical systems, and the dynamic Hopf bifurcation problem has many applications in biology, chemistry, and mathematics. In particular, when the system response jumps between a slowly varying quiescent state and spiking state, it passes over a bifurcation point along a channel, but instead of the spiking state starting to oscillate at the limit point, a hysteresis occurs, which is known as the bifurcation hysteresis phenomenon. Bilinsky ^[32] et al. used analytical and numerical methods to analyze the bifurcation hysteresis phenomenon of the no-turning-point case of the FitzHugh-Nagumo model new hysteresis and memory effect phenomenon. Han ^[33] discussed the periodic bifurcation hysteresis behavior when the parameter excitation passes slowly through the Hopf bifurcation value in a controlled van der Pol system, and concluded that the first bifurcation hysteresis behavior is dependent on the initial conditions, while the bifurcation hysteresis behavior after the first bifurcation hysteresis behavior is not affected by the initial conditions.

At present, studies on the bursting phenomenon of predator-prey systems mainly focus on the existence of limit cycles, approximate solutions of the singular uptake method, experimental analyses, and numerical simulations ^[34,35], with fewer analytical studies on the oscillatory mechanism of its bursting phenomenon. In this paper, we take the Holling-IV predator-prey system as the subject of our research, and consider the environmental change, i.e., the perturbation of environmental capacity, and apply the slow-fast analytical method and asymptotic theory to this slow-fast coupled nonlinear predator-prey system, aiming to reveal the mechanism of the bursting oscillation of this system and deepen the study of the spiking state hysteresis induced by the bifurcation hysteresis. It provides theoretical analyses for the prediction, regulation, and control of predator and prey populations.

2.   System stability and bifurcation analysis

2.1.   System equations

The predator-prey system with Holling-IV functional responses are given by

where $x(t)$ denotes the number of prey populations, $y(t)$ is the number of predator populations, $a$ represents the growth rate of the prey population, $K$ represents environmental capacity, $\frac{x}{\alpha + \beta x+x^2}$ represents the functional response, $d$ represents the death rate of the predator. $a$, $m$, $\alpha$, $\beta$, $\lambda$, and $d$ are all positive constants.

Let $d\theta = \frac{dt}{\alpha + \beta x+x^2}$, and it is also expressed as $dt$, and thus system (2.1) is equivalently transformed into

In ecosystems, environmental capacity does not remain static, but changes with the environment, which in turn affects the dynamics of the behavior of the entire population. For example, decreases in environmental capacity due to water pollution, pesticide use, etc. may have a dramatic effect on the stability of the entire system.

Assuming that there are slow periodic perturbations in the environmental capacity of system (2.2), system (2.2) becomes a non-autonomous system as follows

If $\omega \ll1$, system (2.3) will involve both slow and fast time scales, dividing the system into a fast subsystem controlled by the fast variables $x, y$, and a slow subsystem controlled by the slow variable $\omega t$. Let $F = l\cos(\omega t)$, and obtain the fast subsystem

Taking $F$ as the bifurcation parameter of the autonomous system (2.4), the variation of $F$ will lead to changes in the stability of system (2.4), which in turn will lead to bifurcation behavior. Meanwhile, $K+l \cos(\omega t)$ in the non-autonomous system (2.3) varies slowly and periodically within $[K-l, K+l]$, so the stability and bifurcation of the autonomous system (2.4) will affect the dynamic behavior of the non-autonomous system (2.3).

2.2.   Equilibrium points and its stability

Theorem 1. Let $\Delta = (d \beta -\lambda m)^2-4d^2 \alpha$. When $F > \frac{\lambda m -d \beta \pm \sqrt{\Delta }}{2d}-K$, $\lambda m -d \beta > 0$, the positive equilibrium of system (2.4) exists. The generalized equilibrium of the autonomous system (2.4) and its stability are shown below:

(1) The equilibrium $E_1(0, 0)$ is an unstable saddle point.

(2) The equilibrium $E_2(K+F, 0)$ is a stable node or focus when $\lambda m(K+F) < d(\alpha + \beta (K+F) +(K+F)^2)$ is satisfied, and a saddle point when $\lambda m(K+F) > d(\alpha + \beta (K+F) +(K+F)^2)$ is satisfied.

(3) The equilibrium $E_3(x_3, y_3)$ is $(\frac{\lambda m-d\beta +\sqrt{\Delta}}{2d}, \frac{a}{m}((1-\frac{x_3}{K+F}))(\alpha + \beta x_3+x_3^2))$. $E_3(x_3, y_3)$ is a saddle point.

(4) The equilibrium $E_4(x_4, y_4)$ is $(\frac{\lambda m-d\beta -\sqrt{\Delta}}{2d}, \frac{a}{m}((1-\frac{x_4}{K+F}))(\alpha + \beta x_4+x_4^2))$, $E_4(x_4, y_4)$ is an unstable focus or node when it satisfies $F > \frac{-4d^2\alpha -(3\lambda m-d\beta)[d\beta -\lambda m+\sqrt{(d \beta -\lambda m)^2-4d^2 \alpha}]}{2d[\lambda m-\sqrt{(d \beta -\lambda m)^2-4d^2 \alpha}]}-K$, and a stable focus or node when it satisfies $F < \frac{-4d^2\alpha -(3\lambda m-d\beta)[d\beta -\lambda m+\sqrt{(d \beta -\lambda m)^2-4d^2 \alpha}]}{2d[\lambda m-\sqrt{(d \beta -\lambda m)^2-4d^2 \alpha}]}-K$.

(5) The equilibrium $E_5(x_5, y_5)$ is$(\frac{\lambda m-d\beta}{2d}, \frac{a}{m}(1-\frac{x_5}{K+F})(\alpha + \beta x_5+x_5^2))$. $E_5(x_5, y_5)$ is a higher-order singularity, half-saddle, and half-node.

The proof is easy and is omitted here.

2.3.   Bifurcation analysis

^[20] analyzes the local bifurcation and global bifurcation of this system, including saddle-node bifurcation, Hopf bifurcation, homoclinic bifurcation, and bifurcation of cusp-type with codimension two (Bogdanov-Takens bifurcation). Based on this paper, we analyze autonomous system (2.4) Hopf bifurcation in detail.

Theorem 2. Taking $F$ as bifurcation parameter, the following critical conditions need to be satisfied to generate a Hopf bifurcation

Proof. The Jacobi matrix corresponding to the system (2.4) at the equilibrium $E_4(x_4, y_4)$ is

Because of $m_{22} = \lambda m x_4-d(\alpha +\beta x_4+ x_4^2) = 0$, the corresponding characteristic equation is

where

If $m_{11} = 0$, the equation $m_{11} = 0$ can be solved for Eq (2.5). When $m_{11} = 0$, the characteristic root of the system is a pair of pure imaginary roots, it is denoted as $\Lambda = \pm i \omega_0$.

Taking the derivative of $F$ on both sides of the characteristic equation yields the following formula,

The simplified expression is as follows

Substituting $\Lambda = \pm i \omega_0$ into the above formula, we get

Therefore, the Hopf bifurcation occurs when the system satisfies Eq (2.5).

Fixed parameters $a = 0.2$, $K = 2$, $m = 0.6$, $\alpha = 0.35$, $\beta = 0.35$, $\lambda = 1$, and $d = 0.1$, where there exist equilibrium points $E_1(0, 0)$, $E_2(K+F, 0)$, $E_3(x_3, y_3)$, and $E_4(x_4, y_4)$. Figure 1 gives the bifurcation diagram of the variable $x$ in system (2.4) about $F$. The equilibrium point corresponding to the solid line of the figure is stable, and the dotted line is unstable. The equilibrium line corresponding to the unstable equilibrium point $E_1$ is segment $AB$. The equilibrium $E_2$ corresponds to the $AQ$ segment. The line of equilibrium corresponding to equilibrium $E_3$ is unstable. The equilibrium line corresponding to equilibrium $E_4$ is the segment $CD$, where $PH$ segments are stable nodes or focus, and unstable focus on $HD$. The values of the parameters of the critical points are $F_P = -1.9374$, and $F_H = -1.1466$. $F_P$ is the intersection of equilibrium $E_2$ and equilibrium $E_4$. $F_H$ denotes the critical point of the Hopf bifurcation.

3.   Oscillatory behavior at different disturbance amplitudes

The case of system trajectory near the Hopf bifurcation is detailed below, and the corresponding biological explanation is given. Due to the change in the environmental capacity being slow, let the perturbation frequency be $\omega = 0.001$. The other parameters are $a = 0.2$, $K = 2$, $m = 0.6$, $\alpha = 0.35$, $\beta = 0.35$, $\lambda = 1$, and $d = 0.1$. In the following, the oscillatory behavior of the different modes of the non-autonomous system (2.3) will be investigated.

3.1.   Periodic oscillation

When there is no perturbation in the environmental capacity, system (2.3) exhibits the state of periodic oscillation. Obviously, system (2.3) is autonomous and there exists a stable limit cycle attractor, as shown in Figure 2, which moves in a counterclockwise direction as shown by the arrow in Figure 2(a). In segment $A_1A_2$, the predator becomes extinct and the prey continues to increase. In stage $A_2A_3$, the more abundant the food supply is, the easier it is for the predator to catch the prey, which in turn promotes the reproduction of the predator. But, on the other hand, the increase in the predator population leads to an insufficient supply of prey population, which in turn leads to less and less prey. In stage $A_3A_1$, the predator is unable to obtain food due to the extinction of prey, which in turn leads to a decrease in the number of predator populations.

3.2.   Quasi-periodic oscillation

When the perturbation amplitude changes, system (2.3) shows diverse oscillatory behavior. When the perturbation amplitude $l = 0.5$, the non-autonomous system (2.3) shows quasi-periodic oscillations; the phase diagram is shown in Figure 3(a), and the time history diagram is shown in Figure 3(b). Slow periodic variations in environmental capacity break the regular periodic oscillations of the original system. When $F$ is a fixed value, the phase trajectory line corresponds to a closed curve in one cycle, and when $F$ changes, the phase trajectory line enters into another closed trajectory line and does not return to its original state. The larger the value of $F$, the larger the range of variation of the close-track line. Although there are small variations in the system trajectory line from cycle to cycle, the system as a whole is in a long-term stable state as time tends to infinity.

To better explain the generation mechanism of quasi-periodic oscillations, we study the oscillatory behavior of the non-autonomous system (2.3) with the help of the stability and bifurcation of the autonomous system (2.4). We give the transformed phase diagram about the variable $x$ and slow-varying processes $F = 0.5\cos(0.001t)$ in Figure 4(a), which are superimposed with the bifurcation diagram Figure 1 to obtain Figure 4(b).

To clearly present the operation of the system trajectory, we elaborate on the motion of the trajectory in Figure 4(b) during a perturbation period. Figure 5(a) gives the operation of the trajectory of system (2.3) as it moves to the right with $F$ increasing, and Figure 5(b) gives the operation as it moves to the left with $F$ decreasing. Assuming that the system trajectory starts from point $P_1$, since point $P_1$ is located on the unstable equilibrium line corresponding to the equilibrium point $E_4$, then point $P_1$ is attracted by the stable limit cycle of the autonomous system and generates large oscillations in the spiking state. The oscillation amplitude remains the same as the amplitude of the limit cycle attractor of the autonomous system moves to the right until it reaches the maximum value of the slow variable process $F = 0.5\cos(0.001t)$, i.e., the point $T_1$. The system trajectory changes direction and moves along the stable limit cycle corresponding to the unstable focus $E_4$ in the direction of decreasing $F$. At this time, it is attracted by the stable limit cycle and remains in the spiking state. When $F$ decreases to the minimum value of $F = -0.5$, the direction of motion of the trajectory changes again to the direction of increasing $F$, and so on.

3.3.   Single-Hopf bursting oscillation and generation mechanism

Figure 6(a) and Figure 6(b) show the phase diagram and time history diagram of system (2.3) when the perturbation amplitude is taken as $l = 1.6$. The trajectory of system (2.3) trajectory exhibits the typical bursting oscillation behavior of large-amplitude spiking state oscillations (SP) coupled with micro-amplitude quiescent state oscillations (QS).

The transformed phase diagram is shown in Figure 7(a), and the superimposed with the bifurcation diagram Figure 1 is shown in Figure 7(b). Figure 8(a) shows the operation of the trajectory of system (2.3) from left to right as $F$ increases, where the region to the right of the grey dotted line is the spiking state. Figure 8(b) gives the operation from right to left along the decrease of $F$.

Assuming that the trajectory starts from point $P_2$, since point $P_2$ is on the stable equilibrium line corresponding to the equilibrium point $E_4$ of the autonomous system, it is attracted by the stable equilibrium line. It moves steadily to the right, maintaining the quiescent state. The trajectory moves to the point $H$, and due to the Hopf bifurcation of system (2.4), the stable focus becomes the unstable focus, and the stable limit cycle is generated, which leads to the breaking of the system trajectory quiescent state, and the trajectory shows a slight oscillation and slowly enters the spiking state. The trajectory moves to the maximum value of the slow variable process $F = 1.6\cos(0.001t)$, i.e., the point $T_2$. The trajectory changes direction, is attracted by the attractor of the stable limit cycle, and maintains an oscillatory state consistent with the amplitude of the limit cycle of the autonomous system in the direction of decreasing $F$. As it passes through the Hopf bifurcation point $H$ again, the stable limit cycle attractor becomes the stable focus attractor, which gradually reduces the amplitude when $F$ continues to decrease to the minimum value of $F = -1.6$, i.e., the point $P_2$. The direction of motion of the trajectory again changes to motion in the direction of increasing $F$, and so on. The presence of the Hopf bifurcation makes the system have spiking states and quiescent states in one single period of motion, and hence it is called the single-Hopf bursting phenomenon.

It should be noted here that the non-autonomous system has an obvious hysteresis behavior, which is closely related to the attraction of the attractor. By calculating the eigenvectors corresponding to the unstable focus on $HT_2$, e.g., when $F = 0.5$, the imaginary part of the eigenvectors is $0.050$, indicating that the system trajectory is relatively weakly dispersed, and thus the system trajectory moves to the right due to the increase of $F$ almost before it is dispersed.

4.   Spiking state hysteresis phenomenon and its generation mechanism

4.1.   Bifurcation hysteresis behavior in quasi-stationary process

In Figure 8(a), we find that the spiking state does not start to oscillate at the bifurcation point, but bifurcation hysteresis occurs when the perturbation amplitude $F$ increases. This bifurcation hysteresis is shown as follows: when $l\cos(\omega t)$ increases, the trajectory passes through the Hopf bifurcation and does not enter the spiking state at once, but remains in the quiescent state for a certain period of time, and then enters the spiking state slowly. In the following, the stability of the quasi-stationary process is analyzed by asymptotic theory, and the mechanism of the hysteresis phenomenon is explained.

Theorem 3. Bifurcation hysteresis occurs in the quiescent state process and is related to the stability of the quasi-stationary solution.

Proof. Let $\tau = \omega t$. Substituting into system (2.4) gives

The quiescent state process of system (4.1) is called the quasi-stationary solution, and the quasi-stationary solution can be regarded as a perturbed form of the steady state solution, which is denoted by $(\bar{x}, \bar{y})$, and therefore $(\bar{x}, \bar{y})$ can be written in the form

Substituting (4.2) into (2.3) and comparing the $\omega$ zero power coefficients yields the steady state solution of the equation as $(x_{00}, y_{00})$, where

Thus, the quasi-stationary solution becomes

Since $\omega\ll1$, terms such as $\omega$, $\omega^2$ can be neglected, i.e. $(\bar{x}(\tau, \omega), \bar{y}(\tau, \omega)) = (x_{00}, y_{00})$ can be used as the steady state solution of system (4.1), but the stability of the quasi-stationary solution with respect to $t$ is uncertain and therefore requires further analysis.

Let $P(x, y) = ax(1-\frac{x}{K+F})(\alpha + \beta x+x^2)-m x y$, $Q(x, y) = \lambda m x y-dy(\alpha + \beta x+x^2)$.

We then introduce a deviation of the form

Substituting this into system (4.1), translating the positive equilibrium to the origin, and linearizing, we obtain

The collation gives

Let

From WKB exponential approximation theory ^[33,36], it is known that there exists a solution to Eq (4.7) of the form

Substituting (4.9) into (4.7) gives

Let

The stability of the quasi-stationary process depends on $\varphi = \varphi(\tau)$, and $Re(\varphi)$ can be expressed as

If Eq (4.11) is less than 0, the quasi-stationary solution is stable; conversely, the quasi-stationary solution is unstable. When $F = F_H$, $Re[\lambda (F)] = 0$. Where $F_i$ is the initial value, $F_j$ is the value of $l\cos\tau$ when the spiking state starts to oscillate. However, for the Hopf bifurcation point, to the left side of the Hopf bifurcation point, the real part of the characteristic root is negative and the quiescent state process is stable. To the right of the Hopf bifurcation, the real part of the characteristic root is positive and the quiescent process is unstable. When $F = l\cos\tau$ slowly crosses the Hopf bifurcation critical point, the sign of Eq (4.11) does not change, so bifurcation hysteresis occurs.

According to approximation theory ^[32,37], the hysteresis time is approximately obtained by solving the least positive solution of Eq (4.11):

Substituting $F = l\cos\tau$ into Eq (4.12), we obtain the integral with respect to $\tau$

where $\tau_0 = 2k\pi \pm \arccos(F_i/l)$, $\tau_1 = 2k\pi \pm \arccos(F_j/l)$. Substituting $Re[\lambda(F)]$ into Eq (4.13) yields

Figure 9(a) gives the time history diagram of $\tau$ versus $x$. Figure 9(b) gives the diagram of $\tau_1$ versus Eq (4.15), noting that the whole of the left side of the equality sign of Eq (4.15) is $R$, and the initial value $\tau_0 = 9.8$, $\tau_1 = 10.56$ can be obtained according to Figure 9(b), and the bifurcation hysteresis time is $\tau_h = \tau_1-\tau_{H_{1.6}} = 10.56-10.197 = 0.363$. Figure 9(c) plots in detail the relationship between $\tau$ and $d\varphi/d\tau$. It is found that the area of region $S_1$ is approximately equal to the area of region $S_2$, i.e.

The bifurcation hysteresis occurs when the slow-variable perturbation crosses the Hopf bifurcation point, leading to the system's spiking state hysteresis behavior. Theoretical analysis shows that the hysteresis behavior arises because the stability of the quasi-stationary solution does not change immediately near the Hopf bifurcation point. However, the real part of the characteristic root changes sign. The exact hysteresis time was calculated using numerical simulation.

The proof of Theorem 3 is complete.

The bifurcation hysteresis effect corresponds to the hysteresis of environmental effects derived from environmental pollution or improvement. The ecological significance of this is that, as the environment improves, the predator and prey population does not change immediately, but rather the prey population increases rapidly when the environmental capacity increases to a certain value. The larger the value of the environmental capacity, the larger the peak in predator-prey population levels.

4.2.   Effect of disturbance amplitude on bifurcation hysteresis

Due to the existence of slow variables, the real bifurcation of the slow-fast coupled non-autonomous system (2.3) often has a hysteresis effect with the theoretical bifurcation point of the autonomous system (2.4). The transformed phase diagrams of the nonautonomous system (2.3) when $l = 1.5$, $l = 1.55$, and $l = 1.62$ are given in Figure 10 respectively, where the region to the right of the grey dotted line is the spiking state. It is shown that the hysteresis time of this hysteresis effect increases with the increase of the perturbation amplitude.

Figure 11 gives the diagram of $\tau_1$ versus Eq. (4.15) for different perturbation amplitudes, noting that the whole of the left side of the equality sign of equation (4.15) is $R$, with a fixed initial value $\tau_0 = 9.8$. According to Figure 11, we get $\tau_1$ equal to $10.4$, $10.48$, and $10.59$, respectively. Then, the bifurcation hysteresis time is $\tau_{h_{1.5}} = \tau_1-\tau_{H_{1.5}} = 10.4-10.125 = 0.275$, $\tau_{h_{1.55}} = \tau_1-\tau_{H_{1.55}} = 10.48-10.163 = 0.317$, and $\tau_{h_{1.62}} = \tau_1-\tau_{H_{1.62}} = 10.59-10.209 = 0.381$, respectively. Corresponding to the spiking state hysteresis behavior of Figure 10, it is found that the larger the perturbation amplitude, the longer the hysteresis time.

5.   Conclusions

The predator-prey system is a complex non-linear system influenced by many factors. Environmental capacity is an important element in the predator-prey system, which can affect the survival and extinction of both predator and prey populations. Therefore, if the perturbation of environmental capacity is added to the predator-prey system, the whole system will be a slow-fast coupled non-autonomous system. The stability of the related autonomous system is analyzed using the slow-variable parameter as the bifurcation parameter, and the condition for Hopf bifurcation is discussed. Single-Hopf bursting oscillations are first found in the slow-fast coupled non-autonomous Holling-IV predator-prey system. It is found that this non-autonomous system oscillation mode is closely related to the stability and bifurcation of the autonomous system. The Hopf bifurcation, stable equilibrium points, and limit cycle attractor of autonomous system result in quasi-periodic oscillation and single-Hopf bursting oscillations generation in non-autonomous systems. Additionally, it is found that there is an obvious spiking state hysteresis behavior of the system. The asymptotic theory was for the first time used to analyze the bifurcation hysteresis mechanism of the predator-prey system. Moreover, the hysteresis behavior arises because the stability of the quasi-stationary solution does not change immediately although the real part of the eigenvalue changes its sign in the vicinity of the Hopf bifurcation point, and the hysteresis length is calculated. It is found that the larger the perturbation amplitude, the longer the hysteresis time.

Use of AI tools declaration

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

Acknowledgments

This paper was supported from National Natural Science Foundation of China (12172233, U1934201, 11672191).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.