Global dynamics for a non-smooth Filippov system with the threshold strategy

In this study, the Leslie-Gower model with functional response is extended into a non-smooth Filippov system by applying IPM strategies. Once the number of pests reaches or surpasses the given economic threshold(ET), spraying pesticides and releasing the natural enemy are implemented simultaneously. In order to maintain the pest population at or below ET, global dynamics of the proposed model are investigated completely, including the existence of sliding mode and various equilibria, sliding dynamics and global stability of equilibria. The result shows that real equilibrium cannot coexist with the unique pseudo-equilibrium. In particular, after excluding the existence of any possible limit cycle, the global stability of equilibria is obtained by employing qualitative and numerical techniques. In the end, the eﬀect of our work on pest control are dis-cussed.


Introduction
It is well-known that insects play an important role in ecological balance, but some of them collectively referred to as pests would have a substantial adverse impact on human life [1][2][3]. This is especially conspicuous in agriculture. Once the number of pests exceeds a certain value, causing infestation, it is bound to bring about an immeasurable loss to crop yield, crop quality, together with social economy. Therefore, pest management has been always a hot spot for many scholars [4][5][6]. During the World War II, some organochlorine insecticides, such as DDT, were found to be effective in controlling pests, but the long-term abuse of these chemical pesticides can lead to insect resistance, environmental pollution and even the resurgence of secondary pests due to the damage of pesticides to natural enemies [7][8][9]. According to the survey report from U.S. department of agriculture [10], although the dosage and concentration of pesticides used in agriculture in the United States increased tenfold respectively between 1940 and 1978, the loss of crop yield went up from 7% to 17% and the quantity of pests rose as well. So it was gradually realized that using chemical methods alone was not feasible, and the concept of IPM(Integrated Pest Management) came up in the meantime. IPM is aimed at integrate suitable measures, such as chemical and biological methods, to manage the number of pests below EIL(Economic Injury Level) instead of killing off them, thereby not only avoiding economic losses but also contributing to the sustainable development of agroecosystem [11][12][13]. Actually, taking measures when the number of pests approaches to EIL cannot prevent the economic toll timely, so the Economic Threshold(ET<EIL) is set as a critical value, that is, once the pest population reaches or exceeds ET, appropriate measures should be applied to reduce the quantity below ET, which can ensure that EIL is never possible for the pest population to reach.
There have been many mathematical models established to investigate the efficiency of IPM strategies, including fixed moment and state-dependent impulsive differential equations [14][15][16]. However, some weaknesses of this class of models can be noticed. For the former, control strategies are implemented at the fixed moment regardless of whether the pest population reaches ET, which can enhance the cost and cause a waste of resources. And for the latter, they assume that the death of pests is an instantaneous process, whereas this process is continuous in reality [17]. Hence, compared with above models, the non-smooth Filippov system is more realistic as it means that once the quantity of pests surpasses ET, continuous control measures will be taken until it falls below ET again. Recently, Filippov systems with the threshold policy have been extensively investigated [18][19][20][21][22][23]. But referring to Filippov predatorprey ecosystems therein, it is interesting to note that very few previous studies have considered carrying capacity of the predator. So in this paper, we will use the Leslie-Gower(LG) model with Holling type I functional response to describe the interaction between pests and natural enemies [24]: where x and y represent the population of pests and natural enemies, r and s are the intrinsic growth rates of the pest and the natural enemy respectively, the functional response mx indicates the consumption of the natural enemy to the pest; although pests and natural enemies both grow logistically, carrying capacity of the pest is a constant K while that of the natural enemy is proportional to the pest population. When the number of pests is less than ET, no control measure is carried out and the model(1.1) holds true. But as long as the pest population reaches or surpasses ET, we will apply IPM strategies to spray insecticide and release natural enemy simultaneously, and then the following extended LG model with IPM strategies is set up: where the killing rate of pests and the release rate of natural enemies are expressed by p 1 and p 2 separately. Consequently, the combination of models (1.1) and (1.2) is the Filippov system which we will investigate in the rest of this paper. Throughout this paper, in order to guarantee the existence of internal equilibria and sliding mode, we assume r > p 1 , which means that the killing rate caused by spraying pesticide is less than the intrinsic growth rate of pests; moreover, K > rET r−p1 is assumed, which also ensures the normal survival of the pest population in the ecosystem and makes the implementation of our control strategies needed and meaningful.
Our main purpose of this study is to probe into the global dynamics of model (1.1) with (1.2). Hence in section 2, we first give a detailed explanation of our Filippov system and present some useful preliminaries. Afterwards, dynamical behaviors of models (1.1) and (1.2) are investigated separately. In section 3, we exhibit sliding mode dynamics and the existence of five types of equilibria. Furthermore, local stability of pseudoequilibrium is analyzed. In section 4, the global dynamics are obtained based on the previous discussion, and the results of numerical simulation are shown to verify theoretical analysis. At last, biological implications are provided in section 5.

Preliminary results
We apply the threshold policy into LG model in order to manage the pest population more appropriately, and then the combination of models (1.1) and (1.2), namely the Filippov system, is acquired as follows: with where p 1 x represents the impact on pest population caused by chemical method, spraying insecticides, and p 2 y is the number of releasing natural enemies. Due to the property of this model, let D = {(x, y) ∈ R|x > 0, y ≥ 0} be the field of definition of system(2.1) with (2.2). Then denote Z=(x,y), a smooth scale function H(Z)=x-ET and Hence, the system(2.1) and (2.2) can be rewritten as a generalized Filippov system shown below: where In particular, the discontinuous boundary between the region S 1 and S 2 is defined as Σ = {Z ∈ D| H(Z) = 0}, so we have D = S 1 ∪Σ ∪S 2 . In the following, we call the system(2.3) defined in the region S 1 or S 2 as subsystem F 1 or subsystem F 2 respectively.
For the discontinuous boundary Σ, it can be segmented into the following regions: (i) Sliding region Σ S : F 1 (Z), H Z (Z) > 0 and F 2 (Z), H Z (Z) < 0 hold true, which means that trajectories that once touch Σ S will stay on it, (ii) Crossing region Σ C : F 1 (Z), H Z (Z) F 2 (Z), H Z (Z) > 0 holds true, which implies that trajectories can pass through Σ C to another region, (iii) Escaping region Σ E : F 1 (Z), H Z (Z) < 0 and F 2 (Z), H Z (Z) > 0 hold true, which suggests that trajectories on Σ E will move towards either S 1 or S 2 , where ·, · represents the standard scalar product and the non-vanishing gradient of H(Z) is expressed by H Z (Z). However, note that the escaping region is impossible in this paper as F 1 (Z), H Z (Z) < 0 and F 2 (Z), H Z (Z) > 0 cannot hold simultaneously.
In order to investigate the sliding mode dynamics for Filippov system(2.3), we need to know exactly about the sliding domain on Σ, so the renowned Filippov's convex method [25] and Utkin's equivalent control method [26] are both feasible. Let then the sliding mode domain can be determined by{Z ∈ D| σ(Z) < 0} due to the absence of the escaping region.
Next, the definitions of several types of equilibria for system(2.3) are introduced.
Definition 2.2. If the point E V is satisfied with the equation F i (E V ) = 0 but located in another region S j , i,j=1,2, i = j, it is a virtual equilibrium for system(2.3). Definition 2.3. If an equilibrium E p is on the sliding mode Σ S of system(2.3), that is, the equation it is a pseudo-equilibrium for system(2.3).
Definition 2.4. If a point E b is an equilibrium of any subsystem(F 1 or F 2 ) but belongs to the separating boundary Σ, it is called a boundary equilibrium.
In other words, a boundary equilibrium E b satisfies 3) is the one which belongs to the separating boundary Σ, i.e., H(E T ) = 0 and satisfies For the subsystem F 1 , since Huang [24] has analyzed this class of models, we only present main results from that paper here. This model has an equilibrium E 10 = (K, 0), which is a saddle and unstable all the time, and an unique positive equilibrium , which is locally aymptotically stable. And in region S 2 , the subsystem F 2 exists two equilibria as well, E 20 = K 1 − p 1 r , 0 , and an unique internal equilibrium Proof. For the function f 1 , we can get and for another function f 2 , we have Therefore, Z(t) is positive for all t ∈ [0, T ) as long as the initial solution satisfying x(0) > 0 and y(0) > 0. Proof. For the first function in system(2.1), we have which means that dx dt ≤ 0 as long as x = K. And from the above equation, we can get Similarly, combining the above result, for the second function f 2 in system(2.1), we can obtain thus is a positively invariant set and is attracting likewise.
Theorem 2.1. The equilibrium E 1 is globally asymptotically stable in region S 1 .
Proof. Here we choose a Dulac function B(x, y) = 1 xy , we have Therefore, by the Bendixson-Dulac criterion, there are no limit cycles in region S 1 . Based on the local stability of E 1 , hence E 1 is a globally asymptotically stable equilibrium.
Before analyzing properties of system(2.3), we give some assumptions as follows:

Lemma 2.3
There are no closed orbits existing in region S 2 .
Proof. We use the same Dulac function B(x, y) = 1 xy , and here we can get Thus, according to the Bendixson-Dulac criterion, there are not any closed orbits in region S 2 .
holds, E 20 is globally asymptotically stable.
Proof. Consider the Jacobian matrix of system (2.3) at E 20 : then we can get eigenvalues straightforward: where λ 1 < 0 as r > p 1 . Therefore, if λ 2 > 0, two eigenvalues have opposite signs, that is, E 20 is a saddle. And if λ 2 < 0, both eigenvalues are negative, so E 20 is globally asymptotically stable based on the Lemma 2.3. Theorem 2.3. The internal equilibrium E 2 is globally asymptotically stable in region S 2 .
Proof. For simplicity, we denote that q = hrs+mK(s+ . By calculation, the Jacobian matrix of system(2.3) is obtained as follows: from which we get the characteristic equation where, Notice that B > 0 when So here we substitute values of q, q i , i = 1, 2, into the left-hand side of the above inequality and simplify it: Obviously this inequality holds, that is, B > 0 is true.
Similarly, substitute values of q, q i , i = 1, 2, into the formula of C: > 0, which means that C > 0 is established. Hence, according to Routh-Hurwitz criteria, the internal equilibrium E 2 is locally asymptotically stable, and combined with Lemma 2.3, E 2 is globally asymptotically stable in region S 2 .

Dynamics of sliding mode and equilibria
In this section, we initially pay attention to the existence of the sliding domain, and then we study sliding mode dynamics. In the end, different types of equilibria are discussed.
In view of preliminaries, the sliding mode are defined as As a result, σ(Z) < 0 is equivalent to By calculation, we get y min < y < y max , and y min = r m Hence the sliding segment of Filippov system(2.3) can be expressed by Note. There is no escaping region for Filippov system(2.3), as it is impossible that the two inequalities H Z (Z), F S1 (Z) < 0 and H Z (Z), F S2 (Z) > 0 are satisfied in the meantime. Now, we use Utkin's equivalent control method [26] to probe for the dynamics of the sliding segment.
In light of this method, we can obtain the differential equation for sliding dynamics on the segment Σ S . It follows from H = 0 and the first equation of Filippov system(2.3) that Solving the above equation respect to γ produces Substituting the formula of γ into the second equation in model(2.1) generates which determines the dynamics on Σ S .
The following is to explore different kinds of equilibria. Based on the analysis of section 2, there is an unique positive equilibrium in each subsystem, that is, The discussion of these two equilibria types is shown below: (i) Since hrK mK+hr > hsK(r−p1) hrs+mK(s+p2) , it is unlikely that hrK mK+hr < ET < hsK(r−p1) hrs+mK(s+p2) holds. In other words, E 1 and E 2 cannot become real equilibria simultaneously.
(ii) When hrK mK+hr < ET , E 1 is a real equilibrium denoted by E 1 R , while E 2 is a virtual equilibrium as hrs+mK(s+p2) > ET , E 2 is a real equilibrium denoted by E 2 R , and E 1 is a virtual equilibrium expressed by E 1 V . (iv) Due to the fact that hsK(r−p1) hrs+mK(s+p2) < ET < hrK mK+hr is true, in this case, both equilibria are virtual recorded as E 1 V , E 2 V . Now, we probe for the pseudo-equilibrium. Solving the equation ψ(y) = 0 yields As y 1 < y min , we omit this solution. In order to make sure that y 2 ∈ (y min , y max ), it must have: during which we can obtain hsK(r − p 1 ) hrs + mK(s + p 2 ) < ET < hrK mK + hr , where both of internal equilibria happen to exist as two virtual equilibria E 1 V , E 2 V . Therefore, the pseudoequilibrium E p = (ET, y 2 ) cannot coexist with any regular equilibrium.
According to definition 2.4, the boundary equilibria satisfy x = ET and equations: Hence we can get two boundary equilibria below: .
Correspondingly, tangent points satisfy x = ET and equations: so we get two tangent points as shown below: which are endpoints of sliding segment Σ S .
Lemma 3.1. The pseudo-equilibrium E p is locally asymptotically stable respect to Σ S .
Proof. For the dynamic equation of Σ S , namely ψ(y), we can take the derivative with regard to y at E p : Therefore, on the basis of stability theory for ODEs, E p is locally asymptotically stable on Σ S .

Global dynamics and numerical simulation
In this section, we concentrate on the global dynamical behaviors that the Filippov system(2.3) can exhibit. Initially, we rule out the existence of three kinds of possible limit cycles, and then with combining properties of equilibria, the global stability of equilibria can be obtained.
According to the analysis in section 2, we can make sure that there is no limit cycle totally contained within the region S 1 and S 2 . Then we have the following lemmas.
Lemma 4.1. There is no closed orbit completely contained within the region S i (i=1,2).
Next, we give proofs of non-existence of the other limit cycles which respectively contain a part of the sliding segment and the entire sliding mode.

Lemma 4.2.
There exists no closed orbit which contains partial sliding domain.
Proof. We prove it by reduction to absurdity. There are two cases in all, namely, the pseudo-equilibrium E p exists or not. So when E p exists on the sliding segment, it is stable such that all solutions nearby will converge to this point finally. However, suppose that there is a limit cycle including part of the sliding domain, some trajectories that once touch the sliding domain converge to this limit cycle, which contradicts the stability of the pseudo-equilibrium.
On the other hand, when E p does not exist, there are a real equilibrium and a virtual equilibrium coexisting in a certain region. Without loss of generality, assume that this kind of limit cycle denoted by L contains a real equilibrium E 1 in the region S 1 . This means that any trajectory outside L will converge to itself and cannot reach E 1 anymore, which contradicts the global stability of the real equilibrium E 1 in the region S 1 . Hence, there does not exist any closed orbit including part of the sliding domain.  Proof. Suppose that there is a closed orbit Γ = Γ 1 +Γ 2 around the sliding segment (shown as Fig1), where Γ i = Γ S i , i = 1, 2. And Γ intersects with the switching line x = ET at points P and Q. Denote the region enclosed by Γ as D and D i = D S 1 , i = 1, 2. Then we have two auxiliary lines x = ET − θ and x = ET + θ, which intersect with Γ at points P 1 , Q 1 and P 2 , Q 2 respectively. Meanwhile, P 1 = P − a 1 (θ), Q 1 = Q − b 1 (θ), P 2 = P + a 2 (θ), Q 2 = Q + b 2 (θ), where a i (θ) and b i (θ) (i = 1, 2) are continuous with respect to θ and lim θ→0 a i (θ) = lim θ→0 b i (θ) = 0, i = 1, 2. This means that P 1 Q 1 and P 2 Q 2 tends to P Q with θ → 0. Moreover, describe two regions delimited by Γ 1 , P 1 Q 1 and Γ 2 , P 2 Q 2 as › D 1 and › D 2 respectively, which severally tend to D 1 and D 2 when θ → 0.
Let the Dulac function be B = 1 xy , then we can obtain the following results by using Green's theorem.
After the above analysis, it is easy to get that that is, we denote the result of the above equation by α.
By calculations, we have which means that α < 0 is right. Then, according to Green's theorem, we have: as Q is greater than P along the vertical axis. This causes a contradiction. Hence we can obtain the conclusion as Lemma 4.3.
In what follows, we investigate the global dynamics of Filippov system(2.3), and then the global stability of internal equilibria and pseudo-equilibrium is obtained with the support of numerical simulation.
Theorem 4.1. The pseudo-equilibrium E p is globally asymptotically stable as long as it exists. Proof. From the analysis in section 3, we know that two internal equilibria E 1 , E 2 are both virtual denoted as E 1 V and E 2 V when E p exists. This implies that the internal equilibrium E 1 belonging to the region S 1 is in the region S 2 , and E 2 is an internal equilibrium of the region S 2 but is in the region S 1 . Therefore, all trajectories move towards the opposite domain to reach their own equilibrium and they are bound to touch the switching line x=ET. Since E p is LAS on the sliding segment and there does not exist any limit cycle in terms of Lemmas 4.1-4.3, some trajectories which hit the sliding mode Σ S in the switching line will slide to E p , while other trajectories which pass through the switching line to the opposite area will return and ultimately approach to E p after colliding with the sliding segment. So it follows that E p is globally asymptotically stable when it is present. The result of numerical simulation is shown in the figure below.  Proof. In the light of the discussion in section 3, when hrK mK+hr < ET , E p does not exist, and E 1 is a real equilibrium expressed by E 1 R located in the region S 1 and E 2 is a virtual equilibrium described as E 2 V situated at the region S 1 as well. Meanwhile, on the basis of Theorem 2.1, E 1 is GAS in the region S 1 , and there is no closed orbit for this Filippov system(2.3) by Lemmas 4.1-4.3.
Consequently, as Fig 3 shows, trajectories starting from the region S 1 will either approach to E 1 R straightforward or intersect with the sliding segment Σ S on the switching line and slide to the top endpoint along this segment, then get into the region S 1 again to reach E 1 R . On the other hand, trajectories whose initial points are in the region S 2 will firstly move towards S 1 to reach their own equilibrium E 2 . Hence these trajectories will either pass through the switching line to S 1 and go towards E 1 R or collide with Σ S and slide up to the end, then approach to E 1 R in the region S 1 . In view of this, all trajectories will reach E 1 R ultimately, that is, E 1 R is GAS for system(2.3). Theorem 4.3. If hsK(r−p1) hrs+mK(s+p2) > ET holds, the internal equilibrium E 2 is globally asymptotically stable for the Filippov system(2.3).
Proof. Combining the analysis in section 3, we have that when hsK(r−p1) hrs+mK(s+p2) > ET holds, there are two equilibria in total, namely, E 2 R and E 1 V both in the region S 2 . Moreover, according to Theorem 2.3, E 2 is GAS in the region S 2 , and no closed orbit exists for system(2.3) by Lemmas 4.1-4.3.
Thus, as Fig 4 shows, trajectories initiating from the region S 2 will tend to E 2 R directly, or once they hit the sliding segment, they will swipe down along this segment to the endpoint and then enter the region S 2 to reach E 2 R . Furthermore, trajectories starting from the region S 1 are attracted by E 1 in the opposite region. So they will either go across the switching line and move towards E 2 R or collide with the sliding segment and slide down to approach to E 2 R in the region S 2 . It then follows that all solutions will converge to E 2 R finally, in other words, E 2 R is a GAS node for the system(2.3).
Finally, we find that the parameter ET plays an important role in the pest control. Denote rhK mK+rh and hsK(r−p1) hrs+mK(s+p2) as A and B respectively, and A > B holds based on the analysis in section 3. By employing numerical simulation, the tendency of solutions with different ET are exhibited in Fig 5. We can see clearly that only if ET > B holds, all trajectories would approach to either the pseudo-equilibrium on the switching line Σ or the internal equilibrium in the region S 1 , which implies that the pest population is controlled and eventually stabilize at or below ET, in other words, our control tactics work.

Biological implications and conclusions
Nowadays, the Filippov system has been found to be useful to describe the real-world problems and investigated in many fields [27][28][29]. Compared with the impulsive differential equations, the non-smooth Filippov system is more close to the reality in describing pest control measures. And despite many papers have studied Leslie-Gower model, none of them have considered the effect of control strategies. Hence we expand the LG model with Holling I functional response into a piecewise smooth system by applying IPM strategies, which can be called as a Filippov system with the threshold policy as well. This means that we do not take any action when the number of pests is below the economic threshold(ET), but in case the quantity reaches or exceeds ET, we will spray pesticides and release the natural enemy in the meantime. In order to minimize the economic toll and maximize the crop yield, our goal in this paper is to get the conditions under which the pest population can stabilize at or below the value given in advance. So we make use of Filippov theories and qualitative techniques with numerical simulations to investigate dynamical behaviors of proposed system in detail, including global dynamics of subsystems, the existence of sliding mode and different types of equilibria, sliding mode dynamics and the global stability of equilibria.
It is worth mentioning that no possible limit cycle exists for our proposed system and all trajectories are able to converge to the certain equilibrium finally under some conditions rather than exhibiting divergence or oscillation, which means that our strategies are feasible. Moreover, the results of global stability are concentrated in Theorems 4.1-4.3 and Figs 2-4. Especially, combining the outcome of numerical simulation in Fig  5, it is realized that the selection of ET plays a strong part in the pest control, as this offers the condition of choosing ET. Consequently, our findings are valuable for how to draw up strategies effectively and when to take measures.