Next Article in Journal
An Algorithm for Calculating the Parameter Selection Area of a Doubly-Fed Induction Generator Based on the Guardian Map Method
Next Article in Special Issue
Influence of the Effective Reproduction Number on the SIR Model with a Dynamic Transmission Rate
Previous Article in Journal
Generalized n-Polynomial p-Convexity and Related Inequalities
Previous Article in Special Issue
Ultimate Dynamics of the Two-Phenotype Cancer Model: Attracting Sets and Global Cancer Eradication Conditions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Impulsive Effects and Complexity Dynamics in the Anti-Predator Model with IPM Strategies

Department of Mathematics, Yunnan Minzu University, Kunming 650500, China
*
Author to whom correspondence should be addressed.
Mathematics 2024, 12(7), 1043; https://doi.org/10.3390/math12071043
Submission received: 26 February 2024 / Revised: 25 March 2024 / Accepted: 28 March 2024 / Published: 30 March 2024
(This article belongs to the Special Issue Advances in Mathematical Biology and Applications)

Abstract

:
When confronted with the imminent threat of predation, the prey instinctively employ strategies to avoid being consumed. These anti-predator tactics involve individuals acting collectively to intimidate predators and reduce potential harm during an attack. In the present work, we propose a state-dependent feedback control predator-prey model that incorporates a nonmonotonic functional response, taking into account the anti-predator behavior observed in pest-natural enemy ecosystems within the agricultural context. The qualitative analysis of this model is presented utilizing the principles of impulsive semi-dynamical systems. Firstly, the stability conditions of the equilibria are derived by employing pertinent properties of planar systems. The precise domain of the impulsive set and phase set is determined by considering the phase portrait of the system. Secondly, a Poincaré map is constructed by utilizing the sequence of impulsive points within the phase set. The stability of the order-1 periodic solution at the boundary is subsequently analyzed by an analog of the Poincaré criterion. Additionally, this article presents various threshold conditions that determine both the existence and stability of an order-1 periodic solution. Furthermore, it investigates the existence of order-k ( k 2 ) periodic solutions. Finally, the article explores the complex dynamics of the model, encompassing multiple bifurcation phenomena and chaos, through computational simulations.

1. Introduction

Pest control is a serious challenge facing agricultural production today, and it is also a key factor affecting agricultural production efficiency, increasing income, ensuring the quality of agricultural products, and the sustainable development of agriculture. A more effective control strategy is to use the biological control method of releasing natural enemies. Consequently, studying predator-prey systems [1] is essential to mitigate or manage the types and quantities of pests that pose a threat to human development, thereby reducing economic losses. However, in reality, pest outbreaks result in countless crop losses each year. These outbreaks not only cause significant harm to agricultural production but also facilitate the spread of diseases. Consequently, pest control has become a pressing concern for agricultural and disease control authorities.
Integrated pest management (IPM) is currently recognized as one of the most effective strategies for pest control [2,3,4,5]. It capitalizes on the predatory behavior of natural enemies to regulate pest populations. In specific agricultural ecosystems, the natural enemy populations often fall short in comparison to pest numbers. Consequently, pest populations on farmlands tend to exhibit a natural growth trend due to the abundance of habitat and food resources. Considering the economic significance of farmland, it becomes necessary to artificially introduce natural enemies. Recognizing that complete eradication of pests is typically unattainable. Instead, efforts are focused on limiting pest populations below the economic threshold, the point at which they cause significant economic harm [6,7,8,9].
Artificial release of natural enemies is one of the most effective pest control methods, but it also introduces complexities to the ecosystem. As pests are preyed upon by natural enemies, they exhibit anti-predator behavior in response to these threats [10,11,12,13]. However, natural enemies are not always capable of completely suppressing pest populations. The outcome depends on the density ratio of the two and the effectiveness of anti-predator behavior.
One prevalent form of anti-predator behavior involves role reversals between predator and prey, where juvenile predators may be attacked by adult prey, thus reducing the survival rate of the predators. In this process, prey also incur losses as attacking juvenile predators increases their own risk of predation. For instance, when left alone, Crossaster papposus, which feeds on sea urchins, is attacked and consumed by multiple sea urchins. The specific process of anti-predator behavior is as follows: a leading sea urchin attempts to approach Crossaster papposus and attaches itself to one of its tentacles. However, this process carries risks, as the leading sea urchin may become food for Crossaster papposus if it fails. Once the first sea urchin starts chewing on the tentacle, other sea urchins also attach themselves to the remaining tentacles, gradually consuming Crossaster papposus bit by bit. When prey invests excessively in anti-predator behavior, the costs for both predator and prey become significant.
During the long-term evolutionary process, the emergence of anti-predator behavior can exacerbate the complexity of predation relationships, which poses challenges to the study of pest-natural enemy ecosystems in farmland [14]. To capture the complexity of these systems, a two-dimensional differential equation approach, known as the ODE model, is employed to depict the trend of changes in the number of two populations. This mathematical modeling technique for pest management decisions [15,16] builds upon the classic predator-prey model, which utilizes differential equations to represent the average growth rates of natural enemies and pests.
Furthermore, it is found that the relative growth rate described by linear function is not adaptive; considerable attention and research have been devoted to predator-prey models that incorporate various functional response functions [17,18,19,20,21,22]. Scholars have recognized the significance of considering the influence of anti-predator behavior. In this regard, Ives and Dobson [13] proposed the following predator-prey model, aiming to account for the impact of anti-predator behavior
d x ( t ) d t = λ x ( t ) 1 x ( t ) K v e ε v q x ( t ) y ( t ) 1 + a x ( t ) , d y ( t ) d t = c e ε v q x ( t ) y ( t ) 1 + a x ( t ) m y ( t ) .
Here, the variables x ( t ) and y ( t ) represent the density of the prey and predator populations at time t, respectively. The parameter λ represents the intrinsic growth rate of the pest, while m represents the death rate of the predator. The carrying capacity of the pest population is denoted by K. Additionally, the nonnegative constant v signifies the extent of the pest’s investment in anti-predation behavior, with higher values of v resulting in lower predation rates. The efficiency of the anti-predator behavior is represented by ε ; c represents the conversion rate of prey to predators; q x e ε v q x ( 1 + a x ) ( 1 + a x ) corresponds to the Holling II functional response.
In pursuit of optimal economic development and ecological preservation, impulsive control techniques are employed to reduce pest populations below a specific economic threshold. State-dependent feedback control strategies have recently emerged as a valuable approach in integrated pest management. In many practical scenarios, the decision to apply impulsive control is often contingent upon the state of the pest population, indicating that control tactics are implemented only when the model reaches a certain state. When pest populations reach the economic threshold, measures are taken to impulsively control their numbers to a certain extent. Subsequently, pest populations begin to grow again until they reach the economic threshold once more. This cyclical process of applying impulsive control is repeated to consistently keep pest populations below the economic threshold. Consequently, both economic threshold and state feedback control are considered to offer a rational framework for describing the dynamics of predation systems.
To gain a better understanding of the effects of anti-predator behavior on the dynamics of two populations [23,24,25], we propose an impulsive semi-dynamical system [26,27,28]
d x ( t ) d t = λ x ( t ) 1 x ( t ) K v e ε v q x ( t ) y ( t ) 1 + a x ( t ) , d y ( t ) d t = c e ε v q x ( t ) y ( t ) 1 + a x ( t ) m y ( t ) , x < E T , x ( t + ) = ( 1 p ) x ( t ) , y ( t + ) = y ( t ) + τ , x = E T .
Here, E T represents the economic threshold. The parameter p denotes the killing rate of pests when the pesticide is applied, with 0 < p < 1 . Additionally, τ represents the number of released predators. In the event of immediate application of control tactics, x ( t ) and y ( t ) instantly becomes ( 1 p ) x ( t ) and y ( t ) + τ , respectively. In addition, to maintain the inhibitory effect of natural enemies on pests, it is necessary to retain the number of natural enemies as much as possible. Therefore, insecticides with no or minimal impact on natural enemies were selected.
Based on model (2), the IPM strategy or intermittent control strategy is as follows: When pest density increases to the threshold E T (i.e., x = E T ), chemical and biological control measures are implemented to target both populations. On the other hand, when the pest density is updated to a value less than the prescribed threshold E T (i.e., x < E T ), no control policies are implemented, and model (2) reverts to the free model (1). The purpose of this article is to analyze the factors affecting pest outbreaks by studying the population dynamics of model (2) and to use reasonable pulse control to suppress the number of pests below the economic threshold ET. Its innovation is reflected in the impact of anti-predator behavior on both pests and natural enemies, as well as the combination of anti-predator behavior, functional response, and pulse control.
The remaining sections of the paper are organized as follows: In Section 2, the main properties of an ODE model (1) will be introduced. Section 3 focuses on the exact domain of the impulsive set and phase set of model (2). Additionally, we provide the definition of the Poincaré map. Section 4 and Section 5 are dedicated to the investigation of the boundary and internal order-1 periodic solutions, as well as the order-k periodic solution ( k 2 ) of the proposed model. Furthermore, we explore the complex dynamics of the model, including multiple bifurcation phenomena and chaos, through computational simulations. Section 6 conclude the paper by discussing our findings and their implications for pest control strategies.

2. Basic Properties of an ODE Model (1)

