Dynamic analysis of a predator-prey state-dependent impulsive model with fear e ﬀ ect in which action threshold depending on the prey density and its changing rate

: In ecology, the impact of predators goes beyond killing prey, the mere presence of predators reduces the ability of prey to reproduce. In this study, we extend the predator-prey model with fear e ﬀ ect by introducing the state-dependent control with a nonlinear action threshold depending on the combination of the density of prey and its changing rate. We initially deﬁned the Poincar´e map of the proposed model and studied its fundamental properties. Utilizing the properties of the Poincar´e map, periodic solution of the model is further investigated, including the existence and stability of the order-1 periodic solution and the existence of the order-k ( k ≥ 2) periodic solutions. In addition, the inﬂuence of the fear e ﬀ ect on the system’s dynamics is explored through numerical simulations. The action threshold used in this paper is more consistent with the actual growth of the population than in earlier linear threshold studies, and the results show that the control objectives are better achieved using the action threshold strategy. The analytical approach used in this study provided several novel methods for analyzing the complex dynamics that rely on state-dependent impulsive.


Introduction
Ecosystem modeling and analysis is a fascinating and active research field in biology, in which the study of the interaction between predator and prey is a central topic in ecology and evolutionary biology [1,2].With the progress of ecological research, the views on how predators interact with the prey and affect the prey population have also changed greatly.Predators have a dual impact on the structure of the ecological system [3].On the one hand, predators physically attack and kill the prey [4], thereby reducing the number of prey.On the other hand, the existence of predators alone can make the prey individual in a state of psychological pressure, resulting in many changes in the behavior of prey species [5].This is a common anti-predator response, also known as the 'fear effect', which eventually slows the growth of prey populations [6].
In attempts to better understand the influence of fear effect on the predator-prey system, scholars developed lots of mathematical models of predator-prey system with fear effect.For example, Wang et al. [7] investigated the impacts of fear effect on the stability of the system, and obtained that the fear effect can lead to the bifurcation of limit cycle under appropriate parameter values.Further, Wang and Zou et al. [8] inspired by the experimental work by Zanette and his colleagues [9], expanded their previous study [7] by incorporating age structures into prey populations (juvenile and adult stages) and allowing adult prey to adaptively avoid predators.Das and Samanta [10] investigated the impacts by introducing an exponential form of fear effect into a stochastic predator-prey system when extra meals were offered to predators.Sahoo and Samanta [11] explored a two-prey and one-predator model by incorporating fear factors into switching mechanisms during prey reproduction and predation.Das et al. [12] studied a predator-prey model that incorporates fear factors into including birth and mortality in prey populations with a functional Holling type-II response.See references [13][14][15][16] for more fear factor influences in predator-prey model dynamics.
In addition to the above interactions, the predator-prey system is also disturbed by various external factors, which can cause huge changes in the number of species in the short term.This change is called an impulse in the modeling process, and the impulse differential equation can well describe this change [17].Correspondingly, there are two types impulsive control strategies in impulsive differential equations, namely, fixed-time impulse and state-dependent impulse [17][18][19][20].
State-dependent feedback control methods are commonly represented by impulsive semi-dynamical systems, which can well characterize threshold control strategies.That is, control interventions will only take place when a target species size reaches a preset threshold density.In recent years, statedependent feedback control has received extensive attention from researchers and has been applied in many fields and sciences, for example, fishery harvesting [21][22][23], integral pest management [24][25][26], the effects between biological populations [27][28][29], virus control [30][31][32], etc.
In controlling reality biological species, we focus on choosing a proper time to implement the control interventions and designing a practical strategy in proceeding the control.A basic assumption of many existing models is that when the biological population density reaches a fixed threshold, then the control is carried out, which is subject to the above-mentioned state-dependent feedback control.But from a biological point of view, a fixed threshold control policy usually ignores many important aspects in the real world.There are two practical situations: low population density but high rate of change; high population density but low rate of change.In the above two situations, we may not need to trigger the control interventions even though the population density is high, but the growth rate is low.Thus, the threshold policy needs to include both population density and its changing rate, which is actually a nonlinear threshold control strategy.Note that, the nonlinear threshold policy was recently studied in many works [17,31,33,34], and rich dynamic behaviors were observed in these studies.In the control of biological populations, another assumption in the impulsive model is that the changing rate is proportional to the population density due to the use of control interventions.However, there should definitely be a limitation for this term.Therefore, the saturated form using a nonlinear function can be more realistic.To the best of our knowledge, no study has tried to introduce the nonlinear threshold policy and the nonlinear control function into the predator-prey system with fear effect, and the analysis of the dynamic behaviors remains challenging.
The major goal of this research is to introduce the nonlinear threshold policy and the nonlinear function of the control effect into an impulsive predator-prey system with fear effect, propose a state-dependent impulsive system with nonlinear threshold control, and focus on studying its dynamic behaviors.In Section 2, a predator-prey state-dependent impulsive model with fear effect is proposed, and a ratio-dependent nonlinear action threshold is adopted.In Section 3, the corresponding Poincaré map is constructed by giving the range of impulsive set and phase set, and some fundamental properties of Poincaré map are discussed.In Section 4, the conditions for the existence and stability of order-1 periodic solution are given, and the existence of order-k (k ≥ 2) periodic solutions are discussed.The Section 5 is numerical simulation, which verifies our theoretical results and discusses the effect of fear factors on system dynamics.Finally, the relevant biological significance are discussed and conclusions are given.

