Predicting effectiveness of anti-VEGF injection through self-supervised learning in OCT images

1.   Background and motivation

One of the most significant trends in global agricultural development is the ecological management of pests. From the perspective of ecosystem integrity, reducing and controlling pests through biological and ecological control are of great significance for the construction of ecological civilization. Biological and ecological control can reduce management cost, maintain ecological stability, and avoid environmental pollution and damage to biodiversity. As a large agricultural country, China places a premium on green prevention and control within its agricultural sector and proposed the National Strategic Plan for Quality Agriculture (2018–2022), which proposes to implement green prevention and control actions instead of chemical control and achieve a coverage rate of more than 50$\%$ for green prevention and control of major crop pests. The Crop Pests Regulations on the Prevention and Control of Crop Pests prioritizes the endorsement and support of green prevention and control technologies such as ecological management, fosters the widespread application of information technology and biotechnology, and propels the advancement of intelligent, specialized, and green prevention and control efforts ^[1]. Therefore, the simulation of pest dynamic behavior and the research of control strategies are helpful for more scientific and reasonable pest management.

In a natural ecosystem, the predator-prey relationship is one of the most important relationships, and has become a main topic in ecological research and widely studied by scholars in recent years. Depending on the problem under consideration and the biological background, related research can be divided into two forms: ordinary differential ^[2,3,4,5] and partial differential ^{[6,7,8,9,10,11]}. The earliest work on the mathematical modeling of predation relationships dates back to the twentieth century, named as the Lotka-Volterra model ^[12,13]. Subsequently, scholars have extended the Lotka-Volterra model in different directions such as introducing different types of growth functions ^[14,15,16] and different forms of functional response ^{[17,18,19,20]}. The Gompertz model ^[14] is one of the most frequently used sigmoid models fitted to growth data and other. Scholars have fitted the Gompertz model to everything from plant growth, bird growth, fish growth, and growth of other animals ^[21,22,23]. Compared with the logistic model, it is more suitable for pest or disease curve fitting with S-shaped curve asymmetry, and fast development at first and slow development later. In addition, the Holling-Ⅱ functional response function is the most commonly employed one, in which the searching rate is considered as a constant. Nevertheless, in the real world, the density of the prey and the predator's searching environment can affect the predator's searching speed. Consequently, Hassell et al. ^[24] proposed a saturated searching rate. Guo et al. ^[25] introduced a fishery model with the Smith growth rate and the Holling-Ⅱ functional response with a variable searching rate. In this work, a pest-natural enemy model with the Gomportz growth rate and a variable searching rate is investigated.

To prevent the spread of pests, effective control action should be implemented before the pests cause a certain amount of damage to the environment and crops. One way is to slow down the spreading speed of pests by setting a warning threshold, and when the density of pests exceeds this threshold level, an integrated control measure is imposed on the system. This kind of control system can be modeled by a Filippov system, which has been recognized by scholars and widely used in the study of concrete models with one threshold ^{[26,27,28,29,30,31]}, a ratio-dependent threshold ^[32], or two thresholds ^[33]. In this study, we will also focus on Filippov predation models with dual thresholds. In addition, considering the instantaneous behavior of the control, an integrated pest-management strategy with threshold control is adopted, which is an instantaneous intervention imposed on the system and always taken as a practical approach for pest management. In recent years, there has been a lot of research and application of impulsive differential equations (IDEs) in population dynamics to model the instantaneous intervention activities. There are mainly six types of models involved in the research: periodic ^{[34,35,36,37]}, prey-dependent ^{[38,39,40,41,42,43]}, predator-dependent ^[44], ratio-dependent ^[45], nonlinear prey-dependent ^[46], and combined prey-predator dependent ^{[47,48,49,50,51]}. In the context of integrated pest management, setting a threshold for pest population density to control its spread is crucial. Therefore, in this study, we introduce a pest economic threshold: when the pest population density exceeds this threshold, we will intervene manually, which includes not only spraying pesticides but also releasing natural enemies.

The article is organized in the following way: In Section 2, an integrated pest-management model with a variable searching rate based on double-threshold control is proposed. In Section 3, a dynamical analysis of the continuous system is performed, including the positivity and boundedness of the solutions, the existence and local stability of equilibrium points, and the dynamic behavior of the Filippov pest-management model with double thresholds. In Section 4, the complex dynamic behavior of the system induced by the economic threshold feedback control is focused on. In Section 5, numerical simulations are carried out to illustrate the main results of the above two sections step by step and to illustrate the practical implications. Finally, a summary of the research work is presented, and future research directions are discussed.

2.   Model formulation and preliminaries

2.1.   Pest control model

A pest-natural enemy Gomportz model with a variable searching rate and Holling-Ⅱ functional response is considered:

where $x$ ($y$) represents the pest's (natural enemies) density, respectively; $r$ represents the pest's intrinsic growth rate; $K$ represents the pest's environmental carrying capacity; $b(x) = bx/(x+g)$ represents the variable searching rate ^[24,25] with maximum searching rate $b$ and saturated constant $g$; $e$ represents the conversion efficiency; and $d$ represents the predator's natural mortality. All parameters are positive, and $b$, $e$ and $d$ are less than one. In addition, it requires that $eb-dh > 0$, i.e., the natural enemy species can survive when pests are abundant.

To prevent the rapid spread of pests, two control methods are adapted: one is the continuous control with two thresholds, that is, when the pest density is below the pest warning threshold ${{x}_{ET}}$, no control measures need to be taken, when the pest and natural enemy densities satisfy $x > {{x}_{ET}}$ and $y < {{y}_{ET}}$, the control action by spraying pesticides and releasing a part (${{q}_{1}}$) of the natural enemies is taken, which causes the death of pests (${{p}_{1}}$) and natural enemies (${{q}_{2}}$), when their densities satisfy $x > {{x}_{ET}}, y > {{y}_{ET}}$, only spraying pesticides is adapted. Based on the above control strategy, the impulsive control system can be formulated as follows:

where

Another one is an intermittent control with an economic threshold, that is, when the pest's density is below an economic threshold, no control action is implemented. Once the pest's density reaches the economic threshold, the control action by spraying pesticides and releasing a nonlinear volume $\frac{\tau }{1+ly}$ of natural enemies is taken, which causes the death of pests (${{p}_{1}}$) and natural enemies (${{q}_{2}}$), where $\tau$ and $l > 0$ are the formal parameters of the maximum volume of predators, respectively. Based on this control strategy, we can formulate the impulsive control system as follows:

The aim of this study focuses on analyzing the effects of different control measures on the dynamics of Models (2.2) and (2.4), respectively.

2.2.   Basic knowledge

2.2.1.   Filippov system

Consider a piecewise-continuous system

where

and discontinuous demarcation is

Let ${\bf F}_{i} H = \langle\nabla H, {\bf F}_{i}\rangle$, where $\langle\cdot, \cdot\rangle$ is the standard scalar product. Then ${\bf F}_{i}^{m} H = \langle\nabla\left({\bf F}_{i}^{m-1} H\right), {\bf F}_{i}\rangle$. Thus the discontinuous demarcation $\Sigma$ can be distinguished into three regions: 1) sliding region: $\Sigma_{s} = \left\{(x, y) \in \Sigma: {\bf F}_{1} H < 0\right. \mbox{and} \left.{\bf F}_{2} H > 0\right\}$; 2) crossing region: $\Sigma_{c} = \left\{(x, y) \in \Sigma: {\bf F}_{1} H \cdot {\bf F}_{2} H > 0\right\}$; 3) escaping region: $\Sigma_{e} = \left\{(x, y) \in \Sigma: {\bf F}_{1} H > 0\right. \mbox{and} \left.{\bf F}_{2} H < 0\right\}$.

The dynamics of system (2.5) along $\Sigma_{s}$ is determined by

where ${\bf F}_{s} = \lambda {\bf F}_{1}+(1-\lambda) {\bf F}_{2}$ with $\lambda = \frac{{\bf F}_{2} H}{{\bf F}_{2} H-{\bf F}_{1} H} \in(0.1)$.

Definition 1 (^[24]). For system (2.5), $E^{*}$ is a real equilibrium if $\exists i \in\{1, 2\}$ so that ${\bf F}_{i}\left(E^{*}\right) = 0$, $E^{*} \in S_{i}$; $E^{*}$ is a virtual equilibrium if $\exists i, j \in\{1, 2\}, i \neq j$, so that ${\bf F}_{i}\left(E^{*}\right) = 0$, $E^{*} \in S_{j}$; and $E^{*}$ is a pseudo-equilibrium if ${\bf F}_{s}\left(E^{*}\right) = \lambda {\bf F}_{1}\left(E^{*}\right)+(1-\lambda) {\bf F}_{2}\left(E^{*}\right) = 0, H\left(E^{*}\right) = 0$, and $\lambda = \frac{{\bf F}_{2} H}{{\bf F}_{2} H-{\bf F}_{1} H} \in(0, 1)$.

2.3.   Impulsive semi-continuous system

For the given planar model

we have:

Definition 2 (Order-$k$ periodic solution ^[50,51]). The solution $\tilde{{\bf z}}(t) = (\tilde{x}(t), \tilde{y}(t))$ is called periodic if there exists $n (\geqslant 1)$ satisfying $\tilde{{\bf z}}_{n} = \tilde{{\bf z}}_{0}$. Furthermore, $\tilde{{\bf z}}$ is an order-$k$ T-periodic solution with $k\triangleq\min\{j|1\leq j\leq n, \tilde{{\bf z}}_{j} = \tilde{{\bf z}}_{0}\}$.

Lemma 1 (Stability criterion ^[50,51]). The order-$k$ $T$-periodic solution ${\bf z}(t) = (\xi(t), \eta(t))^T$ is orbitally asymptotically stable if $\vert \mu _q \vert < 1$, where

with

$\chi_1^+ = \chi_1(\xi(\theta_j^ +), \eta(\theta_j^ +))$, $\chi_2^+ = \chi_2(\xi(\theta_j^+), \eta(\theta_j^+))$, and $\chi_1$, $\chi_2$, $\frac{\partial {I_1} }{\partial x}$, $\frac{\partial {I_1} }{\partial y}$, $\frac{\partial {I_2} }{\partial x}$, $\frac{\partial {I_2} }{\partial y}$, $\frac{\partial \omega}{\partial x}$, $\frac{\partial \omega}{\partial y}$ are calculated at $(\xi(\theta_j), \eta(\theta_j))$.

3.   Dynamic properties of the Filippov Model (2.1)

For convenience, denote

Since

then all solutions $(x(t), y(t))$ of Model (2.1) with $x(0) > 0$ and $y(0) > 0$ are positive in the region $\mathbb{D} = \{(x(t), y(t))|0 < x\le K, y\ge 0\}$.

Theorem 1. For Model (2.1), the solutions are ultimately bounded and uniform in the region $\mathbb{D}_1$.

Proof. Define $\iota (x(t), y(t))\triangleq x(t)+y(t)$. Then

Take $0 < \theta \le \min \left\{r, d \right\}$, and there is

Obviously, ${\sigma}'(x) = r(\ln K-r\ln x-1)-\theta$. If $0 < x < K{{\text{e}}^{\frac{\theta }{r}-1}}$, then ${\sigma }'(x) > 0$. If $x > K{{\text{e}}^{\frac{\theta }{r}-1}}$, then ${\sigma }'(x) < 0$. Then $\sigma (x)$ has a maximum ${{\sigma }^{*}}$. Thus $\frac{\text{d}}{\text{d}t}(\iota -\frac{{{\sigma }^{*}}}{\theta })\le -\theta (\iota -\frac{{{\sigma }^{*}}}{\theta })$, and then

For $t\to \infty$, there is $0\le \iota \left(x(t), y(t) \right)\le \frac{{{\sigma }^{*}}}{\theta }$. Therefore, the solutions of Model (2.1) are uniformly bounded in the region

For Model (2.1), the boundary equilibrium ${{E}_{K}}(K, 0)$ always exists. Define

Theorem 2. For Model (2.1), if $b < \bar{b}(d; 0)$, then ${{E}_{B}}(K, 0)$ is locally asymptotically stable. If $b > \bar{b}(d; 0)$, there exists a coexistence equilibrium, denoted as $E_{^{1}}^{*} = \left(x_{^{1}}^{*}, y_{^{1}}^{*} \right)$, which is locally asymptotically stable if $U(x_{1}^{*}) > 0$, where

Proof. For Model (2.1), we have

1) For ${{E}_{K}}(K, 0)$, we have

Then ${{\lambda }_{1}} = -r < 0$ and ${{\lambda }_{2}} = \frac{eb{{K}^{2}}}{(K+g)(1+hx)}-d$. Therefore, ${{E}_{B}}(K, 0)$ is locally asymptotically stable if $b < b(d; 0)$.

2) Since

then for $E_{1}^{*}$, we have

If $U(x_{1}^{*}) > 0$ holds, then ${{\lambda }_{1}}{{\lambda }_{2}} > 0, {{\lambda }_{1}}+{{\lambda }_{2}} < 0$, i.e., $E_{1}^{*}$ is locally asymptotically stable.

3.1.   Complex dynamics of Model (2.2)

Let

Then systems (2.2) and (2.3) can be described as

where

The switching boundaries are, respectively,

Let ${\bf {{n}_{1}}} = (1, 0)$ and ${\bf {{n}_{2}}} = (0, 1)$ be the normal vector for ${{\Sigma }_{1}}$ and ${{\Sigma }_{3}}$. If $\exists {{\Sigma }_{ij}}\subset {{\Sigma }_{i}}$ such that the trajectory of ${\bf {F}_{i}}(x, y)$ approaches or moves away from ${{\Sigma}_{i}}$ ($i\in \{1, 2, 3\}$) on both sides, then a sliding domain exists, and the dynamics on ${{\Sigma }_{i}}$ can be determined by means of the Filippov convex method.

3.1.1.   Dynamic behavior of the sub-models

The dynamic behavior of the model in ${{G}_{1}}$ can be referred to Section 3.2. The model in ${{G}_{2}}$ is described as follows:

Theorem 3. Model (3.2) always has an equilibrium ${E}_{{\overline{B}}}(K{{\mathit{\text{e}}}^{-\frac{{{p}_{1}}}{r}}}, 0)$. If ${{q}_{1}} < {{q}_{2}}+d$ and $b < b(d-q_1+q_2;p_1)$, then ${E}_{\overline{B}}(K{{\mathit{\text{e}}}^{-\frac{{{p}_{1}}}{r}}}, 0)$ is locally asymptotically stable. If ${{q}_{1}} < {{q}_{2}}+d$ and $b > \bar{b}(d+q_2-q_1, p_1)$, Model (3.2) has a coexistence equilibrium, denoted as $E_{^{2}}^{*} = \left(x_{^{2}}^{*}, y_{^{2}}^{*} \right)$, which is locally asymptotically stable if $U(x_{2}^{*}) > 0$, where

Similarly, the model in ${{G}_{3}}$ is described as follows:

Theorem 4. Model (3.3) always has an equilibrium ${{E}_{{\overline{B}}}}(K{{\mathit{\text{e}}}^{-\frac{{{p}_{1}}}{r}}}, 0)$. If $b < b(d+q_2;p_1)$, then ${{E}_{{\overline{B}}}}$ is locally asymptotically stable. If $b > b(d+q_2;p_1)$, Model (3.3) has a coexistence equilibrium, denoted as $E_{^{3}}^{*} = (x_{^{3}}^{*}, y_{^{3}}^{*})$, which is locally asymptotically stable if $U(x_{3}^{*}) > 0$, where

3.1.2.   Dynamic behavior of Model (2.2)

It is assumed that

(H1) ${{p}_{1}} < r$;

(H2) $\frac{d(1+gh)+\sqrt{\Delta(d)}}{2(eb-dh)} < K$;

(H3) $q_1-q_2-d < 0, \frac{(d-q_1+q_2)(1+gh)+\sqrt{\Delta(d-q_1+q_2)}}{2[eb+h(q_1-q_2-d)]} < K{{\text{e}}^{-\frac{{{p}_{1}}}{r}}}$;

(H4) $\frac{(d+{{q}_{2}})(1+gh)+\sqrt{\Delta(d+q_2)}}{2[eb-h(d+{{q}_{2}})]} < K{{\text{e}}^{-\frac{{{p}_{1}}}{r}}}$.

For Model (2.2), we have $x_1^* < x_2^* < x_3^*$ when $q_1 < q_2$ and $x_2^* < x_1^* < x_3^*$ when $q_2 < q_1 < q_2+d$.

Define

where ${{y}_{ET2}} > 0$ and ${{y}_{ET1}} < {{y}_{ET2}}$.

First, we will discuss the sliding mode domain on ${{\Sigma }_{1}}$ and the corresponding dynamics. Since

then the sliding mode domain on ${{\Sigma }_{1}}$ does not exist if ${{y}_{ET}} < {{y}_{ET1}}$. When ${{y}_{ET}} > {{y}_{ET1}}$, we have

Next, the Filippov convex method is used, i.e.,

where

and the sliding mode dynamics of Eq (3.1) along ${{\Sigma }_{11}}$ is determined by the following system:

Let ${{\varsigma }_{1}} = eb{{x}_{ET}}^{2}+(1+h{{x}_{ET}})({{x}_{ET}}+g)\left[\frac{r(q_1-q_2)}{{{p}_{1}}}(\ln K-\ln {{x}_{ET}})-d \right]$. Then a positive equilibrium ${{E}_{a1}}({{x}_{ET}}, {{y}_{a1}})$ exists, where ${{y}_{a1}} = \frac{{{p}_{1}}{{\varsigma }_{1}}}{b{{x}_{ET}}(q_1-q_2)} > 0$. Therefore

If $x_{1}^{*} < {{x}_{ET}}$, then ${{y}_{a1}} > {{y}_{ET2}}$, i.e., ${{E}_{a1}}$ is not located in ${{\Sigma }_{11}}$, and then ${{E}_{a1}}$ is not a pseudo-equilibrium. If $x_{1}^{*} > {{x}_{ET}}$, then ${{y}_{a1}} < {{y}_{ET2}}$.

Similarly, we have

If $x_{2}^{*} > {{x}_{ET}}$, then ${{y}_{a1}} < {{y}_{ET1}}$, i.e., ${{E}_{a1}}$ is not located in ${{\Sigma }_{11}}$, and then ${{E}_{a1}}$ is not a pseudo-equilibrium. If $x_{2}^{*} < {{x}_{ET}}$, then ${{y}_{a1}} > {{y}_{ET1}}$. Therefore, ${{y}_{ET1}} < {{y}_{a1}} < {{y}_{ET2}}$. When ${{y}_{a1}} < {{y}_{ET}}$, ${{E}_{a1}}$ is the pseudo-equilibrium.

Second, we will discuss the sliding mode domain on ${{\Sigma }_{2}}$ and the dynamic characteristics on the sliding mode. Since

then the sliding mode domain on ${{\Sigma }_{2}}$ does not exist if ${{y}_{ET}} > {{y}_{ET2}}$; When ${{y}_{ET}} < {{y}_{ET2}}$, we have

Therefore, when ${{y}_{ET}} > {{y}_{ET2}}$, there is no sliding mode domain on ${{\Sigma }_{2}}$. When ${{y}_{ET}} < {{y}_{ET2}}$, the sliding mode domain of the system (3.1) on ${{\Sigma }_{2}}$ can be expressed as Eq (3.9).

According to the Filippov convex method, we have

where

and the sliding mode dynamics of equation (3.1) along ${{\Sigma }_{22}}$ is determined by the following system:

Let ${{\varsigma }_{2}} = (1+h{{x}_{ET}})({{x}_{ET}}+g)\left[\frac{r{{q}_{2}}}{{{p}_{1}}}(\ln K-\ln {{x}_{ET}})+d \right]-eb{{x}_{ET}}^{2}$. Then the system (3.11) has a positive equilibrium ${{E}_{a2}}({{x}_{ET}}, {{y}_{a2}})$, where ${{y}_{a2}} = \frac{{{p}_{1}}{{\varsigma }_{2}}}{{{q}_{2}}b{{x}_{ET}}} > 0$. Obviously,

If $x_{1}^{*} > {{x}_{ET}}$, then ${{y}_{a2}} > {{y}_{ET2}}$, i.e., ${{E}_{a2}}$ is not located in ${{\Sigma }_{22}}$. If $x_{1}^{*} < {{x}_{ET}}$, then ${{y}_{a2}} < {{y}_{ET2}}$. Similarly, we have

If $x_{3}^{*} < {{x}_{ET}}$, then ${{y}_{a2}} < {{y}_{ET1}}$, i.e, ${{E}_{a2}}$ is not located in ${{\Sigma }_{22}}$. If $x_{3}^{*} > {{x}_{ET}}$, then ${{y}_{a2}} > {{y}_{ET1}}$. Therefore ${{y}_{ET1}} < {{y}_{a2}} < {{y}_{ET2}}$. If ${{y}_{ET1}} < {{y}_{ET}} < {{y}_{a2}}$ or ${{y}_{ET}} < {{y}_{ET1}}$, then ${{E}_{a2}}$ is located in ${{\Sigma }_{22}}$ and is a pseudo-equilibrium.

Finally, we will discuss the sliding mode domain on ${{\Sigma }_{3}}$ and the dynamic characteristics of the sliding mode. We have

According to Eq (3.12), ${{\left. < {\bf {F}_{2}}, {\bf {{n}_{2}}} > \right|}_{(x, y)\in {{\Sigma }_{3}}}} > {{\left. < {\bf {F}_{3}}, {\bf {{n}_{2}}} > \right|}_{(x, y)\in {{\Sigma }_{3}}}}$. If $x_{3}^{*} < {{x}_{ET}}$, then the system (3.1) does not have a sliding mode domain on ${{\Sigma }_{3}}$. If $x_{2}^{*} < {{x}_{ET}} < x_{3}^{*}$, it is found through Eq (3.12) that the system (3.1) can be expressed in the sliding mode domain on ${{\Sigma }_{3}}$ as

According to the Filippov convex method, we have

where

The sliding mode dynamics of Eq (3.1) along ${{\Sigma }_{11}}$ is determined by the following system:

Then

If the root $x = {{x}_{b}} > 0$ of Eq (3.16) satisfies Eq (3.13), then ${{E}_{b}}({{x}_{b}}, {{y}_{ET}})$ of the system (3.15) is a pseudo-equilibrium, and if it does not satisfy Eq (3.13), then ${{E}_{b}}$ is not a pseudo-equilibrium.

4.   Dynamic analysis of the impulse pest-management model

For Model (2.4), let

The curve $y = \widehat{y}(x)$ intersects with $x = {{x}_{ET}}$ and $x = (1-p_1){x}_{ET}$ at $P({{x}_{ET}}, {{y}_{P}})$ $({{y}_{P}} = \widehat{y}({{x}_{ET}}))$ and ${{R}_{0}}((1-{{p}_{1}}){{x}_{ET}}, {{y}_{{{R}_{0}}}})$. The trajectory passing through $P$ is denoted by ${{\gamma }_{1}}$, and it goes backward and intersects $y = \widehat{y}(x)$ at $H({{x}_{H}}, {{y}_{H}})({{y}_{H}} = \widehat{y}({{x}_{H}}))$. If ${{x}_{H}} < (1-{{p}_{1}}){{x}_{ET}}$, then denote ${{Q}_{1}}((1-{{p}_{1}}){{x}_{ET}}, {{y}_{{{Q}_{1}}}}), {{Q}_{2}}((1-{{p}_{1}}){{x}_{ET}}, {{y}_{{{Q}_{2}}}})$ as the intersection points between ${{\gamma }_{1}}$ and $x = (1-{{p}_{1}}){{x}_{ET}}$ with ${{y}_{{{Q}_{1}}}} < {{y}_{{{Q}_{2}}}}$. The trajectory passing through ${{R}_{0}}$ is denoted by ${{\gamma }_{2}}$. If ${{\gamma }_{2}}\cap \{x = x_{ET}\}\ne \varnothing$, then denote ${{R}_{1}}({{x}_{ET}}, {{y}_{{{R}_{1}}}})$ as the intersection point between ${{\gamma }_{2}}$ and $x = {{x}_{ET}}$. The curve ${{\gamma }_{2}}$ defines a function $y = y(x, {{y}_{{{R}_{0}}}})$ on the interval $[(1-p_{1}){{x}_{ET}}, {{x}_{ET}}]$ with

which takes the form

4.1.   Accurate domains of $\mathcal{M}$ and $\mathcal{N}$

For Model (2.4), we have $\mathcal{M} = \{(x, y)\mid x = {{x}_{ET}}, y > 0\}$. The trajectory of the system (2.4) with ${{x}_{0}} < {{x}_{ET}}$ can reach ${\mathcal{M}_{1}} = \{(x, y)| x = {{x}_{ET}}, 0\le y\le {{y}_{P}}\}\subset {\mathcal{M}}$, which is called the effective impulse set, denoted by $\mathcal{M}_{eff}$. The corresponding effective phase set is denoted by $\mathcal{N}_{eff}$. Moreover, define ${\mathcal{M}_{2}} = \left\{ (x, y)\mid x = {{x}_{ET}}, 0\le y\le {{y}_{{{R}_{1}}}} \right\}\subset {\mathcal{M}_{1}}$.

Since $\Delta y = -{{q}_{2}}y+\frac{\tau }{1+ly}$, then define

Obviously, the function $\rho(y)$ reaches a minimum at $y = \overset{ \frown}{y}$, where $\overset{ \frown}{y}\triangleq \frac{\sqrt{\tau l(1-{{q}_{2}})}-(1-{{q}_{2}})}{l(1-{{q}_{2}})}$. Denote $R({{x}_{ET}}, \overset{ \frown}{y})\in \mathcal{M}$, and its phase point is ${{R}^{+}}\left((1-{{p}_{1}}){{x}_{ET}}, \rho (\overset{ \frown}{y}) \right)$.

Define

Denote

The exact domains of $\mathcal{M}$ and $\mathcal{N}$ can be determined by ${\rm{sign}}({\rho}'(y))$ and ${\rm{sign}}(\overset{ \frown}{y})$, which will be discussed in the following two situations:

Case Ⅰ: $x_{ET}^{1} < x_{ET}\le x_{ET}^{2}$.

For this situation, ${\mathcal{M}_{eff}} = {\mathcal{M}_{1}}$. To determine $\mathcal{N}_{eff}$, we are required to judge the magnitude between $\overset{ \frown}{y}$ and ${{y}_{P}}$. Denote $\Lambda = [0, {{y}_{{{Q}_{1}}}}]\bigcup [{{y}_{{{Q}_{2}}}}, +\infty)$.

ⅰ) $\tau \ge \tau_{1}$, then $\overset{ \frown}{y}\le 0$. For $\forall y\in [0, {{y}_{P}}]$, ${\rho }'\ge 0$ holds, and then $\tau \le \rho (y)\le \rho ({{y}_{P}})$ for $y\in[0, {{y}_{P}}]$. Denote ${{\Lambda }_{11}} = [\tau, \rho ({{y}_{P}})]$, ${{\Lambda }_{1}} = \Lambda \bigcap {{\Lambda }_{11}}$, and ${\mathcal{N}_{eff}} = {\mathcal{N}_{1}} = \left\{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in {{\Lambda }_{1}} \right\}$.

ⅱ) $\tau_{1} < \tau < \tau_{2}$, then $0 < \overset{ \frown}{y} < {{y}_{P}}$. For $\forall y\in [0, \overset{ \frown}{y}]$, ${\rho }'\le 0$ holds, and then $\rho (\overset{ \frown}{y})\le \rho (y)\le \tau$ for $y\in [0, \overset{ \frown}{y}]$. Denote ${{\Lambda }_{21}} = [\rho (\overset{ \frown}{y}), \tau]$, $\Lambda _{21}^{*} = \Lambda \bigcap {{\Lambda }_{21}}$, and ${\mathcal{N}_{21}} = \left\{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in \Lambda _{21}^{*} \right\}$. Similarly, for $\forall y\in (\overset{ \frown}{y}, {{y}_{P}}]$, ${\rho }' > 0$ holds, i.e., $\rho (\overset{ \frown}{y}) < \rho (y)\le \rho ({{y}_{P}})$. Denote ${{\Lambda }_{22}} = (\rho (\overset{ \frown}{y}), \rho ({{y}_{P}})]$, $\Lambda _{22}^{*} = \Lambda \bigcap {{\Lambda }_{22}}$, and ${\mathcal{N}_{22}} = \left\{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in \Lambda _{22}^{*} \right\}$. Thus, we have ${\mathcal{N}_{eff}} = {\mathcal{N}_{2}} = {\mathcal{N}_{21}}\bigcup \mathcal{N}_{22}$.

ⅲ) $\tau \le \tau_{2}$, then $\overset{ \frown}{y}\ge {{y}_{P}}$. For $\forall y\in [0, {{y}_{P}}]$, ${\rho }'\le 0$ holds, i.e., $\rho ({{y}_{P}})\le \rho (y)\le \tau$. Denote ${{\Lambda }_{33}} = [\rho ({{y}_{P}}), \tau]$ and ${{\Lambda }_{3}} = \Lambda \bigcap {{\Lambda }_{33}}$. Then ${\mathcal{N}_{eff}} = {\mathcal{N}_{3}} = \{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in {{\Lambda }_{3}}\}$.

Case Ⅱ: $x_{ET}\ge x_{ET}^{1}$.

For this situation, ${\mathcal{M}_{eff}} = {\mathcal{M}_{2}}$. Similar to the discussion in case Ⅰ, we have

ⅰ) $\tau \ge \tau_{1}$. Then ${\mathcal{N}_{eff}} = {\mathcal{N}_{4}} = \left\{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in {{\Lambda }_{4}} \right\}$, where ${{\Lambda }_{4}} = [\tau, \rho ({{y}_{{{R}_{1}}}})]$.

ⅱ) $\tau_{1} < \tau < \tau_{3}$. For $\forall y\in [0, \overset{ \frown}{y}]$, we have ${\mathcal{N}_{51}} = \left\{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in {{\Lambda }_{21}} \right\}$. Similarly, for $\forall y\in (\overset{ \frown}{y}, {{y}_{{{R}_{1}}}}]$, denote ${{\Lambda }_{52}} = (\rho (\overset{ \frown}{y}), \rho ({{y}_{{{R}_{1}}}})]$, and then ${\mathcal{N}_{52}} = \{({{x}^{+}}, {{y}^{+}})|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in {{\Lambda }_{52}}\}$. Therefore, $\mathcal{N}_{eff} = {\mathcal{N}_{5}} = {\mathcal{N}_{51}}\bigcup {\mathcal{N}_{52}}$.

ⅲ) $\tau \le \tau_{3}$. Then $\mathcal{N}_{eff} = {\mathcal{N}_{6}} = \left\{ \left. ({{x}^{+}}, {{y}^{+}}) \right|{{x}^{+}} = (1-{{p}_{1}}){{x}_{ET}}, {{y}^{+}}\in {{\Lambda }_{6}} \right\}$ with ${{\Lambda }_{6}} = [\rho ({{y}_{{{R}_{1}}}}), \tau]$.

4.2.   Poincaré map

Denote ${{G}_{i}}({{x}_{ET}}, {{y}_{i}})\in \mathcal{M}$, $G_{i}^{+}((1-{{p}_{1}}){{x}_{ET}}, y_{i}^{+})\in \mathcal{N}$, $i = 0, 1, 2, ...$, where $G_{i}^{+} = \textbf{I}({{G}_{i}})$. Since $G_{i}^{+}$ and ${{G}_{i+1}}$ lie on the same trajectory ${{\gamma}_{G_{i}^{+}}}$, then we have ${{y}_{i+1}} = \varpi (y_{i}^{+})$ and $y_{i+1}^{+} = \psi (y_{i}^{+})$, where

If $\exists \hat{y}\in \mathcal{N}$ such that $\psi(\hat{y}) = \hat{y}$, then Model (2.4) admits an order-1 periodic trajectory. Next, we will investigate the monotonicity of $\psi (y)$ with $\tau > 0$ for situations Ⅰ and Ⅱ.

Case Ⅰ: $x_{ET}^{1} < x_{ET}\le x_{ET}^{2}$.

ⅰ) $\tau \ge \tau_{1}$. $\rho (y)$ monotonically increases on $[0, {{y}_{P}}]$, and then the map $\psi (y)$ monotonically increases on $[0, {{y}_{{{Q}_{1}}}}]$ and monotonically decreases on $[{{y}_{{{Q}_{2}}}}, +\infty)$.

ⅱ) $\tau_{1} < \tau < \tau_{2}$. we have $\rho'(y) < 0$ for $y\in [0, \overset{ \frown}{y}]$ and $\rho'(y) > 0$ for $y\in (\overset{ \frown}{y}, {{y}_{P}}]$. Denote ${{y}_{{{S}_{1}}}} = \min \{y:\psi (y) = {{y}_{{{R}^{+}}}}\}$, ${{y}_{{{S}_{2}}}} = \max \{y:\psi (y) = {{y}_{{{R}^{+}}}}\}$. Then $\psi (y)$ monotonically increases on $[{{y}_{{{S}_{1}}}}, {{y}_{{{Q}_{1}}}}]$, $[{{y}_{{{S}_{2}}}}, +\infty)$ and monotonically decreases on the interval $[0, {{y}_{{{S}_{1}}}}]$, $[{{y}_{{{Q}_{2}}}}, {{y}_{{{S}_{2}}}}]$, respectively.

ⅲ) $\tau \le \tau_{2}$. $\rho (y)$ monotonically decreases on $[0, {{y}_{P}}]$, and then the map $\psi (y)$ monotonically decreases on $[0, {{y}_{{{Q}_{1}}}}]$ and monotonically increases on $[{{y}_{{{Q}_{2}}}}, +\infty)$.

Case Ⅱ: $x_{ET}\ge x_{ET}^{1}$.

ⅰ) $\tau \ge \tau_{1}$. Then $\psi'(y) > 0$ for $y\in[0, {{y}_{{{R}_{0}}}}]$ and $\psi'(y) < 0$ for $y\in[{{y}_{{{R}_{0}}}}, +\infty)$.

ⅱ) $\tau_{1} < \tau < \tau_{3}$. Denote ${{y}_{{{V}_{1}}}} = \min \{y:\psi (y) = {{y}_{{{R}^{+}}}}\}$ and ${{y}_{{{V}_{2}}}} = \max \{y:\psi (y) = {{y}_{{{R}^{+}}}}\}$. Then $\psi (y)$ monotonically decreases on $[0, {{y}_{{{V}_{1}}}}]$ and $[{{y}_{{{R}_{0}}}}, {{y}_{{{V}_{2}}}}]$, and monotonically increases on $[{y}_{{V}_{1}}, {y}_{{R}_{0}}]$ and $[{y}_{{V}_{2}}, +\infty)$.

ⅲ) $\tau \le \tau_{3}$. The map $\psi (y)$ monotonically decreases on $[0, {{y}_{{{R}_{0}}}}]$ and monotonically increases on $[{{y}_{{{R}_{0}}}}, +\infty)$.

4.3.   Order-1 periodic trajectory

4.3.1.   Natural enemy extinction periodic trajectory

For Model (2.4) with $\tau = 0$, if ${{y}_{0}}\equiv 0$, then $y\equiv 0$ holds. Thus Model (2.4) is degenerated to

Let $x = \widetilde{x}(t)$ be the solution of equation

with initial value $\widetilde{x}(0) = {{x}_{0}}\triangleq (1-p_{1}){{x}_{ET}}$. Define

We have $\tilde{x}({{T}_{0}}) = {{x}_{ET}}$ and $\tilde{x}(T_{0}^{+}) = {{x}_{0}}$. Thus, $\textbf{z}(t) = (\widetilde{x}(t), 0)$ ($(k-1){{T}_{0}} < t\le k{{T}_{0}}$, $k\in {{N}_{+}}$) is a natural enemy extinction periodic trajectory.

Theorem 5. The natural enemy extinction period trajectory $\textbf{z}(t) = (\widetilde{x}(t), 0)$ ($(k-1){{T}_{0}} < t\le k{{T}_{0}}$, $k\in {{N}_{+}}$) is orbitally asymptotically stable if $q_{2} > \hat{q}$, where

Proof. For Model (4.1), we have

Then

Through calculation, we have

and

Thus,

Therefore, if $q_{2} > \hat{q}$, we have $\hat{\rho } < 1$, and by Lemma 1, $\textbf{z}(t) = (\widetilde{x}(t), 0)$ ($(k-1){{T}_{0}} < t\le k{{T}_{0}}$, $k\in {{N}_{+}}$) is orbitally asymptotically stable.

4.3.2.   Coexisting order-1 periodic trajectory

Denote that the points $P$, ${{R}_{1}}$ are mapped to the points ${{P}^{+}}\left((1-{{p}_{1}}){{x}_{ET}}, (1-{{q}_{2}}){{y}_{P}}+\frac{\tau }{1+l{{y}_{P}}} \right)$ and $R_{1}^{+}\left((1-{{p}_{1}}){{x}_{ET}}, (1-{{q}_{2}}){{y}_{{{R}_{1}}}}+\frac{\tau }{1+l{{y}_{{{R}_{1}}}}} \right)$, respectively, after a single impulse. Denote $W((1-{{p}_{1}}){{x}_{ET}}, \tau)$.

Case Ⅰ: $x_{ET}^{1} < x_{ET}\le x_{ET}^{2}$.

Define

Obviously, $\tau_4 < \tau_2$ and for $q_2\leq q_2^*$, we have $l > \max\{l^*, 0\}$.

1) For $\tau = \tilde{\tau}_2$, we have $\psi (y_{Q_2}) = y_{P^+} = y_{Q_2}$.

2) For $\tau > \tilde{\tau}_2$, we have $\psi (y_{Q_2}) = y_{P^+} > y_{Q_2}$. Then

● 2-a) for $\tau \ge \tau_2$, $W$is the highest after the pulse, while $P^+$ is the lowest after the pulse. Then $\psi (\tau) < \tau$, $\psi (y_{P^+}) > y_{P^+}$, and thus $\exists y'\in ({{y}_{P^{+}}}, \tau)$ such that $\psi(y') = y'$.

● 2-b) for $\tilde{\tau}_2 < \tau < \tau_2$, if $\tau\geq \hat{\tau}_2$, then $\rho(\overset{ \frown}{y})\ge y_{Q_2}$. Since the point $R^+$ is the lowest point after the pulse, then $\psi(y_{R^{+}}) > y_{R^{+}}$. If $\tau > \tau_4$, we have $\tau > y_{P^{+}}$. Then $W$ is the highest point after the impulse, i.e., $\psi (\tau) < \tau$. If $\tau\leq \tau_4$, we have $\tau\leq y_{P^{+}}$. Then $P^+$ is the highest point after the impulse, i.e., $\psi (y_{P^+}) < y_{P^+}$. Combine the above two aspects and it can be concluded that $\exists y''\in ({{y}_{R^{+}}}, \max\{\tau, y_{P^+}\})$ such that $\psi (y'') = y''$. While for $\tau < \hat{\tau}_2$, we have $\rho(\overset{ \frown}{y}) < y_{Q_2}$. In such a case, $\psi$ is not defined on $[\rho(\overset{ \frown}{y}), y_{Q_2})$ and it is uncertain whether a fixed point of $\psi$ exists or not.

3) When $\tilde{\tau}_1 < \tau < \tilde{\tau}_2$, then $y_{Q_1} < y_{P^+} < y_{Q_2}$. If $\tau\geq \hat{\tau}_1$, then $\rho(\overset{ \frown}{y})\ge y_{Q_1}$, i.e., $\psi$ does not have a fixed point. While for $\tau < \hat{\tau}_1$, we have $\rho(\overset{ \frown}{y}) < y_{Q_1}$. In such a case, it is uncertain whether a fixed point of $\psi$ exists or not.

4) When $0 < \tau < \tilde{\tau}_1$, and then $\psi (y_{Q_1}) = y_{P^+} < y_{Q_1}$. If $\tau > \tau_1$, the point $R^+$ is the lowest point after the pulse, then $\psi(y_{R^{+}}) > y_{R^{+}}$, i.e., $\psi (y)$ has a fixed point on $(y_{R^{+}}, {{y}_{{Q_1}}})$. If $\tau\leq \tau_1$, the point $W$ is the lowest point after the pulse, and then $\psi(\tau) > \tau$, i.e., $\psi (y)$ has a fixed point on $(\tau, {{y}_{{Q_1}}})$.

Case Ⅱ: $x_{ET}\ge x_{ET}^{1}$.

Define $\tilde{\tau}\triangleq (1+ly_{R_1})(y_{R_0}-(1-q_2)y_{R_1})$.

1) For $\tau = \tilde{\tau}$, we have $\psi (y_{R_0}) = y_{R_{1}^{+}} = y_{R_0}$.

2) For $\tau > \tilde{\tau}$, we have $\psi (y_{R_0}) = y_{R_{1}^{+}} > y_{{R}_{0}}$. Then

● 2-a) for $\tau \ge \tau_3$, since $R^+_1$ and $W$ are the lowest and highest points after the pulse, then we have $\psi ({{y}_{R_{1}^{+}}}) > {{y}_{R_{1}^{+}}}$, $\psi (\tau) < \tau$, and thus $\exists y'\in ({{y}_{R_{1}^{+}}}, \tau)$ such that $\psi (y') = y'$;

● 2-b) for $\tilde{\tau} < \tau < \tau_3$ and if $\tau > y_{R_{1}^{+}}$, then $\psi (\tau) < \tau$; if $\tau\leq y_{R_{1}^{+}}$, then $\psi ({{y}_{R_{1}^{+}}}) > {{y}_{R_{1}^{+}}}$. On the other hand, take the point $A$ in a small neighborhood near the point ${{R}_{0}}$, i.e., $A\in \bigcup ({{R}_{0}}, \delta)$. $A$ is above ${{R}_{0}}$. By the continuity of the impulse function and the Poincaré map, we have $\psi (y_A) > y_A$. Therefore, the map $\psi (y)$ has a fixed point on $(y_A, \max\{y_{R_1^+}, \tau\})$.

3) When $\tau < \tilde{\tau}$, then $\psi (y_{R_0}) = y_{R_{1}^{+}} < y_{{R}_{0}}$. If $\tau > \tau_1$, we have $\psi ({{y}_{{{R}^{+}}}}) > {{y}_{{{R}^{+}}}}$. If $\tau\leq \tau_1$, we have $\psi(\tau) > \tau$. Combine the above two aspects and it can be concluded that $\exists y''\in (\max\{y_{R^+}, \tau\}, y_{R_0})$ such that $\psi(y'') = y''$.

To sum up, we have:

Theorem 6. For the situation of $x_{ET}\ge x_{ET}^{1}$, Model (2.4) admits an order-1 periodic trajectory. While for the situation of $x_{ET}^{1} < x_{ET}\le x_{ET}^{2}$, Model (2.4) admits an order-1 periodic trajectory if $\tau < \tilde{\tau}_1$ or $\tau\geq\max\{\tilde{\tau}_2, \hat{\tau}_2\}$.

Let $\widetilde{{\bf z}}(t) = (\xi (t), {\eta }(t))$ ($(k-1){{T}} < t\le k{{T}}$, $k\in {{N}_{+}}$) be the $T$-periodic trajectory of the system (2.4) with initial values ${{A}_{0}}((1-p_{1}){{x}_{ET}}, {{y}_{{{A}_{0}}}})$. The trajectory intersects $\mathcal{M}$ at $A_{_{0}}^{-}(\xi (T), {\eta }(T))$, where $\xi (T) = {{x}_{ET}}, {\eta }(T) = {{y}_{0}}$, and then it is pulsed to $\mathcal{N}$ at $A_{_{0}}^{+}(\xi ({{T}^{+}}), {\eta }({{T}^{+}}))$. Thus, $\xi ({{T}^{+}}) = (1-{{p}_{1}}){{x}_{ET}}, {{\eta }_{1}}({{T}^{+}}) = (1-{{q}_{2}}){{y}_{0}}+\frac{\tau }{1+l{{y}_{0}}} = {{y}_{{{A}_{0}}}}.$

Theorem 7. The $T$-periodic trajectory $\widetilde{{\bf z}}(t) = (\xi (t), {\eta }(t))$ ($(k-1){{T}} < t\le k{{T}}$, $k\in {{N}_{+}}$) with initial values $\left((1-{{p}_{1}}){{x}_{ET}}, (1-{{q}_{2}}){{y}_{0}}+\frac{\tau }{1+l{{y}_{0}}}\right)$ is orbitally asymptotically stable if

where

Proof. The proof can be referred to that in Theorem 5 and is, therefore omitted.

5.   Computer simulations

For the purpose of simulation, it is assumed that $r = 1.5$, $K = 120$, $b = 0.595$, $h = 0.92$, $e = 0.8$, $d = 0.5$ and $g = 0.8$.

5.1.   Verification of Model (2.1)

When $b = 0.595$, the interior equilibrium $E_{1}^{*} = (54.7,103.1)$ is locally asymptotically stable, as presented in Figure 1. When $b$ increases to $0.61$, a limit cycle occurs, as presented in Figure 2. The effect of the maximum search rate on pests and natural enemies in the coexistence steady state is presented in Figure 3, and it is obvious that $x_{1}^{*}$ decreases with increasing $b$, while $y_{1}^{*}$ increases and then decreases with increasing $b$. Therefore, increasing the search rate for pests helps to reduce the number of pests.

5.2.   Verification of Model (2.2)

When ${{p}_{1}} = 0.25$, ${{q}_{1}} = 0.5$ and ${{q}_{2}} = 0.007$, the positive equilibrium of the ${{G}_{1}}$ region is $E_{1}^{*} = (54.707,103.1337)$, the positive equilibrium of the ${{G}_{2}}$ region is $E_{2}^{*} = (0.1229,141.532)$, and the positive equilibrium of the ${{G}_{3}}$ region is $E_{3}^{*} = (92.5247, 20.443)$. When ${{x}_{ET}} = 100$ and ${{y}_{ET}} = 110$, there is $x_{2}^{*} < x_{3}^{*} < {{x}_{ET}}$, $x_{1}^{*} < {{x}_{ET}}$, ${{y}_{ET1}} = 3.6997$ and ${{y}_{ET2}} = 43.0879$. The sliding mode domain of Model (3.1) on ${{\Sigma }_{1}}$ can be represented as ${{\Sigma }_{11}} = \{(x, y)\in {{\Sigma }_{1}} | 3.6997 < y < 43.0879 \}$, and there is no pseudo-equilibrium on ${{\Sigma }_{1}}$, as illustrated in Figure 4(a). When ${{x}_{ET}} = 80$ and ${{y}_{ET}} = 110$, there is $x_{2}^{*} < {{x}_{ET}} < x_{3}^{*}$, $x_{1}^{*} < {{x}_{ET}}$, ${{y}_{ET1}} = 45.3593$ and ${{y}_{ET2}} = 77.0172$. The sliding mode domain of Model (3.1) on ${{\Sigma }_{1}}$ can be represented as ${{\Sigma }_{11}} = \{ \left. (x, y)\in {{\Sigma }_{1}} \right|45.3593 < y < 77.0172\}$, and there is no pseudo-equilibrium on ${{\Sigma }_{1}}$, as presented in Figure 4(b).

5.3.   Verification of Model (2.4)

When ${{q}_{2}} = 0.11$, for ${{x}_{ET}} = 62.667$, the periodic trajectory of Model (2.4) is presented by changing the killing rate ${{p}_{1}}$ of the prey, the amount of predator released $\tau$, and the value of the parameter $l$. When $\tau = 0$, ${{p}_{1}} = 0.5$ and $l = 0.02$, the natural enemy extinction periodic trajectory is orbitally asymptotically stable (Figure 5(a)). To prevent the extinction of natural enemies, we are required to release natural enemies in an appropriate amount. When $\tau = 5$, the natural enemy extinction periodic trajectory loses its stability and an order-1 periodic trajectory occurs (Figure 5(b)).

Next, the accurate domains of $\mathcal{M}$ and $\mathcal{N}$ for different cases are presented as well as the order-1 periodic trajectory (Figure 6). The accurate domains of $\mathcal{M}$ and $\mathcal{N}$ are marked in red and blue solid lines, respectively. When $p_1 = 0.5$, $H$ lies on the left side of the phase set. The schematic diagram of the exact domain of the phase set and pulse set, and the order-1 periodic trajectories for different cases are presented in subfigures Figure 6(a)–(c). When $p_1 = 0.6$, $H$ lies on the right side of the phase set. The schematic diagram of the accurate domain of $\mathcal{M}$ and $\mathcal{N}$ and the order-1 periodic trajectories for different cases are presented in subfigures Figure 6(d)–(f).

When ${{p}_{1}} = 0.5$, $l = 0.02$, $\tau = 10$, $b = 0.595$ and ${{x}_{ET}} = 50$, $E^*$ is locally asymptotically stable, and Model (2.4) admits an order-1 periodic trajectory for $x_{ET} < x_{ET}^{1}$, as presented in Figure 7.

Finally, order-$n$ periodic solutions are presented for different $\tau$ and $l$. When $b = 0.595$, $E^*$ is locally asymptotically stable. For control parameters $l = 0.03$, $\tau = 120$, $x_{ET} = 62.667$ or $l = 0.002$, $\tau = 40$, $x_{ET} = 62.667$, Model (2.2) admits an order-$k$ periodic trajectory, as presented in subfigures 8(a) and (b). When $b = 0.61$, Model (2.1) admits a limit cycle. For ${{x}_{ET}} = 50$, Model (2.2) admits an order-$k$ periodic trajectory, as presented in subfigures 8(c)–8(f).

6.   Conclusions and discussion

Pests are important factors that harm agricultural production. In order to effectively control the spread of pests, a pest-natural enemy model with a variable search rate and threshold dependent feedback control was proposed. The dynamic properties such as the existence, positivity, and boundedness of solutions for continuous systems were discussed, and the results show that pests and natural enemies will not increase indefinitely due to system constraints (Theorem 1). In addition, it is shown that the natural enemy's searching rate $b$ plays an important role in determining the dynamics of the system, i.e., when $b$ is smaller than the level $\bar{b}(d; 0)$, the predators in the system will go to extinction and when $b$ is greater than $\bar{b}(d; 0)$, there exists a steady state $E^*$ at which the natural enemies and the pests in the system keep a balance. Moreover, the steady state is locally asymptotically stable as long as $U(x^*) > 0$ (Theorems 2–4, Figure 1). When $U(x^*) < 0$, the stability is lost and a limit cycle surrounding $E^*$ is obtained (Figure 2). The relationship between the number of pests (natural enemies) and the maximum search rate at the steady state was presented in Figure 3.

To prevent the spread of pests, two different types of control strategies were adopted. The first is a non-smooth control and the model is described by a Filippov system with two warning thresholds. By analyzing the sliding dynamics, we discussed the existence of pseudo-equilibrium $E_p$ (Figure 4). The pseudo-equilibrium $E_p$ is a new state of the control system at which the pests and the natural enemies keep a balance and the pest populations can be controlled at appropriate levels, which in turn indicates the effectiveness of the control. The second is an intermittent control with an economic threshold. When the pests reach the economic threshold, manual intervention is carried out by spraying pesticides and releasing a certain amount of natural enemies. For the control model, the accurate domain of the phase set was presented and the Poincaré map was constructed, through which the conditions for the existence of the order-1 periodic trajectories were presented (Theorems 5 and 6 and Figures 5–7). The order-1 periodic solution provides a possibility for periodic pest control, thus avoiding the need and difficulty of implementing pest population monitoring. The stability of the order-1 periodic trajectory was also verified (Theorems 5 and 7). This ensures the robustness of the control, and even if there is a condition monitoring error, it can still converge to the periodic solution of the system, thus providing a guarantee for the periodic control. We also presented the order-$k$ periodic solutions in numerical simulations (Figure 8), which further explain the complexity of the control system and the necessity of maintaining the stability of the system. The results illustrate the complex dynamics of the proposed models, which can serve as a valuable reference for the advancement of sustainable agricultural practices and the control of pests.

Use of AI tools declaration

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

Acknowledgement

The work was supported by the National Natural Science Foundation of China (No. 11401068).

Conflict of interest

The authors declare that there is no known competing financial interests to influence the work in this paper.