Effects of additional food availability and pulse control on the dynamics of a Holling-($p$+1) type pest-natural enemy model

1.   Introduction

Bemisia tabaci is widely distributed in more than 90 countries and regions. It was recorded in China as far back as the 1940s and has become a major pest of fruits, vegetables and garden flowers ^[1,2,3]. Biologists have made many attempts to eliminate the effect of Bemisia tabaci on plants, including the use of pesticides, but the negative effect of pesticide is also evident. As an alternative way, the idea of biological control was put forward, and the key was to find the right natural enemy species of Bemisia tabaci. In view of the predatory relationship between Neoseiseius barkeri and Bemisia tabaci, Neoseiseius barkeri was regarded as one of the best biological control species of Bemisia tabaci due to its short development period, low death rate and high spawning rate. Moreover, Neoseiulus barkeri is commonly found in plants such as papaya, rice, mango and artemisia annua, and it has a number of food sources, feeding on phylloides, tarsus mites and gall mites, thrips, aphids, moths and whiteflies ^[4,5,6,7]. Plant pollen and insect honeydew are also complementary foods of Neoseiulus barkeri when prey is scarce. Therefore, Neoseiulus barkeri has a high survival rate in the wild and can provide good and sustainable control of plant pests ^[8,9,10].

In natural ecosystems, predator-prey interactions are key determinants of species behavior, community species composition, and their dynamics. There exists a dynamic balance of predation and mutual selection between the two species during the long evolutionary process. In addition, the existence of predation relationships also affects and restricts the development of populations to a certain extent. Since mathematical models can provide an effective way for in-depth understanding of the mechanism, existence and stability of population growth, to reveal the effects of predation relationships between populations and their group effects on ecosystem stability, scholars have conducted extensive studies on the behavior of predator-prey models with different effects in recent years ^{[11,12,13,14,15]}. It should be pointed out that Lotka ^[16] and Volterra ^[17] initiated the earliest research work and built the classic Lotka-Volterra model from the perspective aspect of molecular chemistry and ocean ecological system respectively. For different application scenarios and real problems, in subsequent studies, many scholars extended the Lotka-Volterra model including proposals of logistic growth function ^[18], different functional responses ^[19,20,21] and so on ^{[22,23,24,25,26]}. The generalized Gause-type predator-prey model takes the following form

where $\varphi(x)$ characterizes the functional response of predator on prey. It was shown in experiments and analysis that Holling type functional response functions are basically applicable to different organisms, thus attracting the attention and research of many subsequent scholars, such as Holling type Ⅱ ^{[11,21,23,27]} and Holling type Ⅲ ^[28,29,30]. Yang et al. ^[31] introduced an extension form of the uptake function called Holling-($p$+1) ($p\geq 2$), and discussed the dynamics of the proposed model by using $p$ as an extended parameter. In this paper, Holling-type functional responses between Bemisia tabaci and Neoseiulus barkeri conforming to the extended form will also be considered.

From the perspective of species protection and pest management point of view, providing additional food for predators in the predatory system has been regarded in the literature as one of the effective biological ways ^[32,33]. To further investigate the effect of additional food, scholars have carried out related researches on this topic ^[34,35,36]. For species like Neoseiulus barkeri, plant pollen and insect honeydew are always their complementary foods, so it is vital and practical to study the Bemisia tabaci and Neoseiulus barkeri system with additional food provision. To suppress the proliferation and spread of Bemisia tabacti, effective control strategies must be adopted. The concept of integrated pest management (IPM) was widely used in the mid to late 20th century ^[37,38], and it emphasizes an integration of different methods (biological, chemical and others) to minimize the use of harmful pesticides and other undesirable measures to control pests. To describe the instantaneous changes in biological populations caused by IPM strategies, impulsive differential equations (IDEs) can be used as a powerful tool to model this phenomenon ^[39]. Tang and Chen ^[40] developed the Lotka-Volterra model by introducing periodic impulsive control. Liu et al. ^[41] analyzed a Holling type Ⅰ pest management model with periodic impulsive control. Zhang et al. ^[42] and Pei ^[43] discussed the mathematical models concerning continuous and impulsive pest control strategies, respectively. Song and Li ^[44] discussed a Holling type Ⅱ Leslie-Gower model with periodic impulsive effect. Wang et al. ^[45] proposed a Watt-type prey-predator model with periodic impulsive effect. Qian et al. ^[46] analyzed a prey-predator model with competition of predator and periodic impulsive control. Tang et al. ^[47,48] and Pei et al. ^[49] discussed the optimum timing control problem. Besides, scholars had also introduced IDEs to model different situations, such as pest control ^{[50,51,52,53]}, ecological system ^[54,55,56] and fishery exploitation ^{[57,58,59,60]}. In this work, we also applied an impulsive control strategy into the system to suppress the spread of Bemisia tabacti.

The article is organized in the following way. In Section 2, a Bemisia tabaci and Neoseiulus barkeri model with additional food provision is established, and the dynamical properties such as the existence and local stability of equilibria as well as the existence of the limit cycle are analyzed. In Section 3, the Bemisia tabaci and Neoseiulus barkeri model with periodic impulsive control is put forward, and it is shown that the Bemisia tabaci can be eradicated as long as the period of control is less than a certain threshold, while for a larger control period, persistence of the system can be achieved. Using the bifurcation theory, it is shown that a supercritical coexistence periodic solution can be bifurcated from the pest-extinction one when the control period exceeds a certain time threshold. In Section 4, numerical simulations are implemented with the purpose of verifying the main results. At the end of the paper, a brief conclusion with practical significance is presented.

2.   Pest management model

Bemisia tabaci is a major pest of fruits, vegetables and garden flowers. Neoseiulus barkeri is considered an appropriate natural enemy against Bemisia tabaci, which is polyphagous and can maintain survival on plant pollen in a short time whenever Bemisia tabaci is scarce ^[8]. Based on this phenomenon, it is hypothesized that Bemisia tabaci and Neoseiulus barkeri follow the Holling-(p+1) functional response associated with additional food supply, i.e.,

where

$\bullet$ $x$ represents the density of Bemisia tabaci,

$\bullet$ $y$ represents the density of Neoseiulus barkeri,

$\bullet$ $r$ represents the Bemisia tabaci's intrinsic growth rate,

$\bullet$ $K$ represents the Bemisia tabaci's environmental carrying capacity,

$\bullet$ $p$ denotes the Holling index of the functional response, $p\geq 2$,

$\bullet$ $\gamma$ denotes Neoseiulus barkeri's average capture rate of Bemisia tabaci,

$\bullet$ $a$ denotes the ratio of Neoseiulus barkeri's handling time per unit of extra food to that per unit of prey,

$\bullet$ $\eta$ denotes ratio of Neoseiulus barkeri's ability to find additional food to Neoseiulus barkeri's ability to find prey,

$\bullet$ $A$ denotes the amount of extra food,

$\bullet$ $d$ denotes the natural mortality rate of Neoseiulus barkeri,

$\bullet$ $\upsilon$ denotes the conversion coefficient.

Since the additional food is only a secondary option, it cannot sustain the long-term survival of Neoseiulus barkeri, and the following hypotheses are assumed:

A1) Without considering additional food sources, Neoseiulus barkeri will not become extinct due to the existence of Bemisia tabaci, i.e., $\upsilon \gamma -d > 0$ holds in (2.1).

A2) In absence of Bemisia tabaci, the additional food resource cannot maintain the survival of Bemisia tabaci, i.e., $\eta A \leq \overline{\eta A}\triangleq d q/(\upsilon \gamma- a d)$ holds.

For Bemisia tabacti, it is impossible to achieve the expected control effect by only relying on predator predation, so an additional control measure is required. Considering the seasonality and periodicity of the Bemisia tabaci occurrence, a periodic impulsive control strategy is introduced into the system, which can be formulated as follows:

The illustration of (2.2) is as follows: At fixed times $t = j \tau\left(j \in \mathbb{N}_{+}\right)$, integrated management measures (i.e., release a certain amount of Neoseiulus barkeri ($\delta$) and simultaneously kill a certain proportion of Bemisia tabaci ($\kappa_1$), also result in a certain proportion of Neoseiulus barkeri ($\kappa_2$) killed) are taken, which causes the densities of Bemisia tabaci $x$ and Neoseiulus barkeri $y$ to be immediately changed to $(1-\kappa_1)x$ and $(1-\kappa_2)y+\delta$.

3.   Main work

3.1.   Dynamic analysis of (2.1)

Define

For $\upsilon\gamma > d$, denote

3.1.1.   Equilibria and Local stability

Theorem 1. Model (2.1) always has a saddle point $E_0(0, 0)$ and a boundary equilibrium $E_K(K, 0)$. $E_{K}$ is globally asymptotically stable if $\eta A < \underline{\eta A}\triangleq (dq-(\upsilon \gamma -d) K^{p})/(\upsilon \gamma -a d)$. The coexistence equilibrium $E^{*} (x^*, y^*)$ exists in case of $\underline{\eta A} < \eta A < \overline{\eta A}$ and is locally asymptotically stable if one of the constraints holds: 1) $p\geq\overline{p}$; 2) $p < \overline{p}$ and $K < \overline{K}$, where

Proof. The Jacobi matrix is

1) For $E_{0}$, there is

By Assumption 2, there is ${\upsilon \gamma \eta A}/{(q+a \eta A)}-d < 0$, then $E_{0}$ is a saddle and unstable.

2) For $E_{K}$, there is

the eigenvalues are $\lambda_{1} = -r < 0$, $\lambda_{2} = (\upsilon \gamma \left(K^{p}+\eta A\right))/(q+K^{p}+a \eta A)-d$. Thus, $E_{K}$ is locally asymptotically stable if $\eta A < \underline{\eta A} = (dq-(\upsilon \gamma -d) K^{p})/(\upsilon \gamma -a d)$.

3) In case of $\underline{\eta A} < \eta A < \overline{\eta A}$, $E_K$ is unstable. Since

then for $E^{*}$, there is

Since $f'_{y}\left(x^{*}, y^{*}\right) < 0$ an $g'\left(x^{*}\right) > 0$, then $\lambda_1\lambda_2 > 0$. If one of the conditions holds: 1) $p\geq\overline{p}$; 2) $p < \overline{p}$ and $K < \overline{K}$ holds, then $f'_x(x^*, y^*) < 0$, that is $\lambda_1+\lambda_2 < 0$, thus, $E^{*}$ is locally asymptotically stable.

3.1.2.   Limit cycle

When the sign of inequality (3.1) is reversed, the coexistence equilibrium $E^*$ becomes unstable. In this case, we have:

Theorem 2. In case of $p < \overline{p}$ and $K > \overline{K}$, there exists at least a limit cycle surrounding $E^{*}$.

Proof. If $p < \overline{p}$ and $K > \overline{K}$ hold, then $E^{*}\left(x^{*}, y^{*}\right)$ exists but is unstable.

Next, we constructed a closed region $G$ containing $E^*$, such that all solutions of (2.1) are bounded in $G$.

1) Since

then the trajectory of (2.1) will pass the line $x = K$ from the right to the left.

2) For the line $l \triangleq \upsilon x+y-k = 0$, then $y = k-\upsilon x$. Since

then

Obviously, $h(x)$ is a continuous bounded function and has a maximum for $0 \leq x \leq K$. Denote $h_{\max} = \max\{h(x)\mid 0\leq x\leq K\}$, so it is only to set $k > h_{\max}/m$, then $\frac{{\rm d } l}{{\rm d } t}\mid_{l = 0} < 0$. So the trajectory of (2.1) will pass the line $l$ from the top to the bottom. In addition, to ensure $E^{*}$ lies in the region, there must be $k > k_{1}\triangleq y^{*}+\upsilon x^{*}$, so as to $k > \max \left\{{h_{\max}}/{m}, k_{1}\right\}$.

The straight lines $x = 0$ and $y = 0$ are both trajectories of (2.1), so $x = K$, $y = k-\upsilon x$, the $x$-axis and the $y$-axis encloses a bounded region $G$. The well known Poincaré-Bendixson Theorem implies that there exists at least a limit cycle surrounding $E^{*}$ in $G$.

Remark 1. Theorem 2 gives the existence of limit cycles, but its uniqueness cannot be determined. Therefore, the stability of the limit cycle can not be determined, which remains an open problem. If the limit cycle is unique, then it is stable from two sides.

3.2.   Complex dynamics of (2.2)

3.2.1.   Pest-eradication periodic solution

In absence of Bemisia tabaci, a reduced subsystem is obtained

The solution of (3.2) is

with

Theorem 3. For (2.2), there exists a pest-eradication periodic solution $\mathbf{\bar{z}}(t) = (0, \overline{y}(t))\; ((j-1)\tau\leq t < j\tau)$. Moreover, for any solution $y(t)$ of (3.2), there is $y(t) \to \overline{y}(t)\; (t \to \infty)$.

Proof. Clearly,

with $\mathbf{\bar{z}}(0^{+}) = (0, \overline{y}\left(0^{+}\right))$ is a pest-eradication periodic solution of (2.2).

For $0 < t \leq \tau$, there is

then $y\left(\tau^{+}\right) = \left(1-\kappa_2\right) y(\tau)+\delta$.

For $\tau < t \leq 2\tau$, there is

and $y\left(2 \tau^{+}\right) = \left(1-\kappa_2\right) y(2 \tau)+\delta$.

Similarly, for $(j-1)\tau < t \leq j\tau$, there is

and $y\left({(j-1)\tau}^{+}\right) = \left(1-\kappa_2\right) y(j\tau)+\delta$.

Then we have

Thus, for $j\tau < t\le (j+1)\tau$, there is

due to the assumption that $qd + ad\eta A - \upsilon \gamma \eta A > 0$.

Theorem 4. For (2.2), $\mathbf{\overline{z}}(t) = (0, \overline{y}(t))$ is locally asymptotically stable and globally attractive if $\tau < \tau_0\triangleq -\ln (1-\kappa_1)/r.$

Proof. Denote $\sigma_1(t) \triangleq x(t)$, $\sigma_2(t)\triangleq y(t)-\overline{y}(t)$. Then, $(\sigma_1, \sigma_2)$ satisfies

where $\Phi$ is the fundamental solution matrix, i.e.,

with $\Phi (0) = \textbf{I}$.

Since

then

Undoubtedly,

According to Floquet theory, if $\tau < \tau_0 = -\ln (1-\kappa_1)/r$, then $\mathbf{\overline{z}}(t) = (0, \overline{y}(t))$ is locally asymptotically stable.

Choose $\epsilon > 0$ such that

Denote $\rho \triangleq\left(1-\kappa_1\right) \exp \left\{\int_{0}^{\tau}\left(r\left(1-\frac{x}{K}\right)-\frac{\gamma x^{p-1}(\overline{y}(t)-\epsilon)}{q+x^{p}+a \eta A}\right) {\rm d } t\right\}$. Since $\dot{y} > -d y(t)$, then for system

there is $y(t)\geq v(t)\rightarrow \overline{y}(t)$ as $t \rightarrow \infty$. Thus,

holds when $t$ is sufficiently large. Here it is assumed that (3.3) always holds. Since

then $x((j+1)\tau) \leq \rho x(j \tau)$, which means that $x(j \tau) \leq x\left(0^{+}\right) \rho^{j}$ and $x(j\tau) \rightarrow 0$ as $n \rightarrow \infty$ due to $\rho < 1$. Since $0 < x(t) \leq x(j \tau)\left(1-\kappa_1\right) e^{r \tau}$ for $j \tau < t \leq(j+1)\tau$, then $x(t) \rightarrow 0$ as $n \rightarrow \infty$.

Next, for $0 < \varepsilon < [(d q+a d \eta A-\upsilon \gamma\eta A)/(\upsilon \gamma-d)]^{\frac{1}{p}}$, $\exists T_0 > 0$, $0 < x(t) < \varepsilon\; (t \geq T_0)$. Here, we suppose that $0 < x(t) < \varepsilon\; \; (t \geq 0)$. Since

then there is $\overline{y}(t)\leftarrow v_{1}(t) \leq y(t) \leq v_{2}(t)\rightarrow \overline{v}_{2}(t)\; (t\rightarrow \infty)$, where $v_{1}(t)$, $v_{2}(t)$ satisfies the following two equations, respectively

and

Since

then $\forall\varepsilon_{1} > 0$, $\exists T_1 > 0$, $\overline{y}(t)-\varepsilon_{1} < y(t) < \overline{v}_{2}(t)+\varepsilon_{1}\; (t > T_1)$, and for $\varepsilon \rightarrow 0$, there is $\overline{y}(t)-\varepsilon_{1} < y(t) < \overline{y}(t)+\varepsilon_{1}$ for $t > T_1$, which means that $y(t) \rightarrow \overline{y}(t)$ when $t\rightarrow +\infty$.

3.2.2.   Persistence analysis

Theorem 4 indicates that if the control period is less than $\tau_0$, then Bemisia tabaci will be completely eradicated from the system. However, from practical and economic aspects of point of view, it is impossible to take controls frequently. Therefore, it is necessary to consider that $\tau > \tau_0$.

Lemma 1. System (2.2) is bounded, i.e., $\exists M > 0$, there is $\max\{x(t), y(t)\}\leq M$ holds for sufficiently large $t$.

Proof. Denote $L(t) = y(t)+\upsilon x(t)$. Then for $t\neq j\tau$,

Thus

Define $m^{*}\triangleq\min \left\{d, \frac{a d-\upsilon \gamma }{a}\right\}$. Clearly, for $m = m^{*}/2$, there is

i.e., (3.4) is bounded. Then

For the system

there is

Therefore, $L(t)$ is bounded, which also means that $\exists M = \frac{M_{0}}{m^{*}/2}+\delta \frac{e^{m^{*}\tau/2}}{e^{m^{*}\tau/2}-1} > 0$, $\max\{x(t), y(t)\} < L(t)\leq M$.

For system (2.2), if $\exists M_1, M_2 > 0$ and $\overline{T} > 0$ such that $M_1\leq x_1(t) \leq M_2, M_1 \leq x_2(t) \leq M_2$ for all $t \geq \overline{T}$, then it is said to be persistent.

Theorem 5. System (2.2) is persistent in case of $\tau > \tau_0\triangleq -\ln (1-\kappa_1)/r$.

Proof. By Lemma 1, there is

By (3.3), there is

Next, we show that $x(t)\geq x_{\min} > 0$ for sufficiently large $t$, which is given below in two steps:

Step 1: Choose $x_{M} > 0$ and $\varepsilon > 0$, such that

and

It can be shown that $x(t) < x_{M}$ cannot be true for all $t > 0$. Otherwise, due to

thus

implies that $y(t) \leq u(t)\rightarrow \overline{u}(t)$ as $t \rightarrow \infty$ by comparison theorem, where

satisfies

Therefore, $\exists T_2 > 0$, $y(t)\leq \overline{u}(t)+\varepsilon$ when $t > T_2$. This also means that

Select $\overline{j}_{1} \in \mathbb{N}_{+}$ with $\overline{j}_{1} \tau \geq T_2$. Then,

so that $x((\overline{j}_{1}+k) \tau) \geq x(\overline{j}_{1} \tau)\; \varrho_0^k \rightarrow \infty (k \rightarrow \infty)$, which contradicts to Lemma 1. Therefore, $\exists t_1 > 0$, $x\left(t_1\right) \geq x_{M}$.

Step 2: If $x(t)\geq x_{M}\; (\forall t \geq t_1)$ holds, the proof is completed. Define ${\overline{t}}^{+}\triangleq \inf _{t > t_{1}}\left\{x(t) < x_{M}\right\}$, then the value of $\overline{t}$ has two possible cases:

1) $\overline{t} = j\tau, \; j\in \mathbb{N}_{+}$. There is $x(t) \geq x_{M}\; (t \in\left[t_{1}, \overline{t}\right])$ and $x\left(\overline{t}^{+}\right) = \left(1-\kappa_1\right) x\left(\overline{t}\right) < x_{M}$. Since $\varrho_{1} = r-\frac{r x_{M}}{K}-\frac{\gamma x_{M}{ }^{p-1} \mathrm{M}}{q+x_{M}{ }^{p}+a \eta A} < 0$, choose $\overline{j}_{2}, \overline{j}_{3} \in \mathbb{N}_{+}$ such that

Let $T_3\triangleq (\overline{j}_{2}+\overline{j}_{3}) \tau$, then it can be deduced that $\exists t_{2} \in\left(\overline{t}, \overline{t}+T_3\right]$ satisfies $x\left(t_{2}\right) > x_{M}$. Otherwise, $x(t) \leq x_{M}\; (\forall t \in\left(\overline{t}, \overline{t}+T_3\right])$ holds. By (3.5),

Then by (3.6),

which means that $y(t) \leq \overline{u}(t)+\varepsilon$, $\forall t \in\left[\overline{t}+\overline{j}_{2} \tau, \overline{t}+T_3\right]$. In addition, take the same discussion on $\left[\overline{t}+\overline{j}_{2} \tau, \overline{t}+T_3\right]$, and we get $x\left(\overline{t}+T_3\right) \geq x\left(\overline{t}+\overline{j}_{2} \tau\right) \varrho_0^{\overline{j}_{3}}$.

Since

for $t \in\left(\overline{t}, \overline{t}+\overline{j}_{2} \tau\right]$, then there is $x\left(\overline{t}+\tau\right) \geq\left(1-\kappa_1\right) x\left(\overline{t}\right) \exp \left\{\varrho_{1} \tau\right\}$ for $t \in\left(\overline{t}, \overline{t}+\tau\right]$; $x\left(\overline{t}+2 \tau\right) \geq\left(1-\kappa_1\right)^{2} x\left(\overline{t}\right) \exp \left\{2 \varrho_{1} \tau\right\}$ for $t \in\left(\overline{t}+\tau, \overline{t}+2 \tau\right]$. Similarly, $x\left(\overline{t}+\overline{j}_{2} \tau\right) \geq\left(1-\kappa_1\right)^{\overline{j}_{2}} x\left(\overline{t}\right) \exp \left\{\overline{j}_{2} \varrho_{1} \tau\right\}$ for $t \in\left(\overline{t}+\left(\overline{j}_{2}-1\right) \tau, \overline{t}+\overline{j}_{2} \tau\right]$. Then for $t \in\left(\overline{t}, \overline{t}+\overline{j}_{2} \tau\right]$, there is

By (3.6), there is

which contradicts to $x(t) \leq x_{M}$ on $\left(\overline{t}, \overline{t}+T_3\right]$. Thus, $\exists t_{2} \in\left(\overline{t}, \overline{t}+T_3\right]$ satisfies $x\left(t_{2}\right) > x_{M}$.

Denote $\widetilde{t}\triangleq\inf _{t > \overline{t}}\left\{x(t) > x_{M}\right\}$. For $t\in\left(\overline{t}, \widetilde{t}\right)$, $\exists k \leq \overline{j}_{2}+\overline{j}_{3}$ such that $t \in\left(\overline{t}+(k-1) \tau, \overline{t}+k \tau\right]\subset\left(\overline{t}, \overline{t}+T_3\right]$. By (3.7), if $t \in\left(\overline{t}, \overline{t}+\tau\right]$, there is $x\left(\overline{t}+T\right) \geq x\left(\overline{t}^{+}\right) \exp \left\{\varrho_{1} \tau\right\}$; if $t \in\left(\overline{t}+\tau, \overline{t}+2 \tau\right]$,

Similarly, if $t\in\left(\overline{t}+(k-1)\tau, \overline{t}+k \tau\right]$, there is

Define $x'_{\min}\triangleq x_{M}\left(1-\kappa_1\right)^{\overline{j}_{2}+\overline{j}_{3}} \exp \left\{\left(\overline{j}_{2}+\overline{j}_{3}\right) \varrho_{1} \tau\right\}$. Then there is $x(t) \geq x'_{\min}$ for $t\in\left(\overline{t}, \widetilde{t}\right)$. For $t > \widetilde{T}$, we can find a lower bound $x'_{\min}$ of $x$ in a similar procedure.

2) $\overline{t} \neq j \tau, j \in \mathbb{N}_{+}$. Then $x\left(\overline{t}\right) = x_{M}$ and $x(t) > x_{M}$ for $t \in\left[t_{1}, \overline{t}\right)$. $\exists \overline{j}_{3} \in \mathbb{N}_{+}$ such that $\overline{t}\in\left(\overline{j}_{3}\tau, \left(\overline{j}_{3}+1\right) \tau\right)$. For $t \in\left(\overline{t}, \left(\overline{j}_{3}+1\right) \tau\right)$, there are two subcases:

2-ⅰ) $x(t) \leq x_{M}$. In this case, we can claim that $\exists t'_{2} \in\left[\left(\overline{j}_{3}+1\right) \tau, \left(\overline{j}_{3}+1\right) \tau+T_3\right]$ such that $x\left(t'_{2}\right) > x_{M}$. Otherwise, there must be $x(t) \leq x_{M}$, $\forall t \in \left[\left(\overline{j}_{3}+1\right) \tau, \left(\overline{j}_{3}+1\right) \tau+T_3\right]$. By (3.5), $u((\overline{j}_{3}+1)T^{+}) = y\left(\left(\overline{j}_{3}+1\right) \tau^{+}\right)$, thus

for $t \in(j \tau, (j+1)\tau]$, where $\overline{j}_{3}+1 < j < \overline{j}_{3}+1+\overline{j}_{2}+\overline{j}_{3}$. Then, we have

and $y(t) \leq \overline{u}(t)+\varepsilon$ for $t \in\left[\left(\overline{j}_{3}+1\right) \tau+\overline{j}_{2} \tau, \left(\overline{j}_{3}+1\right) \tau+T_3\right]$. Similar to the discussion in the first step, there is $x\left(\left(\overline{j}_{3}+1\right) \tau+T_3\right) \geq x\left(\left(\overline{j}_{3}+1\right) \tau+\overline{j}_{2} \tau\right) \varrho_0^{\overline{j}_{3}}$. In addition, for $t \in\left(\overline{t}, \left(\overline{j}_{3}+1\right) \tau\right)$, there is $x(t) \leq x_{M}$. Integrating (3.7) for $t \in\left(\overline{t}, \left(\overline{j}_{3}+\right.\right. 1) \left.\tau+\overline{j}_{2} \tau\right]$, if $t \in\left(\left(\overline{j}_{3}+1\right) \tau, \left(\overline{j}_{3}+2\right) \tau\right]$,

If $t \in\left(\left(\overline{j}_{3}+2\right) \tau, \left(\overline{j}_{3}+3\right) \tau\right]$, then

Similarly, if $t \in\left(\left(\overline{j}_{3}+1\right) \tau+\left(\overline{j}_{2}-1\right) \tau, \left(\overline{j}_{3}+1\right) \tau+\overline{j}_{2}\tau\right]$, then

Since $\overline{t} \in\left(\overline{j}_{3} \tau, \left(\overline{j}_{3}+1\right) \tau\right)$, then for $t \in\left(\overline{t}, \left(\overline{j}_{3}+1\right) \tau\right]$, there is $x\left(\left(\overline{j}_{3}+1\right) \tau\right) \geq x\left(\overline{t}\right) \exp \left\{\varrho_{1}\tau\right\}$, thus

By (3.6), there is

which contradicts to the assumption that $x(t) \leq x_{M}$ on $\left(\overline{t}, \left(\overline{j}_{3}+1\right) \tau+T_3\right]$. Therefore, $\exists t'_{2} \in\left[\left(\overline{j}_{3}+1\right)\tau, \left(\overline{j}_{3}+1\right) \tau+T_3\right]$ such that $x\left(t'_{2}\right) > x_{M}$. Define $\bar{t}\triangleq \inf _{t > \overline{t}}\left\{x(t) > x_{M}\right\}$, then $x(t) \leq x_{M}$ for $t \in\left(\overline{t}, \bar{t}\right)$. For $t\in\left(\overline{t}, \bar{t}\right)$, $\exists k^{\prime} \leq \overline{j}_{2}+\overline{j}_{3}$ such that $t \in\left(\left(\overline{j}_{3}+1\right) \tau+\left(k^{\prime}-1\right) \tau, \left(\overline{j}_{3}+1\right) \tau+k^{\prime}\tau\right]$, then by (3.7), there are

Thus, for $t \in\left(\left(\overline{j}_{3}+1\right) \tau+\left(k^{\prime}-1\right) \tau, \left(\overline{j}_{3}+1\right) \tau+k^{\prime}\tau\right]$, there is

Define $x_{\min}\triangleq x_{M}\left(1-\kappa_1\right)^{\overline{j}_{2}+\overline{j}_{3}} \exp \left\{(\overline{j}_{2}+\overline{j}_{3}+1)\varrho_{1} \tau\right\} < x'_{\min}$. Then, $x(t) \geq x_{\min}$ for $t \in (\overline{t}, \bar{t})$, while for $t > \bar{t}$, a similar procedure can be adopted.

2-ⅱ) $\exists t \in (\overline{t}, (\overline{j}_{3}+1)\tau)$ such that $x(t) > x_{M}$. Define $\widehat{t}\triangleq \inf _{t > \overline{t}}\left\{x(t) > x_{M}\right\}$, then, $x(t) \leq x_{M}$ for $t \in (\overline{t}, \widehat{t})$. Similarly, (3.7) holds. The integration of (3.7) over $(\overline{t}, \widehat{t})$ gives that $x(t) \geq x (\overline{t}) \exp \left\{\varrho_{1}(t-\overline{t})\right\} \geq x_{M} \exp \left\{\varrho_{1} \tau\right\} > x_{\min}$, while for $t > \widehat{t}$, a similar procedure can be done.

To sum up, $\exists x_{\min} > 0$, $x(t) \geq x_{\min}$ when $t$ is sufficiently large. Therefore, the persistence is proved.

3.2.3.   Coexistence periodic solution

Theorem 6. The system (2.2) can bifurcate a supercritical coexistent periodic solution when $\tau > \tau_0\triangleq -\ln (1-\kappa_1)/r$.

Proof. In order to be consistent with the branching theory in Appendix, let us make the transformation $x_{1}(t) = y(t)$, $x_{2}(t) = x(t)$, which leads to the following system

Thus, we have

where $F_{1}, F_{2}, \theta_{1}, \theta_{2}$ are sufficiently smooth functions, $\theta_{1}, \theta_{2}$ are strictly positive and $F_{2}\left(x_{1}, 0\right) = \theta_{2}\left(x_{1}, 0\right) \equiv 0$. The corresponding subsystem of (3.8) in the case of $x_{2}(t) = 0$ is

which has a stable $\tau_{0}$ periodic solution denoted as

Thus, $\xi = \left(x_{s}, 0\right)$ is the boundary period solution of (3.8). Let $x_{0} = x_{s}(0)$, $\xi(0) = \left(x_{0}, 0\right)$. Then, we obtain

and thus, we have

Therefore, $B C < 0$, which means that (2.2) can bifurcate from the boundary periodic solution to a supercritical coexistent periodic solution.

4.   Numerical simulations

To illustrate the main results and verify the effectiveness of the control method, numerical simulations are implemented in MATLAB programs. The biological parameters in the models are chosen as follows: $r = 0.25$, $K = 200$, $\gamma = 0.5$, $p = 2$, $q = 1.2$, $a = 0.3$, $\eta = 0.2$, $A = 2$, $d = 0.44998$.

4.1.   Verification of (2.1)

By Theorem 1, it can be concluded that the dynamic properties of (2.1) depend on the magnitude of parameter $\upsilon$ for given $\eta A$, as illustrated in Figure 1.

For $\upsilon = 0.8$, there is $\upsilon\gamma < d$, and in this case, $E_K$ is globally asymptotically stable, as illustrated in the left subplot of Figure 1. In such situations, predator species will go extinct, which is not the result we want. To ensure the survival of predators, the conversion rate of predators has to be increased by themselves, i.e., $\upsilon > d/\gamma$.

For $\upsilon = 0.9$, there is $x^* < K$, then $E^*(143.87, 20.19)$ exists and is locally asymptotically stable due to $K < \overline{K} = 287.8$, as illustrated in the middle subplot of Figure 1. It can be easily checked that $\eta A > \underline{\eta A}$, i.e., the additional food availability induces the coexistence of the two species.

For $\upsilon = 0.90004$, there is $K > \overline{K} = 86.8$, then $E^*(43.4, 17)$ is unstable and a limit cycle exists, as illustrated in the right subplot of Figure 1. It can be observed that the limit cycle is unique, so it is orbitally asymptotically stable.

Figure 2 shows the impact of additional food resource on the pest and natural enemy in the coexistent steady state. It can be observed that $x^*$ decreases as $\eta A$ increases, and $y^*$ varies with $\eta A$ and is influenced by $\upsilon$. For $\upsilon = 0.9$, $y^*$ increases as $\eta A$ increases; for $\upsilon = 0.90004$, $y^*$ increases first, and then decreases as $\eta A$ increases; while for $\upsilon = 0.90008$, $y^*$ decreases as $\eta A$ increases.

4.2.   Verification of (2.2)

For (2.2), it is assumed that $\kappa_1 = 0.5$, $\kappa_2 = 0.3$ and $\delta = 1$. When $\tau = 0.5 < \tau_0 = 0.8926$, the pest-eradication periodic solution is locally asymptotically stable and globally attractive, as presented in Figure 3. Although the pest is eradicated in this situation, it requires a frequent control, which is neither necessary nor desirable for practical systems.

When $\tau = 5 > \tau_0$, the pest-eradication periodic solution is unstable, then a coexistence periodic solution exists, as shown in Figure 4.

Figure 5 presents the coexistence periodic solution for different $\tau$. It can be observed that the density of pests at control time increases as $\tau$ increases. This implies that the control period should be neither too small nor too large. This needs to be determined according to the actual situation, to ensure that the amount of pests cannot exceed its allowable threshold.

5.   Conclusions

Agricultural pests are an important factor that endanger agricultural production. Common pest control methods include chemical, biological, and IPM. In order to effectively solve the problem of pest control, we proposed a novel pest-natural enemy model with additional food sources and a generalized Holling type functional response of predators based on the consideration of multiple food sources of natural enemy species. The results showed that the natural enemy has a certain restraining effect on the hostile pests (Theorem 1, Figure 1) and the additional food availability induces a direct impact on the positive equilibrium (Theorem 1, Figure 2). For $\eta A < \underline{\eta A} = (dq-(\upsilon \gamma -d) K^{p})/(\upsilon \gamma -a d)$, $E_{K}$ is globally asymptotically stable. In such situations, predator species go extinct. This is not the result we want. To ensure the survival of predators, there has to be enough extra food available (i.e., $\eta A > \underline{\eta A}$), or the conversion rate of predators has to be increased by themselves.

In order to achieve the effect of rapid pest management, we introduced a periodic impulsive control into the system and established a pest management model. It was shown that the pest-eradication periodic solution exists and its stability depends on the control period (Figure 3) and when the control period is less than a given time threshold (i.e., $\tau < \tau_0\triangleq -\ln (1-\kappa_1)/r$), the pest-eradication periodic solution is globally asymptotically stable. Although the pest is eradicated, it requires a frequent control, which is neither necessary nor desirable for practical systems. In fact, it is not necessary to remove all pests; as long as the amount of pests in the system does not exceed a certain threshold, it will not cause environmental and ecological harm. Thus, we magnified the period over which the control was imposed. The result showed that the system was persistent and a coexistence periodic solution exists as a supercritical branch when the control period exceeds the time threshold.

Simulations indicated that the parameter $\upsilon$ directly impacts the local stability of the coexistence equilibrium. When $\upsilon = 0.9$, $E^*$ is locally asymptotically stable, while when $\upsilon = 0.90004$, it loses the stability and a limit cycle occurs. For the control system, a pest-eradication periodic solution exists when $\tau < 0.893$ (Figure 3), or a coexistence periodic solution exists when $\tau > 0.893$ (Figures 4 and 5). It can be observed that the density of pest at control time increases as $\tau$ increases, so the control period should be neither too small nor too large. This needs to be determined according to the actual situation, to ensure that the amount of pests cannot exceed its allowable threshold. This study not only enriched the related content of population dynamics, but also provided certain reference for the management of plant pests.

Use of AI tools declaration

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

Acknowledgments

We express our warmest thanks to Shasha Song, Handling Editor, for communicating the paper and her useful suggestions. We also express our warmest thanks to the referees for their interest in our work, providing their valued time to review the manuscript carefully and their valuable comments for improving the paper.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

In this section, some notation, definitions, and lemmas to facilitate understanding of the main results are presented ^[41,44]. For an impulse differential equation

Denote $\mathbb{R}_{+} = [0, +\infty)$ and $\mathbb{R}^2_{+} = \{(x_1, x_2)\mid x_1\geq 0, x_2\geq 0\}$. Let $\mathbf{z} = (x_1, x_2)$, $\mathbf{F} = (F_1, F_2)^{T}$. A map $\pi$: $\mathbb{R}_{+}\times \mathbb{R}^2_{+}\rightarrow \mathbb{R}_{+}$ is said to belong $\Pi_0$ if it satisfies

1) $\pi$ is continuous on $(j\tau, (j+1)\tau]\times \mathbb{R}^2_{+}$, $\lim_{(t, \mathbf{z}')\rightarrow (nT^{+}, \mathbf{z})} = \pi(nT^{+}, \mathbf{z})$ exists;

2) $\pi$ is locally Lipschitzian for $\mathbf{z}$.

Definition 1. Let $\pi\in \Pi_0$, $(t, \mathbf{z})\in (j\tau, (j+1)\tau]\times \mathbb{R}^2_{+}$. The upper right derivatives of $\pi(t, (j\tau, (j+1)\tau]\times \mathbb{R}^2_{+})$ with respect to model (A1) is defined as

Lemma 2 (Comparison Theorem). Suppose that $\omega \in \mathrm{PC}^{1}\left(\mathbb{R}_{+}, \mathbb{R}\right)$ and

Then for $t > 0$, there is

Similarly, if the group of non-equations (A2) is reversed, then for $t > 0$ we have

Lemma 3. For $t > 0$, suppose that $\omega \in \mathrm{PC}\left(\mathbb{R}_{+}, \mathbb{R}_{+}\right)$ and

where $f\in \mathrm{PC}\left(\mathbb{R}_{+}, \mathbb{R}_{+}\right), \gamma _{k} \geqslant 0$ and $\omega_{0}$ is constant. Then

Denote $\textbf{z}(0) = \textbf{z}_{0} = (x_{10}, x_{20})^{T}$, $\textbf{z}(t) = (x_1(t), x_2(t))^{T} = \mathbf{\Phi}(t, \mathbf{z}_{0}) = (\Phi_{1}(t, \mathbf{z}_{0}), \Phi_{2}(t, \mathbf{z}_{0}))$. Define

Thus, we have the following lemma.

Lemma 4. In the case of $\mid 1-a_{0}^{\prime}\mid < 1, d_{0}^{\prime} = 0$, there are: a) If $B C \neq 0$, then the system (A1) can branch out from the boundary period solution to a nontrivial periodic solution. Moreover, in the case of $BC < 0$, the periodic solution is a supercritical branching one; in the case of $BC > 0$, the periodic solution is a subcritical branching one; b) If $BC = 0$, it cannot be determined.