Modelling and analysis of periodic impulsive releases of the Nilaparvata lugens infected with wStri-Wolbachia

Abstract

In this paper, we formulate a population suppression model and a population replacement model with periodic impulsive releases of Nilaparvata lugens infected with wStri. The conditions for the stability of wild-$N.lugens$-eradication periodic solution of two systems are obtained by applying the Floquet theorem and comparison theorem. And the sufficient conditions for the persistence in the mean of wild $N.lugens$ are also given. In addition, the sufficient conditions for the extinction and persistence of the wild $N.lugens$ in the subsystem without wLug are also obtained. Finally, we give numerical analysis which shows that increasing the release amount or decreasing the release period are beneficial for controlling the wild $N.lugens$, and the efficiency of population replacement strategy in controlling wild populations is higher than that of population suppression strategy under the same release conditions.

Keywords:

• Nilaparvata lugens
• wStri-Wolbachia
• cytoplasmic incompatibility
• periodic release
• global stability

Mathematics Subject Classifications:

• 00A71
• 92B05
• 34A37

1. Introduction

Nilaparvata lugens (N. lugens) is a monophagous pest, which can only feed and reproduce on rice and common wild rice. It is the most destructive pest on rice in many Asian countries by sucking rice phloem sap and transmitting rice ragged stunt virus ($\mathbf{R}\mathbf{R}\mathbf{S}\mathbf{V}$) [1,2]. Few feasible control strategies are obtainable because of the evolution of high levels of insecticide resistance, and new environmentally friendly methods are urgently needed [3]. In recent studies, disease control methods based on artificial Wolbachia infection to inhibit mosquito vector-borne pathogens have been applied [4,5]. Scientists are currently working on similar methods for controlling agricultural pests [6].

Nilaparvata lugens can be naturally infected by Wolbachia strain wLug, which lacks the ability to induce cytoplasmic incompatibility ($\mathbf{C}\mathbf{I}$) [7]. Many experimental results show that the $N.lugens$ infected with wLug has stronger reproductive ability than the uninfected $N.lugens$ , and show imperfect maternal transmission characteristics [6,8,9]. Fortunately, Gong et al. [6] successfully developed a stable artificial Wolbachia infection of $N.lugens$ by introducing the Wolbachia strain wStri from $L.striatellus$ host into $N.lugens$ . The results [6] showed that $N.lugens$ infected with wStri maintained perfect maternal transmission and induced moderately high levels of CI. When the wStri-infected males mated with either uninfected or the wLug-infected females, the mean hatch rates per female were $32.5\mathrm{\%}$ and $27.7\mathrm{\%}$, respectively. The results of Gong et al. [6] lay a foundation for future experiments on paddy fields.

In order to theoretically study the interaction mechanism between $N.lugens$ infected with wStri and wild $N.lugens$ population, Liu and Zhou [10] proposed and studied a Wolbachia spreading dynamics model in $N.lugens$ with two strains, and obtained sufficient conditions for the $N.lugens$ infected with wStri to invade wild $N.lugens$ successfully. But the authors did not address when and how many $N.lugens$ infected with wStri is released to suppress or replace wild $N.lugens$ , so studying the release of wStri-infected $N.lugens$ to control wild $N.lugens$ for future field trials with wStri is biologically significant. Currently, there are few models to study how to release the $N.lugens$ infected Wolbachia, but many mosquito population dynamics models have been proposed and studied [11–24]. Cai et al. [11] considered three strategies for continuous release of sterile mosquitoes: constant release rate, release rate proportional to the wild mosquitoes and proportional release rate with saturation, and developed corresponding continuous mathematical models to study the influences of the three release strategies on the interactive dynamics of mosquitoes. In practice, continuous release of sterile mosquitoes is difficult to realize in practical applications. Huang et al. [16] formulated and investigated two mathematical models with impulsive releases of sterile mosquitoes. The first model considered the periodic pulse release of sterile mosquitoes strategy and obtained a sufficient condition for the stability of the wild mosquito-eradication periodic solution. The second model considered the state-feedback pulse release strategy and proved the existence of the first-order periodic solution. In studies [17–20], the authors adopted a new modelling idea that only those sexually active sterile mosquitoes were considered in the modelling process. Yu et al. [19] developed and analyzed a population suppression model considering that the release period was longer than the sexual life of sterile mosquitoes. They obtained sufficient conditions for the global asymptotic stability of positive periodic solutions. Later, Zheng et al. [17] considered the following situation that the release period was shorter than the sexual lifespan of sterile male mosquitoes. Li and Ai [18] incorporated the maturation process of mosquito larvae to adults into their model and employed time delay to describe the maturation period of the larvae, the results showed that the delay affects the control of wild mosquitoes. In addition, some scholars considered the interaction between mosquitoes and disease transmission (such as $Denguefever$, Zika, etc.), and established some mathematical models to control mosquito-borne diseases [23,24]. For example, Taghikahani et al. [23] formulated a new two-sex mathematical model for the population ecology of dengue fever and Wolbachia-infected mosquitoes, and used it to evaluate the impact of periodic release the Wolbachian-infected mosquitoes on the population-level. However, we notice from studies [11–24] that most models were established around the release of male mosquitoes infected with Wolbachia (i.e. population suppression strategy), and the research on the impulsive release strategy of female mosquitoes infected with Wolbachia mostly adopts numerical simulation method.

All the models mentioned above about the release strategies of mosquitoes infected with Wolbachia also provide good help for the release of $N.lugens$ infected with wStri. Compared with the transmission characteristics of mosquitoes infected with Wolbachia, the $N.lugens$ infected with Wolbachia has many differences. For example, the cytoplasmic incompatibility induced by the mating of male $N.lugens$ infected with wStri with wild uninfected female $N.lugens$ or wild female $N.lugens$ infected with wlug is incomplete [6]. To determine when is the best time to release $N.lugens$ infected with wStri and the number of $N.lugens$ infected with wStri per release. In 2023, Liu et al. [25] established and discussed two semi-continuous models with state-feedback impulsive releases of $N.lugens$ infected with wStri. But the authors only considered the scenario where all wild $N.lugens$ populations were uninfected $N.lugens$ , and there are few studies of wild $N.lugens$ including both the wild $N.lugens$ uninfected and infected with wLug. In this paper, we adapt the modelling idea in [10] and consider the periodic impulsive release of male or female $N.lugens$ infected with wStri into the field. And then according to the transmission characteristics of two Wolbachia strains in $N.lugens$ , we establish two impulsive mathematical models with periodic impulsive release of $N.lugens$ infected with wStri: population suppression model and population replacement model. We theoretically discuss the stability of the wild-$N.lugens$ -eradication periodic solution for both models and the corresponding numerical analyses are also carried out. Furthermore, the sufficient conditions of persistence in the mean of the wild $N.lugens$ are also established. Finally, we compare the degree of control of two periodic impulsive release strategies on the wild $N.lugens$ within a short time.

This paper is organized as follows: In Section 2, population suppression and population replacement models with periodic impulsive release are proposed. In Section 3, we study the stability of wild-$N.lugens$ -extinction periodic solution of two models, respectively, and discuss the persistence of wild $N.lugens$ . Then we numerically analyze the influences of release period and release amount on the control of wild $N.lugens$ in Section 4. Finally, we present our conclusions in Section 5.

2. Model formulation

Let $U_f\left(t\right)$ and $U_m\left(t\right)$ be the densities of the wild female $N.lugens$ and male $N.lugens$ uninfected by wLug and wStri at time t, respectively. $L_f\left(t\right)$ and $L_m\left(t\right)$ denote the densities of the wild female $N.lugens$ and male $N.lugens$ infected with wLug at time t, respectively. $S_f\left(t\right)$ and $S_m\left(t\right)$ represent the densities of the female $N.lugens$ and male $N.lugens$ infected with wStri at time t.

In order to establish the mathematical model more conveniently, we give a chart to show all possible mating patterns of $N.lugens$ according to the research conclusions of studies [10,25]. The fourth line of Figure 1 shows that the males $N.lugens$ infected with wStri $S_m$ can induce higher CI level when they mate with $U_f$ and $L_f$, $0<c_1<1$ and $0<c_2<1$ represent the CI intensity of $S_m$ against $U_f$ and $L_f$, respectively. The third and fourth columns of Figure 1 indicate that the $N.lugens$ infected with wStri has perfect maternal transmission, while the wild $N.lugens$ infected with wLug has imperfect maternal transmission. $0<\theta <1$ is the percentage of uninfected progeny produced by a wLug-infected mother. Ju et al. [8] shows that Wolbachia of the wlug strain has a reproductive promoting effect on its natural host $N.lugens$ , but wStri does not. Therefore, we assume that the uninfected $N.lugens$ female $U_f$ and the wStri-infected $N.lugens$ female $S_f$ have the same birth rate, while the birth rate of wLug-infected $N.lugens$ female $L_f$ is greater than $U_f$ and $S_f$. Let a>0 be the birth rate of the uninfected $N.lugens$ , and $\rho a$ is the birth rate of the wLug-infected $N.lugens$ , where $\rho >1$. Moreover, we adhere to the conventional method by substituting the natural death of $N.lugens$ population with a Logistic-like density dependent term [10,25]. Let $d_1>0$, $d_2>0$ and $d_3>0$ denote the decay rate constants of the uninfected $N.lugens$ , wLug-infected $N.lugens$ and wStri-infected $N.lugens$ , respectively. Similar to Refs. [10,26,27], we always assume that the proportion of individuals born male is equal to the proportion of individuals born female in this paper, that is, $U_f=U_m$, $L_f=L_m$ and $S_f=S_m$.

Figure 1. All possible mating patterns of $N.lugens$. The wild $N.lugens$ infected with wLug has imperfect maternal transmission.

From Figure 1, we give all the mating patterns of $N.lugens$ , such as $U_f$ come from the mating patterns $U_f\times U_m$, $U_f\times L_m$ and $U_f\times S_m$. Then the expression $\frac{U_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}$ denotes the probability of a female $N.lugens$ mates with an uninfected male $N.lugens$, $\frac{L_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}$ represents the probability of a female $N.lugens$ mates with a wLug-infected male $N.lugens$ , and $\frac{S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}$ gives the probability of incompatible crossing that a wild female $N.lugens$ mates with a wStri-infected male $N.lugens$ . According to the above assumptions and mating patterns of $N.lugens$ , we propose the following two mathematical models.

2.1. Population suppression model with periodic impulsive release