To investigate the dynamics of the state-dependent feedback control model (2), it is essential to examine the basic properties of model (1). The following are the two isoclines of model (1)
L 1 : y = ( 1 + a x ) ( λ x 2 + K λ x K v ) K b q x , L 2 : c b q x 1 + a x m = 0 ,
here we define b = e ε v . When y = 0 , the natural enemies become extinct, and x satisfies the equation λ x 2 + λ K x K v = 0 . Considering the condition λ K 4 v > 0 , it indicates that the prey investment v satisfies v < λ K 4 ; we obtain two boundary equilibriums: E 1 ( x 1 , 0 ) and E 2 ( x 2 , 0 ) . The values of x 1 and x 2 are
x 1 = λ K λ 2 K 2 4 λ K v 2 λ , x 2 = λ K + λ 2 K 2 4 λ K v 2 λ ,
it is evident that 0 < x 1 < x 2 < K , and x 1 + x 2 = K .
When y > 0 , model (1) admits an internal equilibrium E 3 ( x 3 , y 3 ) with
x 3 = m c b q a m , y 3 = ( 1 + a x 3 ) ( λ x 3 2 + K λ x 3 K v ) K b q x 3 ,
in order to ensure that the internal equilibrium E 3 ( x 3 , y 3 ) is positive, we need the inequality λ x 3 2 + K λ x 3 K v > 0 to hold. The necessary and sufficient condition for this inequality to have a solution is x 1 < x 3 < x 2 and v < λ K 4 . Thus, when x 3 ( x 1 , x 2 ) or when the investment v > λ K 4 , the system (1) does not possess an internal equilibrium. Additionally, x 3 > 0 is equivalent to the condition c b q > a m . To simplify calculations, we define c b q a m = h . Therefore, the internal equilibrium E 3 ( x 3 , y 3 ) can also be expressed as follows:
x 3 = m h , y 3 = c ( λ m 2 + K λ m h K v h 2 ) K m h 2 .
In particular, when v = 0 , model (1) exhibits a zero equilibrium ( 0 , 0 ) and a boundary equilibrium ( K , 0 ) . This demonstrates that without human interference, as time approaches infinity ( t + ) , the predator population gradually tends towards extinction. Unless stated otherwise, in order to guarantee the existence of a positive internal equilibrium, we assume that v < λ K 4 .
To perform a qualitative analysis of the dynamics of model (2) under varying values of v, we can define it as a hybrid dynamic system that combines discrete events and continuous system interactions. Specifically, when v = 0 , the model corresponds to a classical Lotka-Volterra predator-prey model [29] with a Holling II functional response. This variant exhibits rich dynamics depending on the relative positions of the equilibria and the choice of economic threshold. These dynamics may include multiple globally stable equilibria, limit cycles, and transcritical bifurcations, among others.
On the other hand, when v > 0 , the predation rate decreases due to the prey’s investment in anti-predator behavior. This investment indirectly affects the prey’s per capita growth rate, either positively or negatively. Based on these considerations, this paper aims to theoretically and numerically examine the dynamics of model (2).
If v = 0 , we have reformulated model (1) as
d x ( t ) d t = λ x ( t ) 1 x ( t ) K q x ( t ) y ( t ) 1 + a x ( t ) , d y ( t ) d t = c q x ( t ) y ( t ) 1 + a x ( t ) m y ( t ) ,
The properties of model (3) are presented in [30]. In this paper, we primarily focus on analyzing the dynamics of model (1) when v > 0 . By examining the relationship between stability and eigenvalues at equilibriums in a planar system, we obtain the following theorem.
Theorem 1.
(i) If x 3 > x 2 , model (1) admits no positive internal equilibrium. In this case, E 1 ( x 1 , 0 ) is unstable, and E 2 ( x 2 , 0 ) is stable; (ii) If x 1 < x 3 < x 2 and β < 0 and Δ < 0 , then E 3 ( x 3 , y 3 ) is a stable focus; (iii) If x 1 < x 3 < x 2 and β < 0 , Δ > 0 , then E 3 ( x 3 , y 3 ) is a stable node. Here we denote
β = λ 2 λ m K h λ m 2 + K λ m h K v h 2 c b q K m ,
γ = λ m 2 + K λ m h K v h 2 c b q K ,
Δ = β 2 4 γ .
Proof. 
Consider equations d x ( t ) d t and d y ( t ) d t in the model (1) as P and Q, and take the derivative of x and y to obtain the Jacobian matrix, which is defined as
J = P x P y Q x Q y = λ 2 λ K x b q y ( 1 + a x ) 2 b q x 1 + a x c b q y ( 1 + a x ) 2 c b q x 1 + a x m .
Substituting E 1 into the Jacobian matrix obtains
J | E 1 = λ 2 λ K x 1 b q x 1 1 + a x 1 0 c b q x 1 1 + a x 1 m ,
it is easy to get that λ 2 λ K x 1 = λ 2 K 2 4 λ K v K > 0 . If c b q x 1 1 + a x 1 m > 0 , then t r J | E 1 > 0 and det J | E 1 > 0 , which indicates that E 1 is an unstable focus or node. If c b q x 1 1 + a x 1 m < 0 , then det J | E 1 < 0 , which indicates that E 1 is a saddle. Then substituting E 2 into the Jacobian matrix gives
J | E 2 = λ 2 λ K x 2 b q x 2 1 + a x 2 0 c b q x 2 1 + a x 2 m ,
It is evident that λ 2 λ K x 2 = λ 2 K 2 4 λ K v K < 0 . If c b q x 2 1 + a x 2 m > 0 , then det J | E 2 < 0 , indicating that E 2 is a saddle. On the other hand, if c b q x 2 1 + a x 2 m < 0 , then t r J | E 2 < 0 and det J | E 2 > 0 , implying that E 2 is stable, as shown in Figure 1A. Analogously, the internal equilibrium E 3 is substituted into the Jacobian matrix as follows
J | E 3 = λ 2 λ m K h λ m 2 + K λ m h K v h 2 c b q K m m c c ( λ m 2 + K λ m h K v h 2 ) c b q K m 0 ,
Here, we have γ > 0 . Therefore, when x 1 < x 3 < x 2 and β < 0 and Δ < 0 , E 3 ( x 3 , y 3 ) is a stable focus, as illustrated in Figure 2. On the other hand, when x 1 < x 3 < x 2 and β < 0 and Δ > 0 , E 3 ( x 3 , y 3 ) is a stable node, as shown in Figure 3. □
In addition, for the case of x 3 < x 1 , as shown in Figure 1B, due to the small value of x 1 , the trajectory of the system is not completely in the first quadrant, which is clearly not in line with the actual biological significance, so it will not be discussed.

3. Establishment of Poincaré Map

