An eco-epidemic model for assessing the application of integrated pest management strategies

1.   Introduction

Controlling insects and other arthropods plays a crucial role in agriculture due to the potential ecological and economic consequences of pest outbreaks ^[1,2,3,4,5]. Integrated pest management (IPM) offers a highly effective and safer approach by employing a combination of pest control techniques including biological, cultural and chemical methods. This strategy involves carefully selecting and implementing tactics that are specifically tailored to the unique characteristics of the cropping systems, pest complexes and biological environments involved ^[6,7,8,9].

The main objective of IPM is to keep pest populations below the economic injury level (EIL) ^[10] rather than eradicating them entirely. The selection of suitable tactics is based on the pest densities reaching the economic threshold (ET) ^[11], which minimizes the negative impact of insecticides on non-target pests and helps maintain environmental quality. The effectiveness of IPM has been demonstrated in various applications ^{[12,13,14,15]}.

Biological control, employing microbial pathogens that target pest species, is a crucial pest control technique that has received significant attention in both experimental and theoretical studies. Drawing inspiration from the classical SIR (susceptible-infected-recovered) model introduced by Kermack and McKendrick ^[16], the pest population is divided into susceptible and infected pests. As a result, epidemiological models have attracted substantial interest from researchers in the field ^{[17,18,19,20,21,22]}.

However, conventional ecological models primarily focus on two-species interactions, with many studies emphasizing predator-prey dynamics using Lotka-Volterra type functional responses. Limited attention has been given to IPM models involving a three-species food web. Gakkhar and Naji ^[23] investigated the dynamics of a tri-trophic food web comprising a generalist-specialist-prey system, incorporating a modified Holling type II functional response. Apreutesei ^[24] explored the necessary optimality conditions for a tri-trophic food web system. In 2013, Priyadarshi and Gakkhar ^[25] analyzed the dynamics of a tri-trophic food web system that included a Leslie-Gower type generalist predator. In 2016, Liang et al. ^[26] proposed an eco-epidemiological model for a plant-pest-predator system, considering the impulsive effect. These studies highlight the importance of considering three-species interactions in IPM models and their ecological implications.

To achieve effective pest control, it is vital to have a thorough understanding of the ecological dynamics within the cropping system, encompassing the pests, their predators, the environment, and their intricate interrelationships. However, certain control strategies, such as pesticide spraying, natural enemy releases and infected pest releases are implemented instantaneously or over relatively short periods. Hence, determining the optimal timing and frequency of these actions becomes crucial. In this study, we employ impulsive differential equations to model human interventions and develop a novel mathematical model that captures the dynamics of these actions and their impact on pest control.

The paper is organized as follows. In the following section, we will present a novel three-species mathematical model along with its main biological assumptions. Additionally, we will explore the stability of the susceptible pest-eradication periodic solution in Section 3. Section 4 presents an updated eco-epidemic model that incorporates various frequencies of pesticide sprays, releases of infected pests and releases of natural enemies. The threshold value for achieving the susceptible pest-eradication periodic solution is provided for different scenarios, showcasing the global attractivity of this solution. Subsequently, a concise discussion of the obtained results is presented, followed by potential directions for future research.

2.   Eco-epidemiological model with IPM strategies at a fixed moment

Recently, there has been a growing interest in investigating IPM models using a three-species food web, which includes a predator (the natural enemy of the pest), the pest itself and a basal producer (a renewable biotic resource such as plants or vegetation). Several papers have delved into this topic ^{[23,24,25,26,27,28]}. In this context, let $R$ represent the state of the basal producer, while $P$ and $N$ denote the densities of the pest and natural enemy, respectively. In 2006, Apreutesei ^[24] investigated the following pest-natural enemy-plant model