To begin with, we consider the case that the wStri-infected males $N.lugens$ are released into the paddy field. It is easy to see from Figure 1 that the mating pattern of the third column will not appear in this ecosystem. According to the mating patterns of $N.lugens$, we first propose the following population suppression model with periodic release: (1) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}{U}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_f\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{U}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_m\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{L}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_f\left(t\right),\\ \frac{\mathrm{d}{L}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_m\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }{U}_f\left(t\right)& =0,\mathrm{\Delta }{U}_m\left(t\right)=0,\mathrm{\Delta }{L}_f\left(t\right)=0,\mathrm{\Delta }{L}_m\left(t\right)=0,\\ \mathrm{\Delta }{S}_m\left(t\right)& =b,\end{array}t=n\tau ,n\in Z_{+},\end{array}\end{array}$(1) where $b\ge 0$ represents the release amount of male $N.lugens$ infected with wStri each time. τ is the release period, and $\mathrm{\Delta }{U}_f\left(t\right)=U_f\left(t^{+}\right)-U_f\left(t\right),\mathrm{\Delta }{U}_m\left(t\right)=U_m\left(t^{+}\right)-U_m\left(t\right)$, $\mathrm{\Delta }{L}_f\left(t\right)=L_f\left(t^{+}\right)-L_f\left(t\right)$, $\mathrm{\Delta }{L}_m\left(t\right)=L_m\left(t^{+}\right)-L_m\left(t\right)$, $\mathrm{\Delta }{S}_m\left(t\right)=S_m\left(t^{+}\right)-S_m\left(t\right)$. $N\left(t\right)=U_f\left(t\right)+U_m\left(t\right)+L_f\left(t\right)+L_m\left(t\right)+S_m\left(t\right)$.

Denote $U=U_m+U_f$ and $L=L_m+L_f$. Then model (1(1) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}{U}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_f\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{U}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_m\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{L}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_f\left(t\right),\\ \frac{\mathrm{d}{L}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_m\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }{U}_f\left(t\right)& =0,\mathrm{\Delta }{U}_m\left(t\right)=0,\mathrm{\Delta }{L}_f\left(t\right)=0,\mathrm{\Delta }{L}_m\left(t\right)=0,\\ \mathrm{\Delta }{S}_m\left(t\right)& =b,\end{array}t=n\tau ,n\in Z_{+},\end{array}\end{array}$(1) ) can be simplified to the following model: (2) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aU\left(t\right)(1-\frac{2c_1{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)})-d_{1}N\left(t\right)U\left(t\right)\\ & +\frac{1}{2}\rho a\theta L\left(t\right)(1-\frac{2c_2{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)}),\\ \frac{\mathrm{d}L\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L\left(t\right)(1-\frac{2c_2{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)})-d_{2}N\left(t\right)L\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }U\left(t\right)& =0,\mathrm{\Delta }L\left(t\right)=0,\mathrm{\Delta }{S}_m\left(t\right)=b,t=n\tau ,n\in Z_{+}.\end{array}\end{array}\end{array}$(2) Using a linear scaling for system (2(2) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aU\left(t\right)(1-\frac{2c_1{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)})-d_{1}N\left(t\right)U\left(t\right)\\ & +\frac{1}{2}\rho a\theta L\left(t\right)(1-\frac{2c_2{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)}),\\ \frac{\mathrm{d}L\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L\left(t\right)(1-\frac{2c_2{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)})-d_{2}N\left(t\right)L\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }U\left(t\right)& =0,\mathrm{\Delta }L\left(t\right)=0,\mathrm{\Delta }{S}_m\left(t\right)=b,t=n\tau ,n\in Z_{+}.\end{array}\end{array}\end{array}$(2) ), $\begin{array}{rl}s& =\frac{a}{2}t,x=\frac{2d_3}{a}U,y=\frac{2d_3}{a}L,z=\frac{2d_3}{a}{S}_m,{\delta }_1=\frac{d_1}{d_3},{\delta }_2=\frac{d_2}{d_3},\\ \text{β}& =\frac{2d_{3}b}{a},T=\frac{a\tau }{2},\end{array}$and rewriting $d./ds$ as $d./\mathrm{d}t$, system (2(2) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aU\left(t\right)(1-\frac{2c_1{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)})-d_{1}N\left(t\right)U\left(t\right)\\ & +\frac{1}{2}\rho a\theta L\left(t\right)(1-\frac{2c_2{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)}),\\ \frac{\mathrm{d}L\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L\left(t\right)(1-\frac{2c_2{S}_m\left(t\right)}{U\left(t\right)+L\left(t\right)+2S_m\left(t\right)})-d_{2}N\left(t\right)L\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }U\left(t\right)& =0,\mathrm{\Delta }L\left(t\right)=0,\mathrm{\Delta }{S}_m\left(t\right)=b,t=n\tau ,n\in Z_{+}.\end{array}\end{array}\end{array}$(2) ) can be transmuted into (3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3)

2.2. Population replacement model with periodic impulsive release

We also establish a population replacement model with periodic release. The females $N.lugens$ infected with wStri and males $N.lugens$ infected with wStri are released into the farmland ecosystem, then the farmland ecosystem has uninfected $N.lugens$ , wLug-infected $N.lugens$ and wStri-infected $N.lugens$ population. According to the mating patterns of $N.lugens$ , a population replacement model with periodic impulsive release is given as follows: (4) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}{U}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_f\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{U}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_m\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{L}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_f\left(t\right),\\ \frac{\mathrm{d}{L}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_m\left(t\right),\\ \frac{\mathrm{d}{S}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{S}_f\left(t\right)-d_{3}N\left(t\right)S_f\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{S}_f\left(t\right)-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }{U}_f\left(t\right)& =0,\mathrm{\Delta }{U}_m\left(t\right)=0,\\ \mathrm{\Delta }{L}_f\left(t\right)& =0,\mathrm{\Delta }{L}_m\left(t\right)=0,\\ \mathrm{\Delta }{S}_f\left(t\right)& =\frac{1}{2}b,\mathrm{\Delta }{S}_m\left(t\right)=\frac{1}{2}b,\end{array}\}t=n\tau ,\end{array}\end{array}$(4) where $N\left(t\right)=U_f\left(t\right)+U_m\left(t\right)+L_f\left(t\right)+L_m\left(t\right)+S_f\left(t\right)+S_m\left(t\right)$.

Denote $U=U_m+U_f$, $L=L_m+L_f$ and $S=S_m+S_f$. Then model (4(4) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}{U}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_f\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{U}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_m\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{L}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_f\left(t\right),\\ \frac{\mathrm{d}{L}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_m\left(t\right),\\ \frac{\mathrm{d}{S}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{S}_f\left(t\right)-d_{3}N\left(t\right)S_f\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{S}_f\left(t\right)-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }{U}_f\left(t\right)& =0,\mathrm{\Delta }{U}_m\left(t\right)=0,\\ \mathrm{\Delta }{L}_f\left(t\right)& =0,\mathrm{\Delta }{L}_m\left(t\right)=0,\\ \mathrm{\Delta }{S}_f\left(t\right)& =\frac{1}{2}b,\mathrm{\Delta }{S}_m\left(t\right)=\frac{1}{2}b,\end{array}\}t=n\tau ,\end{array}\end{array}$(4) ) can be simplified to the following model: (5) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aU\left(t\right)(1-\frac{c_{1}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)})-d_{1}N\left(t\right)U\left(t\right)\\ & +\frac{1}{2}\rho a\theta L\left(t\right)(1-\frac{c_{2}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)}),\\ \frac{\mathrm{d}L\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L\left(t\right)(1-\frac{c_{2}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)})-d_{2}N\left(t\right)L\left(t\right),\\ \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aS\left(t\right)-d_{3}N\left(t\right)S\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }U\left(t\right)& =0,\mathrm{\Delta }L\left(t\right)=0,\mathrm{\Delta }S\left(t\right)=b,\end{array}t=n\tau .\end{array}\end{array}$(5) Using a linear scaling for system (5(5) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aU\left(t\right)(1-\frac{c_{1}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)})-d_{1}N\left(t\right)U\left(t\right)\\ & +\frac{1}{2}\rho a\theta L\left(t\right)(1-\frac{c_{2}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)}),\\ \frac{\mathrm{d}L\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L\left(t\right)(1-\frac{c_{2}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)})-d_{2}N\left(t\right)L\left(t\right),\\ \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aS\left(t\right)-d_{3}N\left(t\right)S\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }U\left(t\right)& =0,\mathrm{\Delta }L\left(t\right)=0,\mathrm{\Delta }S\left(t\right)=b,\end{array}t=n\tau .\end{array}\end{array}$(5) ), $\begin{array}{rl}s& =\frac{a}{2}t,X=\frac{2d_3}{a}U,Y=\frac{2d_3}{a}L,Z=\frac{2d_3}{a}S,{\delta }_1=\frac{d_1}{d_3},{\delta }_2=\frac{d_2}{d_3},\\ \text{β}& =\frac{2d_{3}b}{a},T=\frac{a\tau }{2},\end{array}$and rewriting $d./ds$ as $d./\mathrm{d}t$, then we transmute system (5(5) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aU\left(t\right)(1-\frac{c_{1}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)})-d_{1}N\left(t\right)U\left(t\right)\\ & +\frac{1}{2}\rho a\theta L\left(t\right)(1-\frac{c_{2}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)}),\\ \frac{\mathrm{d}L\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L\left(t\right)(1-\frac{c_{2}S\left(t\right)}{U\left(t\right)+L\left(t\right)+S\left(t\right)})-d_{2}N\left(t\right)L\left(t\right),\\ \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}aS\left(t\right)-d_{3}N\left(t\right)S\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }U\left(t\right)& =0,\mathrm{\Delta }L\left(t\right)=0,\mathrm{\Delta }S\left(t\right)=b,\end{array}t=n\tau .\end{array}\end{array}$(5) ) into (6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6)

3. Main results

We will study the dynamics of models (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) and (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) in this section. For simplicity, denote ${\delta }_0=max\{{\delta }_1,{\delta }_2\},\delta =min\{{\delta }_1,{\delta }_2\},d=\frac{\rho }{\delta }.$

Definition 3.1

1. The population χ is said to be extinct if $\underset{t\to +\mathrm{\infty }}{lim}\chi \left(t\right)=0$.

2. The population χ is said to be strongly persistent in the mean if $\underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^t\chi \left(s\right)ds>0$.

3.1. The dynamics of model (3)

Let $\phi \left(t\right)=\left(x\right(t),y(t),z(t){)}^T$ be the solution of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) with initial value $x\left(0\right)>0$, $y\left(0\right)\ge 0$ and $z\left(0\right)\ge 0$. Clearly, $\phi :R_{+}\to R_{+}^3$ is a piecewise continuous function, where $R_{+}=\{x:x>0\},R_{+}^3=\left\{\right(x,y,z):x\ge 0,y\ge 0,z\ge 0\}$. From [28], the global existence and uniqueness of solutions of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is guaranteed by the smoothness properties of $f=(f_1,f_2,f_3{)}^T$, which denotes the mapping defined by the right-hand side of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ). Denote $\begin{array}{rl}{R}_{m1}& =\frac{1-c_1}{{\delta }_1}-\frac{1}{T}\ln (1+z^{\ast }T),R_{m2}=\frac{\rho (1-\theta )(1-c_2)}{{\delta }_2}-\frac{1}{T}\ln (1+z^{\ast }T),\\ R_{m3}& =T+\frac{2c_1{\delta }_1}{2\mathrm{d}{\delta }_1-1}\ln \left(\frac{d+g^{\ast }\left(1-e^{-\mathrm{d}t}\right)+2{\delta }_1{g}^{\ast }\mathrm{d}{e}^{-\mathrm{d}t}}{d+2{\delta }_1{g}^{\ast }d}\right)-{\delta }_1\ln \left(1+\frac{g^{\ast }\left(1-e^{-\mathrm{d}t}\right)}{d}\right),\\ R_{m4}& =\rho (1-\theta )(T+\frac{2c_2}{d}\ln \frac{d+g^{\ast }(1+e^{-\mathrm{d}t})}{d+2g^{\ast }})-{\delta }_2\ln (1+\frac{g^{\ast }(1-e^{-\mathrm{d}T})}{d}),\\ R_{m5}& =min\{1-c_1,\rho (1-c_2\left)\right\}-{\delta }_0\frac{1}{T}\ln (1+z^{\ast }T),\end{array}$where $z^{\ast }$ and $g^{\ast }$ will be given in the following Lemmas 3.1 and 3.2.