Model formulation and preliminaries
In the literature [35], the authors considered the impulse control guided by a linear threshold of the density of prey based on a predator-prey model with fear effect, and partially studied the dynamics of the proposed model.In the current study, we extend the model in [35] by introducing a nonlinear threshold control policy, taking the combination of the density of prey and its changing rate as the action threshold.The model is given by: where x(t) and y(t) are the sizes of the prey and predator populations, respectively.b, d and c represent the birth rate, natural mortality rate, and density-dependent decay rate caused by intraspecific competition for prey, respectively, k represents the level of fear.px(t) 1+wx(t) indicates the Holling II functional response, e denotes the ratio of biomass conversion (satisfying the restriction 0 < e < p), m is the death rate of the predator.The parameters u 1 , v 1 and ET are all positive constants with u 1 + v 1 = 1.
Another major assumption in existing research is that the prey capture rate is a linear function dependent on population density, and the predator release amount is a positive number.Nevertheless, in real nature, natural resources are finite, and implementing of control policies should be determined by the present size of the population.Both the capture rates of prey and the release amount of predator should be based on population saturation functions.Consequently, we consider the finiteness of natural resources in our model is very meaningful.We apply the following nonlinear impulsive functions relevant to biological population density and capture rate.
that is, when the prey population density reaches the action threshold u 1 x(t) + v 1 dx(t) dt = ET , the control measures will be taken to update the number of prey and predator to x(t) − δx(t) 2 x(t)+µ and y(t) + τ 1+ηy(t) , respectively.Here, δ and µ denote maximum capture rates and the half-saturation constant for the prey, τ represents the maximum number of predators released, η is the morphological coefficient.The release of predators will not surpass the τ due to realistic factors, such as definite resources.For simplicity, we denote x(t + ) = x + and y(t + ) = y + .
Without the feedback control, the ODE system With the dynamics of system (2.2) has been investigated in [35], and here we beriefly recall its dynamics that are useful in the present research.The two nullclines of system (2.2) are noted by L 1 and L 2 , where .
System (2.2) always exists a trivial equilibrium O(0, 0), which is a saddle point, and a boundary equilibrium .
) is unstable and there exists a unique stable limit cycle, where .