where $\alpha_P$ is the intrinsic growth rate for the pest, $\alpha_N$ is the death rate for the natural enemy and $\alpha_R$ is the growth rate for the plant. $\beta_N$ describes the search rate of the natural enemy, $\beta_P$ denotes the efficiency rate with which captured pests are converted to the new natural enemy, $\beta_R$ describes the search rate of the pest and $\gamma_R$ denotes the efficiency rate. Assume that all the coefficients are positive constants.

A novel eco-epidemiological model is proposed, which incorporates the utilization of microbial control through pathogenic agents and the release of natural enemies in integrated pest management. Building upon the research conducted by Xiao et al. ^[17,18,19] and Apreutesei ^[24], the pest species $P$ is classified into two categories: susceptible pests and infected pests, with the introduction of a pest pathogen acting as a biotic insecticide within the pest population. Additionally, the following model takes into account the presence of a food source for the pests

the parameter definitions for model (2.2) are summarized in Table 1. Specifically, in this model, $R(t)$ represents a basal producer, which corresponds to a renewable biotic resource such as plants or vegetation ^[26]. Additionally, the non-monotonic (sigmoidal) saturation function is defined as follows

This function is a simplified form of the Monod-Haldane or Holling IV function proposed by Sokol and Howell ^[29] and it describes the collective defense behavior of the prey population when their density is relatively high. Here, $a_X > 0$ denotes the half-saturation constant of species $X$. It is important to note that only the susceptible pests have the ability to feed on the plants with a Holling type-IV response function, while the infected pests are assumed to be too weak to hunt for food.

Building upon the pest-natural enemy model and the implementation of IPM strategies, the combined approach of periodic pesticide spraying, release of infected pests and release of natural enemies is employed to effectively control the pest population. The reformulation of model (2.2) is as follows

The parameter definitions for model (2.3) are summarized in Table 1. It is important to emphasize here that model (2.3) incorporates the impulsive effect at fixed time intervals $T$, representing the simultaneous application of periodic pesticide sprays, release of infected pests and release of natural enemies to control the pest population. At time points $t = nT$, the susceptible pest population experiences an instantaneous killing rate $p_S$ due to pesticide spraying. Importantly, following pesticide spraying, the susceptible pests become weakened and more susceptible to predation by their natural enemies.

3.   The stability of susceptible pest-eradication periodic solution

Since it is well-known that the infective pest does not cause harm to the plant, our main objective in this section is to examine the stability of the periodic solution for eradicating the susceptible pest and establish the conditions for its global attractivity.

In the following, we investigate the existence and stability of of the susceptible pest-eradication solution, denoted as $S(t)\doteq 0$ for all $t\geq 0$. We derive three subsystems from model (2.3):

and

It follows from the first equation of model (3.1) that $I(t) = \exp(-\delta_I(t-nT))I(nT^+)$ for $t\in (nT, (n + 1)T]$. Consequently, we can conclude that $I((n + 1)T) = \exp(-\delta_IT)I(nT^+)$.

Now, combining the second equation of model (3.1) with the previously mentioned equations results in the following

By denoting $I((n+1)T^+)$ as $I(n+1)$, we can rewrite the above equation as the following difference equation

Equation (3.4) represents the stroboscopic map of model (3.1). It captures the relationship between the density of infected pests at two consecutive pulse points. Specifically, the presence of a positive steady state in Eq (3.4) indicates the existence of a positive periodic solution in model (3.1). Thus, our initial investigation focuses on determining the positive steady state of the stroboscopic map (3.4). The derivative of $F(I(n))$ with respect to $I(n)$ is

According to the stability theory of differential equations, it is necessary to ensure that the condition $|F'(I(n))| < 1$ holds true, which is equivalent to the inequality

Furthermore, we continue to investigate the non-negative fixed point of the stroboscopic map (3.4). For simplicity, we denote the fixed point as $\tilde{I}$. It can be deduced from Eq (3.4) that $\tilde{I} = (1+q_I) \exp(-\delta T)\tilde{I} + R_I$. This implies that the stroboscopic map (3.4) possesses a single non-negative fixed point, given by