If $x\left(t\right)=0$ and $y\left(t\right)=0$, the subsystem of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is presented by impulsive differential equations (7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }z\left(t\right)& =\text{β},t=nT.\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(7)

Lemma 3.1

System (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }z\left(t\right)& =\text{β},t=nT.\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(7) ) has a unique periodic solution $\tilde{z}\left(t\right)$ with period T, and for any solution $z\left(t\right)$ of system (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }z\left(t\right)& =\text{β},t=nT.\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(7) ) with $z\left(0^{+}\right)\ge 0$, we have $|z\left(t\right)-\tilde{z}\left(t\right)|\to 0$ as $t\to +\mathrm{\infty }$, where $\begin{array}{rl}\tilde{z}\left(t\right)& =\frac{z^{\ast }}{1+z^{\ast }(t-nT)},t\in (nT,(n+1\left)T\right],\\ \tilde{z}\left(0^{+}\right)& =z^{\ast }>0,\end{array}$and $z^{\ast }=\frac{\text{β}+\sqrt{{\text{β}}^2+\frac{4\text{β}}{T}}}{2}$.

Proof.

Integrating the first equation of (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }z\left(t\right)& =\text{β},t=nT.\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(7) ) between $[nT,(n+1\left)T\right]$, we have $\begin{array}{r}z\left(t\right)=\frac{z\left(n{T}^{+}\right)}{1+(t-nT)z\left(n{T}^{+}\right)},t\in (nT,(n+1\left)T\right].\end{array}$By the second equation of (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }z\left(t\right)& =\text{β},t=nT.\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(7) ), we can obtain the stroboscopic map: $\begin{array}{r}{z}_{n+1}=z\left(\right(n+1\left)T^{+}\right)=\frac{z\left(n{T}^{+}\right)}{1+Tz\left(n{T}^{+}\right)}+\text{β}=\frac{z_n}{1+T{z}_n}+\text{β}=f\left(z_n\right).\end{array}$It is easy to calculate that the above discrete system has a unique positive fixed point $\begin{array}{r}{z}^{\ast }=\frac{\text{β}+\sqrt{{\text{β}}^2+\frac{4\text{β}}{T}}}{2}.\end{array}$Because $|\frac{\mathrm{\partial }f}{\mathrm{\partial }z}{|}_{z=z^{\ast }}|=\frac{1}{(1+T{z}^{\ast }{)}^2}<1$, the unique positive fixed point $z^{\ast }$ is locally asymptotically stable. In addition, $\begin{array}{rl}|z_{n+1}-z^{\ast }|& =|\frac{z_n}{1+T{z}_n}-\frac{z^{\ast }}{1+T{z}^{\ast }}|=\frac{|z_n-z^{\ast }|}{(1+T{z}_n)(1+T{z}^{\ast })}\le \frac{1}{1+T{z}^{\ast }}|z_n-z^{\ast }|\to 0,\\ & \mathrm{a}\mathrm{s}n\to +\mathrm{\infty }.\end{array}$So fixed point $z^{\ast }$ is globally asymptotically stable. Further, the positive periodic solution of system (7(7) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }z\left(t\right)& =\text{β},t=nT.\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(7) ) $\begin{array}{rl}\tilde{z}\left(t\right)& =\frac{z^{\ast }}{1+z^{\ast }(t-nT)},t\in (nT,(n+1\left)T\right],\\ \tilde{z}\left(0^{+}\right)& =z^{\ast }>0,\end{array}$is also locally asymptotically stable.

Furthermore, $\begin{array}{r}|z\left(t\right)-\tilde{z}\left(t\right)|=\left|\frac{z\left(t\right)\tilde{z}\left(t\right)}{z_n{z}^{\ast }}\right||z_n-z^{\ast }|<|z_n-z^{\ast }|\to 0,\mathrm{a}\mathrm{s}t\to +\mathrm{\infty }.\end{array}$The proof is completed.

Consider the following system: (8) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}& =-\mathrm{d}g\left(t\right)-g^2\left(t\right),t\ne nT,\\ g\left(t^{+}\right)& =g\left(t\right)+\text{β},t=nT,\\ g\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(8) Then we give the following lemma according to Lemma 3.2 in Huang et al. [16].

Lemma 3.2

See [16]

System (8(8) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}& =-\mathrm{d}g\left(t\right)-g^2\left(t\right),t\ne nT,\\ g\left(t^{+}\right)& =g\left(t\right)+\text{β},t=nT,\\ g\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(8) ) has a unique periodic solution $\widetilde{g\left(t\right)}$ with period T, and for any solution $g\left(t\right)$ of system (8(8) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}g\left(t\right)}{\mathrm{d}t}& =-\mathrm{d}g\left(t\right)-g^2\left(t\right),t\ne nT,\\ g\left(t^{+}\right)& =g\left(t\right)+\text{β},t=nT,\\ g\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(8) ) with $g\left(0^{+}\right)\ge 0$, $|g\left(t\right)-\tilde{g}\left(t\right)|\to 0$ as $t\to +\mathrm{\infty }$, where $\begin{array}{rl}\tilde{g}\left(t\right)& =\frac{g^{\ast }\mathrm{d}{e}^{-d(t-nT)}}{d+g^{\ast }(1-e^{-d(t-nT)})},t\in (nT,(n+1\left)T\right],\\ \tilde{g}\left(0^{+}\right)& =g^{\ast }>0,\end{array}$and $g^{\ast }$ is the positive root of the following equation $\begin{array}{r}{g}^2+(d-\text{β})g-\frac{\text{β}d}{(1-e^{-\mathrm{d}T})}=0.\end{array}$

In the following, we will discuss the stability of wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{z}(t\left)\right)$. We first show the periodic solution $(0,0,\tilde{z}(t\left)\right)$ is locally asymptotically stable and then prove it is also a global attractor.

Theorem 3.1

The wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{z}(t\left)\right)$ of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is locally asymptotically stable if $R_{m1}<0$ and $R_{m2}<0$.

Proof.

The local stability of periodic solution $(0,0,\tilde{z}(t\left)\right)$ can be determined by considering the small-amplitude perturbation of the solution.

Define $u\left(t\right)=x\left(t\right),v\left(t\right)=y\left(t\right),w\left(t\right)=z\left(t\right)-\tilde{z}\left(t\right).$ The linearized system of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) in $(nT,(n+1\left)T\right]$ is obtained as follows: $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}u\left(t\right)}{\mathrm{d}t}& =\left(\right(1-c_1)-{\delta }_1\tilde{z}(t\left)\right)u\left(t\right)+\rho \theta (1-c_2)v\left(t\right),\\ \frac{\mathrm{d}v\left(t\right)}{\mathrm{d}t}& =\left(\rho \right(1-\theta \left)\right(1-c_2)-{\delta }_2\tilde{z}(t\left)\right)v\left(t\right),\\ \frac{\mathrm{d}w\left(t\right)}{\mathrm{d}t}& =-\tilde{z}\left(t\right)u\left(t\right)-\widetilde{z\left(t\right)}v\left(t\right)-2\tilde{z}\left(t\right)w\left(t\right).\end{array}\end{array}$By simply calculating, the fundamental solution matrix in interval $(nT,(n+1\left)T\right]$ can be given by $\begin{array}{r}\varphi \left(t\right)=\left(\begin{array}{ccc}{e}^{{\int }_{nT}^t(1-c_1)-{\delta }_1\tilde{z}\left(p\right)\mathrm{d}p}& {\mathrm{♯}}_1& 0\\ 0& e^{{\int }_{nT}^t\rho (1-\theta )(1-c_2)-{\delta }_2\tilde{z}\left(p\right)\mathrm{d}p}& 0\\ {\mathrm{♯}}_2& {\mathrm{♯}}_3& e^{{\int }_{nT}^t-2\tilde{z}\left(p\right)\mathrm{d}p}\end{array}\right),\end{array}$where the exact expression of function ${\mathrm{♯}}_i$ (i=1,2,3) is not presented because it is not used in the calculation of the eigenvalues of matrix Φ.

It follows from the fourth equations of (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) that $\begin{array}{r}\left(\begin{array}{c}u\left(\right(n+1\left)T^{+}\right)\\ v\left(\right(n+1\left)T^{+}\right)\\ w\left(\right(n+1\left)T^{+}\right)\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}u\left(\right(n+1\left)T\right)\\ v\left(\right(n+1\left)T\right)\\ w\left(\right(n+1\left)T\right)\end{array}\right).\end{array}$Further from the Floquet theory, we obtain that the periodic solution $(0,0,\tilde{z}(t\left)\right)$ is locally asymptotically stable, which can be determined by the absolute values of all eigenvalues of matrix $\mathrm{\Phi }=\varphi \left(T\right)$ are less than 1. The eigenvalues of matrix Φ are $\begin{array}{rl}{\lambda }_1& =e^{{\int }_0^T(1-c_1)-{\delta }_1\tilde{z}\left(t\right)\mathrm{d}t},{\lambda }_2=e^{{\int }_0^T\rho (1-\theta )(1-c_2)-{\delta }_2\tilde{z}\left(t\right)\mathrm{d}t},\\ {\lambda }_3& =e^{-{\int }_0^{T}2\tilde{z}\left(t\right))\mathrm{d}t}<1.\end{array}$According to the conditions given in Theorem 3.1, we have ${\lambda }_1<1$ and ${\lambda }_2<1$, thus the periodic solution $(0,0,\tilde{z}(t\left)\right)$ of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is locally asymptotically stable.

We next study the global asymptotical stability of the wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{z}(t\left)\right)$.

Theorem 3.2

If $R_{m3}<0$ and $R_{m4}<0$, the wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{z}(t\left)\right)$ of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is globally asymptotically stable.

Proof.