Define R + 2 = { ( x , y ) | x 0 , y 0 } . Considering its biological significance, the subsequent research will focus on the values of x and y within this set. Moreover, given the condition E T < K , we define an open set Λ = { ( x , y ) | 0 < x < E T , y > 0 } R + 2 . Now, based on the third and fourth equations of model (2), we can proceed with the following definitions
L 3 : x = ( 1 p ) E T , L 4 : x = E T .
By combining L 1 and L 3 , we can find their intersection point, denoted as Q 1 = ( ( 1 p ) E T , y Q 1 ) , here
y Q 1 = 1 + a ( 1 p ) E T λ ( 1 p ) E T 2 + K λ ( 1 p ) E T K v K b q ( 1 p ) E T .
Analogously, we can obtain the intersection point of L 1 and L 4 , denoted as Q 2 = ( E T , y Q 2 ) , here
y Q 2 = ( 1 + a E T ) ( λ E T 2 + K λ E T K v ) K b q E T .
The impulsive set
M = { ( x , y ) R + 2 | x = E T , 0 y y Q 2 } ,
M is a closed subset of R + 2 ; we can define I : ( x , y ) M ( x + , y + ) = ( ( 1 p ) E T , y + τ ) based on model (2). Consequently, the phase set can be defined as follows
N = I ( M ) = { ( x + , y + ) Λ | x + = ( 1 p ) E T , y + Y 0 } ,
where Y 0 = [ τ , y Q 2 + τ ] ; an impulsive point of ( x , y ) is denoted as ( x + , y + ) . We assume that the initial value ( x 0 + , y 0 + ) belongs to N in this paper.
Next, we discuss the precise domain of the pulse and phase sets; this can be achieved by excluding the points in the impulsive set M that do not experience impulsive effects. Since the impulsive set M represents the largest interval in the ordinate, we can obtain the exact domain of the impulsive set by removing those points.
In Case (i) of Theorem 1, model (2) possesses a unique stable boundary equilibrium, denoted as E 2 ( x 2 , 0 ) , as illustrated in Figure 4A. If E T < x 2 < K , a trajectory Γ 2 within the system touches the line x = ( 1 p ) E T , and the point of tangency is denoted as Q 1 ( ( 1 p ) E T , y Q 1 ) , which represents the intersection of L 1 and L 3 . It is evident that trajectory Γ 2 will also intersect the line x = E T at a certain point, denoted as T ( E T , y T ) . It is worth noting that 0 < y T < y Q 2 is satisfied. Based on these considerations, the impulsive set is
M 1 = ( x , y ) R + 2 x = E T , 0 y y T ,
and the phase set is
N 1 = ( x + , y + ) R + 2 x + = ( 1 p ) E T , y + Y 1 ,
where Y 1 = [ τ , y T + τ ] . After such a definition, any solutions of model (2) that originate from the interior of trajectory Γ 2 will either not reach the impulsive set or intersect with it, leading to a mapping to the phase set N after a single impulsive effect. If x 2 < E T , the solution starting from the phase set may directly approach the stable equilibrium E 2 ( x 2 , 0 ) without reaching the impulsive set. In this case, the exact domain of the impulsive set and phase set cannot be determined.
In Case (ii) of Theorem 1, model (2) possesses an unique stable focus, we consider two cases
C 1 : x 3 < E T ; C 2 : x 3 > E T .
Case C 1 , as shown in Figure 4B, the trajectory Γ 1 is tangent to L 4 , as observed in the phase portrait. It can be seen that the point of tangency Q 2 ( E T , y Q 2 ) represents the intersection of L 1 and L 4 . The intersection of trajectory Γ 1 and L 1 on the left side is denoted as E 4 ( x 4 , y 4 ) . If ( 1 p ) E T < x 4 , employing the same analysis method as in case (i), we can define M 1 and N 1 . On the other hand, if ( 1 p ) E T > x 4 , the trajectory Γ 1 must intersect the line x = ( 1 p ) E T at two points, marked as P 1 = ( ( 1 p ) E T , y P 1 ) and P 2 = ( ( 1 p ) E T , y P 2 ) respectively. When y P 2 < y + < y P 1 , any solutions of the system originating from ( x + , y + ) N will not reach the impulsive set. Therefore, in this scenario, the impulsive set can be defined as M and
N 2 = ( x + , y + ) R + 2 x + = ( 1 p ) E T , y + Y 2 ,
where Y 2 = { [ 0 , y P 2 ] [ y P 1 , + ) } Y 0 . For Case C 2 , similar to Case (i), we can derive the impulsive set M 1 and the phase set N 1 .
In Case (iii) of Theorem 1, it is also discussed in two cases, C 1 and C 2 . In case C 1 , if ( 1 p ) E T < x 3 < E T , it is impossible to determine whether the solution of the system will reach the impulsive set. Therefore, the impulsive set and phase set cannot be defined in this case. In case C 2 , if ( 1 p ) E T < E T < x 3 , then we derive the impulsive set M 1 and the phase set N 1 .
Based on the above analysis, we have determined the exact domains of the impulsive set and phase set for model (2) in all cases. The results are summarized in Table 1. In the table, the cases of the same pulse set and phase set are defined as A i ( i = 1 , 2 , 3 ) . Note that since the pulse set of case A 3 cannot be determined, this case will not be discussed in the following sections.
Consider an arbitrary solution Ψ ( t , t 0 , ( 1 p ) E T , y 0 + ) with an initial value ( x 0 + , y 0 + ) N . This solution undergoes a single impulsive effect at time t = t 1 , a second impulsive effect at time t = t 2 , and so on, until the kth impulsive effect at time t = t k . We denote the points associated with these impulsive effects as p i = ( E T , y i ) M and p i + = ( ( 1 p ) E T , y i + ) N . Thus, we can express p i + as a continuous function: p i + = I ( p i ) . Consequently, we have y i + = y i + τ . Since p i + and p i + 1 lie on the same trajectory, p i + 1 can be uniquely determined by p i + , and y i + 1 can be uniquely determined by y i + . Specifically, we can express y i + 1 as y i + 1 = f ( y i + ) . According to the Cauchy-Lipschitz theorem, the function f is continuously differentiable in cases A 1 and A 2 .
In order to determine the domain of the Poincaré map, we need to establish the condition under which a solution initiating from p 0 + N j ( j = 1 , 2 ) does not experience impulsive effects. In case A 1 , it is straightforward to determine the domain and range of the Poincaré map as [ 0 , + ) and Y 1 = [ τ , y T + τ ) , respectively. In Case A 2 , we observe that when y P 2 < y 0 + < y P 1 , the initial point lies within the trajectory Γ 1 . In this scenario, the trajectory originating from the initial point p 0 + will not reach the impulsive set. Consequently, the domain is [ 0 , y P 2 ) ( y P 1 , + ) , and the range is [ τ , y Q 2 + τ ) . To ensure the continuous occurrence of impulsive effects, it suffices to select the initial point p i + from N 2 .
Therefore, except for the condition y P 2 < y 0 + < y P 1 in case A 2 , any solutions of model (2) initiating from p i + N j will reach p i + 1 M after a single impulsive effect. We can therefore obtain the domain and range of the Poincaré map as follows:
y i + 1 + = f ( y i + ) + τ = . φ ( y i + ) , i = 1 , 2 , , k .
The lemma regarding the Poincaré map within the phase set domain for cases A 1 and A 2 is presented.
Lemma 1.
The Poincaré map generated by the sequence of impulsive points in model (2) can be defined as follows:
y i + 1 + = φ ( y i + ) , y i + Y 1 o r Y 2 ,
where φ ( y i + ) is defined by Formula (4).
For Cases A 1 and A 2 , since the function f is continuously differentiable with y i + , φ ( y i + ) is also continuously differentiable with y i + . Therefore, the fixed point of φ ( y i + ) corresponds to an order-1 periodic solution within model (2).

4. Existence and Stability of Boundary Order-1 Periodic Solution