The Poincaré map and its properties
In a biological sense, our discussion is limited to the region R 2 + = {(x, y) : x ≥ 0, y ≥ 0}.The impulsive set and phase set become two curves because system (2.1) adopts the action threshold control strategy of prey density and its changing rate.According to u 1 x(t) + v 1 dx(t) dt = ET and the first equation of system (2.1), we can obtain denoted by L M , where where denoted by L N .Next, unless otherwise specified, we always take an initial point S + 0 (x + 0 , y + 0 ) from the curve L N .The trajectory starting from L N may not reach L M due to different initial conditions.Now, define the range of the impulsive set and phase set based on the different positions of the trajectories of system (2.1) with the equilibrium point (see Figure 1).
Case (i) , there exists a curve Γ T 1 such that the L N tangent to this cure at the T 1 (x T 1 , y T 1 ).Clearly, the curve Γ T 1 intersects the L M at point T 2 (x T 2 , y T 2 ) (see Figure 1 (a)).Therefore, we denote the impulsive set and phase set as: where (x, y) is a point on the curve L M and (x + , y + ) is a point on the curve L N .
Case (ii) x * ≤ x * M When x * ≤ x * M , there must have a point A(x A , y A ) on impulsive set such that the trajectory Γ A is tangent to L M at that point and intersects L 2 at point E 1 (x E 1 , y E 1 ).Suppose the horizontal coordinate of curve L N at y = y E 1 is x N 1 .Therefore, we consider the following two different cases: the situation is the same as Case (i), and the trajectory Γ A does not intersect L N .There exists a trajectory Γ B 1 tangent to L N at point B 1 (x B 1 , y B 1 ), and after a period of time the trajectory Γ B 1 will intersect the impulsive set L M at point B 2 (x B 2 , y B 2 ) (see Figure 1 (b)).For this situation, the impulsive set is defined by and the corresponding range of phase set is defined by , where y C 1 > y C 2 .In this case, it is easy to observe that none of the solution trajectory from the interior of segment C 1 C 2 reach the curve L M , that is, it is not affected by the impulsive effect (see Figure 1 (c)).So we define the impulsive set by and the phase set is Based on the above discussion, the construction of Poincaré map is given below.Given an initial point S + 0 (x + 0 , y + 0 ) ∈ N, we define the trajectory starting from S + 0 as Γ(t, t 0 , S + 0 ) = Γ(x(t, t 0 , (x + 0 , y + 0 )), y(t, t 0 , (x + 0 , y + 0 )).For any S + i (x + i , y + i ) ∈ N, the trajectory from S + i reaches the L N at point S i+1 (x i+1 , y i+1 ) in a finite time t.We express this process as where y i+1 = y + ( t, x + i , y + i ).From the Cauchy-Lipschitz Theorem we know that y i+1 can only be represented by y + i .Therefore, we have Consider the scalar differential equation from system (2.1): The function g(x, y) is continuously differentiable.Further, we denote x + 0 = X, y + 0 = Z, and here the value of x is between the L N and the L M .According to (3.2) have Thus, from Eqs. (3.1) and (3.4), the Poincaré map can be expressed as Next, we provide some fundamental properties of the Poincaré map Q(Z) in the following theorem.(ii) Q(Z) is continuously differentiable in [0, +∞), and has a unique fixed point.
When x * ≤ x * M and x E 1 > x N 1 , the properties of Poincaré map is similar to Case (i).In what follows, we consider the case where x E 1 < x N 1 .
Theorem 3.2.Suppose x * ≤ x * M and x E 1 < x N 1 .The Poincaré map satisfies (Figure 3): we know the trajectory Γ A is tangent to the impulsive set L M and intersect with the phase set L N at two points C 1 and C 2 .For ∀S , the trajectory starting from the S + i cannot reach the L M .Meanwhile, the Q(Z) is well defined on (0, The trajectory from these two points will reach the impulse set after time t, i.e., Γ( t, x + i , y + i ) = Γ( t, x i+1 , y i+1 ).Since the uniqueness of the solution we have y k 1 +1 < y k 2 +1 .Hence, according to the expression of Poincaré map The trajectory starts from point y + k 1 , y + k 2 cross the nullcline L 1 and intersects the L N at the points S k 1 (x k 1 , y k 1 ) and S k 2 (x k 2 , y k 2 ).We can easily obtain y k 1 > y k 2 , and from the uniqueness of the solution, we have (ii) From equation (3.2) we can easily determine that g(x, y) is a continuously differentiable function.Therefore, it follows from the continuity and differentiability theorem for solution of ordinary differential equation that the Poincaré map Q(Z) is also continuously differentiable on (0, y C 2 ] ∪ [y C 1 , +∞).
(iii) The curve Γ A from C 1 reaches the L M at point A, and then arrives in L N at point By the existence of zeros theorem and the continuous differentiability of Q(Z), there exists a point ỹ ∈ (y C 1 , y + A ) and satisfies Q(ỹ) = ỹ, which implies that Q(Z) has a unique fixed point on [y C 1 , +∞).

Order-k periodic solution
From the discussion in the above section, we explore the existence of an equilibrium point in the system, which implies that there is an order-1 periodic solution to system (2.1).Next, we will explore more properties of order-k (k ≥ 1) periodic solutions.(i) If Q(y T 1 ) < y T 1 , the order-1 periodic solution of system (2.1) is globally asymptotically stable; (ii) If Q(y T 1 ) > y T 1 , then system (2.1) has a globally stable order-1 periodic solution if and only if Therefore, the trajectory from S + 0 (x + 0 , y + 0 ) ∈ N reaches the point S 1 (x 1 , y 1 ) on the impulsive set after a period of time, and will arrive at point S + 1 (x + 1 , y + 1 ) ∈ N through the impulsive effect.This process can be defined as For different value ranges of initial point y + 0 , we will have the following three situations: Case a.If ỹ < y + 0 < y T 1 , since Q(y T 1 ) < y T 1 and Q(Z) is monotonously increasing on [0, y T 1 ], so we have , by mathematical induction we have That is Q k (y + 0 ) increases monotonically and lim this is the situation a.; when 0 < Q(y + 0 ) < ỹ, this is the same as situation b.
In summary, we always can obtain (ii) Sufficient condition: When Q(y T 1 ) > y T 1 , Q(Z) has unique fixed point ỹ ∈ [y T 1 , +∞) by Theorem 3.1 and monotonically decreasing on [y T 1 , +∞).Thus, if Q 2 (y + 0 ) > y + 0 for any y + 0 ∈ [y T 1 , ỹ], then have Further more, we derive that Thus by the monotone bounded theorem we can get lim = ỹ, then order-1 periodic solution is globally stable.
Necessary condition: Suppose ỹ is an order-1 periodic solution of system (2.1), i.e., By the global stability of ỹ it follows that there exists ỹ+ , where ε is sufficiently small.It also follows from the continuity of Poincaré map that there must exist a ỹ ∈ (ỹ + 2 , ỹ+ 1 ) satisfying Q 2 (ỹ ) = ỹ , which shows that system (2.1) has an order-2 periodic solution and contradicts the conditions we know.Thus, if there exists a globally stable order-1 periodic solution for system (2.1), then Q 2 (y + 0 ) > y + 0 for any then system (2.1) exists an order-1 periodic solution which is globally stable.Proof.By Theorem 3.2, when Q(y C 1 ) > y C 1 , we obtain that the Poincaré map Q(Z) has a unique fixed point on the interval [y C 1 , +∞), which implies that system (2.1) has an order-1 periodic solution on the interval [y C 1 , +∞).
For any By further deduction we can get Therefore, we have lim 1) has a stable order-1 or order-2 periodic solution.
Proof.For any initial point S + 0 (x + 0 , y + 0 ) in the phase set, where x + 0 , y + 0 > 0, there exists a positive integer k such that y + k = Q k y + 0 holds after the pulse.When y + 0 < y T 1 , it follows from Theorem 3.1 that Q(Z) has no fixed point in the interval [0, y T 1 ] and is monotonically increasing.Therefore, there exists a integer q such that y + q > y T 1 and y + q−1 < y T 1 , then can get y + q = Q(y + q−1 ) < Q y T 1 , i.e., y T 1 < y + q < Q y T 1 .

Numerical simulations
In this subsection, we will design some numerical examples to verify the theoretical results.In addition, the influence of fear factor k on the dynamics of system (2.1) is also discussed.
We set the parameters as k = 0.04, d = 0.02, c = 0.01, w = 1, e = 0.44, m = 0.2, δ = 0.2, µ = 1, τ = 1.8, η = 5, u 1 = 0.6, v 1 = 0.4, ET = 0.7.For Case (i), we set b = 0.6, p = 1, it can be seen from Theorem 4.1 that when x * ≥ x * M , system (2.1) has a stable order-1 periodic solution.Observation of Figure 4(a) shows that system (2.1) forms a stable order-1 periodic solution after the impulse.For Case (ii)(II), let b = 0.7, p = 0.5, when Q(y C 1 ) > y C 1 , system (2.1) has a stable order-1 periodic solution, as is illustrated in Figure 4(d).Let b = 2.3, p = 0.5, w = 0.1 and k = 0.5, when Q(y C 1 ) > y C 1 , Figure 4(g) indicates that system (2.1) does not have periodic solution, after multiple impulses, the system is stable at the equilibrium E * (0.998, 3.313).Figure 5(a) shows the trajectory of the order-2 periodic solution of system (2.1). Figure 5(b) and Figure 5(c) are the solutions of the prey and predator corresponding to the order-2 periodic solution, respectively.By observing Figure 5, we can find that the existence of order-2 periodic solution of system (2.1) indicates that reasonable predation behaviour can keep the system in a stable ecological equilibrium.As the parameter b increases, the existence of the order-1 periodic solution of system (2.1) cannot be guaranteed, for example, when b = 0.8, system (2.1) has an order-3 periodic solution (Figure 6), and when b = 1.4,system (2.1) has an order-4 periodic solution (Figure 7).
In the action threshold, when the weighted coefficients u 1 = 1 and v 1 = 0, the impulsive set and phase set become two straight lines (Figure 8(a)).This reduces to the linear impulsive systems that have been studied in most previous, in which case the action threshold only depends on the population density of the prey.When changing the weighted coefficients of the action threshold, here v 1 0, we can observe that the impulsive set and phase set become two curves and that the curvature of the curves varies with u 1 and v 1 , see Figure 8.When we use an action threshold control strategy determined by the density of the prey and its rate of change, system (2.1) still has a stable order-1 periodic solution (Figure 8(a)-(c)).On the other hand, if the comprehensive control strategy is more dependent on the change rate of the prey, the control goal can be successfully achieved after a limited number of control measures, making the prey population density is below a given action threshold, as can be seen in Figure 8(d).
To explore the impact of the fear coefficient k on the dynamics of system (2.1), we choose k as the bifurcation parameter and perform a one-dimensional bifurcation analysis on system (2.1).It can be seen from Figure 9 that the internal equilibrium E * of system (2.1) is unstable when the parameter set satisfies k < k * .The bifurcation diagram of k (0 < k ≤ 0.8) shows that there is a very complex dynamic phenomenon in system (2.1).In particularly, the order-3 periodic solutions exist in the relatively small scope of k (0.06 < k < 0.08), such as Figure 9(a), which further verifies the validity of the results shown in Theorem 4.4.Comparing the bifurcation diagrams at different birth rates b, we conclude that changing b can significantly alter the dynamics of system (2.1), and an interesting phenomenon is that period doubling and period halving branches occur as b increases, as shown in Figures 9(b

Conclusion
State-dependent impulsive semi-dynamic systems are among the most discontinuous non-smooth systems and have been investigated in many applications like integrated pest managememt, virus dynamical systems, and diabetes treatment [17,[36][37][38].In the last few years, substantial progress has been made in the study of such systems in terms of analytical techniques, qualitative analysis, and applied research [39,40].In particular, with the help of the mathematical tool of Poincaré map, a comprehensive and in-depth branch study of impulsive systems with nonlinear impulsive functions obtained are more realistic, as the various factors affecting the population are taken into account in the modelling.

Figure 6 .Figure 7 .
Figure 6.(a) An order-3 periodic solution of system (2.1).The time series of prey population x(t) (b) and predator population y(t) (c).The parameter b = 0.8, other parameter values are same as those shown in Figure 5.
Figure5(a) shows the trajectory of the order-2 periodic solution of system (2.1).Figure5(b) and Figure5(c) are the solutions of the prey and predator corresponding to the order-2 periodic solution, respectively.By observing Figure5, we can find that the existence of order-2 periodic solution of system (2.1) indicates that reasonable predation behaviour can keep the system in a stable ecological equilibrium.As the parameter b increases, the existence of the order-1 periodic solution of system (2.1) cannot be guaranteed, for example, when b = 0.8, system (2.1) has an order-3 periodic solution (Figure6), and when b = 1.4,system (2.1) has an order-4 periodic solution (Figure7).In the action threshold, when the weighted coefficients u 1 = 1 and v 1 = 0, the impulsive set and phase set become two straight lines (Figure8(a)).This reduces to the linear impulsive systems that have been studied in most previous, in which case the action threshold only depends on the population density of the prey.When changing the weighted coefficients of the action threshold, here v 1 0, we can observe that the impulsive set and phase set become two curves and that the curvature of the curves varies with u 1 and v 1 , see Figure8.When we use an action threshold control strategy determined by the density of the prey and its rate of change, system (2.1) still has a stable order-1 periodic solution (Figure8(a)-(c)).On the other hand, if the comprehensive control strategy is more dependent on the change rate of the prey, the control goal can be successfully achieved after a limited number of control measures, making the prey population density is below a given action threshold, as can be seen in Figure8(d).To explore the impact of the fear coefficient k on the dynamics of system (2.1), we choose k as the bifurcation parameter and perform a one-dimensional bifurcation analysis on system (2.1).It can be seen from Figure9that the internal equilibrium E * of system (2.1) is unstable when the parameter set satisfies k < k * .The bifurcation diagram of k (0 < k ≤ 0.8) shows that there is a very complex dynamic phenomenon in system (2.1).In particularly, the order-3 periodic solutions exist in the relatively small scope of k (0.06 < k < 0.08), such as Figure9(a), which further verifies the validity of the results shown in Theorem 4.4.Comparing the bifurcation diagrams at different birth rates b, we conclude that changing b can significantly alter the dynamics of system (2.1), and an interesting phenomenon is that period doubling and period halving branches occur as b increases, as shown in Figures9(b) and 9(c).The island in the middle of Figure9(b) indicates that the density of prey and predators shows a complex pattern when k lies in the interval around (0.06, 0.19).