From system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), we have $\begin{array}{r}\frac{\mathrm{d}\left(x\right(t)+y(t\left)\right)}{\mathrm{d}t}\le \left(x\right(t)+y(t\left)\right)(\rho -\delta (x\left(t\right)+y\left(t\right)\left)\right),\end{array}$and $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& \le -z^2\left(t\right),t\ne nT,\\ z\left(t^{+}\right)& =z\left(t\right)+\text{β},t=nT,\\ z\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$According to the comparison theorem and Lemma 3.1, we obtain that for any $ϵ>0$, there always exists a $n_1\in Z_{+}$ such that (9) $\begin{array}{r}x\left(t\right)+y\left(t\right)\le W\left(t\right)<d+ϵ,z\left(t\right)<\tilde{z}\left(t\right)+ϵ\end{array}$(9) for $t\ge n_{1}T$, where $W\left(t\right)$ is the solution of the following system: $\begin{array}{r}\frac{\mathrm{d}W\left(t\right)}{\mathrm{d}t}=W\left(t\right)(\rho -\delta W(t\left)\right),W\left(0\right)=x\left(0\right)+y\left(0\right)>0.\end{array}$From the third equation of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), for any $t>n_{1}T$, we have (10) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& \ge -(d+ϵ)z\left(t\right)-z^2\left(t\right),t\ne nT,\\ z\left(t^{+}\right)& =z\left(t\right)+\text{β},t=nT.\end{array}\end{array}$(10) Then it follows from Lemma 3.2 and the comparison theorem of impulsive differential equation [28] that, if $ϵ>0$ small enough, there exists a positive integer $n_2>n_1$ such that (11) $\begin{array}{r}z\left(t\right)\ge g\left(t\right)>\tilde{g}\left(t\right)-ϵ\ge \tilde{g}\left(T\right)-ϵ,t\ge n_{2}T.\end{array}$(11) Substituting $\tilde{g}\left(t\right)-ϵ$ into the second equation of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), we have (12) $\begin{array}{r}\frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}\le y\left(t\right)\left(\rho \right(1-\theta )(1-\frac{2c_2\left(\tilde{g}\right(t)-ϵ)}{d+2\tilde{g}\left(t\right)-ϵ})-{\delta }_2(\tilde{g}\left(t\right)-ϵ\left)\right),t\ge n_{2}T.\end{array}$(12) If $\begin{array}{rl}{R}_{m4}& ={\int }_0^T\rho (1-\theta )(1-\frac{2c_2\tilde{g}\left(t\right)}{d+2\tilde{g}\left(t\right)})-{\delta }_2\tilde{g}\left(t\right)\mathrm{d}t\\ & =\rho (1-\theta )[T+\frac{2c_2}{d}\ln \frac{d+g^{\ast }(1+e^{-\mathrm{d}T})}{d+2g^{\ast }}]-{\delta }_2\ln [1+\frac{g^{\ast }(1-e^{-\mathrm{d}T})}{d}]<0,\end{array}$we can select a ϵ small enough such that $\begin{array}{r}\lambda ={\int }_0^T\rho (1-\theta )(1-\frac{2c_2\left(\tilde{g}\right(t)-ϵ)}{d+2\tilde{g}\left(t\right)-ϵ})-{\delta }_2\left(\tilde{g}\right(t)-ϵ)\mathrm{d}t<0.\end{array}$It follows from (12(12) $\begin{array}{r}\frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}\le y\left(t\right)\left(\rho \right(1-\theta )(1-\frac{2c_2\left(\tilde{g}\right(t)-ϵ)}{d+2\tilde{g}\left(t\right)-ϵ})-{\delta }_2(\tilde{g}\left(t\right)-ϵ\left)\right),t\ge n_{2}T.\end{array}$(12) ) that $\begin{array}{r}y\left(\right(n+1\left)T\right)\le e^{\lambda (n-n_2)T}y\left(n_{2}T\right),n>n_2,\end{array}$which implies that $y\left(t\right)\to 0$ as $t\to +\mathrm{\infty }$ ($n\to +\mathrm{\infty }$). Thus, for any $ϵ>0$ small enough, there must be a positive integer $n_3>n_2$ such that $y\left(t\right)<ϵ$ for $t>n_{3}T$. By the first equation of (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), we derive (13) $\begin{array}{r}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}\le \rho \theta ϵ+x\left(t\right)(1-{\delta }_{1}x(t\left)\right),t\ge n_{3}T,\end{array}$(13) then it follows from (13(13) $\begin{array}{r}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}\le \rho \theta ϵ+x\left(t\right)(1-{\delta }_{1}x(t\left)\right),t\ge n_{3}T,\end{array}$(13) ) that for sufficiently small ${ϵ}_1>0$, there is a $n_4>n_3$ such that $x\left(t\right)<\frac{1}{{\delta }_1}+\sqrt{\frac{\rho \theta ϵ}{{\delta }_1}}+{ϵ}_1$ for $t>n_{4}T$, and then we discuss the first equation of (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) again, we have (14) $\begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& \le \rho \theta ϵ+x\left(t\right)((1-\frac{2c_1{\delta }_1\left(\tilde{g}\right(t)-ϵ)}{1+\sqrt{{\delta }_1\rho \theta ϵ}+{\delta }_1{ϵ}_1+2{\delta }_1\tilde{g}\left(t\right)-{\delta }_1ϵ})-{\delta }_1(\tilde{g}\left(t\right)-ϵ\left)\right),\\ & t\ge n_{4}T.\end{array}$(14) If conditions $R_{m3}<0$ and $R_{m4}<0$ hold, we can deduce from (14(14) $\begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& \le \rho \theta ϵ+x\left(t\right)((1-\frac{2c_1{\delta }_1\left(\tilde{g}\right(t)-ϵ)}{1+\sqrt{{\delta }_1\rho \theta ϵ}+{\delta }_1{ϵ}_1+2{\delta }_1\tilde{g}\left(t\right)-{\delta }_1ϵ})-{\delta }_1(\tilde{g}\left(t\right)-ϵ\left)\right),\\ & t\ge n_{4}T.\end{array}$(14) ) that $x\left(t\right)\to 0$ as $t\to +\mathrm{\infty }$.

Substituting ϵ into the third equation of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) for $x\left(t\right)$ and $y\left(t\right)$, we have $\begin{array}{r}\{\begin{array}{rl}& \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}\ge -2ϵz\left(t\right)-z^2\left(t\right),t\ne nT,\\ & z\left(t^{+}\right)=z\left(t\right)+\text{β},t=nT.\end{array}\end{array}$From Lemma 3.2 and the comparison theorem of impulsive differential equation [28], we obtain that for sufficiently small $ϵ>0$, there is a positive integer $n_5(n_5>n_4)$ such that $z\left(t\right)>\tilde{z}\left(t\right)-ϵ$ for $t>n_{5}T$.

According to the above discussion, if the conditions of Theorem 3.2 are satisfied, for $ϵ>0$ small enough, we have $\begin{array}{r}0<x\left(t\right)<ϵ,0<y\left(t\right)<ϵ,and\tilde{z}\left(t\right)-ϵ<z\left(t\right)<\tilde{z}\left(t\right)+ϵ,t\ge n_{5}T.\end{array}$Letting $ϵ\to 0$, we obtain $x\left(t\right)\to 0$, $y\left(t\right)\to 0$ and $z\left(t\right)\to \tilde{z}\left(t\right)$ as $t\to +\mathrm{\infty }$, which means that the periodic solution $(0,0,\tilde{z}(t\left)\right)$ of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is a global attractor. Combining Theorem 3.1, we obtain that the periodic solution $(0,0,\tilde{z}(t\left)\right)$ is globally asymptotically stable when conditions $R_{m3}<0$ and $R_{m4}<0$ hold. The proof is completed.

Theorem 3.3

If $R_{m1}>0$ and $R_{m4}<0$, $x\left(t\right)$ is strongly persistent in the mean and $y\left(t\right)$ goes to extinction.

Proof.

From the proof of Theorem 3.2, we obtain that $\underset{t\to +\mathrm{\infty }}{lim}y\left(t\right)=0$ when $R_{m4}<0$. Further, for any $ϵ>0$, there exists $N_1>0$ such that (15) $\begin{array}{r}y\left(t\right)<0.5ϵ,t>N_{1}T.\end{array}$(15) By Lemmas 3.1 and 3.2, we can obtain that for $ϵ>0$ sufficiently small, there exists a $N_2>N_1$ such that (16) $\begin{array}{r}\tilde{g}\left(T\right)-0.5ϵ\le \tilde{g}\left(t\right)-0.5ϵ\le z\left(t\right)\le \tilde{z}\left(t\right)+0.5ϵ\le z^{\ast }+0.5ϵ,t>N_{2}T.\end{array}$(16) Substituting inequalities (15(15) $\begin{array}{r}y\left(t\right)<0.5ϵ,t>N_{1}T.\end{array}$(15) ) and (16(16) $\begin{array}{r}\tilde{g}\left(T\right)-0.5ϵ\le \tilde{g}\left(t\right)-0.5ϵ\le z\left(t\right)\le \tilde{z}\left(t\right)+0.5ϵ\le z^{\ast }+0.5ϵ,t>N_{2}T.\end{array}$(16) ) into the first equation of (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), we have (17) $\begin{array}{r}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}\ge \left[\right((1-c_1)-{\delta }_1\left(\tilde{z}\right(t)+ϵ))-{\delta }_{1}x(t\left)\right]x\left(t\right),t>N_{2}T.\end{array}$(17) Dividing (17(17) $\begin{array}{r}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}\ge \left[\right((1-c_1)-{\delta }_1\left(\tilde{z}\right(t)+ϵ))-{\delta }_{1}x(t\left)\right]x\left(t\right),t>N_{2}T.\end{array}$(17) ) by $x\left(t\right)$ and then integrating both sides of the resulted equation on the interval $[N_{2}T,t]$, we can obtain (18) $\begin{array}{r}\frac{\ln \left(x\right(t\left)/x\right(N_{2}T\left)\right)}{t}\ge (1-c_1)-{\delta }_1{t}^{-1}{\int }_{N_{2}T}^t\left(\tilde{z}\right(p)+ϵ)\mathrm{d}p-{\delta }_1{t}^{-1}{\int }_{N_{2}T}^{t}x\left(p\right)\mathrm{d}p,t>N_{2}T.\end{array}$(18) By (9(9) $\begin{array}{r}x\left(t\right)+y\left(t\right)\le W\left(t\right)<d+ϵ,z\left(t\right)<\tilde{z}\left(t\right)+ϵ\end{array}$(9) ), we have $\underset{t\to +\mathrm{\infty }}{lim sup}\frac{\ln \left(x\right(t\left)/x\right(N_{2}T\left)\right)}{t}\le 0$, and $\begin{array}{rl}\frac{(n-N_2)}{(n+1)T}{\int }_0^T\tilde{z}\left(p\right)\mathrm{d}p& \le t^{-1}{\int }_{N_{2}T}^t\tilde{z}\left(p\right)\mathrm{d}p<\frac{(n+1-N_2)}{nT}{\int }_0^T\tilde{z}\left(p\right)\mathrm{d}p,\\ & t\in [nT,(n+1\left)T\right),n>N_2.\end{array}$Then from (18(18) $\begin{array}{r}\frac{\ln \left(x\right(t\left)/x\right(N_{2}T\left)\right)}{t}\ge (1-c_1)-{\delta }_1{t}^{-1}{\int }_{N_{2}T}^t\left(\tilde{z}\right(p)+ϵ)\mathrm{d}p-{\delta }_1{t}^{-1}{\int }_{N_{2}T}^{t}x\left(p\right)\mathrm{d}p,t>N_{2}T.\end{array}$(18) ) and $R_{m1}>0$, we have $\underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^{t}x\left(p\right)\mathrm{d}p\ge \frac{(1-c_1)-{\delta }_1{T}^{-1}{\int }_0^T\tilde{z}\left(p\right)\mathrm{d}p}{{\delta }_1}=m>0.$ The proof is completed.

Theorem 3.4