Analogously, we can obtain the non-negative fixed points of the corresponding stroboscopic map for subsystems (3.2) and (3.3)

given that

for $t\in(nT, (n+1)T]$, respectively. Hence, we can conclude the following result.

Lemma 3.1. Model $\rm (3.1)$ exists a positive periodic solution $I^{*}$. Furthermore, for every positive solution $I(t)$ of model $\rm (3.1)$, it holds that $\lim_{t\to +\infty}|I(t)-I^*(t)| = 0$. Here, $I^*(t)$ is defined as $I^*(t) = \exp(-\delta_I (t-nT))\tilde{I}$, where

for $t\in(nT, (n+1)T)$ and $n\in \mathbb{Z}$.

Similarly, we can deduce that

and

for all $n\in \mathbb{Z}$.

In order to analyze the global asymptotic stability of the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ of model (2.3), we establish the following two results.

Theorem 3.2. If $\min\{\lambda_S, \lambda_I, \lambda_R, \lambda_N\} < 1$, then the susceptible pest-eradication periodic solution $(0, I^*(t), N^*(t), R^*(t))$ of model (2.3) is locally asymptotically stable, where

Theorem 3.3. The periodic solution $(0, I^*(t), N^*(t), R^*(t))$ of model (2.3) is globally asymptotically stable provided that

The proofs for Theorems 3.2 and 3.3 can be found in Appendix A. It is worth noting that in practical pest control, the main objective is to manage the population size of pests within a certain range rather than completely eradicating them. This approach helps maintain the balance of the ecosystem and aligns with the fundamental principles of IPM. Therefore, the pest-eradication periodic solution only exists in theory for the system.

4.   Eco-epidemiological model with optimum pulse timing of IPM strategies at different fixed moment

We now extend model (2.3) by incorporating periodic spraying of pesticides and the release of infected pests and natural enemies at specific fixed time intervals, as suggested by Tang et al. ^[30]. Consequently, we obtain the following revised model.

We assume that the application of pesticides has an impact on the populations of pests and natural enemies. Some parameters remain the same as in model (2.3), other paramaters see Table 1 for details.

Due to the implementation of different strategies for insecticide applications and releases of infected pests in model (4.1)^[30,31], we need to consider two cases regarding the timing of IPM applications.

Case 1. Spraying pesticides more frequently than releasing infected pests and natural enemies.

Case 2. Releases of infected pests and natural enemies more frequent than spraying pesticides.

4.1.   Threshold conditions of Case 1

To facilitate the analysis, we assume that the infected pests and natural enemies are released periodically, with a period denoted as $T_N$. We define $\lambda_{m+1}-\lambda_m = T_N$, where $m \in \mathbb{Z}$ represents an integer. Additionally, we have $\tau_{n+k_p} = T_N + \tau_n$.

Now, we define some variables: $\Delta_i = \tau_{i+1} - \tau_i$ $E_i^h = (hT_N+\tau_i, hT_N+\tau_{i+1}]$, where $i = 0, 1, 2, \cdots, k_p$ $\Delta_0 = \tau_1$, $\Delta_{k_p} = T_N - \tau_{k_p}$, $\tau_0 = 0$, $\tau_{k_p+1} = T_N$ $h$ represents a non-negative number. Furthermore, we analyze model (4.1) by dividing it into three subsystems where $S(t) = 0$ as follows

and

From the first equation of model (4.2) within the interval $E_0^h$, we can deduce that

and

Furthermore, we get

At the moment $hT_N + \tau_i$, where $i = 1, 2, \cdots, k_p$, corresponding to the application of a pesticide, we observe

and

Denote $I[(hT_N)^+] = I(h)$, we can get the following difference equation

It is evident that the linear difference equation mentioned above possesses a unique equilibrium

This equilibrium is considered stable globally, given that $\mathcal{\lambda}_I^1 = (1+q_I)(1-p_I)^{k_p}\exp(-\delta_IT_N) < 1$.