When τ = 0 , φ ( y i + ) has only one fixed point, which is y = 0 . This fixed point corresponds to the unique boundary order-1 periodic solution of model (2) initiated from ( ( 1 p ) E T , 0 ) . In fact, substituting y ( t ) = 0 and τ = 0 into model (2), we obtain the subsystem
d x d t = λ x ( t ) 1 x ( t ) K v , x < E T , x ( t + ) = ( 1 p ) x ( t ) , x = E T .
It follows from model (5) with x ( 0 + ) = ( 1 p ) E T that we can find the solution with period T. In order for the value of x to increase from ( 1 p ) E T to E T , we require d x d t > 0 . Thus, we have λ x 2 + λ k x k v > 0 , which implies v < λ K 4 . Additionally, the value of x must satisfy x 1 < x < x 2 . To ensure the existence of the period T, we need E T to satisfy x 1 < ( 1 p ) E T < E T < x 2 .
After performing the calculation, we have determined that the solution to model (5) is
x = x 1 + x 2 C e h 1 t 1 + C e h 1 t ,
where h 1 = λ K ( x 2 x 1 ) and C = x 1 ( 1 p ) E T ( 1 p ) E T x 2 .
Referring to Equation (6), let’s define E T = x 1 + x 2 C e h 1 T 1 + C e h 1 T . Now, if we consider T as a variable, we obtain
T = 1 h 1 I n ( E T x 1 ) x 2 ( 1 p ) E T ( 1 p ) E T x 1 ( x 2 E T ) .
Therefore, model (2) admits the boundary order-1 periodic solution with a period of T, which is given by
( x T ( t ) , 0 ) = x 1 + x 2 C e h 1 t 1 + C e h 1 t , 0 .
In the following, the local stability of solution ( x T ( t ) , 0 ) is analyzed using the analogue of the Poincaré criterion. This criterion indicates that if the Floquet multiplier u 2 satisfies condition u 2 < 1 , then ( x T ( t ) , y T ( t ) ) is orbitally asymptotically stable, where
u 2 = Δ 1 exp 0 T P x x T ( t ) , y T ( t ) + Q y x T ( t ) , y T ( t ) d t ,
Δ 1 = P + β y ϕ x β x ϕ y + ϕ x + Q + α x ϕ y α y ϕ x + ϕ y P ϕ x + Q ϕ y .
Here y T ( t ) = 0 , P + = P ( x T ( T + ) , y T ( T + ) ) , Q + = Q ( x T ( T + ) , y T ( T + ) ) and P , Q , α x , α y , β x , β y , ϕ x , ϕ y is calculated at point ( x T ( t ) , y T ( t ) ) .
Theorem 2.
The solution ( x T ( t ) , 0 ) of model (2) is asymptotically stable if and only if
r < 0 , θ r < ω + C ω + C θ K c b q λ ( 1 + a x 1 ) ( 1 + a x 2 ) ,
where r = 1 x 2 x 1 c b q x 1 K ( 1 + a x 1 ) λ m K λ , ω = 1 + a x 1 1 + a x 2 , θ = ( E T x 1 ) [ x 2 ( 1 p ) E T ] [ ( 1 p ) E T x 1 ] ( x 2 E T ) .
Proof. 
For τ = 0 , we define
P ( x , y ) = λ x 1 x K v e ε v q x y 1 + a x = λ K ( x x 1 ) ( x x 2 ) e ε v q x y 1 + a x ,
Q ( x , y ) = c e ε v q x y 1 + a x m y , α ( x , y ) = p x , β ( x , y ) = τ , ϕ ( x , y ) = x E T ,
x T ( T ) , y T ( T ) = ( E T , 0 ) , x T ( T + ) , y T ( T + ) = ( 1 p ) E T , 0 .
Since b = e ε v , we can determine the value by performing the calculation
P x = λ 2 λ K x b q y ( 1 + a x ) 2 , Q y = c b q x 1 + a x m ,
α x = p , ϕ x = 1 , α y = β x = β y = ϕ y = 0 ,
then
Δ 1 = P + β y ϕ x β x ϕ y + ϕ x + Q + α x ϕ y α y ϕ x + ϕ y P ϕ x + Q ϕ y = P x T ( T + ) , y T ( T + ) P x T ( T ) , y T ( T ) = ( 1 p ) E T x 1 ( 1 p ) E T x 2 ( E T x 1 ) ( E T x 2 ) .
Additionally, let’s discuss
0 T P x x T ( t ) , y T ( t ) + Q y x T ( t ) , y T ( t ) d t = 0 T λ m 2 λ K x T ( t ) + c b q x T ( t ) 1 + a x T ( t ) d t = ( λ m ) T 0 T 2 λ K x T ( t ) d t + 0 T c b q x T ( t ) 1 + a x T ( t ) d t ,
where
0 T 2 λ K x T ( t ) d t = 2 λ x 1 K 0 T 1 1 + C exp ( h 1 t ) d t + 2 λ x 2 K 0 T C exp ( h 1 t ) 1 + C exp ( h 1 t ) d t = 2 λ x 1 K 0 T 1 C exp ( h 1 t ) 1 + C exp ( h 1 t ) d t + 2 λ x 2 K 0 T C exp ( h 1 t ) 1 + C exp ( h 1 t ) d t = 2 λ x 1 K T 1 h 1 I n 1 + C exp ( h 1 t ) | 0 T + 2 λ x 2 K 1 h 1 I n 1 + C exp ( h 1 t ) | 0 T = 2 λ x 1 K T + 2 λ K h 1 ( x 2 x 1 ) I n 1 + C exp ( h 1 T ) 1 + C
and
0 T c b q x T ( t ) 1 + a x T ( t ) d t = I 1 + I 2 ,
let I 1 and I 2 be denoted as
I 1 = c b q x 1 0 T 1 ( 1 + a x 1 ) + ( 1 + a x 2 ) C exp ( h 1 t ) d t = c b q x 1 1 + a x 2 0 T 1 ( 1 + a x 1 ) 1 + a x 1 ) ( 1 + a x 2 ) ( 1 + a x 2 ) + C exp ( h 1 t ) d t = c b q x 1 1 + a x 2 1 + a x 2 1 + a x 1 T 1 h 1 I n ( 1 + a x 1 ) 1 + a x 1 ) ( 1 + a x 2 ) ( 1 + a x 2 ) + C exp ( h 1 t ) | 0 T = c b q x 1 1 + a x 1 T 1 h 1 I n ( 1 + a x 1 ) 1 + a x 1 ) ( 1 + a x 2 ) ( 1 + a x 2 ) + C exp ( h 1 T ) ( 1 + a x 1 ) ( 1 + a x 1 ) ( 1 + a x 2 ( 1 + a x 2 ) + C
and
I 2 = c b q x 2 0 T C exp ( h 1 t ) ( 1 + a x 1 ) + ( 1 + a x 2 ) C exp ( h 1 t ) d t = c b q x 2 1 + a x 2 0 T C exp ( h 1 t ) 1 + a x 1 1 + a x 1 1 + a x 2 1 + a x 2 + C exp ( h 1 t ) d t = c b q x 2 1 + a x 2 1 h 1 I n 1 + a x 1 1 + a x 1 1 + a x 2 1 + a x 2 + C exp ( h 1 t ) | 0 T = c b q x 2 1 + a x 2 h 1 I n 1 + a x 1 1 + a x 1 1 + a x 2 1 + a x 2 + C exp ( h 1 T ) 1 + a x 1 1 + a x 1 1 + a x 2 + C 1 + a x 2 + C ,
respectively. Let θ = ( E T x 1 ) [ x 2 ( 1 p ) E T ] [ ( 1 p ) E T x 1 ] ( x 2 E T ) , ω = 1 + a x 1 1 + a x 2 ; this yields
exp 0 T P x x T ( t ) , y T ( t ) + Q y x T ( t ) , y T ( t ) d t = θ 1 x 2 x 1 x 2 x 1 m K λ + c b q x 1 K ( 1 + a x 1 ) λ 1 + C 1 + C θ 2 ω + C θ ω + C K c b q λ ( 1 + a x 1 ) ( 1 + a x 2 ) ,
Thus, the Floquet multiplier
u 2 = Δ 1 exp 0 T P x x T ( t ) , y T ( t ) + Q y x T ( t ) , y T ( t ) d t ,
after calculation
Δ 1 θ 1 + C 1 + C θ 2 = 1 ,
thus
u 2 = θ 1 x 2 x 1 c b q x 1 K ( 1 + a x 1 ) λ m K λ ω + C θ ω + C K c b q λ ( 1 + a x 1 ) ( 1 + a x 2 ) .
To determine the domain of u 2 , we consider θ as a function of p. Since 0 p < 1 , when p = 0 , we have θ = 1 and u 2 = 1 . When 0 < p < 1 , based on the inequality x 1 < ( 1 p ) E T < E T < x 2 , we can deduce that θ > 1 and ω + C θ ω + C > 1 . Since K c b q λ ( 1 + a x 1 ) ( 1 + a x 2 ) > 0 , we have ω + C θ ω + C K c b q λ ( 1 + a x 1 ) ( 1 + a x 2 ) > 1 . In order to satisfy u 2 < 1 , we need r < 0 and θ r < ω + C ω + C θ K c b q λ ( 1 + a x 1 ) ( 1 + a x 2 ) . Therefore, for all p ( 0 , 1 ) , we have 0 < u 2 < 1 . Under the conditions of Theorem 2, any solution trajectory of model (2) will eventually approach solution ( x T ( t ) , 0 ) infinitely, indicating that solution ( x T ( t ) , 0 ) is asymptotically stable. □
For the global attractivity of solution ( x T ( t ) , 0 ) , when x 1 < x < x 3 , the sequence of impulsive point y k + of any solutions originating from the ( x + , y + ) N monotonically decreases with increasing k, so we have d y d t < 0 . Additionally, Figure 5 shows that lim k y k + = y * = 0 , indicating the global attractivity of the solution ( x T ( t ) , 0 ) . On the side, when x 3 < x < x 2 , the sequence of impulsive points y k + monotonically increases with increasing k, leading to the solution ( x T ( t ) , 0 ) becoming unstable, as depicted in Figure 6.

5. Existence and Stability of Order-k Positive Periodic Solution

The previous section has systematically described the boundary solution of the model (2), and our primary focus of this section is on examining the properties of the internal solution when τ > 0 , especially the existence and stability of the order-1 periodic solution as well as the existence of the order-k periodic solution. This aim can be achieved by demonstrating the presence and stability of the fixed point of the Poincaré map φ ( y i + ) .
Firstly, there is a threshold condition for the existence of order-1 periodic solution in case A 1 .
Theorem 3.
In case A 1 , φ ( y i + ) possesses a fixed point. As a result, the model (2) exhibits an order-1 periodic solution.
Proof. 
For this particular case, the trajectory Γ 2 starting from Q 1 ( 1 p ) E T , y Q 1 crosses x = E T , and the point where they intersect is labeled as T E T , y T 1 . Where Q 1 is the tangent point between curve Γ 2 and line x = ( 1 p ) E T . Since both point Q 1 and point T lie on the same curve, it can be concluded that y T = f ( y Q 1 ) .
When a single impulsive effect occurs, T is mapped to Q 1 + ( ( 1 p ) E T , y Q 1 + ) , and obtain y Q 1 + = y T + τ from the pulse function. Accordingly, since y Q 1 + is determined by τ , we assume the existence of τ 0 such that y T + τ 0 = y Q 1 . Then, when τ = τ 0 , the point Q 1 and the point Q 1 + coincide, establishing a fixed point for φ ( y i + ) . Consequently, the trajectory Q 1 T ^ is treated as an order-1 periodic solution for model (2).
If τ > τ 0 , as illustrated in Figure 7A, the position of the impulsive point Q 1 + is higher than that of point Q 1 . Consequently, we have φ ( y Q 1 ) y Q 1 = y Q 1 + y Q 1 > 0 . This leads to the following inequality:
φ ( y Q 1 ) > y Q 1 .
The trajectory originating from Q 1 + crosses x = E T , and the point where they intersect is labeled as T 1 ( E T , y T 1 ) . Since any two trajectories of the plane system cannot intersect, it follows that the position of point T 1 is lower than that of point T. After undergoing a pulse effect, T 1 is mapped to T 1 + ( ( 1 p ) E T , y T 1 + ) . As y T 1 < y T , we have y T 1 + = y T 1 + τ < y T + τ = y Q 1 + . Consequently, the position of point T 1 + is lower than that of point Q 1 + . Hence, φ ( y Q 1 + ) y Q 1 + = y T 1 + y Q 1 < 0 . This leads to the following inequality:
φ ( y Q 1 + ) < y Q 1 + .
By combining (8) and (9), we conclude that there possesses a fixed point within the interval ( y Q 1 , y Q 1 + ) for φ ( y i + ) .
If τ < τ 0 , the position of the impulsive point Q 1 + is lower than that of point Q 1 . Consequently, we have φ ( y Q 1 ) y Q 1 = y Q 1 + y Q 1 < 0 . This leads to the following inequality:
φ ( y Q 1 ) < y Q 1 .
Select a point H + ( 1 p ) E T , y H + on the line x = ( 1 p ) E T that satisfies the condition 0 < y H + < τ . The trajectory originating from H + crosses x = E T , and the point where they intersect is labeled as H 1 ( E T , y H 1 ) , where y H 1 = f ( y H + ) . After undergoing a pulse effect, H 1 is mapped to H 1 + ( 1 p ) E T , f ( y H + ) + τ . Now, considering φ ( y H + ) y H + = f ( y H + ) + τ y H + > 0 , we obtain the following inequality:
φ ( y H + ) > y H + .
Combining (10) and (11), we conclude that there possesses a fixed point within the interval ( y H + , y Q 1 ) for φ ( y i + ) .
The above process illustrates that no matter what value τ takes, the system (2) always has an order-1 periodic solution. □
Next, explain the conditions for the existence of an order-1 periodic solution in case A 2 . Before that, define a key point Q 2 + ( x Q 2 + , y Q 2 + ) , which is the impulsive point that appears on the trajectory and initiates from point P 1 ( ( 1 p ) E T , y P 1 ) after a single impulsive effect.
Theorem 4.
In case A 2 , if y Q 2 + y P 1 or y Q 2 + y P 2 , it can be observed that the Poincaré map φ ( y i + ) possesses a fixed point. As a consequence, model (2) exhibits an order-1 periodic solution.
Proof. 
For this situation, let’s mark E 5 ( x 5 , y 5 ) as the second intersection of Γ 1 and L 2 . If ( 1 p ) E T x 5 , as depicted in Figure 7B, since y ( t ) is increasing on the right area of L 2 , the position of Q 2 is higher than that of P 2 . In addition, Q 2 is mapped to Q 2 + ( ( 1 p ) E T , y Q 2 + ) after undergoing a pulse effect. Considering τ > 0 , it follows that the position of point Q 2 + is also higher than that of point P 2 .
If y Q 2 + = y P 1 , then the trajectory P 1 Q 2 ^ forms an order-1 periodic solution of the system (2). If y Q 2 + > y P 1 , we have φ ( y P 1 ) y P 1 = y Q 2 + y P 1 > 0 . This leads to the following inequality:
φ ( y P 1 ) > y P 1 .
Furthermore, the trajectory originating from Q 2 + crosses x = E T , and the point where they intersect is labeled as Q 3 E T , y Q 3 . By virtue of the uniqueness of the solution, we deduce that the position of Q 3 is lower than that of Q 2 . After undergoing a pulse effect, Q 3 is mapped to Q 3 + ( 1 p E T , y Q 3 + ) . Clearly, the position of Q 3 + is also lower than that of Q 2 + , hence φ ( y Q 2 + ) y Q 2 + = y Q 3 + y Q 2 + < 0 . This yields
φ ( y Q 2 + ) < y Q 2 + .
Combining (12) and (13), we conclude that there possesses a fixed point within the interval ( y P 1 , y Q 2 + ) for φ ( y i + ) .
If x 4 < ( 1 p ) E T < x 5 , two cases arise based on the relative positions of points P 2 and Q 2 . In the first case, when y Q 2 > y P 2 , the point Q 2 is positioned above P 2 . As demonstrated in the previous proof, φ ( y i + ) possesses a fixed point in this case.
In the second case, when y Q 2 y P 2 , two sub-cases can be considered. If y Q 2 + y P 1 , the conclusion is evidently valid. On the other hand, if y Q 2 + y P 2 , the proof process from case τ τ 0 in Theorem 3, indicates φ ( y i + ) possesses a fixed point for this case.
In summary, the order-1 periodic solution of model (2) always exists. □
Theorems 3 and 4 have confirmed the existence of the order-1 periodic solution in the system (2), and then we will prove the stability of this periodic solution.
Define ξ ( t ) , η ( t ) to represent an order-1 periodic solution with a period of T. The points where this periodic solution intersects with the lines L 3 and line L 4 are denoted as follows:
N + ξ ( T + ) , η ( T + ) = ( 1 p ) E T , y N + τ , N ξ ( T ) , η ( T ) = E T , y N ,
then it can be calculated.
Δ 1 = P + 1 p E T , y N + τ P E T , y N = λ 1 p E T 1 1 p E T K v e ε v q 1 p E T y N + τ 1 + a ( 1 p ) E T λ E T 1 E T K v e ε v q E T y N 1 + a E T ,
and
0 T P x + Q y d t = 0 T λ m 2 λ ξ t K b q η t 1 + a ξ t 2 + c b q ξ t 1 + a ξ t = 0 T H t d t .
Hence, the Poincaré criterion serves as a valuable tool to analyze the stability of the order-1 periodic solution ξ ( t ) , η ( t ) in model (2).
Theorem 5.
The solution ξ ( t ) , η ( t ) of model (2) is asymptotically stable when
Δ 1 exp 0 T H t d t < 1 .
Proof. 
For case A 1 , considering the proof methodology of Theorem 2, u 2 can be expressed as
u 2 = Δ 1 exp 0 T P x + Q y d t = Δ 1 exp 0 T H t d t .
Looking back at Theorem 3, it is evident when τ = τ 0 , y T + τ 0 = y Q 1 holds, which implies that the points Q 1 + and Q 1 coincide. Consequently, we obtain u 2 = 0 < 1 , indicating that ξ ( t ) , η ( t ) is stable in this case.
When τ > τ 0 , point Q 1 + lies above Q 1 , and it is evident that point N + is also positioned above Q 1 . Moreover, y Q 1 < y N + τ . Thus, for τ > τ 0 , we have λ 1 p E T 1 1 p E T K v e ε v q 1 p E T y N + τ 1 + a ( 1 p ) E T < 0 and exp 0 T H t d t > 0 . Thus, it holds true that u 2 < 0 . Similarly, when τ < τ 0 , we have u 2 > 0 .
In summary, when Formula (14) is satisfied, i.e., u 2 < 1 . Regardless of the initial value, the solution of system (2) will approach solution ξ ( t ) , η ( t ) , so the solution ξ ( t ) , η ( t ) is asymptotically stable. Similarly, under the conditions of Theorem 4, ξ ( t ) , η ( t ) is also asymptotically stable for case A 2 . □
Next, we will investigate the existence of order k ( k 2 ) periodic solutions for model (2).
Theorem 6.
In Case A 1 , if y Q 1 + y Q 1 , then model (2) only possesses an order-1 periodic solution. However, if y Q 1 + > y Q 1 and y T 1 + y Q 1 , then the system (2) can have either an order-1 periodic solution or an order-2 periodic solution.
Proof. 
If y Q 1 + y Q 1 , according to Theorem 3, a fixed point must exist in φ ( y i + ) . The trajectory originating from Q 1 + crosses x = E T , and the point where they intersect is labeled as T 1 . Since any two trajectories cannot intersect, the position of T 1 is lower than that of T, and T 1 will be mapped to T 1 + ( ( 1 P ) E T , y T 1 + ) after undergoing a pulse effect. Since y T 1 < y T , the impulsive point T 1 + of point T 1 satisfies y T 1 + = y T 1 + τ < y T + τ = y Q 1 + , and thus, the position of point T 1 + is lower than that of point Q 1 + . After undergoing another pulse effect, the next impulsive point, T 2 + , is obtained, and the position of T 2 + is lower than that of T 1 + . By analogy, the position of point T i + 1 + is lower than that of point T i + ( i = 1 , 2 , , n ). Therefore, the sequence of impulsive points monotonically decreases and gradually approaches the point N + , yielding
y N + τ < y T n + < < y T i + < < y T 1 + < y Q 1 + .
Furthermore, since each trajectory does not intersect with others, the system (2) only possesses an order-1 periodic solution and does not have an order-2 periodic solution.
If y Q 1 + > y Q 1 and y T 1 + = y Q 1 , it is evident that the conclusion is valid. If y T 1 + > y Q 1 , The trajectory originating from T 1 + crosses x = E T , and the point where they intersect is labeled as T 2 , satisfying the inequality y T 1 < y T 2 < y T . The point T 2 is mapped to the point T 2 + after undergoing a pulse effect and satisfies the inequality y T 1 + < y T 2 + < y Q 1 + . Repeat the above steps as the control strategy is implemented multiple times. It can be observed that the position of T i is between T i 1 and T i 2 , and the position of T i + is between T i 1 + and T i 2 + ( i = 3 , 4 , , k ). If the series of impulsive points is divided into odd and even series, two possibilities arise.
( a ) y Q 1 < y T 1 + < y T 3 + < < y T 2 n 1 + < y T 2 n + 1 + < < y Q 1 + , ( b ) y Q 1 < < y T 2 n + < y T 2 n 2 + < < y T 4 + < y T 2 + < y Q 1 + .
Based on a and b , the odd series y T 2 n + 1 + and even series y T 2 n + are increasing and decreasing, respectively, on the interval ( y Q 1 , y Q 1 + ) . Both series are also convergent. Furthermore, since all impulsive segments do not intersect each other, either has a fixed point or exhibits a stable period two-point cycle for φ ( y i + ) . This corresponds to the order-1 or order-2 periodic solutions of the model (2), respectively. □
According to the aforementioned analysis, we can also conclude that for case A 2 , if y Q 2 + < y P 2 , then model (2) only possesses an order-1 periodic solution. If y Q 2 + > y P 1 and y Q 3 + y P 1 , then model (2) has either an order-1 periodic solution or an order-2 periodic solution.
Now, we will proceed to discuss the conditions under which the order-3 periodic solution exists.
Theorem 7.
Let y m = min y + N 1 φ y + = y Q 1 . In case A 1 , if y Q 1 + > y Q 1 and y T 1 + < y m , then the order-3 periodic solution of model (2) exists.
Proof. 
According to Theorem 3, when y Q 1 + > y Q 1 , the φ ( y i + ) possesses a fixed point on ( y Q 1 , y Q 1 + ) . To prove the model (2) admits an order-3 periodic solution, it should be demonstrated that there is a y on interval N 1 makes φ y y and φ 3 y = y . According to the requirements stipulated, yield
φ 3 0 = φ 2 τ > 0 ,
and
φ 3 y m = φ 2 y Q 1 = φ ( y Q 1 + ) = y T 1 + < y m .
Considering the property of φ ( y i + ) being continuous, we can say that there is a number y ¯ in the interval 0 , y m that makes φ 3 y ¯ = y ¯ . Hence, model (2) has an order-3 periodic solution. □
Theorem 7 indicates that a discrete mapping φ ( y i + ) has a periodic point with a minimum period of 3. According to the Period Three theorem and the Sarkovskii theorem, it is shown that mapping φ ( y i + ) also has periodic points with a minimum period of any positive integer. The reason is as follows: From the reasoning process of the Sarkovsskii theorem, it can be inferred that there is such a sequence rearranged by all positive integers:
3 , 5 , 7 , 9 , 11 , 3 × 2 , 5 × 2 , 7 × 2 , 9 × 2 , 11 × 2 , 3 × 2 2 , 5 × 2 2 , 7 × 2 2 , 9 × 2 2 , 11 × 2 2 , 3 × 2 3 , 5 × 2 3 , 7 × 2 3 , 9 × 2 3 , 11 × 2 3 , 2 5 , 2 4 , 2 3 , 2 2 , 2 1 , 2 0 .
This sequence is called the Sarkovsskii sequence; periodic points will operate in the order of the sequence. The Sarkovsskii theorem states that when m is sorted before n in the Sarkovsskii sequence, if there are periodic points with period m, then there must be periodic points with period n, and 3 is the first digit of the Sarkovsskii sequence, before all positive integers. According to the Sarkovsskii theorem, if there is a period of 3 in the mapping φ ( y i + ) during the iteration process, there must be a period point with a period of any positive integer.
This means the existence of order-k periodic solutions ( k 3 ) in model (2). This observation suggests the presence of chaos in model (2). Compared with reference [13], our study not only emphasizes the phase diagrams and sets of each equilibrium point but also highlights the existence and stability of periodic solutions under pulse control.
For cases A 1 and A 2 , we have established the existence and stability of the order-1 periodic solution, as well as the conditions for the existence of the order-k periodic solution ( k 2 ), in model (2). However, to enhance our understanding of model dynamics, numerical simulations are performed to validate the accuracy of our theoretical findings.
As shown in Figure 8, we observe a bifurcation behavior of model (2) by selecting parameter p as bifurcation parameter. It is evident that as p increases, the system undergoes abrupt transitions from an order- ( k + 1 ) periodic solution to an order-k periodic solution ( k 2 ) through period-decreasing bifurcations accompanied by chaotic regions. Furthermore, when p > 0.61 , period-halving bifurcations occur, leading to the emergence of an order-1 periodic solution. This demonstrates the significant influence of the key parameter p on the dynamics of model (2). Moreover, it indicates that the growth patterns and quantity trends of pests are dependent on the initial quantity ( x 0 + , y 0 + ) .
The bifurcation diagram of the parameter τ reveals the intricate dynamics of model (2), as illustrated in Figure 9. When 0 < τ < 0.6582 , it does not show the solution of model (2). However, at τ = 0.6582 , a transcritical bifurcation occurs, giving rise to an order-1 periodic solution. In addition, this solution will become unstable in pace as the parameter τ increases further, leading to a period-doubling bifurcation and the subsequent generation of an order-2 periodic solution in model (2). As there is a continued increase in τ , the period-doubling bifurcations drive the system (2) into a chaotic regime. Notably, the system (2) demonstrates distinct transitions from an order-k periodic solution to an order- ( k + 1 ) periodic solution (where k 2 ) through period-adding bifurcations involving chaotic regions. These solutions validate that both populations can coexist in different states, thus further corroborating the conclusions of the above theorems.
Figure 9 clearly illustrates the transition from regular dynamics to chaos, where a period-doubling bifurcation occurs as the bifurcation parameter values increase. The originally stable periodic orbit constantly becomes unstable, and new stable orbits appear, but the emergence of new stable orbits does not mean that the original orbit disappears. It can only be said that the original orbit exists in an unstable form in the new stable state. As the number of iterations increases, a blurry area eventually appears, indicating that the system has entered a completely chaotic state. And Figure 8 can be seen as the inverse process of the period-doubling bifurcation.
Next, the correctness of the above theorems is verified, and schematic diagrams of the order-1, order-2, order-3, and even order-4 periodic solutions of the system (2) are presented in conjunction with the bifurcation diagram. Here, we mainly discuss case A 2 . As shown in Figure 10, Figure 11, Figure 12 and Figure 13. In addition, in order to further understand the population dynamics of prey and predators, we have also provided the phase portrait and time series. This also once again demonstrates the complex, dynamic behavior of the system. Here, we choose the order-2 and order-4 periodic solutions as examples to illustrate.
For the given parameter values, it can be observed from Figure 11 that the order-2 periodic solution exhibits distinct characteristics. It consists of an active period with short trajectories and a quiescent period with long trajectories. During the long stretches of trajectory where both populations reside, the pest population undergoes a period of rapid growth until an outbreak occurs, as the number of predators is relatively low. However, with the introduction of predators, a quiescent period comes to an end, and the active period begins following a single impulsive effect. The population density of the pest immediately adjusts to ( 1 p ) E T and experiences rapid growth until reaching E T on short notice. At this point, control tactics are implemented, causing the quantities of both populations to return to their long stretches of trajectory, thus initiating a new cycle.
Similarly, by varying the values of the key parameters p and τ , we can obtain the order-4 periodic solution of model (2), as illustrated in Figure 13. The phase portrait and time series exhibit a more intricate pattern, featuring two quiescent periods and two active periods with long trajectories and short trajectories, respectively. These active and quiescent periods alternate, giving rise to highly complex dynamics. It is evident that random perturbations play a vital role in the occurrence of pest outbreaks. Small changes in the key parameters can directly impact the population dynamics and outbreak patterns of the pests. This further emphasizes the sensitivity of the system to variations in parameter values, highlighting the need for effective control measures to mitigate and manage pest outbreaks.
The above descriptions of bifurcation and chaos all emphasize that the complex dynamics of the model pose challenges to the comprehensive control strategy of pests. This also indicates that analyzing biological problems and explaining biological phenomena and related conclusions from a mathematical perspective is still extremely difficult. In addition to using predatory models to describe and analyze the common biological activity of pest control, other biological activities such as virus outbreaks and vaccination can also be explained using SIR mathematical models. For example, in an improved SIR model in [31], the author and his team linked the classical SIR model with pulse vaccination, analyzed the dynamic properties of disease-free periodic solution and stable positive T-periodic solution, and described the bifurcation phenomenon and chaos in the pulse model. In addition, there are also models in [32,33] that combine mathematics with biology, which reveal that biological systems can be simplified into computable models, thereby helping us understand and predict biological processes.

6. Conclusions

To enhance the efficacy of pest control, impulsive tactics such as releasing natural enemies and using insecticides are employed when the pest density increases to E T . In light of this, these control strategies dependent on pest status prove more suitable for practical implementation. However, in previous models, the detrimental impacts of pests’ anti-predator behavior on themselves were often disregarded. Recognizing the significance of this aspect, the Lotka-Volterra model incorporates the effect of anti-predator behavior on both parties, and its primary objective is to offer more rational pest control strategies and illustrate the dynamics of pest-natural enemy ecosystems more accurately.
The phase portrait of different equilibria was examined as the parameters m and v underwent changes. Subsequently, the relationship between the equilibrium E 3 and the line x = E T was utilized to discuss the precise domain of the phase set and impulsive set in cases ( i ) to ( i i i ) . So as to facilitate the analysis of the periodic solution in subsequent studies, these cases where the domain is identical are defined as A i ( i = 1 , 2 , 3 ) , as illustrated in Table 1. Building upon the threshold control process outlined in the impulsive semi-dynamical system, it was observed that the ordinates of two adjacent impulsive points satisfy y i + 1 + = f ( y i + ) + τ , leading to the establishment of the corresponding Poincaré map.
The stability condition about the boundary order-1 periodic solution ( x T ( t ) , 0 ) is determined by Theorem 2. For the given parameters, an illustration depicting the global stability of solution ( x T ( t ) , 0 ) when τ = 0 is provided. In this scenario, the number of pests will exhibit periodic oscillations, while predators will gradually approach 0. This behavior is depicted in Figure 5. As the parameters a and p are altered, ( x T ( t ) , 0 ) is no longer stable, leading to the coexistence of pests and natural enemies in a stable focus; see Figure 6 for details.
Where τ 0 , we provide the threshold conditions sequentially for the existence of order-k ( k = 1 , 2 , 3 ) periodic solutions; this suggests the presence of chaos within the system. Notably, these conditions also indicate that pests and natural enemies can coexist under specific parameter values. Numerical results demonstrate that the system exhibits diverse and intricate dynamic properties, including transcritical bifurcations, period-decreasing bifurcations, period-adding bifurcations, chaotic regions, and more.
One of the highlights of this article is the bifurcation and chaos of the model (2). When related to other biological models [14], it can be found that we all express the stability and threshold conditions of the periodic solution and also elaborate on the switching process of the periodic solution under different bifurcations. However, the innovation of this article is to discuss chaos under different sensitivity parameters and provide periodic solutions to various orders of the system in combination with chaos.
Significant development has been achieved in the theoretical exploration of predator-prey systems [34,35,36,37] over recent years. However, this paper introduces several innovative contributions that set it apart from previous research. The key distinctions are emphasized in the following aspects: (1) Consideration of the reciprocal influence of anti-predator behavior on both populations. By acknowledging the strategic actions taken by pests to minimize predation rates, the dynamics of both populations are affected, establishing a symbiotic relationship between the species; (2) Precise delineation of the domains of the phase set and impulsive set based on distinct threshold conditions. This enables a more accurate understanding of when impulsive control strategies should be employed and how the system behaves accordingly; (3) Provision of multiple bifurcation diagrams for essential parameters. These diagrams provide a comprehensive depiction of the intricate and diverse dynamics exhibited by system (2), facilitating a deeper comprehension of its behavior across varying parameter values; (4) Illustration of order-k periodic solutions, specifically focusing on order-2 and order-4 solutions. These solutions shed light on the periodic behavior exhibited by the system, contributing to a comprehensive understanding of its temporal patterns.
Collectively, these distinct characteristics presented in this article contribute to the progress in our comprehension of predator-prey systems and provide valuable insights into the intricate dynamics of model (2). However, there are still some areas for improvement in the article, such as assuming that harvesting has a meaning. In such a situation, if the harvesting can be either a constant number, a proportional number, or a periodic number of the density of the prey over a given time interval, what can change in the system’s behavior? Additionally, this study did not account for the impact of limited resources, and the utilization of linear pulse control fails to optimize the economic conditions of farmland. Exploring the impact of resource constraints and investigating alternative control strategies to enhance the economic viability of agricultural settings will be key research objectives in the future.

Author Contributions

Idea, W.Q.; Writing—original draft, Z.D.; Writing—review and editing, L.H.; Funding acquisition, W.Q. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by the National Natural Science Foundation of China (Nos. 12261104, 12361104), the Youth Talent Program of Xingdian Talent Support Plan (No. XDYC-QNRC-2022-0708), and the Yunnan Provincial Basic Research Program Project (Nos. 202401AT070036, 202301AT070016).

Data Availability Statement

No data were used for the research described in the article.

Acknowledgments

The authors would like to thank the editor and the anonymous reviewers for their constructive comments and suggestions to improve the quality of the paper.

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

References

  1. Doran, C.; Bierbach, D.; Lukas, J.; Klamser, P.; Landgraf, T.; Klenz, H.; Habedank, M.; Arias-Rodriguez, L.; Krause, S.; Romanczuk, P.; et al. Fish waves as emergent collective antipredator behavior. Curr. Biol. 2022, 32, 708–714. [Google Scholar] [CrossRef] [PubMed]
  2. Khan, M.; Mahmood, H.Z.; Damalas, C.A. Pesticide use and risk perceptions among farmers in the cotton belt of punjab, pakistan. Crop. Prot. 2015, 67, 184–190. [Google Scholar] [CrossRef]
  3. Tian, Y.; Gao, Y.; Sun, K. A fishery predator-prey model with anti-predator behavior and complex dynamics induced by weighted fishing strategies. Math. Biosci. Eng. 2023, 20, 1558–1579. [Google Scholar] [CrossRef] [PubMed]
  4. Marten, A.L.; Moore, C.C. An options based bioeconomic model for biological and chemical control of invasive species. Ecol. Econ. 2011, 70, 2050–2061. [Google Scholar] [CrossRef]
  5. Tan, X.; Qin, W.; Liu, X.; Yang, J.; Jiang, S. Sliding bifurcation analysis and global dynamics for a filippov predator-prey system. J. Nonlinear Sci. Appl. 2016, 9, 3948–3961. [Google Scholar] [CrossRef]
  6. Tan, X.; Tang, S.; Chen, X.; Xiong, L.; Liu, X. A stochastic differential equation model for pest management. Adv. Differ. Equ. 2017, 2017, 197. [Google Scholar] [CrossRef]
  7. Swierk, L.; Petrula, M.; Esquete, P. Diving behavior in a neotropical spider (Trechalea extensa) as a potential antipredator tactic. Ethology 2022, 128, 508–512. [Google Scholar] [CrossRef]
  8. Burger, A.L.; Fennessy, J.; Fennessy, S.; Dierkes, P.W. Nightly selection of resting sites and group behavior reveal antipredator strategies in giraffe. Ecol. Evol. 2020, 10, 2917–2927. [Google Scholar] [CrossRef]
  9. Udayagiri, S.; Norton, A.P.; Welter, S.C. Integrating pesticide effects with inundative biological control: Interpretation of pesticide toxicity curves for anaphes iole in strawberries. Entomol. Exp. Appl. 2000, 95, 87–95. [Google Scholar] [CrossRef]
  10. Tang, S.; Tang, B.; Wang, A.; Xiao, Y. Holling ii predator–prey impulsive semi-dynamic model with complex poincare map. Nonlinear Dyn. 2015, 81, 1575–1596. [Google Scholar] [CrossRef]
  11. Tang, S.; Cheke, R.A. State-dependent impulsive models of integrated pest management (ipm) strategies and their dynamic consequences. J. Math. Biol. 2005, 50, 257–292. [Google Scholar] [CrossRef]
  12. Huang, M.; Li, J.; Song, X.; Guo, H. Modeling impulsive injections of insulin: Towards artificial pancreas. SIAM J. Appl. Math. 2012, 72, 1524–1548. [Google Scholar] [CrossRef]
  13. Ives, A.R.; Dobson, A.P. Antipredator behavior and the population dynamics of simple predator-prey systems. Am. Nat. 1987, 130, 431–447. [Google Scholar] [CrossRef]
  14. Tian, Y.; Tang, S.; Cheke, R.A. Dynamic complexity of a predator-prey model for ipm with nonlinear impulsive control incorporating a regulatory factor for predator releases. Math. Model. Anal. 2019, 24, 134–154. [Google Scholar]
  15. Huang, M.; Yang, A.; Yuan, S.; Zhang, T. Stochastic sensitivity analysis and feedback control of noiseinduced transitions in a predator-prey model with anti-predator behavior. Math. Biosci. Eng. 2023, 20, 4219–4242. [Google Scholar] [CrossRef]
  16. Qin, W.; Tang, S.; Cheke, R.A. The effects of resource limitation on a predator-prey model with control measures as nonlinear pulses. Math. Probl. Eng. 2014, 2014, 450935. [Google Scholar] [CrossRef]
  17. Lcg, A.; Rs, B.; Sd, C.; Cdv, C.; Dsja, B. Short-term predation risk and habitat complexity influence cheetah antipredator behaviours. Anim. Behav. 2021, 178, 175–184. [Google Scholar]
  18. Gao, Y.; Yang, F. Persistence and extinction of a modified leslie-gower holling-type ii two-predator one-prey model with levy jumps. J. Biol. Dyn. 2022, 16, 117–143. [Google Scholar] [CrossRef]
  19. Li, S.; Liu, W. A delayed holling type iii functional response predator-prey system with impulsive perturbation on the prey. Adv. Differ. Equ. 2016, 2016, 42. [Google Scholar] [CrossRef]
  20. Tang, S.; Pang, W.; Cheke, R.A.; Wu, J. Global dynamics of a state-dependent feedback control system. Adv. Differ. Equ. 2015, 2015, 322. [Google Scholar] [CrossRef]
  21. Zhang, J.; Chen, X. Dynamic behaviors of a discrete commensal symbiosis model with holling type functional response. IAENG Int. J. Appl. Math. 2023, 53, 277. [Google Scholar]
  22. Li, G.; Yang, Q. Dynamics of a stochastic holling ii predator-prey model with levy jumps and habitat complexity. Int. J. Biomath. 2021, 14, 2150077. [Google Scholar] [CrossRef]
  23. Tian, Y.; Tang, S.; Cheke, R.A. Nonlinear state-dependent feedback control of a pest-natural enemy system. Nonlinear Dyn. 2018, 94, 2243–2263. [Google Scholar] [CrossRef]
  24. Tang, G.; Qin, W.; Tang, S. Complex dynamics and switching transients in periodically forced filippov prey–predator system. Chaos Solitons Fractals 2014, 61, 13–23. [Google Scholar] [CrossRef]
  25. Qin, W.; Tan, X.; Tosato, M.; Liu, X. Threshold control strategy for a non-smooth filippov ecosystem with group defense. Appl. Math. Comput. 2019, 362, 124532. [Google Scholar] [CrossRef]
  26. Dong, T.; Zhu, H. Filippov fitzhugh-nagumo neuron model with membrane potential threshold control policy. Neural Process. Lett. 2021, 53, 3801–3824. [Google Scholar] [CrossRef]
  27. Zhang, Q.; Tang, S.; Zou, X. Rich dynamics of a predator-prey system with state-dependent impulsive controls switching between two means. J. Differ. Equ. 2023, 364, 336–377. [Google Scholar] [CrossRef]
  28. Liu, J.; Qi, Q.; Liu, B.; Gao, S. Pest control switching models with instantaneous and non-instantaneous impulsive effects. Math. Comput. Simul. 2023, 205, 926–938. [Google Scholar] [CrossRef]
  29. Yadav, S.; Kumar, V. A prey–predator model and control of a nematodes pest using control in banana: Mathematical modeling and qualitative analysis. Int. J. Biomath. 2022, 15, 2150089. [Google Scholar] [CrossRef]
  30. Yang, J.; Tang, S. Holling type ii predator–prey model with nonlinear pulse as state-dependent feedback control. J. Comput. Appl. Math. 2016, 291, 225–241. [Google Scholar] [CrossRef]
  31. de Carvalho, J.P.M.; Rodrigues, A.A. Pulse vaccination in a modified sir model: Global dynamics, bifurcations and seasonality. arXiv 2023, arXiv:2312.02306. [Google Scholar]
  32. de Carvalho, J.P.M.; Rodrigues, A.A. Sir model with vaccination: Bifurcation analysis. Qual. Theory Dyn. Syst. 2023, 22, 105. [Google Scholar] [CrossRef]
  33. Kanchanarat, S.; Nudee, K.; Chinviriyasit, S.; Chinviriyasit, W. Mathematical analysis of pulse vaccination in controlling the dynamics of measles transmission. Infect. Dis. Model. 2023, 8, 964–979. [Google Scholar] [CrossRef]
  34. Tang, S.; Liang, J.; Tan, Y.; Cheke, R.A. Threshold conditions for integrated pest management models with pesticides that have residual effects. J. Math. Biol. 2013, 66, 1–35. [Google Scholar] [CrossRef]
  35. Tang, S.; Pang, W. On the continuity of the function describing the times of meeting impulsive set and its application. Math. Biosci. Eng. 2017, 14, 1399–1406. [Google Scholar] [CrossRef] [PubMed]
  36. Tang, B.; Xiao, Y. Bifurcation analysis of a predator–prey model with anti-predator behaviour. Chaos Solitons Fractals 2015, 70, 58–68. [Google Scholar] [CrossRef]
  37. Tan, Y.; Cai, Y.; Yao, R.; Hu, M.; Wang, W. Complex dynamics in an eco-epidemiological model with the cost of anti-predator behaviors. Nonlinear Dyn. 2022, 107, 3127–3141. [Google Scholar] [CrossRef]
Figure 1. The stable boundary equilibrium of model (1). Parameters are λ = 2 , K = 50 , q = 0.2 , a = 0.01 ,   c = 0.5 , ε = 0.8 . (A): v = 0.65 , m = 2.12 ; (B): v = 5 , m = 0.1 .
Figure 1. The stable boundary equilibrium of model (1). Parameters are λ = 2 , K = 50 , q = 0.2 , a = 0.01 ,   c = 0.5 , ε = 0.8 . (A): v = 0.65 , m = 2.12 ; (B): v = 5 , m = 0.1 .
Mathematics 12 01043 g001aMathematics 12 01043 g001b
Figure 2. The stable focus of model (1) with same parameters as in Figure 1 except v = 0.24 , m = 1 .
Figure 2. The stable focus of model (1) with same parameters as in Figure 1 except v = 0.24 , m = 1 .
Mathematics 12 01043 g002
Figure 3. The stable node of model (1) with same parameters as in Figure 1 except v = 0.5 , m = 2 .
Figure 3. The stable node of model (1) with same parameters as in Figure 1 except v = 0.5 , m = 2 .
Mathematics 12 01043 g003
Figure 4. Impulsive set and phase set. (A) Case (i) of Theorem 1, i.e., there is no interior equilibrium of system (1). All parameters are identical to those in Figure 1; (B) Case (ii) of Theorem 1, i.e., there is a stable focus of system (1). All parameters are identical to those in Figure 2.
Figure 4. Impulsive set and phase set. (A) Case (i) of Theorem 1, i.e., there is no interior equilibrium of system (1). All parameters are identical to those in Figure 1; (B) Case (ii) of Theorem 1, i.e., there is a stable focus of system (1). All parameters are identical to those in Figure 2.
Mathematics 12 01043 g004
Figure 5. Stability of boundary order-1 periodic solution ( x T ( t ) , 0 ) of model (2) with parameters λ = 1, K = 100, q = 0.2, c = 0.5, v = 0.434, ε = 0.8, m = 1, ET = 20, a = 0.34, p = 0.6.
Figure 5. Stability of boundary order-1 periodic solution ( x T ( t ) , 0 ) of model (2) with parameters λ = 1, K = 100, q = 0.2, c = 0.5, v = 0.434, ε = 0.8, m = 1, ET = 20, a = 0.34, p = 0.6.
Mathematics 12 01043 g005
Figure 6. Instability of boundary order-1 periodic solution ( x T ( t ) , 0 ) with a = 0.01 , p = 0.2 . All other parameters are identical to those in Figure 5.
Figure 6. Instability of boundary order-1 periodic solution ( x T ( t ) , 0 ) with a = 0.01 , p = 0.2 . All other parameters are identical to those in Figure 5.
Mathematics 12 01043 g006
Figure 7. (A) Location of the periodic solution of system (2) for case A 1 ; (B) Location of the periodic solution of system (2) for case A 2 . All parameter values are fixed as follows: λ = 2 , K = 50 , q = 0.2 , a = 0.04 , c = 0.5 , v = 0.12 , ε = 0.8 , m = 1 .
Figure 7. (A) Location of the periodic solution of system (2) for case A 1 ; (B) Location of the periodic solution of system (2) for case A 2 . All parameter values are fixed as follows: λ = 2 , K = 50 , q = 0.2 , a = 0.04 , c = 0.5 , v = 0.12 , ε = 0.8 , m = 1 .
Mathematics 12 01043 g007
Figure 8. Bifurcation diagram of model (2) with p. Here the parameters are: λ = 1 , K = 50 , a = 0.19 , c = 0.45 , m = 0.36 , v = 0.03 , ε = 3.513 , q = 0.211 , E T = 25 , τ = 4.1 .
Figure 8. Bifurcation diagram of model (2) with p. Here the parameters are: λ = 1 , K = 50 , a = 0.19 , c = 0.45 , m = 0.36 , v = 0.03 , ε = 3.513 , q = 0.211 , E T = 25 , τ = 4.1 .
Mathematics 12 01043 g008
Figure 9. Bifurcation diagram of model (2) with τ . Here the parameters are: λ = 1 , K = 50 , a = 0.19 , c = 0.45 , m = 0.36 , v = 0.03 , ε = 3.513 , q = 0.211 , E T = 25 , p = 0.054 .
Figure 9. Bifurcation diagram of model (2) with τ . Here the parameters are: λ = 1 , K = 50 , a = 0.19 , c = 0.45 , m = 0.36 , v = 0.03 , ε = 3.513 , q = 0.211 , E T = 25 , p = 0.054 .
Mathematics 12 01043 g009
Figure 10. An order-1 periodic solution of model (2) with p = 0.7 ,   τ = 4.7 ,   E T = 17 ; the other parameters are identical to those in Figure 9.
Figure 10. An order-1 periodic solution of model (2) with p = 0.7 ,   τ = 4.7 ,   E T = 17 ; the other parameters are identical to those in Figure 9.
Mathematics 12 01043 g010
Figure 11. An order-2 periodic solution of model (2) with p = 0.2 ,   τ = 8.5 ,   E T = 12 ; the other parameters are identical to those in Figure 9.
Figure 11. An order-2 periodic solution of model (2) with p = 0.2 ,   τ = 8.5 ,   E T = 12 ; the other parameters are identical to those in Figure 9.
Mathematics 12 01043 g011
Figure 12. An order-3 periodic solution of model (2) with p = 0.3 ,   τ = 5.8 ,   E T = 19 ; the other parameters are identical to those in Figure 9.
Figure 12. An order-3 periodic solution of model (2) with p = 0.3 ,   τ = 5.8 ,   E T = 19 ; the other parameters are identical to those in Figure 9.
Mathematics 12 01043 g012
Figure 13. An order-4 periodic solution of model (2) with p = 0.13 ,   τ = 4 ,   E T = 12 ; the other parameters are identical to those in Figure 9.
Figure 13. An order-4 periodic solution of model (2) with p = 0.13 ,   τ = 4 ,   E T = 12 ; the other parameters are identical to those in Figure 9.
Mathematics 12 01043 g013
Table 1. Impulsive and phase sets of model (2).
Table 1. Impulsive and phase sets of model (2).
CaseCase E T ( 1 p ) E T M s N s
A 1 (i) E T < x 2 < K ( 1 p ) E T < E T M 1 N 1
A 3 (i) x 2 < E T ( 1 p ) E T < x 2 ××
A 1 (ii) C 1 ( 1 p ) E T < x 4 M 1 N 1
A 2 (ii) C 1 ( 1 p ) E T > x 4 M N 2
A 1 (ii) C 2 ( 1 p ) E T < E T M 1 N 1
A 3 (iii) C 1 ( 1 p ) E T < x 3 ××
A 1 (iii) C 2 ( 1 p ) E T < E T M 1 N 1
Here, × denotes there exists no pulse set or phase.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Qin, W.; Dong, Z.; Huang, L. Impulsive Effects and Complexity Dynamics in the Anti-Predator Model with IPM Strategies. Mathematics 2024, 12, 1043. https://doi.org/10.3390/math12071043

AMA Style

Qin W, Dong Z, Huang L. Impulsive Effects and Complexity Dynamics in the Anti-Predator Model with IPM Strategies. Mathematics. 2024; 12(7):1043. https://doi.org/10.3390/math12071043

Chicago/Turabian Style

Qin, Wenjie, Zhengjun Dong, and Lidong Huang. 2024. "Impulsive Effects and Complexity Dynamics in the Anti-Predator Model with IPM Strategies" Mathematics 12, no. 7: 1043. https://doi.org/10.3390/math12071043

APA Style

Qin, W., Dong, Z., & Huang, L. (2024). Impulsive Effects and Complexity Dynamics in the Anti-Predator Model with IPM Strategies. Mathematics, 12(7), 1043. https://doi.org/10.3390/math12071043

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop