Threshold control strategy for a Filippov model with group defense of pests and a constant-rate release of natural enemies

1.   Introduction

Since ancient times, the struggle between human beings and pests has never been interrupted, which has brought serious impacts on agricultural production, so pest control has always been a hot issue. Particularly, rice planthoppers, armyworms, aphids, cotton bollworms, meadow moths, corn borers, tea silkworms and locusts have become a common and widespread problem in recent years ^[1,2], with part of them exhibiting group defense. Group defence is a term used to describe the phenomenon whereby predation is decreased, or even prevented altogether, due to the increased ability of the prey to better defend or disguise themselves when their density is large enough ^[3,4]. This phenomenon described by Holmes and Bethel involves certain insect populations ^[5]. For example, large swarms of the locusts can release a lot of toxic substances, which helps them escape natural enemies, and then the locusts can reproduce in large numbers, causing pest population outbreaks, which are extremely destructive to crops. In early 2020, the worst outbreak of desert locusts in 70 years destroyed tens of thousands of hectares of farmland in East Africa, posing a huge threat to global food security ^[6,7]. Obviously, the ability of the group defense of pests has brought a great challenge to pest control.

When controlling pests, it is necessary to take into account integrated pest management strategies ^{[8,9,10,11,12]} in order to reduce the density of a pest population to a tolerable level, which involves choosing appropriate measures such as spraying pesticides and releasing natural enemies. One of the main purposes of an integrated pest management strategy is to take control measures when the pest population density exceeds the economic threshold (ET) ^[13,14] so as to prevent the pest population density from reaching the economic injury level ^[14,15]. Therefore, a semi-continuous differential dynamic system, namely, a differential dynamic system with impulsive state feedback control has been proposed ^{[16,17,18,19,20]}. This can not only avoid unnecessary environmental pollution and resource waste caused by the abuse of pesticides and release of natural enemies, but it also protects the diversity and stability of the ecosystem. The basic premise of the above model is that the control measures are completed in an instant. However, in the practical application process, once the pest population density exceeds the ET, people will continue to use pesticides for a period of time or continuously release natural enemies to make the pest population density equal to the ET or fall below the ET; otherwise, no control measures will be taken. This control strategy is called the threshold control strategy or the switching strategy, which is characterized by the switched Filippov differential dynamical system ^{[21,22,23,24]}. In recent years, piecewise smooth Filippov systems have been widely used in many fields, such as control engineering, biology, medicine, physics and economics ^{[25,26,27,28,29,30,31,32,33,34]}.

Based on the predator-prey model with group defense proposed by Xiao et al. ^[35], Hou et al. ^[36] established and studied a Filippov system to simulate the process of integrated pest management in which pest populations had group defense behaviors and natural enemies were released at a linear rate. The results showed that, although there were no closed orbits in the subsystem, a stable periodic solution may exist for a Filippov system after switching-induced perturbation. The linear release rate of natural enemies in ^[36] requires us to monitor the density of natural enemies to control pests. However, in the real world, it is difficult to do that. In this paper, we study an integrated pest management Filippov model with group defense and a constant-rate release of natural enemies to simulate the pest control process. Our main purpose of this paper is to study the effects of a threshold control strategy in a Filippov system on the dynamics of subsystems and the pest control.

In the following section, a Filippov pest control model with group defense behavior is proposed. Then, relevant definitions and preliminaries of the Filippov system are given in Section 3. By using Filippov's theory and numerical simulations, the comprehensive dynamics of the Filippov system are analyzed in Section 4. A brief discussion is given in the last section.

2.   Model formulation

Based on the model proposed by Xiao et al. ^[35], it is assumed that the pest population exhibits group defense behavior against natural enemies, which is described by the following model:

Once the pest population density exceeds the ET, the integrated pest management strategy is implemented, that is, a constant rate of natural enemies and insecticide spraying is applied at the same time; otherwise no intervention measures will be taken. Therefore, when $x(t) > ET$, the following model is established.

where $x(t)$ and $y(t)$ represent the densities of the pest and natural enemy populations at time $t$, respectively; $K > 0$ is the carrying capacity of the pest population and $r > 0$ represents its intrinsic growth rate; $\alpha{x(t)}e^{-\beta{x(t)}} ({\alpha} > 0, {\beta} > 0)$ represents the predation rate of natural enemies against pests with group defense behavior, which is a non-monotonic saturating response function; $d > 0$ represents the death rate of natural enemies; $k > 0$ represents the conversion rate of natural enemies preying on pests; $p_1 > 0$ and $p_2 > 0$ represent the reduction rates of the densities of pests and natural enemies due to the effects of pesticides, respectively; $\tau > 0$ represents the constant release rate of natural enemies, and as time passes, natural enemies are continuously and linearly released. The constant $\tau$ represents the degree of control of pests. In other words, the larger the constant, the more significant the reduction in the pest population through the release of natural enemies. Therefore, by studying and adjusting this constant, we can better understand and optimize pest management strategies to improve their effectiveness and sustainability. It is assumed that the effect of pesticides on natural enemies is ignored, so the value of ${p_2}$ is approximately $0$.

Then, models (2.1) and (2.2) can be rewritten as the following generalized non-smooth Filippov system ^[22,24]:

where $Z = (x, y)^T, F_{{M}_{1}} = (F_{11}, F_{12})^T, F_{{M}_{2}} = (F_{21}, F_{22})^T, M_1 = \{Z\in{R_+^2}|H(Z) < 0\}, M_2 = \{Z\in{R_+^2}|H(Z) > 0\}, H(Z) = x(t)-ET$ and $R_+^2 = \{Z = (x, y)^T|x\geq0, y\geq0\}.$ Note that $\Sigma = \{(x, y)\in{R_+^2}|H(Z) = 0\}$ represents the dividing line of $M_1$ and $M_2$.

3.   Preliminaries

First, we give the sliding mode of the Filippov system (2.3), which can be defined as

where $\langle \cdot \rangle$ represents the standard scalar product, $H(Z)$ is a smooth scalar function and the gradient of $H(Z)$ on $\Sigma$ is ${H_Z} = (1, 0)$. Note that a sliding mode is stable for $\langle {{H_Z}, {F_{{M_1}}}(Z)} \rangle > 0$ and $\langle {{H_Z}, {F_{{M_2}}}(Z)} \rangle < 0$; correspondingly, when $\langle {{H_Z}, {F_{{M_1}}}(Z)} \rangle < 0$ and $\langle {{H_Z}, {F_{{M_2}}}(Z)} \rangle > 0$, the sliding mode is unstable ^[24]. Obviously, based on the right-hand expressions of system (2.1) and system (2.2), it can be concluded that there is no type of unstable sliding mode. Therefore, we only discuss stable cases.

A sliding mode may occur if there exist subsets $\Sigma_s$ of the manifold $\Sigma$ such that the vectors ${F_{{M_1}}}(Z)$ and ${F_{{M_2}}}(Z)$ are directed toward each other. By Filippov's convex theory ^[24], if the sliding mode is smooth on a one-dimensional sliding interval $\Sigma_s$, then the sliding mode dynamics can be determined by

where $0 < \lambda {\rm{ = }}\dfrac{{\langle {{H_Z}, {F_{{M_2}}}(Z)} \rangle }}{{\langle {{H_Z}, {F_{{M_2}}}(Z) - {F_{{M_1}}}(Z)} \rangle }} < 1$ and $Z = {(x, y)^T} \in {\Sigma _s}$.

In a Filippov system, there exist different types of equilibria, which are defined as follows:

Definition 1 ^[24] (1) Let $F_{{M}_{i}}(E) = 0 (i = 1, 2);$ then, a point $E$ is called a regular equilibrium of system (2.3) if it belongs to $M_i$, and a virtual equilibrium if it belongs to $M_j, j\neq{i}$, as denoted by $E^r, E^v,$ respectively.

(2) A point $E$ is called a pseudo-equilibrium of system (2.3) if it is an equilibrium of the sliding mode $\Sigma_s,$ i.e., $\lambda{F_{M_1}}(E)+(1-\lambda){F_{M_2}}(E) = 0, \; H(E) = 0$ and $0 < \lambda < 1$, as denoted by $E_p.$

(3) A point $E$ is called a tangent point of system (2.3) if $E\in\Sigma_s, \langle{H_Z(E), F_{{M}_{1}}(E)}\rangle = 0$ or $\langle{H_Z(E), F_{{M}_{2}}(E)}\rangle = 0,$ as denoted by $E_t.$

(4) A point $E$ is called a boundary equilibrium of system (2.3) if $E\in\Sigma_s, F_{{M}_{1}}(E) = 0$ or $F_{{M}_{2}}(E) = 0,$ as denoted by $E_b.$

The global sliding bifurcations of a Filippov system can be divided into the following types:

Definition 2 ^[24] (1) If the standard limit circle of the Filippov system collides with the discontinuity boundary, it is called grazing bifurcation.

(2) If a stable sliding limit cycle of the Filippov system collides with another boundary point of the system, and then contains the entire sliding mode, it is called the buckling bifurcation.

(3) If the sliding limit cycle becomes a crossing limit cycle in the Filippov system with the parameters changing, it is called crossing bifurcation.

4.   Dynamics of Filippov system

4.1.   Dynamics of subsystems

The Lambert W function ^[37] plays an important role in solving the equilibria of subsystem (2.1). Now, we introduce the following definition:

Definition 3 The multivalued inverse of the function $z\rightarrow ze^{z}$ is called the Lambert W function, satisfying

It can be calculated that the derivative of the Lambert W function with respect to the variable $z$ is

The inverse function of $ze^z$ is denoted as $W(z)\doteq W(0, z)$, as defined on the interval $[-e^{-1}, +\propto)$ for $z\in [-1, +\propto);$ Similarly, the inverse function of $ze^z$ is denoted as $W(z)\doteq W(-1, z)$, as defined on the interval $[-e^{-1}, 0)$ for $z\in [-\propto, -1),$as shown in Figure 1. For further details about the background and properties of the Lambert W function, see ^[37].

When $x(t) < ET$, the subsystem (2.1) has a trivial equilibrium $(0, 0)$ and a boundary equilibrium $(K, 0)$. According to the definition of the Lambert W function, if $e\beta d < k\alpha$, the equation $k\alpha x(t)e^{-\beta x(t)}-d = 0$ has two positive roots, that is, the subsystem (2.1) has two positive equilibria ${E_{11}}({x_{11}}, {y_{11}}), {E_{12}}({x_{12}}, {y_{12}}),$ where

The stability of the equilibria and the existence of limit cycles in subsystem 2.1 are discussed in detail from ^[36], which will not be elaborated here.

When $x(t) > ET$, subsystem (2.2) has a boundary equilibrium, that is, the pest extinction equilibrium $(0, \frac{\tau }{d}).$ It is easy to know that the vertical and horizontal isoclines of subsystem (2.2) are $y = \frac{{r(1 - \frac{x}{K}) - {p_1}}}{{\alpha {e^{ - \beta x}}}}$ and $y = \frac{\tau }{{d - k\alpha x{e^{ - \beta x}}}}$, respectively. Combining the two formulas, we know that the positive equilibrium of subsystem (2.2) satisfies the following equation:

Lemma 4.1 If $\frac{\tau}{r-p_1} < \frac{d}{\alpha} < \frac{k}{\beta}$, then Eq (4.1) has at least one positive real solution.

Proof. Equation (4.1) can be transformed into

and then it is further transformed into $Ax+B = \frac{1}{e^{\beta x}-Cx}$, where $A = \frac{-rd}{\tau \alpha k} < 0$, $B = \frac{d(r-p_1)}{\tau \alpha}$ and $C = \frac{\alpha k}{d}$.

We first discuss the existence of a positive real solution of equation $e^{\beta x}-Cx = 0$. Let $g(x) = e^{\beta x}-Cx$; therefore, the derivative of $g(x)$ with respect to $x$ is $g'(x) = \beta e^{\beta x}-C.$ From the condition $\frac{d}{\alpha} < \frac{k}{\beta}$, one can derive the equation $g'(x) = 0$ to have the unique positive real root $x = \frac{1}{\beta} \ln \frac{C}{\beta}.$ Thus we conclude that $g'(x) < 0$ on the interval $[0, \frac{1}{\beta} \ln \frac{C}{\beta})$, which implies that $g(x)$ is decreasing on $[0, \frac{1}{\beta} \ln \frac{C}{\beta})$, and that $g'(x) > 0$ on the interval $[\frac{1}{\beta} \ln \frac{C}{\beta}, + \infty)$, which implies that $g(x)$ is increasing on $[\frac{1}{\beta} \ln \frac{C}{\beta}, + \infty)$. Furthermore, the minimum value of the function $g(x)$ is $g(\frac{1}{\beta} \ln \frac{C}{\beta}) = \frac{C}{\beta}(1-\ln \frac{C}{\beta})$. Thus, the existence of the positive real solution of the equation $g(x) = 0$ can be divided into three cases:

$(A_1)$ if $1-\ln \frac{C}{\beta} > 0$, that is $\frac{C}{\beta} = \frac{\alpha k}{\beta d} < e$, then the equation has no positive roots;

$(A_2)$ if $1-\ln \frac{C}{\beta} = 0$, that is $\frac{\alpha k}{\beta d} = e$, then the equation has one positive root;

$(A_3)$ if $1-\ln \frac{C}{\beta} < 0$, that is $\frac{\alpha k}{\beta d} > e$, then the equation has two positive roots.

For convenience, denote $h(x) = Ax+B$ and $f(x) = \frac{1}{g(x)}$, whereby the y-intercept of $h(x)$ is $B$ and the y-intercept of the function $f(x)$ is $1$. Now, we discuss the intersection of $h(x)$ and $f(x)$ in different cases. For case $(A_1)$, $f(x)$ is a continuous function on the interval $[0, +\infty)$, and the x axis is the asymptote in the first quadrant. According to $\frac{\tau}{r-p_1} < \frac{d}{\alpha}$, we know that $B > 1$. Since $A < 0$, $h(x)$ and $f(x)$ must intersect (see Figure 2(a)). For case $(A_2)$, the $f(x)$ has one discontinuity point. For the sake of convenience, we establish that $x = \alpha$ is the vertical asymptote. When $0 \leq x < \alpha$, $f'(x) > 0$, and when $x\rightarrow \alpha^{-}$, $f(x) \rightarrow + \infty$. Thus, if $A < 0$ and $B > 1$, that is, $\frac{\tau}{r-p_1} < \frac{d}{\alpha} < \frac{k}{\beta}$, then $h(x)$ and $f(x)$ must also intersect (see Figure 2(b)). For case $A_3$, we can discuss it in a similar manner as case $A_2$, and we will not elaborate on it here. This completes the proof.

Remark 4.1: As the intercept of function $h(x)$ with the x-axis increases, the number of intersection points between $h(x)$ and $f(x)$ will increase. As shown in Figure 2, it can be seen that there can be, at most, three intersection points.

According to the above analysis, the existence of positive equilibria of subsystem (2.2) has been proved. When the subsystem (2.2) has only one positive equilibrium, we can analyze the stability of this equilibrium via the vector field of subsystem (2.2). As shown in Figure 3, the isolines $x'(t) = 0$ and $y'(t) = 0$ divide the phase plane into four regions: $G_1$, $G_2$, $G_3$ and $G_4$. In the region $G_1$, $x'(t) < 0$ and $y'(t) > 0$; this means that the trajectories move upward to the left and are tangent to $y'(t) = 0$ into region $G_2$. In the region $G_2$, $x'(t) < 0$ and $y'(t) < 0$; this means that the trajectories go down to the left and are tangent to $x'(t) = 0$ into region $G_3$. In the region $G_3$, $x'(t) > 0$ and $y'(t) < 0$; this means that the trajectories move to the right and are tangent to $y'(t) = 0$ into region $G_4$. In the region $G_4$, $x'(t) > 0$ and $y'(t) > 0$; this means that the trajectories move to the upper right and are tangent to $x'(t) = 0$ into region $G_1$ again. According to the direction of the track, as shown in Figure 3, it can be inferred that the trajectories tend to the equilibrium. Therefore, the equilibrium is an asymptotically stable focus.

4.2.   Sliding mode dynamics and equilibria of Filippov system (2.3)

According to the definition of the sliding mode in the Filippov system, we have

For convenience, we denote ${y_{{c_1}}}\buildrel.\over = \frac{{r{e^{\beta ET}}(1 - \frac{{ET}}{K})}}{\alpha } - \frac{{{p_1}{e^{\beta ET}}}}{\alpha }, \; {y_{{c_2}}}\buildrel.\over = \frac{{r{e^{\beta ET}}(1 - \frac{{ET}}{K})}}{\alpha }.$ Therefore, the sliding mode of system (2.3) can be denoted as ${\Sigma_s} = \left\{ {{{(x, y)}^T} \in \Sigma|{y_{{c_1}}} \le y \le {y_{{c_2}}}} \right\}.$ From the Filippov convex theory, the dynamics of the sliding mode system are determined by the following equation:

Let $\psi (y){\rm{ = }}0$; we get ${y_{{Q_1}}} = \frac{{r\tau (K - ET){e^{\beta ET}}}}{{\alpha K(\tau - {p_1}kET) + {p_1}dK{e^{\beta ET}}}}$. If ${y_{{Q_1}}} \in ({y_{{c_1}}}, {y_{{c_2}}})$, that is

then Filippov system (2.3) has a unique pseudo-equilibrium ${E_p}(ET, {y_{{Q_1}}})$. Differentiate $\psi (y)$ with respect to $y$; then, we have

Therefore, if $\psi '(y) < 0$, that is, if the condition

holds, then the pseudo-equilibrium ${E_p}(ET, {y_{{Q_1}}})$ is locally asymptotically stable.

If $x(t) < ET$, subsystem (2.1) may have two internal equilibria ${E_{11}}({x_{11}}, {y_{11}})$ and ${E_{12}}({x_{12}}, {y_{12}})$, where ${x_{11}} < {x_{12}}.$ We have the following results about Filippov system (2.3).

(1) If $ET < {x_{11}}$, ${E_{11}}({x_{11}}, {y_{11}})$ and ${E_{12}}({x_{12}}, {y_{12}})$ are virtual equilibria, denoted by $E_{11}^v$ and $E_{12}^v$, respectively.

(2) If ${x_{11}} < ET < {x_{12}}$, ${E_{11}}({x_{11}}, {y_{11}})$ is a regular equilibrium and ${E_{12}}({x_{12}}, {y_{12}})$ is a virtual equilibrium, denoted by $E_{11}^r$ and $E_{12}^v$, respectively.

(3) If $ET > {x_{12}}$, ${E_{11}}({x_{11}}, {y_{11}})$ and ${E_{12}}({x_{12}}, {y_{12}})$ are regular equilibria, denoted by $E_{11}^r$ and $E_{12}^r$, respectively.

If $x(t) > ET$, assume that subsystem (2.2) have three positive equilibria, ${E_{21}}({x_{21}}, {y_{21}}),$ ${E_{22}}({x_{22}}, {y_{22}})$ and ${E_{23}}({x_{23}}, {y_{23}})$, where ${x_{21}} < {x_{22}} < {x_{23}}$. We have the following results about Filippov system (2.3).

(1) If $ET < {x_{21}}$, ${E_{21}}({x_{21}}, {y_{21}})$, ${E_{22}}({x_{22}}, {y_{22}})$ and ${E_{23}}({x_{23}}, {y_{23}})$ are regular equilibria, denoted by $E_{21}^r$, $E_{22}^r$ and $E_{23}^r$, respectively.

(2) If ${x_{21}} < ET < {x_{22}}$, ${E_{21}}({x_{21}}, {y_{21}})$ is a virtual equilibrium, and ${E_{22}}({x_{22}}, {y_{22}})$ and ${E_{23}}({x_{23}}, {y_{23}})$ are regular equilibria, denoted by $E_{21}^v$, $E_{22}^r$ and $E_{23}^r$, respectively.

(3) If ${x_{22}} < ET < {x_{23}}$, ${E_{21}}({x_{21}}, {y_{21}})$ and ${E_{22}}({x_{22}}, {y_{22}})$ are virtual equilibria, and ${E_{23}}({x_{23}}, {y_{23}})$ is a regular equilibrium, denoted by $E_{21}^v$, $E_{22}^v$ and $E_{23}^r$, respectively.

(4) If $ET > {x_{23}}$, ${E_{21}}({x_{21}}, {y_{21}})$, ${E_{22}}({x_{22}}, {y_{22}})$ and ${E_{23}}({x_{23}}, {y_{23}})$ are virtual equilibria, denoted by $E_{21}^v$, $E_{22}^v$ and $E_{23}^v$, respectively.

The tangent point of ${\Sigma _s}$ satisfies the following equations:

so, the tangent points of the Filippov system (2.3) are ${E_t}_{_1}(ET, {y_{{c_2}}})$ and ${E_t}_{_2}(ET, {y_{{c_1}}})$.

The boundary equilibrium of ${\Sigma _s}$ satisfies the following equations:

If $d = k\alpha ET{e^{ - \beta ET}}$, Filippov system (2.3) has a boundary equilibrium ${E_b}_{_1}(ET, {y_{{c_2}}})$; if $y = \frac{\tau }{{d - k\alpha ET{e^{ - \beta ET}}}}$, Filippov system (2.3) has a boundary equilibrium ${E_b}_{_2}(ET, {y_{{c_1}}})$.

4.3.   Global stability of Filippov system (2.3)

Select the parameters as those in Figure 4, where subsystem (2.1) has a globally asymptotically stable limit cycle (see Figure 4(a)), and subsystem (2.2) has a global asymptotic stability focus (see Figure 4(b)). When the ET is small, Filippov system (2.3) has a virtual equilibrium $E_{21}^v$ and regular equilibrium $E_{21}^r$, which yields a globally asymptotically stable focus (see Figure 5(a)). When $ET$ increases to 7.1, the equilibrium $E_{21}$, the tangent point ${E_{{t_2}}}$ and the pseudo-equilibrium ${E_p}$ collide at the boundary equilibrium ${E_{{b_2}}}$ and all orbits approach ${E_{{b_2}}}$ (see Figure 5(b)). As the ET continues to increase, both equilibrium points $E_{11}$ and $E_{21}$ become virtual equilibria; in this case, there is a pseudo-equilibrium ${E_p}$, which is globally asymptotically stable (see Figure 5(c)). When $ET = 8.7$, the equilibrium $E_{11}$, tangent point ${E_{{t_1}}}$ and pseudo-equilibrium ${E_p}$ collide at the boundary equilibrium ${E_{{b_1}}}$ all orbits will approach ${E_{{b_1}}}$ (see Figure 5(d)). When $ET$ increases to 9.5, the unstable equilibrium $E_{11}^r$ and $E_{21}^v$ coexist and a stable sliding limit cycle intersecting with the sliding mode ${\Sigma _s}$ appears in Filippov system (2.3) (see Figure 5(e)). With the increase of $ET$, the sliding limit cycle will occur in the order of grazing bifurcation, buckle bifurcation and crossing bifurcation (see Figure 5(e)–(h)). When the ET is sufficiently large, a globally asymptotically stable standard limit cycle is completely located in the region ${M_1}$, and all orbits approach the stable standard limit cycle (2.3) (see Figure 5(i)).

4.4.   Existence region of sliding mode and pseudo equilibrium

From the analysis of the sliding mode ${\Sigma _s}$ and pseudo-equilibrium ${E_p}$, Filippov system (2.3) has rich and complex dynamic properties. The following numerical simulations method is used to discuss the effects of killing rate ${p_1}$ on the sliding mode ${\Sigma _s}$ and pseudo-equilibrium ${E_p}$. It can be seen in Figure 6 that, when the killing rate ${p_1} = 0.25$, with the gradual increase of $ET$, the sliding mode will gradually increase and the pseudo- equilibrium ${E_p}$ always exists. This shows that the density of the pest population is controlled to $ET$, which avoids the large-scale outbreak of the pest population. The pests and natural enemies will always be in a relatively stable coexistence state, which is beneficial to the stability of the ecological environment (see Figure 6(a)). When the killing rate ${p_1} = 0.06$, which is relatively small, with the gradual increase of $ET$, the pseudo-equilibrium changes from existence to non-existence, which may lead to the outbreak of the pest population (see Figure 6(b)).

4.5.   Boundary equilibrium bifurcation

In what follows, $ET$ is used as the bifurcation parameter to analyze the boundary equilibrium bifurcation in Filippov system (2.3). The boundary focus or boundary node bifurcation phenomenon occurs at $ET = {x_{21}}$ and $ET = {x_{11}}$, respectively. At this time, the tangent point ${E_t}$, the pseudo- equilibrium ${E_p}$ and the equilibrium ${E_{11}}$ (or ${E_{21}}$) all collide together at the boundary equilibrium ${E_b}$.

From Section 4.2, Filippov system (2.3) exhibits the boundary focus bifurcation phenomenon at the two endpoints of the sliding mode when $ET = 7.1$ and $ET = 8.7$, respectively (see Figure 5(b), (d)). Selecting the parameters in Figure 7, when $ET = 6.1$ and $ET = 8.7$, Filippov system (2.3) exhibits the boundary node bifurcation phenomenon at the two endpoints of the sliding mode (see Figure 7(b), (d)).

4.6.   Global sliding bifurcation

When $ET = 30$, an unstable regular equilibrium $E_{11}^r$ and a virtual $E_{21}^v$ of Filippov system (2.3) coexist and all orbits tend to be a stable standard limit cycle in $M_1$ (see Figure 8(a)). With the decrease of $ET$, when $ET = 28.8$, the standard limit cycle of Filippov system (2.3) collides with the upper tangent point $E_{t_1}$, that is, the grazing bifurcation phenomenon occurs (see Figure 8(b)). When the parameter $ET$ continues to decrease to $28$, the sliding mode ${\Sigma _s}$ becomes part of the limit cycle (see Figure 8(c)). Specifically, when $ET$ decreases to $25.7$, all orbits of system (2.3) tend to a stable sliding limit cycle, which is located in ${M_1}$ and contains the whole sliding mode ${\Sigma _s}$ (Figure 8(d)). In this case, the pest population density is controlled within the reasonable range of the ET, which is beneficial for pest control.

When $ET = 23$, the sliding limit cycle of system (2.3) is located in ${M_1}$ and ${M_2}$, which contains part of the sliding mode and the upper tangent point $E_{t_1}$ (Figure 8(e)). When $ET$ decreases to $22.2$, the sliding limit cycle is located in ${M_1}$ and ${M_2}$, passes through a tangent point ${E_{{t_1}}}$ and contains the entire sliding mode ${\Sigma _s}$ of Filippov system (2.3) (see Figure 8(f)). When the parameter $ET$ continues to decrease to $20$, Filippov system (2.3) still has an unstable regular equilibrium $E_{11}^r$ and a virtual equilibrium $E_{21}^v$, and it exhibits a stable crossing limit cycle containing the whole sliding mode (see Figure 8(g)). Obviously, Figures 8(c)–(e) depicts the process of buckling bifurcation, and Figures 8(e)–(g) depicts the process of crossing bifurcation.

In conclusion, with $ET$ decreasing from $30$ to $20$, the sliding limit cycle of system (2.3) changes gradually from being completely located in ${M_1}$ to containing part or the whole section of ${\Sigma _s}$, and to being partially located in ${M_2}$, and this process progresses from the grazing bifurcation to the complete buckling bifurcation phenomenon, and then to the complete crossing bifurcation phenomenon. It can be seen from the numerical simulations that, when $ET$ is small, the complex sliding bifurcation phenomenon occurs in system (2.3), which makes the density of the pest population sometimes exceed $ET$, brings uncertainty to pest control and is not conducive to the monitoring and control of pests (see Figure 8(e)–(g)).

5.   Discussion

In this paper, an integrated pest management Filippov model with group defense and a constant-rate release of natural enemies is established, and the complex dynamics of Filippov system (2.3) are studied. Combined with the conclusions of this paper, the establishment of ${x_{21}} > ET$ can be avoided by controlling the dose of sprayed pesticides or the release rate of natural enemies. Otherwise, $E_{21}^r$ becomes a globally asymptotically stable equilibrium. In this case, the density of the pest population exceeds $ET$, which may lead to pest population outbreak and is not conducive to pest control (see Figures 5(a) and 7(a)). With the gradual increase of pest population group defense ability $\beta$, ${x_{21}}$ also gradually increases (see Figure 9). This suggests that the stronger the group defense ability of the pest population, the more difficult it is to control pests. If a regular equilibrium that is globally asymptotically stable in both subsystems and a virtual equilibrium coexist, then the regular equilibrium will become the globally asymptotically stable equilibrium in Filippov system (2.3) after the switching disturbance. If a regular equilibrium that is unstable in both subsystems and a virtual equilibrium coexist, then Filippov system (2.3) may appear in a globally asymptotically stable limit cycle. When two virtual equilibria coexist, Filippov system (2.3) has a globally asymptotically stable pseudo-equilibrium (see Figures 5(c) and 7(c)). When the equilibrium reaches the dividing line $x(t) = ET$, the tangent point ${E_t}$, the pseudo-equilibrium ${E_p}$ and the regular equilibrium collide at the boundary equilibrium ${E_b}$. This collision results in a bifurcation at the boundary, which can take the form of either a focus or a node, depending on whether the equilibrium is itself a focus or a node (see Figures 5(b), (d) and 7(b), (d)). With the change of the ET, the limit cycle will generate complex global sliding bifurcation phenomena such as grazing bifurcation, buckling bifurcation, and crossing bifurcation. If the grazing bifurcation occurs naturally for the Filippov system as the limit cycle collides with the upper tangent point $E_{t_1}$ and the standard limit cycle is completely located in ${M_1}$, then this situation can control the density of the pest population within a reasonable range that does not exceed $ET$, which is conducive to integrated pest management. On the contrary, if there is buckling bifurcation, crossing bifurcation and standard limit cycles located in ${M_2}$, then these situations can cause the density of the pest population to exceed the ET, which is not conducive to pest control (see Figure 8). In this paper, we assume that if $x > ET$, we adopt a continuous threshold control strategy and the killing rate of pests is linear owing to the effects of pesticides. In fact, pesticides are sprayed periodically and have certain residual and delay effects on the pests. In this case, it becomes more difficult to analyze the dynamics of the system with the threshold control strategy. We will leave it for future work.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (12171004).

Conflict of interest

All authors declare no conflicts of interest regarding the publication of this paper.