Therefore, model (4.2) has a globally stable periodic solution denoted as $\tilde{I}_{N_1}(t)$ and

where

particularly, $I[(\tau_0)^+] = I_1^*$.

Analogously, we can obtain a globally stable periodic solution denoted as $\tilde{N}_{N_1}(t)$, for model (4.3). This is achieved under the condition that $\mathcal{\lambda}_N^1 = (1+q_N)(1-p_N)^{k_p}\exp(-\delta_NT_N) < 1$ and

where

and

Model (4.3) has a globally stable periodic solution, denoted as $\tilde{R}_{N_1}(t)$, under the condition that $\mathcal{\lambda}_R^1 = \exp(r_RT_N) < 1$ and

where

To analyze the local and global asymptotic stability of the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ of model (4.1), we derive the following two results.

Theorem 4.1. If $\min\{\lambda_S^1, \lambda_I^1, \lambda_R^1, \lambda_N^1\} < 1$, then the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ is locally asymptotically stable for model (4.1), where

Theorem 4.2. The periodic solution $(0, I^*(t), N^*(t), R^*(t))$ is globally asymptotically stable for model (4.1) provided that

The proofs for Theorems 4.1 and 4.2 can be found in Appendix B.

Next, we study the impact of important and critical paramters on the threshold $\mathcal{R}_{\mathcal S}^1$. To accomplish this, we select the killing rate $p_I$ for the infected pest, and the pesticide application frequency $k_p$ as bifurcation parameters, while keeping all other parameters fixed as shown in Figure 1. In fact, the effects of $k_p$ and $p_I$ on $\mathcal{R}_{\mathcal S}^1$ are depicted in Figure 1. The results indicate that if the pesticide kills the infected pests with a relatively lower killing rate $p_I$ (such as $p_2 = 0.2$). Then the threshold value $\mathcal{R}_{\mathcal S}^1$ is a monotonically increasing function with respect to $k_p$, as illustrated in Figure 1(a). If $p_I$ is increased (in this case, $p_2 = 0.7, 0.9$), Figure 1(b), (c) demonstrates that the threshold value $\mathcal{R}_{\mathcal S}^1$ becomes nonmonotonic with respect to $k_p$. This indicates that there exists an optimal number of pesticide applications within the releasing period $T_N$. It can be observed that in this case, the optimal control strategy is to apply pesticides two or three times over the period $T_N$. In this scenario, the threshold value $\mathcal{R}_{\mathcal S}^1$ consistently remains greater than $1$, as depicted in Figure 1(b). However, Figure 1(c) illustrates that $\mathcal{R}_{\mathcal S}^1$ is less than $1$ if pesticide applications are applied twice. This demonstrates that if pesticide applications lead to significant effects on the infected pests, then the number of pesticide applications should be increased.

Furthermore, we can conduct two-parameter bifurcation analyses for $\mathcal{R}_{\mathcal S}^1$, as illustrated in Figure 2. In each subplot, we vary two key parameters simultaneously and observe their impact on $\mathcal{R}_{\mathcal S}^1$. All simulation results presented in Figure 2 indicate that $\mathcal{R}_{\mathcal S}^1$ is highly sensitive to slight variations in the release constants $\theta_I$, $\theta_S$, $\theta_N$, the releasing period $T_N$ and the killing rates $p_I$, $p_S$ and $p_N$. Moreover, these findings emphasize that for a given releasing period, the number of pesticide applications within the period $T_N$, the releasing constants $\theta_I$, $\theta_S$, $\theta_N$ and the killing rates $p_I$, $p_S$, $p_N$ are crucial for effective pest control. This information can assist pest control experts in designing and determining the optimal timing for spray applications and the appropriate release rates.

Although we mentioned in Section 3 that the pest-eradication periodic solution only exist in theory and are difficult to achieve in practical pest control, we have conducted numerical studies to investigate the existence of the pest-eradication periodic solution, which means that the susceptible pest-eradication periodic solution is globally attractive, as shown in Figure 3. The theoretical existence of periodic solutions is indeed correct. In fact, as long as the parameters of the model satisfy the threshold conditions by changing the control strategy in practical pest control, the corresponding pest-eradication periodic solution of system can be obtained.

4.2.   Threshold conditions of Case 2

We assume that the pesticide is sprayed periodically, with the period recorded as $T_p$. This means that $\tau_{n+1}-\tau_n = T_p$ for every $n$ belonging to the set of integers $\mathbb{Z}$. During each pesticide spraying period of $T_p$, susceptible pests and natural enemies are introduced $k_p$ times. Consequently, we have the relationship

Define

where $\Delta_0 = \lambda_R, \Delta_{k_q} = T_p-\lambda_{k_q}, \tau_0 = 0, \tau_{k_ q+1} = T_p,$ $h$ is non-negative number, based on these definitions, we can conclude

and

It follows from model (4.1) that

where

Denote $y([(h+1)T_p]^+)$ as $y(h+1).$ Then, the above equation can be rewritten as the following difference equation

It is evident that the above equation has a unique global equilibrium given by

This equilibrium is considered unique and globally stable under the condition that

Therefore, model (4.5) has a globally stable periodic solution

where

$k = 1, 2, \cdots, k_p$, particularly $I_2[(\lambda_0)^+] = I_2^*$.

Analogously, model (4.6) has a a globally stable periodic solution

as $\mathcal{\lambda}_{N}^2 = (1-p_N)(1+q_N)^k\exp(-\delta_NT_p) < 1$. Here,

$k = 1, 2, \cdots, k_p$ and

Model (4.7) has a a globally stable periodic solution

as $\mathcal{\lambda}_{R}^2 = (1-q_R)\exp(r_RT_P) < 1.$ Here,

and

Hence, we have determined the local and global stability conditions for the periodic solution of model (4.1).

Theorem 4.3. If $min\{\mathcal{\lambda}_S^2, \mathcal{\lambda}_I^2, \mathcal{\lambda}_R^2, \mathcal{\lambda}_N^2\} < 1$, then the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ of model (4.1) is locally asymptotically stable, where

Theorem 4.4. The periodic solution $(0, I^*(t), N^*(t), R^*(t))$ of model (4.1) is globally asymptotically stable provided that

By employing the same methods as those used in Section 4.1, we can investigate the impact of key parameters on the threshold value $\mathcal{R}_\mathcal{S}^2$ and discuss its biological implications. Figure 4 provides detailed information on how different release constants $\theta_I$ affect the threshold value $\mathcal{R}_\mathcal{S}^2$. Regarding the parameters $\delta_I, \delta_S, r_R$, Figure 4 suggests that even small changes can significantly influence $\mathcal{R}_\mathcal{S}^2$, indicating that a slight increase in the release constant $\theta_I$ of infected pests can considerably reduce the threshold parameter $\mathcal{R}_\mathcal{S}^2$. The larger the release rate $\theta_I$ of infective pests, the smaller the resulting threshold value, thereby benefiting pest control efforts.

5.   Concluding remarks

The objective of this paper is to introduce and analyze an eco-epidemiological pest-natural enemy-plant model that investigates the synergistic combination of chemical and biological control methods. Our model examines the effects of varying frequencies of pesticide sprays and releases of infected pests and natural enemies on the suppression of pest populations. Through the use of impulsive differential equations, we establish threshold conditions for achieving a periodic solution that eradicates the susceptible pest population under different scenarios.

We examine two distinct strategies for releasing natural enemies and two different tactics for pest control. One strategy involves more frequent pesticide sprays, while the other focuses on more frequent releases of susceptible pests. We establish sufficient conditions for the global attraction of the periodic solution that eradicates the susceptible pest population. This finding implies that complete elimination of the pest population is achievable under certain conditions regarding the release quantities of infected pests and natural enemies. In fact, in practical pest control, it is difficult to achieve the complete eradication of pests if the implemented control strategies cannot precisely bring the relevant parameter set to meet the theoretical threshold. In other words, the complete eradication of pests is a challenging goal to achieve.

The theoretical findings presented in Section 4 offer several pest management control strategies. Nevertheless, complete eradication of the pest population is often unattainable in practice. Instead, the primary objective is to prevent the pest population from reaching harmful levels and to maintain its density below the economic threshold. To address this goal, we propose the development of a state-dependent impulsive differential equation, which will serve as a threshold control model. This approach will be explored as a potential future direction for our research.

Use of AI tools declaration

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

Acknowledgements

We would like to thank the editor and the anonymous referees for their valuable comments and suggestions that greatly improved the presentation of this work. This work is supported by National Natural Science Foundation of China (12261104, W. J. Qin) and Youth Talent of Xingdian Talent Support Program (XDYC-QNRC-2022-0708, W. J. Qin). Meanwhile, it is partly supported by National Natural Science Foundation of China (12201086, Y. Yang) and Youth Project of Science and Technology Research Program of Chongqing Municipal Education Commission (KJQN202201209, Y. Yang).

Conflict of interest

No potential conflict of interest was reported by the authors.

Appendix

Appendix A

The proof of Theorem 3.2. To establish the local stability of the periodic solution, we employ the small amplitude perturbation method. Let

where $s(t), i(t), n(t),$ and $r(t)$ represent small perturbations. By employing Taylor expansions and disregarding higher-order terms, we can linearize model (2.3), resulting in the following set of linearized equations

where

Let $\Phi(t)$ be the fundamental solution matrix of (A.1), $\Phi(t)$ must satisfy

$\Phi(0) = I_4$ is the identical matrix. For all $t \in (nT, (n+1)T]$, the fundamental solution matrix is

Furthermore, the last four equations of model (A.1) can be expressed as follows

Hence, the eigenvalues of

are

Since $\min\{\lambda_S, \lambda_I, \lambda_R, \lambda_N\} < 1$, then the periodic solution is locally stable. The proof is complete.

The proof of Theorem 3.3. The periodic solution $(0, I^*(t), N^*(t), R^*(t))$ exhibits local asymptotic stability, and we are primarily concerned with its global attractiveness. Let

where $\varepsilon_1 > 0$. Furthermore, we have $I'(t) = \beta S(t)I(t)-\delta_II(t)\geq -\delta_I I(t)$. By utilizing the comparison principle and Lemma (3.1), we can establish the existence of $n_1$ such that $I(t)\geq I^*(t) -\varepsilon_1$ for all $t \geq n_1T$. Similarly, we can find a $n_2$ (here $n_2 > n_1$) such that $N(t)\geq N^*(t)-\varepsilon_1$ for all $t\geq n_2T$. Additionally, there exists a $n_3$ (here $n_3 > n_2$) such that $R(t)\leq R^*(t)-\varepsilon_1$ for all $t\geq n_3T$. Consequently, we can conclude that

for $t \geq n_3T$. The first equation of model (2.3) can be expressed as follows

By integrating over the interval $(n_3 + k)T$ to $(n_3 + k + 1)T$, we can obtain

where

Since $\mathcal{R}_5 < 1$, we get $S(t) \to 0$ as $k \to \infty$. Thus, for an arbitrary positive constant $\varepsilon_2$ small enough, there exists a $n_3$ ($n_3 > n_2$) such that $S(t) < \varepsilon_2$ for all $t \geq n_3T$.

Next, we will demonstrate that $I(t)$ approaches $I^*(t)$. For $0 < \varepsilon_2 < {\delta_I}/{\beta}$, there exists $0 < S(t) < \varepsilon_1$ such that $0 < S(t) < \varepsilon_1$ for $t > n_4T$ (here $n_4 > n_3)$. By considering the second equation of model (2.3), we can conclude that

Let $I_1(t)$ and $I_2(t)$ be the solutions of the following equations, respectively

and