If $R_{m5}>0$, we have $\begin{array}{r}\underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^t\left(x\right(p)+y(p\left)\right)\mathrm{d}p\ge \frac{min\{1-c_1,\rho (1-c_2\left)\right\}-{\delta }_0{T}^{-1}{\int }_0^T\tilde{z}\left(t\right)\mathrm{d}t}{{\delta }_0},\end{array}$that is, the wild $N.lugens$ is strongly persistent in the mean.

Proof.

From the first and second equations of (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) and (18(18) $\begin{array}{r}\frac{\ln \left(x\right(t\left)/x\right(N_{2}T\left)\right)}{t}\ge (1-c_1)-{\delta }_1{t}^{-1}{\int }_{N_{2}T}^t\left(\tilde{z}\right(p)+ϵ)\mathrm{d}p-{\delta }_1{t}^{-1}{\int }_{N_{2}T}^{t}x\left(p\right)\mathrm{d}p,t>N_{2}T.\end{array}$(18) ), we have $\begin{array}{r}\frac{\mathrm{d}\ln \left(x\right(t)+y(t\left)\right)}{\mathrm{d}t}\ge min\{1-c_1,\rho (1-c_2\left)\right\}-{\delta }_0\left(\tilde{z}\right(t)+ϵ)-{\delta }_0\left(x\right(t)+y(t\left)\right),\end{array}$for $t>N_{2}T$. Similar to the proof of Theorem 3.3, we obtain that $\begin{array}{r}\underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^t\left(x\right(p)+y(p\left)\right)\mathrm{d}p\ge \frac{min\{1-c_1,\rho (1-c_2\left)\right\}-{\delta }_0{T}^{-1}{\int }_0^T\tilde{z}\left(t\right)\mathrm{d}t}{{\delta }_0}>0.\end{array}$The proof is completed.

The uninfected $N.lugens$ and the wLug-infected $N.lugens$ are two types of natural populations. Since the wild $N.lugens$ infected with wLug is imperfectly maternally transmitted, it is unlikely that only the wild $N.lugens$ infected with wLug occur in the agricultural ecosystem (see Figure 1), but it is possible that only uninfected wild $N.lugens$ populations exist in the agricultural ecosystem. Next, we will show the dynamics of the following subsystem which absents the $N.lugens$ infected with wLug. (19) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x\left(t\right)(1-\frac{2c_{1}z\left(t\right)}{x\left(t\right)+2z\left(t\right)})-{\delta }_1\left(x\right(t)+z(t\left)\right)x\left(t\right),\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-\left(x\right(t)+z(t\left)\right)z\left(t\right),\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(19) From Theorems 3.1–3.4, we can give the following results for system (19(19) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x\left(t\right)(1-\frac{2c_{1}z\left(t\right)}{x\left(t\right)+2z\left(t\right)})-{\delta }_1\left(x\right(t)+z(t\left)\right)x\left(t\right),\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-\left(x\right(t)+z(t\left)\right)z\left(t\right),\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(19) ).

Corollary 3.1

(a)	The wild-$N.lugens$ -eliminate periodic solution $(0,\tilde{z}(t\left)\right)$ of system (19(19) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x\left(t\right)(1-\frac{2c_{1}z\left(t\right)}{x\left(t\right)+2z\left(t\right)})-{\delta }_1\left(x\right(t)+z(t\left)\right)x\left(t\right),\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-\left(x\right(t)+z(t\left)\right)z\left(t\right),\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(19) ) is locally asymptotically stable if $R_{m1}<0$.
(b)	The wild-$N.lugens$ -eliminate periodic solution $(0,\tilde{z}(t\left)\right)$ of system (19(19) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x\left(t\right)(1-\frac{2c_{1}z\left(t\right)}{x\left(t\right)+2z\left(t\right)})-{\delta }_1\left(x\right(t)+z(t\left)\right)x\left(t\right),\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-\left(x\right(t)+z(t\left)\right)z\left(t\right),\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(19) ) is globally asymptotically stable if $\begin{array}{r}T+2c_1{\delta }_1\ln \frac{1+{\delta }_1{g}_1^{\ast }(1+e^{-\frac{1}{{\delta }_1}T})}{1+2{\delta }_1{g}_1^{\ast }}-{\delta }_1\ln [1+{\delta }_1{g}_1^{\ast }(1-e^{-\frac{1}{{\delta }_1}T}\left)\right]<0,\end{array}$where $g_1^{\ast }$ is the positive root of the following equation $\begin{array}{r}{g}_1^2+(\frac{1}{{\delta }_1}-\text{β})g_1-\frac{\text{β}}{{\delta }_1(1-e^{-\frac{1}{{\delta }_1}T})}=0.\end{array}$
(c)	System (19(19) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x\left(t\right)(1-\frac{2c_{1}z\left(t\right)}{x\left(t\right)+2z\left(t\right)})-{\delta }_1\left(x\right(t)+z(t\left)\right)x\left(t\right),\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-\left(x\right(t)+z(t\left)\right)z\left(t\right),\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(19) ) is strongly persistent in the mean if $R_{m1}>0$.

3.2. The dynamics of model (6)

Denote $\begin{array}{rl}{R}_{f1}& =\frac{1-c_1}{{\delta }_1}-\frac{1}{T}\ln (1+Z^{\ast }(e^T-1\left)\right),\\ R_{f2}& =\frac{\rho (1-c_2)(1-\theta )}{{\delta }_2}-\frac{1}{T}\ln (1+Z^{\ast }(e^T-1\left)\right),\\ R_{f3}& =T-\frac{c_1{\delta }_1{h}^{\ast }}{(1-d)-h^{\ast }}(e^{(1-d)T}-1)-{\delta }_1\ln (1+h^{\ast }\frac{e^{(1-d)T}-1}{1-d}),\\ R_{f4}& =\rho (1-\theta )(T-c_2\ln (\frac{d}{d+h_1^{\ast }}+\frac{\mathrm{d}{h}^{\ast }}{d+h^{\ast }}\frac{e^{(1-d)T}-1}{1-d}+\frac{h^{\ast }{e}^{(1-d)T}}{d+h^{\ast }}))\\ & -{\delta }_2\ln (1+h^{\ast }\frac{e^{(1-d)T}-1}{1-d}),\\ R_{f5}& =T-\frac{c_1{\delta }_1}{1+{\delta }_1(1-d)}\ln (\frac{1}{{\delta }_1{h}^{\ast }+1}+\frac{h^{\ast }}{{\delta }_1{h}^{\ast }+1}\frac{e^{(1-d)T}-1}{1-d}+\frac{{\delta }_1{h}^{\ast }{e}^{(1-d)T}}{1+{\delta }_1{h}^{\ast }})\\ & -{\delta }_1\ln (1+h^{\ast }\frac{e^{(1-d)T}-1}{1-d}),\\ R_{f6}& =T-c_1{\delta }_1\ln (1+\frac{h_1^{\ast }T}{1+{\delta }_1{h}_1^{\ast }})-{\delta }_1\ln (1+h_1^{\ast }T),\\ R_{f7}& =\rho (1-\theta )(T-c_2\ln (1+\frac{h_1^{\ast }T}{1+h_1^{\ast }}))-{\delta }_2\ln (1+h_1^{\ast }T),\\ R_{f8}& =min\{1-c_1,\rho (1-c_2\left)\right\}-{\delta }_0\frac{1}{T}\ln (1+Z^{\ast }(e^T-1\left)\right),\end{array}$where $Z^{\ast }$, $h^{\ast }$ and $h_1^{\ast }$ are given in the following Lemmas 3.3 and 3.4.

To begin with, we consider the following subsystem of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) when x=0 and y=0: (20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20)

Lemma 3.3

System (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ) has a unique positive periodic solution $\tilde{Z}\left(t\right)$ with period T, and for any solution $Z\left(t\right)$ of system (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ) with $Z\left(0^{+}\right)\ge 0$, we have $|Z\left(t\right)-\tilde{Z}\left(t\right)|\to 0$ as $t\to +\mathrm{\infty }$, where $\begin{array}{rl}\tilde{Z}\left(t\right)& =\frac{Z^{\ast }{e}^{(t-nT)}}{1+Z^{\ast }(e^{(t-nT)}-1)},t\in (nT,(n+1\left)T\right],\\ \tilde{Z}\left(0^{+}\right)& =Z^{\ast }>0,\end{array}$and $Z^{\ast }$ is the positive root of the following equation: $\begin{array}{r}{Z}^2-(\text{β}+1)Z-\frac{\text{β}}{e^T-1}=0.\end{array}$

Proof.

Integrating the first equation in (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ) between pulses, we have $\begin{array}{r}Z\left(t\right)=\frac{Z\left(n{T}^{+}\right)e^{(t-nT)}}{1+Z\left(n{T}^{+}\right)(e^{(t-nT)}-1)}\end{array}$for $t\in (nT,(n+1\left)T\right]$.

By the second equation of (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ), we can obtain the stroboscopic map: (21) $\begin{array}{r}\begin{array}{rl}Z\left(\right(n+1\left)T^{+}\right)& =Z\left(\right(n+1\left)T\right)+\text{β}\\ & =\frac{Z\left(n{T}^{+}\right)e^T}{1+Z\left(n{T}^{+}\right)(e^T-1)}+\text{β}=h\left(Z\right(n{T}^{+}\left)\right).\end{array}\end{array}$(21) Consider the equation $Z=h\left(Z\right)$, namely, $\begin{array}{r}\frac{Z{e}^T}{1+Z(e^T-1)}+\text{β}=Z,\end{array}$which is equivalent to the standardized quadratic equation: $\begin{array}{r}{Z}^2-(\text{β}+1)Z-\frac{\text{β}}{e^T-1}=0.\end{array}$Obviously, system (21(21) $\begin{array}{r}\begin{array}{rl}Z\left(\right(n+1\left)T^{+}\right)& =Z\left(\right(n+1\left)T\right)+\text{β}\\ & =\frac{Z\left(n{T}^{+}\right)e^T}{1+Z\left(n{T}^{+}\right)(e^T-1)}+\text{β}=h\left(Z\right(n{T}^{+}\left)\right).\end{array}\end{array}$(21) ) has unique positive fixed point $Z^{\ast }$. Hence, system (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ) has unique positive periodic solution $\begin{array}{r}\tilde{Z}\left(t\right)=\frac{Z^{\ast }{e}^{(t-nT)}}{1+Z^{\ast }(e^{(t-nT)}-1)}\end{array}$with $\tilde{Z}\left(0^{+}\right)=Z^{\ast }$.

Similar to Lemma 3.1, we can also get that $Z^{\ast }$ is globally asymptotically stable for system (21(21) $\begin{array}{r}\begin{array}{rl}Z\left(\right(n+1\left)T^{+}\right)& =Z\left(\right(n+1\left)T\right)+\text{β}\\ & =\frac{Z\left(n{T}^{+}\right)e^T}{1+Z\left(n{T}^{+}\right)(e^T-1)}+\text{β}=h\left(Z\right(n{T}^{+}\left)\right).\end{array}\end{array}$(21) ), then the corresponding period solution $\tilde{Z}\left(t\right)$ of system (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ) is also globally asymptotically stable. This completes the proof.

Remark 3.1