From here it is easy to show that

It can be easily deduced that $I(t)$ approaches $I^*(t)$ as $t$ tends to infinity.

Similarly, we can prove that $N(t)$ tends to $N^*(t)$ as $t$ approaches infinity. For $R(t)$, the following inequality holds

for $t > n_5T$ (here $n_5 > n_4)$. By combining the solution of model (3.3) with the comparison theorem, we obtain the following result

Therefore, as $\varepsilon_2$ approaches 0, we have $R(t) \to R^*(t)$. In fact, if $\mathcal{R}_5 < 1$, then $I(t) \to I^*(t)$, $N(t) \to N^*(t)$, and $R(t) \to R^*(t)$ as $t\to 0$. This implies that the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ is globally asymptotically stable. Thus, the proof is complete.

Appendix B

The proof of Theorem 4.1. We define

where $s(t), i(t), n(t),$ and $r(t)$ represent small perturbations. By linearizing model (4.1) through Taylor expansions and neglecting higher-order terms, the resulting linearized equations are given by

where

Let $\Phi(t)$ be the fundamental solution matrix of (B.1). In this case, $\Phi(t)$ must satisfy the following

where $\Phi(0) = I_4$ is the identical matrix. For all $t \in E_i^h$, the fundamental matrix is

Furthermore, system (B.1) can be expressed as

where

Given that $\min{\lambda_S^1, \lambda_I^1, \lambda_R^1, \lambda_N^1} < 1$, indicating that all eigenvalues of matrix $M$

thus, the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ is locally asymptotically stable. Hence, the proof is complete.

The proof of Theorem 4.2. We define

where $\varepsilon_1 > 0$. By utilizing the comparison principle, we can establish the existence of $n_1$ such that $I(t)\geq I^*(t) -\varepsilon_1$ for all $t \geq n_1T$. Similarly, there exists a $n_2$ (here $n_2 > n_1$) such that $N(t)\geq N^*(t)-\varepsilon_1$ for all $t\geq n_2T$. Additionally, there is a $n_3$ (here $n_3 > n_2$) such that $R(t)\leq R^*(t)+\varepsilon_1$ for all $t\geq n_3T$. Hence, we can conclude that

for $t \geq n_3T$. It is follows from model (4.1) that

Taking integral on $E_i^h, i = 1, 2, \cdots, k_p$ that

and

We define

If $\lambda_S^1 < 1$, we can observe that $S(t) \to 0$ as $k \to \infty$. Consequently, for any sufficiently small positive constant $\varepsilon_2$, there exists an $n_3$ (here $n_3 > n_2$) such that $S(t) < \varepsilon_2$ for all $t \geq n_3T$.

Next, we will prove that $I(t)$ approaches $I^*(t)$. For $0 < \varepsilon_2 < \frac{\delta_I}{\beta}$, there exists $\varepsilon_1$ such that $0 < S(t) < \varepsilon_1$ for $t > n_4T$ (here $n_4 > n_3$). Consequently, we have

Let $I_1(t)$ and $I_2(t)$ denote the solutions of the following equations, respectivelyd

and

Therefore, we can conclude that

So it gives that $\lim_{\varepsilon_2 \to 0}I_i^{**}(t) = I^*(t)$ from $\lim_{t \to \infty}I(t) = I^*(t)$.

Analogously, we can prove that $N(t) \to N^*(t)$ as $t \to \infty$.

For $R(t)$, it gives that

for $t > n_5T$ (here $n_5 > n_4$), and we have

Therefore, we can conclude that $\lim_{\varepsilon_2 \to 0}R(t) = R^*(t)$. In summary, if $\lambda_S^1 < 1$, it follows that $I(t) \to I^*(t)$, $N(t) \to N^*(t)$, and $R(t) \to R^*(t)$ as $t \to 0$. This implies that the periodic solution $(0, I^*(t), N^*(t), R^*(t))$ is globally asymptotically stable. Hence, the proof is complete.