From the Lemma 3.3, we can easily calculate that $Z^{\ast }\ge 1$, thus there is a special case that $Z^{\ast }=1$ when $\text{β}=0$ and $Z\left(0^{+}\right)>0$, and for any solution $Z\left(t\right)$ of system (20(20) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)-Z^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }Z\left(t\right)& =\text{β},t=nT,\\ Z\left(0^{+}\right)& =Z_0\ge 0.\end{array}\end{array}$(20) ), we have $Z\left(t\right)\to 1$ as $t\to +\mathrm{\infty }$.

For the following system: (22) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}h\left(t\right)}{\mathrm{d}t}& =(1-d)h\left(t\right)-h^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }h\left(t\right)& =\text{β},t=nT,\\ h\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(22) Similar to Lemma 3.3, we obtain the following results:

Lemma 3.4

System (22(22) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}h\left(t\right)}{\mathrm{d}t}& =(1-d)h\left(t\right)-h^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }h\left(t\right)& =\text{β},t=nT,\\ h\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(22) ) has a unique periodic solution $\tilde{h}\left(t\right)$ with period T, and for any solution $h\left(t\right)$ of system (22(22) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}h\left(t\right)}{\mathrm{d}t}& =(1-d)h\left(t\right)-h^2\left(t\right),t\ne nT,\\ \mathrm{\Delta }h\left(t\right)& =\text{β},t=nT,\\ h\left(0^{+}\right)& =z_0\ge 0.\end{array}\end{array}$(22) ) with $h\left(0^{+}\right)\ge 0$, we have $|h\left(t\right)-\tilde{h}\left(t\right)|\to 0$ as $t\to +\mathrm{\infty }$.

When $d\ne 1$, the periodic solution $\tilde{h}\left(t\right)$ is given by $\begin{array}{rl}\tilde{h}\left(t\right)& =\frac{h^{\ast }(1-d)e^{(1-d)(t-nT)}}{(1-d)+h^{\ast }(e^{(1-d)(t-nT)}-1)},t\in (nT,(n+1\left)T\right],\\ \tilde{h}\left(0^{+}\right)& =h^{\ast }>0,\end{array}$where $h^{\ast }$ is the positive root of equation $h^2-\left(\right(1-d)+\text{β})h-\frac{\text{β}(1-d)}{e^{(1-d)T}-1}=0.$

When d=1, the periodic solution $\tilde{h}\left(t\right)$ is given by $\begin{array}{r}\tilde{h}\left(t\right)=\frac{h_1^{\ast }}{1+h_1^{\ast }(t-nT)},t\in (nT,(n+1\left)T\right],\end{array}$and $\tilde{h}\left(0^{+}\right)=h_1^{\ast }=\frac{\text{β}+\sqrt{{\text{β}}^2+\frac{4\text{β}}{T}}}{2}>0.$

Similar to the discussion of Theorems 3.1–3.4, we have the following results:

Theorem 3.5

If $R_{f1}<0$ and $R_{f2}<0$, the wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{Z}(t\left)\right)$ of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) is locally asymptotically stable.

Proof.

Define $u\left(t\right)=X\left(t\right)$, $v\left(t\right)=Y\left(t\right)$, $w\left(t\right)=Z\left(t\right)-\tilde{Z}\left(t\right)$. Then the linearized system of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) in $(nT,(n+1\left)T\right]$ is obtained as follows: $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}u\left(t\right)}{\mathrm{d}t}& =\left(\right(1-c_1)-{\delta }_1\tilde{Z}(t\left)\right)u\left(t\right)+\rho \theta (1-c_2)v\left(t\right),\\ \frac{\mathrm{d}v\left(t\right)}{\mathrm{d}t}& =\left(\rho \right(1-\theta \left)\right(1-c_2)-{\delta }_2\tilde{Z}(t\left)\right)v\left(t\right),\\ \frac{\mathrm{d}w\left(t\right)}{\mathrm{d}t}& =-\tilde{Z}\left(t\right)u\left(t\right)-\tilde{Z}\left(t\right)v\left(t\right)+(1-2\tilde{Z}(t\left)\right)w\left(t\right).\end{array}\end{array}$Clearly, the fundamental solution matrix in interval $(nT,(n+1\left)T\right]$ is $\begin{array}{r}\Phi \left(t\right)=\left(\begin{array}{ccc}{e}^{{\int }_{nT}^t(1-c_1)-{\delta }_1\tilde{Z}\left(p\right)\mathrm{d}p}& {†}_1& 0\\ 0& e^{{\int }_{nT}^t\rho (1-\theta )(1-c_2)-{\delta }_2\tilde{Z}\left(p\right)\mathrm{d}p}& 0\\ {†}_2& {†}_3& e^{{\int }_{nT}^{t}1-2\tilde{Z}\left(p\right)\mathrm{d}p}\end{array}\right),\end{array}$and the expression of function ${†}_i$ (i=1,2,3) does not need to be given.

From the fourth equations of (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ), we have $\begin{array}{r}\left(\begin{array}{c}u\left(\right(n+1\left)T^{+}\right)\\ v\left(\right(n+1\left)T^{+}\right)\\ w\left(\right(n+1\left)T^{+}\right)\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}u\left(\right(n+1\left)T\right)\\ v\left(\right(n+1\left)T\right)\\ w\left(\right(n+1\left)T\right)\end{array}\right).\end{array}$The local stability of the periodic solution $(0,0,\tilde{Z}(t\left)\right)$ is determined by the eigenvalues of $\Phi \left(T\right)$, where they are $\begin{array}{r}{\lambda }_1=e^{{\int }_0^T(1-c_1)-{\delta }_1\tilde{Z}\left(t\right)\mathrm{d}t},{\lambda }_2=e^{{\int }_0^T\rho (1-\theta )(1-c_2)-{\delta }_2\tilde{Z}\left(t\right)\mathrm{d}t},\end{array}$and $\begin{array}{r}{\lambda }_3=e^{{\int }_0^{T}1-2\tilde{Z}\left(t\right)\mathrm{d}t}=\frac{\tilde{Z}\left(T\right)}{Z^{\ast }}{e}^{-{\int }_0^T\tilde{Z}\left(t\right)\mathrm{d}t}=\frac{Z^{\ast }-\text{β}}{Z^{\ast }}{e}^{-{\int }_0^T\tilde{Z}\left(t\right)\mathrm{d}t}<1.\end{array}$According to the Floquet theorem and conditions given in Theorem 3.5, we have ${\lambda }_1<1$ and ${\lambda }_2<1$. Hence, the periodic solution $(0,0,\tilde{Z}(t\left)\right)$ of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) is locally asymptotically stable.

Theorem 3.6

If any of the following conditions is true,

(i)	$1-d\ne 0,(1-d){\delta }_1+1=0,$ $R_{f3}<0,$ and $R_{f4}<0$.
(ii)	$1-d\ne 0,(1-d){\delta }_1+1\ne 0,$ $R_{f5}<0$ and $R_{f4}<0$.
(iii)	1−d=0, $R_{f6}<0$ and $R_{f7}<0$.

Then the wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{Z}(t\left)\right)$ of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) is globally asymptotically stable.

Theorem 3.7

If any of the following conditions is true,

(i)	$1-d\ne 0$, $R_{f4}<0$, and $R_{f1}>0$.
(ii)	1−d=0, $R_{f7}<0$, and $R_{f1}>0$.

Then $X\left(t\right)$ is strongly persistent in the mean and $Y\left(t\right)$ goes to extinction.

Theorem 3.8

If $R_{f8}>0$, we have $\begin{array}{r}\underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^t\left(X\right(p)+Y(p\left)\right)\mathrm{d}p\ge \frac{min\{1-c_1,\rho \theta (1-c_2\left)\right\}-{\delta }_0{T}^{-1}{\int }_0^T\tilde{Z}\left(t\right)\mathrm{d}t}{{\delta }_0}.\end{array}$

The proofs of Theorem 3.6–3.8 are basically similar to those of Theorems 3.2–3.4, therefore, we omit them here.

In the following, we consider the subsystem without wLug of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ). (23) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =X\left(t\right)\frac{X\left(t\right)+(1-c_1)Z\left(t\right)}{X\left(t\right)+Z\left(t\right)}-{\delta }_1\left(X\right(t)+Z(t\left)\right)X\left(t\right),\\ \frac{\mathrm{d}\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)(1-(X\left(t\right)+Z\left(t\right)\left)\right),\end{array}\}t\ne nT.\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(23)

From Theorem 3.5, we can obtain the following results for system (23(23) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =X\left(t\right)\frac{X\left(t\right)+(1-c_1)Z\left(t\right)}{X\left(t\right)+Z\left(t\right)}-{\delta }_1\left(X\right(t)+Z(t\left)\right)X\left(t\right),\\ \frac{\mathrm{d}\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)(1-(X\left(t\right)+Z\left(t\right)\left)\right),\end{array}\}t\ne nT.\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(23) ).

Corollary 3.2

(a)	The wild-$N.lugens$ -eliminate periodic solution $(0,\tilde{Z}(t\left)\right)$ of system (23(23) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =X\left(t\right)\frac{X\left(t\right)+(1-c_1)Z\left(t\right)}{X\left(t\right)+Z\left(t\right)}-{\delta }_1\left(X\right(t)+Z(t\left)\right)X\left(t\right),\\ \frac{\mathrm{d}\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)(1-(X\left(t\right)+Z\left(t\right)\left)\right),\end{array}\}t\ne nT.\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(23) ) is locally asymptotically stable if $R_{f1}<0$.
(b)	If any of the following conditions is true, (i) ${\delta }_1=1$ and $T-c_1\ln (1+\frac{h_1^{\ast }T}{1+h_1^{\ast }})-\ln (1+h_1^{\ast }T)<0$. (ii) ${\delta }_1\ne 1$ and $T-c_1\ln (\frac{1}{{\delta }_1{h}_2^{\ast }+1}+\frac{h_2^{\ast }}{{\delta }_1{h}_2^{\ast }+1}\frac{e^{(1-\frac{1}{{\delta }_1})T}-1}{1-\frac{1}{{\delta }_1}}+\frac{{\delta }_1{h}_2^{\ast }{e}^{(1-\frac{1}{{\delta }_1})T}}{1+{\delta }_1{h}_2^{\ast }})-{\delta }_1\ln (1+h_2^{\ast }\frac{e^{(1-\frac{1}{{\delta }_1})T}-1}{1-\frac{1}{{\delta }_1}})<0$, where $h_2^{\ast }$ is the solution of equation $h^2-(1-\frac{1}{{\delta }_1}+\text{β})h-\frac{\text{β}(1-\frac{1}{{\delta }_1})}{e^{(1-\frac{1}{{\delta }_1})T}-1}=0$. The wild-$N.lugens$ -eliminate periodic solution $(0,\tilde{Z}(t\left)\right)$ of system (23(23) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =X\left(t\right)\frac{X\left(t\right)+(1-c_1)Z\left(t\right)}{X\left(t\right)+Z\left(t\right)}-{\delta }_1\left(X\right(t)+Z(t\left)\right)X\left(t\right),\\ \frac{\mathrm{d}\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)(1-(X\left(t\right)+Z\left(t\right)\left)\right),\end{array}\}t\ne nT.\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(23) ) is globally asymptotically stable.
(c)	System (23(23) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =X\left(t\right)\frac{X\left(t\right)+(1-c_1)Z\left(t\right)}{X\left(t\right)+Z\left(t\right)}-{\delta }_1\left(X\right(t)+Z(t\left)\right)X\left(t\right),\\ \frac{\mathrm{d}\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)(1-(X\left(t\right)+Z\left(t\right)\left)\right),\end{array}\}t\ne nT.\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(23) ) is strongly persistent in the mean if $R_{f1}>0$.

4. Numerical simulation and discussions

In this section, we will verify our results by numerical simulation. From Table 1, we obtain that $\theta =0.85$, $c_1=0.675$, $c_2=0.732$, and $\rho =2$ by taking a=45 and $\rho a=90$ at room temperature ${25}^{\circ }$C. We choose the parameters ${\delta }_1=0.85$ and ${\delta }_2=0.2$. If we do not release the $N.lugens$ infected with wStri, the ecosystem will be transformed into the following model (24) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x+\rho \theta y-{\delta }_1(x+y)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y-{\delta }_2(x+y)y.\end{array}\end{array}$(24) According to Theorem 10 in [10], system (24(24) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x+\rho \theta y-{\delta }_1(x+y)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y-{\delta }_2(x+y)y.\end{array}\end{array}$(24) ) has unique a positive equilibrium point $(1.2911,0.2089)$, and it is globally asymptotically stable. The time series of the wild uninfected $N.lugens$ and the wild $N.lugens$ infected with wLug are shown in Figure 2. In the following, we will discuss the influences of the release rate β and release period T of the $N.lugens$ infected with wStri on the dynamics of systems (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) and (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ).

4.1. Long-time behaviours of population suppression model

In order to evaluate the influences of the release period T and the release amount β of male $N.lugens$ infected with wStri on the dynamics of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ). We first take T=3 and $\text{β}=1$, by calculating, then $R_{m1}=-0.1399<0$ and $R_{m2}=-0.1203<0$, which satisfy the conditions of Theorem 3.1, the wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{z}(t\left)\right)$ of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) is locally asymptotically stable. System (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) may have two steady states coexisting, the wild-$N.lugens$ -eradication period solution where the wild $N.lugens$ population will be extinct, and a positive periodic solution where the $N.lugens$ population oscillates positively and periodically as shown in Figure 3. These results also indicate that when the number of the wild $N.lugens$ population is small, the wild $N.lugens$ population will be controlled by regularly releasing fewer male $N.lugens$ populations infected with wStri into the field.

Figure 2. The solution of system (24(24) $\begin{array}{r}\{\begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x+\rho \theta y-{\delta }_1(x+y)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y-{\delta }_2(x+y)y.\end{array}\end{array}$(24) ) with $\theta =0.85$, $c_1=0.675$, $c_2=0.732$, $\rho =2$, ${\delta }_1=0.85$ and ${\delta }_2=0.2$.

Figure 3. The solutions of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) with T=3 and $\text{β}=1$ for different initial values. Here we take initial values (0.1, 0.01, 1), (0.08, 0.015, 1), (0.4, 0.08, 0) and (0.8, 0.15, 0).

Table 1. Parameter values of the $N.lugens$ in the laboratory.

Parameter	Definition	Lab.	source
a	Number of eggs laid by a female wild $N.lugens$ uninfected with wLug ($da{y}^{-1}$)	(37,55)	[8]
$\rho a$	Number of eggs laid by a female wild $N.lugens$ infected with wLug ($da{y}^{-1}$)	(70,100)	[8]
θ	The maternal transmission leakage rate of wLug	0.85	[6]
$c_1$	The CI intensity of $S_m$ against $U_f$	0.675	[6]
$c_1$	The CI intensity of $S_m$ against $L_f$	0.732	[6]

However, when we take T=0.5 and $\text{β}=11$, then $R_{m3}=-0.2858<0$ and $R_{m4}=-0.0036<0$, from Theorem 3.2, we can get that the periodic solution $(0,0,\tilde{z}(t\left)\right)$ is globally asymptotically stable. It means that whatever the initial value is, the wild $N.lugens$ population goes to extinction for T=0.5 and $\text{β}=11$, see Figure 4. If we take T=6 and $\text{β}=1$, according to Theorem 3.4, we have $R_{m5}=0.1961>0$, and the wild $N.lugens$ population is strongly persistent in the mean as shown in Figure 5. From Figures 2–5, we can also see that increasing the release amount β of male $N.lugens$ infected with wStri or decreasing the release period T are beneficial to suppress the density of wild $N.lugens$ population.

4.2. Long-time behaviours of population replacement model

We first take $\text{β}=0.1$ and T=6, by calculating, we obtain $R_{f1}=-0.6335<0$ and $R_{f2}=-0.6139<0$, there exists a locally asymptotically stable wild $N.lugens$ eradication periodic solution for system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) where the wild $N.lugens$ population goes to extinction. Furthermore, by selecting appropriate initial values, the numerical simulations show that system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) also exists a locally asymptotically stable positive periodic solution (see Figure 6).

In addition, we consider a special case of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) that $\text{β}=0$ and $Z\left(0^{+}\right)>0$. We obtain $R_{f1}=-0.6176<0$ and $R_{f2}=-0.5980<0$ by Theorem 3.5, and the equilibrium infected with wStri (0,0,1) is locally asymptotically stable as shown in Figure 7. It can be seen that we only need to release the $N.lugens$ infected with wStri once to the field, and it is possible to complete the replacement of wild $N.lugens$ populations.

Figure 4. The solutions of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) with T=0.5 and $\text{β}=11$ for different initial values. Here we take initial values (0.6, 0.1, 1), (0.8, 0.15, 0.2), (1, 0.2, 0), (1.4, 0.4, 0) and (1.6, 0.5, 0.5).

Figure 5. The solution of system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) with T=6, $\text{β}=1$ and the initial value (0.1, 0.05, 1).

Figure 6. The solutions of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) with T=6 and $\text{β}=0.1$ for different initial values. Here we take initial values (0.1,0.02,0.9), (0.2, 0.05, 0.7), (0.3, 0.08, 1), (0.6, 0.1, 0.1), (0.9, 0.15, 0.2) and (1.2, 0.2, 0.4).

Figure 7. The solutions of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) with $\text{β}=0$ for different initial values. Here we take initial values (0.1,0.02,0.9), (0.2, 0.05, 0.7), (0.3, 0.08, 1), (0.6, 0.1, 0.1), (0.9, 0.15, 0.2) and (1.2, 0.2, 0.4).

And take $\text{β}=10$ and T=0.5, we can get 1−d=−9, $(1-d){\delta }_1+1=-6.65$, $R_{f4}=-0.0106<0$ and $R_{f5}=-0.2591<0$. It follows from Theorem 3.6 that there exists a globally asymptotically stable wild-$N.lugens$ -eradication periodic solution $(0,0,\tilde{Z}(t\left)\right)$ as shown in Figure 8. Figure 8 also shows that whatever the initial value is, the wild $N.lugens$ population will be replaced by the $N.lugens$ infected with wStri when $\text{β}=10$ and T=0.5.

Figure 8. The solutions of system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) with T=0.5 and $\text{β}=10$ for different initial values. Here we take initial values (0.6,0.1,0), (0.8, 0.2, 0), (1.2, 0.3, 0.1), (1.6, 0.5, 0.1) and (2, 0.8, 0.2).

4.3. Control of the wild $N.lugens$ within a short time

In the first two sub-sections, we fixed parameters ρ, θ, $c_1$, $c_2$, ${\delta }_1$ and ${\delta }_2$, and then showed the long-time behaviours of population suppression model and population replacement model by changing parameters β and T. However, it is very important to control the wild $N.lugens$ population to a low level in a short time in the actual paddy field management. In this subsection, we will study the effects of release amount β and release period T on the control efficiency of wild $N.lugens$ population within a short time under different release strategies. Therefore, we define a concept of control degree as the ratio of the amount of wild $N.lugens$ population suppressed or replaced with a finite time to the initial wild $N.lugens$ amount [29], denoted by e, and the calculation formula is shown by $\begin{array}{r}e\left(t\right)=\frac{w\left(0\right)-w\left(t\right)}{w\left(0\right)}\times 100\mathrm{\%},\end{array}$where $w\left(t\right)$ is the amount of wild $N.lugens$ population at time t, and $w\left(0\right)$ denotes the initial amount of wild $N.lugens$ population. We assume that $x\left(0\right)=1$, $y\left(0\right)=0.3$, $z\left(0^{+}\right)=\text{β}$ and $T_0=20$.

We first consider the effect of the release period T on the control efficiency of wild $N.lugens$ in a finite time $T_0$. Fix parameters $c_1=0.675$, $c_2=0.732$, ${\delta }_1=0.85$, ${\delta }_2=0.2$, $\text{β}=1.5$ and vary T. We take T=1, 1.25, 2 and 2.5, respectively, Both system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) and system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ) have stable wild-$N.lugens$ -eradication boundary periodic solution. From Figure 9, we can see that the control efficiency of wild $N.lugens$ decreases as T increases. For system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), Figure 9 shows that when $T\le 1.25$, the control efficiency of the wild $N.lugens$ reaches more than $90\mathrm{\%}$ within time $T_0$, and when T=2.5, the control efficiency of the uninfected wild $N.lugens$ and the wild $N.lugens$ infected with wLug within time $T_0$ is only $29.29\mathrm{\%}$ and $55.10\mathrm{\%}$, respectively. However, the control efficiency of the wild $N.lugens$ within time $T_0$ is more than $95\mathrm{\%}$ for system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ).

Figure 9. The control degree of the wild $N.lugens$ within a finite time for different release periods. The first four lines of the legend represent the control efficiency of wild $N.lugens$ in system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), and the last four lines of the legend represent the control efficiency of wild $N.lugens$ in system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ). (a) The control degree of the uninfected wild $N.lugens$. (b) The control degree of the wild $N.lugens$ infected with wLug.

And then, we fix parameters $c_1=0.675$, $c_2=0.732$, ${\delta }_1=0.85$, ${\delta }_2=0.2$, T=2 and vary β. Let $\text{β}=1.8,2,2.5$ and 3, there is also a stable wild-$N.lugens$ -eradication boundary periodic solution for both system. Figure 10 shows that the control efficiency of wild $N.lugens$ increases as β increases. If we adopt a population suppression strategy, the control efficiency of uninfected wild $N.lugens$ within time $T_0$ is more than $93\mathrm{\%}$ when $\text{β}\ge 2$, and the control efficiency of wild $N.lugens$ infected with wLug within time $T_0$ is more than $90\mathrm{\%}$ when $\text{β}\ge 3$. However, when we adopt the population replacement strategy, the control efficiency of wild $N.lugens$ within time $T_0$ is more than $98\mathrm{\%}$ when $\text{β}\ge 1.8$.

Figure 10. The control degree of the wild $N.lugens$ within a finite time for different release amounts. The first four lines of the legend represent the control efficiency of wild $N.lugens$ in system (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ), and the last four lines of the legend represent the control efficiency of wild $N.lugens$ in system (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ). (a) The control degree of the uninfected wild $N.lugens$. (b) The control degree of the wild $N.lugens$ infected with wLug.

In addition, we consider an ecosystem in which the population of $N.lugens$ infected with wLug is absent, fix $c_1=0.675$, ${\delta }_1=0.85$, $x\left(0\right)=1$, $T_0=20$, and the values of other parameters are the same as those in Figures 9 and 10, then the time series of the control degree of this wild $N.lugens$ population under different release strategies are shown in Figure 11. It can be seen from Figure 11 that no matter which strategy is adopted, the wild $N.lugens$ population will be controlled to more than $99\mathrm{\%}$ in $T_0$ time. Comparing the results in Figures 9–11, the wild $N.lugens$ infected with wLug can reduce the control efficiency of uninfected wild $N.lugens$.

Figure 11. The control degree of the wild $N.lugens$ within a finite time. The first four lines of the legend represent the control efficiency of wild $N.lugens$ in system (19(19) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x\left(t\right)(1-\frac{2c_{1}z\left(t\right)}{x\left(t\right)+2z\left(t\right)})-{\delta }_1\left(x\right(t)+z(t\left)\right)x\left(t\right),\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-\left(x\right(t)+z(t\left)\right)z\left(t\right),\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(19) ), and the last four lines of the legend represent the control efficiency of wild $N.lugens$ in system (23(23) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =X\left(t\right)\frac{X\left(t\right)+(1-c_1)Z\left(t\right)}{X\left(t\right)+Z\left(t\right)}-{\delta }_1\left(X\right(t)+Z(t\left)\right)X\left(t\right),\\ \frac{\mathrm{d}\left(t\right)}{\mathrm{d}t}& =Z\left(t\right)(1-(X\left(t\right)+Z\left(t\right)\left)\right),\end{array}\}t\ne nT.\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(23) ). (a) The effect of β on control degree of the wild $N.lugens$. (b) The effect of T on control degree of the wild $N.lugens$.

From Figures 9–11, it is easy to see that under the same release amount and release cycle, the efficiency of population replacement strategy in controlling wild $N.lugens$ within a finite time is significantly higher than that of population suppression strategy. Furthermore, these results also suggest that increasing the release amount β or decreasing the release period T are beneficial for controlling the wild $N.lugens$.

Finally, we will simulate the change of control degree of wild $N.lugens$ population under a finite release number of times. Fixed parameters $c_1=0.675$, $c_2=0.732$, ${\delta }_1=0.85$, ${\delta }_2=0.2$, $\text{β}=0.5$, T=4, $x\left(0\right)=1$, $y\left(0\right)=0.3$ and $z\left(0^{+}\right)=\text{β}$. Based on these parameters, we discuss the situation that the number of release times of $N.lugens$ infected with wStri is 1, 2, 3 and 4 respectively, as shown in Figure 12(12) $\begin{array}{r}\frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}\le y\left(t\right)\left(\rho \right(1-\theta )(1-\frac{2c_2\left(\tilde{g}\right(t)-ϵ)}{d+2\tilde{g}\left(t\right)-ϵ})-{\delta }_2(\tilde{g}\left(t\right)-ϵ\left)\right),t\ge n_{2}T.\end{array}$(12) . It can be seen from Figure 12(12) $\begin{array}{r}\frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}\le y\left(t\right)\left(\rho \right(1-\theta )(1-\frac{2c_2\left(\tilde{g}\right(t)-ϵ)}{d+2\tilde{g}\left(t\right)-ϵ})-{\delta }_2(\tilde{g}\left(t\right)-ϵ\left)\right),t\ge n_{2}T.\end{array}$(12) that for the population replacement strategy, when the release number of times of $N.lugens$ infected with wStri exceeds 3, the wild $N.lugens$ population will be controlled. However, for population suppression strategy, it is impossible to achieve complete control of wild $N.lugens$ through limited release times. This indicates that the cost of population replacement strategy in controlling wild population will be less than that of population suppression strategy.

In Section 3, we obtained the conditions for the stability of the wild-$N.lugens$ -eradication periodic solution and verified them numerically in Subsections 4.1–4.2 (see Figures 3–4 and 6–8). This also shows that both release strategies are able to control wild $N.lugens$ populations. In this paper, both release strategies assume that the $N.lugens$ infected with wStri is released in periodic pulses, so the important parameters for control measures are the pulse release amount β and pulse release period T of the $N.lugens$ infected with wStri. For population suppression strategy, if the values of the release period T and the release amount β satisfy the conditions of Theorem 3.2, the wild $N.lugens$ will become extinct. If the values of β and T satisfy the conditions of Theorem 3.1 but not the conditions of Theorem 3.2, the eradication of wild $N.lugens$ depends on the initial values. For population replacement strategy, if the values of T and β satisfy the conditions of Theorem 3.6, the wild $N.lugens$ will be replaced. Since the complexity of the expressions for $R_{mi}$ and $R_{fi}$, we cannot directly compare the two strategies theoretically, but we have employed numerical simulations to compare the two strategies in Subsection 4.3, and the population replacement strategy can control wild $N.lugens$ faster than the population suppression strategy for the same release period T and the same amount of release β (see Figures 9–11). Moreover, the population replacement strategy can achieve the replacement of the wild $N.lugens$ population with a finite number of releases (see Figure 12). Of course, the above scenario occurs mainly because there is perfect maternal transmission of female $N.lugens$ infected with wStri. However, if the rice population density is added to models (1(1) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}{U}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_f\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{U}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_m\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{L}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_f\left(t\right),\\ \frac{\mathrm{d}{L}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_m\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }{U}_f\left(t\right)& =0,\mathrm{\Delta }{U}_m\left(t\right)=0,\mathrm{\Delta }{L}_f\left(t\right)=0,\mathrm{\Delta }{L}_m\left(t\right)=0,\\ \mathrm{\Delta }{S}_m\left(t\right)& =b,\end{array}t=n\tau ,n\in Z_{+},\end{array}\end{array}$(1) ) and (4(4) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}{U}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_f\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{U}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{U}_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_1)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{1}N\left(t\right)U_m\left(t\right)\\ & +\frac{1}{2}\rho a\theta L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)},\\ \frac{\mathrm{d}{L}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a(1-\theta )L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+(1-c_2)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_f\left(t\right),\\ \frac{\mathrm{d}{L}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}\rho a\left(1-\theta \right)L_f\left(t\right)\frac{U_m\left(t\right)+L_m\left(t\right)+\left(1-c_2\right)S_m\left(t\right)}{U_m\left(t\right)+L_m\left(t\right)+S_m\left(t\right)}-d_{2}N\left(t\right)L_m\left(t\right),\\ \frac{\mathrm{d}{S}_f\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{S}_f\left(t\right)-d_{3}N\left(t\right)S_f\left(t\right),\\ \frac{\mathrm{d}{S}_m\left(t\right)}{\mathrm{d}t}& =\frac{1}{2}a{S}_f\left(t\right)-d_{3}N\left(t\right)S_m\left(t\right),\end{array}\}t\ne n\tau ,\\ & \begin{array}{rl}\mathrm{\Delta }{U}_f\left(t\right)& =0,\mathrm{\Delta }{U}_m\left(t\right)=0,\\ \mathrm{\Delta }{L}_f\left(t\right)& =0,\mathrm{\Delta }{L}_m\left(t\right)=0,\\ \mathrm{\Delta }{S}_f\left(t\right)& =\frac{1}{2}b,\mathrm{\Delta }{S}_m\left(t\right)=\frac{1}{2}b,\end{array}\}t=n\tau ,\end{array}\end{array}$(4) ), both models will become more complex, which will be one of our next major works. Furthermore, we provided the conditions for the persistence in the mean of wild $N.lugens$ in Section 3, but Figure 5 shows the existence of a positive periodic solution for model (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ). Unfortunately, we have not found a good mathematical method to solve this problem due to the nonlinearity of the equations, which is also the direction of our future research.

Figure 12. Control of the wild $N.lugens$ under a finite release number of times. The second, fourth, sixth, and eighth lines of the legend represent the control degree curve of wild $N.lugens$ by using population replacement strategy, and the first, third, fifth, and seventh lines of the legend represent the control degree curve of wild $N.lugens$ by using population suppression strategy. (a) The control degree of the uninfected wild $N.lugens$. (b) The control degree of the wild $N.lugens$ infected with wLug.

5. Conclusions

Gong et al. [6] successfully developed a stable artificial Wolbachia infection of $N.lugens$ by introducing the Wolbachia strain wStri from $L.striatellus$ host into $N.lugens$ . The use of $N.lugens$ infected with wStri to control wild $N.lugens$ will be one of the important ways in the future. It is therefore of great interest to study how the release of $N.lugens$ infected with wStri. In this paper, we established two models with periodic impulsive release: the population suppression model (3(3) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =x(1-\frac{2c_{1}z}{x+y+2z})+\rho \theta y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_1(x+y+z)x,\\ \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )y(1-\frac{2c_{2}z}{x+y+2z})-{\delta }_2(x+y+z)y,\\ \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}& =-(x+y+z)z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }x\left(t\right)& =0,\mathrm{\Delta }y\left(t\right)=0,\mathrm{\Delta }z\left(t\right)=\text{β},t=nT,n\in Z_{+}.\end{array}\end{array}\end{array}$(3) ) and the population replacement model (6(6) $\begin{array}{r}\{\begin{array}{rl}& \begin{array}{rl}\frac{\mathrm{d}X\left(t\right)}{\mathrm{d}t}& =X(1-\frac{c_{1}Z}{X+Y+Z})+\rho \theta Y(1-\frac{c_{2}Z}{X+Y+Z})\\ & -{\delta }_1(X+Y+Z)X,\\ \frac{\mathrm{d}Y\left(t\right)}{\mathrm{d}t}& =\rho (1-\theta )Y(1-\frac{c_{2}Z}{X+Y+Z})-{\delta }_2(X+Y+Z)Y,\\ \frac{\mathrm{d}Z\left(t\right)}{\mathrm{d}t}& =Z-(X+Y+Z)Z,\end{array}\}t\ne nT,\\ & \begin{array}{rl}\mathrm{\Delta }X\left(t\right)& =0,\mathrm{\Delta }Y\left(t\right)=0,\mathrm{\Delta }Z\left(t\right)=\text{β},t=nT.\end{array}\end{array}\end{array}$(6) ). Applying Floquet theory and the comparison theorem of impulsive differential equations, we obtained the conditions for the stability of the wild-$N.lugens$ -eradication periodic solution of both models (see Figures 3–4 and 6–8). Meanwhile, sufficient conditions for the persistence in the mean of wild uninfected $N.lugens$ and the extinction of $N.lugens$ infected with wLug were obtained, and the conditions for the persistence in the mean of wild uninfected $N.lugens$ and $N.lugens$ infected with wLug were also given. Numerical simulations were performed to verify our theoretical results. Finally, we compared the control effect of two periodic pulse release strategies on wild $N.lugens$. The population replacement strategy is obviously much better than the population suppression strategy.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by National Natural Scientific Fund of China [Grant Number 12261018], Universities Key Laboratory of Mathematical Modeling and Data Mining in Guizhou Province [No. 2023013], and the Project of High Level Creative Talents in Guizhou Province [Grant Number 20164035].