# On an Impulsive Food Web System with Mutual Interference and Distributed Time Delay

Authors: Zhen Wang; Liwei Liu; Guangwang Su; Yuanfu Shao
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, we present a predator-prey system with mutual interference and distributed time delay and study its dynamical behavior. Based on the existence and universality of mutual interference among species, it is necessary to further study an impulsive food web system. By using stability theory, slight perturbation technique, and comparison theorem, we obtain some theoretical results of the system, such as boundedness and permanence. Moreover, numerical experiments are used to verify the theoretical results and to explore the dynamical behavior of the system, which exhibits rich dynamical behavior such as chaotic oscillation, periodic oscillation, symmetry-breaking bifurcations, chaotic crises, and period bifurcation. Finally, we give some practical guidelines for biological systems based on the theoretical results and numerical experiments of the system.


---

## Body

## 1. Introduction

In the past few decades, food web or food chain predator-prey system has been widely studied by many scholars. Many researchers have obtained some classical theoretical knowledge of the predator-prey system and given some practical control strategies, regulation means, and simulation analysis; all these results can be found in these papers [1–4]. With the evolution and development of the ecological system, the relationships of dynamics between predators and prey will continue to be the significant theme of ecology and mathematics, since its existence and importance.For the predator-prey system, the average consumption rate of a prey by a predator is a hot topic, which is called functional response, and it affects the population relationship between predator and prey. In recent work, for example, these papers [5–7] presented the influence of different functional response functions on the predator-prey system. Moreover, for typical food chain or food web system with an impulsive effect, we can see these papers [8–10]. Specially, some models combined different functional response functions to present the complexity of the system, such as Sahoo and Poria [11] combined Holling type-I and Holling type-II to study the role of alternative food as a disease controller in a disease-induced predator-prey system. Sahoo and Poria [12] combined general Holling type and Holling type-I to present a diseased prey-predator model, which gave an idea for constructing a realistic food chain model through a proper choice of general Holling parameters. Ghosh et al. [13] discussed the effect of a modified form of Holling type-II response function on additional food for a prey-predator model with prey refuge. In the evolution of biological systems, the state of species changes often occurs at a certain moment. These processes are affected by short-term disturbances; compared with process duration, the state change time always can be neglected [14, 15]. So it can be assumed that these perturbations are instantaneous; i.e., these perturbations are the form of impulse. Obviously, the impulsive effect is often manifested in the medical mutation rhythm model, biological control model, economic optimization model, frequency-modulated dynamic model, and pharmacokinetics; we can see the related literature in [16–19]; and the papers [20, 21] pointed out that impulsive differential equations are a good mathematical tool for describing the evolution of these models. In addition, time delay has also been extensively studied in predator-prey systems. In general, two types of time delays have been extensively studied [22]; i.e., continuous/distributed time delay and discrete delay. Some impulsive predator-prey models with distributed time delays are investigated; for example, these papers [23–25] showed some theoretical results and guiding significance of several impulsive systems with distributed time delays. Zhao et al. [26] presented an ecological model with the effect of impulsive and time delay as follows:(1)dx1tdt=rx1tk0−x1tk1−x1t−b3x12t−a3x1tx2t−a1x1tx3tb1+x1t+c1x3t,dx2tdt=d1x2t∫−∞tFt−sx1sds−a2x2tx3tb2+x2t+c2x3t−m1x2t,dx3tdt=e1a1x1tx3tb1+x1t+c1x3t+e2a2x2tx3tb2+x2t+c2x3t−m2x3t,t≠nT,x1t+−x1t=0,x2t+−x2t=0,x3t+−x3t=p,t=nT,where x1t is the lowest level prey’s density, x2t is the midlevel predator’s density, and x3t is the top predator’s density at time t.However, in our best knowledge, mutual interference is rarely considered in the predator-prey system. Actually, the phenomenon of mutual interference between predators and prey has always existed, which was studied by Hassell in 1971, and Hassell detailedly introduced the mutual interference constantm0<m≤1. For example, a Volterra system which contains mutual interference is explored by Hassell as follows [27]:(2)x′=xgx−ψxym,y′=y−d+kψxym−1−qy,where ψt is a functional response. Compared with the predator-prey system without mutual interference, the predator-prey system with mutual interference contains more dynamic behavior. For example, Wang and Zhu [28] proposed the existence of periodic solutions with mutual interference as follows:(3)x′t=xta1t−b1txt−c1txtk+xtymt,y′t=yt−a2t−b2tyt+c1txtk+xtymt,m∈0,1.Guo and Chen [29] explored a predator-prey model with mutual interference as follows:(4)x′t=xta1t−b1txt−k1txta+bxt+cytymt,y′t=yt−a2t−b2tyt+k2txta+bxt+cytymt,0<m≤1.There are few literatures that present the complex dynamic behavior of mutual interference systems [30–32].On the one hand, chaos has a serious impact on the management of dynamic population systems. Schaffer and Kot [33] discussed that chaos is a more complex phenomenon, and controlling chaos is a challenging problem. Sahoo and Poria [34] applied a technique of controlling chaotic predator-prey population dynamics by supplying additional food to the top predator, which study is aimed to introduce a new nonchemical chaos control mechanism in a predator-prey system with the applications in fishery management and biological conservation of prey-predator species. The control and management of chaotic population dynamics are one of the main goals of today’s ecological mathematical modeling. In our work, we also explore the impact of chaos on complex systems.In our work, our main purpose is to explore the effect of mutual interference in an impulsive food web with distributed time delays. The food web system described by differential equations is presented as follows:(5)x′t=rxtk0−xtk1−xt−bx2t−a3xtyt−a1xtztmb1+xt+c1zt,y′t=dyt∫−∞tGt−sxsds−a2ytztmb2+yt+c2zt−m1yt,z′t=e1a1xtztmb1+xt+c1zt+e2a2ytztmb2+yt+c2zt−m2zt,t≠nT,t≠n+l−1T,xt+=1−δ1xt,yt+=1−δ2yt,zt+=1−δ3zt,t=n+l−1T,xt+−xt=0,yt+−yt=0,zt+−zt=p,t=nT,where xt,yt,andzt represent the densities of bottom prey, medium-level predators, and top level predator at time t, respectively. All parameters about the system are positive. Prey’s intrinsic increase rate is denoted by r. Capture rate of predator to prey was expressed by aii=1,2,3. In addition, eii=1,2 represents the transformation rate of predator to prey. m1 and m2 represent the mortality rate of predators y and z, respectively. The fraction of available resources was represented by k1. rk0 represents the carrying capacity of the prey, where 0≤k0/k1≤1. b represents the internal competition rate of prey xt. m is the parameter of mutual interference, and its range of values is 0,1. The species are harvested at time t=n+l−1T, and the fraction of harvest is δ1, δ2, and δ3, respectively, where their range of values are 0,1. Define d as the product of the conversing rate of prey by predators and average predation rate. Predators zt are added to the system at t=nT, and the increased amount is p. In this mathematical model, the impulsive period is represented by T, 0<l<1; n is the nonnegative integer. In fact, it needs some time for predators to grow by consuming prey, and the time delay occurs in the functional response term, where the function Gt=ae−at satisfies ∫0+∞Gsds=1 and a is defined as the parameter of distributed time delay.Now, we can take a chain transformationqt=∫−∞tGt−sxsds to better explore the system (5). Since ∫−∞tGt−sxsds is convergent and(6)∫−∞tGt−sds=limA⟶−∞∫Atae−at−sds=1,we can get(7)Δqt=∫−∞t+Gt−sxsds−∫−∞tGt−sxsds=0,t=nT,t=n+l−1T.Hence, system (5) can be converted to system (8) as follows:(8)x′t=rxtk0−xtk1−xt−bx2t−a3xtyt−a1xtztmb1+xt+c1zt,y′t=dytqt−a2ytztmb2+yt+c2zt−m1yt,z′t=e1a1xtztmb1+xt+c1zt+e2a2ytztmb2+yt+c2zt−m2zt,q′t=axt−aqt,t≠nT,t≠n+l−1T,xt+=1−δ1xt,yt+=1−δ2yt,zt+=1−δ3zt,qt+=qt,t=n+l−1T,xt+=xt,yt+=yt,zt+=zt+p,qt+=qt,t=nT.Based on the above discussion, we can explore the dynamic property of system (5) by studying system (8). Hence, system (8) will be mainly considered in this paper.The remaining sections of our paper are presented as follows. In order to get the theoretical results of this paper, we give some preliminaries in Section2. In Section 3, the main theoretical results of this article are given. We compare this paper with another paper in Section 4. Numerical experiments are used to explore the mutual interference and the other dynamic property of the system in Section 5. In Section 6, a conclusion and a brief discussion of this article are given.


## 2. Preliminaries

Some preliminaries for proving the theoretical results of systems are given in this section, which is applied to Section3.X
t
=
x
t
,
y
t
,
z
t
,
q
t
:
R
+
⟶
R
+
4 is defined as the solution of system (8), which is continuous in interval n−1T,n+l−1T and n+l−1T,nT and R+=0,∞, R+4=xt,yt,zt,qtx,y,z,q≥0, n∈N. Define s as the mapping on the right side of system (8). Letting V:R+×R+4⟶R+, then V∈V0, if(i)
For each solutionXt∈R+4, the two limits limt,y⟶n+l−1T+,XVt,y=Vn+l−1T+,X and limt,y⟶nT+,XVt,y=VnT+,X exist, where V is the continuous in interval n−1T,n+l−1T×R+4 and n+l−1T,nT×R+4,n=1,2…(ii)
V is locally Lipschitzian in XtDefinition 1.
If there are two positive constantsm and M that make m≤xt,yt,zt,qt≤M, then system (8) is permanent.Definition 2.
DefineD+Vt,X=limh⟶0+sup1/hVt+h,X+hft,X−Vt,X as the upper right derivative of Vt,X of system (8), where V∈V0, t,X∈n−1T,n+l−1T×R+4, and n+l−1T,nT×R+4.Remark 1.
X
t: R+⟶R+4, Xt is continuous in n−1T,n+l−1T∪n+l−1T,nT, and Xt0+=limt⟶t0+Xt exists, where t0=n+l−lT+ and nT+. The solution Xt is a piecewise continuous function. At the same time, the mapping f with smoothness properties can ensure the global existence and uniqueness of the model (it can see [35]).Lemma 1 (see [35]).
LetV∈V0, h:R+×R+⟶R, and h is the continuous in interval n−1T,n+l−1T and n+l−1T,nT. For g∈R+2,n∈N, the two limits limt,y⟶nT+,ght,y=hnT+,g and limt,y⟶n+l−1T+,ght,y=hn+l−1T+,g exist. Suppose that ψn1,ψn2:R+⟶R+, which are nondecreasing function for all n∈N. And assume that, when t≠nT,n+l−1T,(9)D+Vt,X≤ht,Vt,X.When t=n+l−1T,(10)Vt,Xt+≤ψn1Vt,X,and when t=nT,(11)Vt,Xt+≤ψn2Vt,X.Suppose thatrt be the maximal solution of the following impulsive differential:(12)g′t=ht,gt,t≠nT,t≠n+l−1T,gt+=ψn1gt,t=n+l−1T,gt+=ψn2gt,t=nT,g0+=g0≥0,which exists on 0,∞. Then, the inequality V0+,X0≤g0 holds, which shows that the inequality Vt,Xt≤rt (all t≥0) holds.Now, we study the following system, which is the subsystem of system (8):(13)z′t=−m2zt,t≠nT,t≠n+l−1T,zt+=1−δ3zt,t=n+l−1T,zt+=zt+p,t=nT,z0+=z0.After some algebraic calculations, we can get positive periodic solution of this subsystem as follows:(14)z∗t=pexp−m2t−n−1T1−1−δ3exp−m2T,n−1T<t≤n+l−1T,p1−δ3exp−m2t−n−1T1−1−δ3exp−m2T,n+l−1T<t≤nT,(15)n∈N,z∗0+=z∗nT+=p1−1−δ3exp−m2T,z∗n+l−1T+=z∗lT+=p1−δ3exp−m2lT1−1−δ3exp−m2T,which is a positive periodic solution of system (13). Since(16)zt=1−δ3n−1z0+−p1−1−δ3exp−m2Texp−m2t+z∗t,n−1T<t≤n+l−1T,1−δ3nz0+−p1−1−δ3exp−m2Texp−m2t+z∗t,n+l−1T<t≤nT,is the solution of system (13) with initial value z0+≥0, we can get following lemma.Lemma 2 (see [35]).
Supposez∗t is a positive periodic solution of system (13), for each solution zt with z0+≥0, it follows that zt−z∗t⟶0 as t⟶∞.Lemma 3 (see [36]).
Suppose that functionVt∈PC1R+,R satisfies the following equation:(17)dVdt≤wtVt+ht,t≠tk,t>0,Vtk+≤jkVtk+lk,t=tk>0,g0+≤g0,where w,h∈CR+,R, jk≥0, and lk are constants k=1,2,…. For t>0, we have(18)Vt≤V0∏t0<tk<tjkexp∫t0twsds+∫t0t∏t0<tk<tjkexp∫stwγdγhsds+∑t0<tk<t∏tk<tj<tjjexp∫tktwsdslk.
Next, we explore boundedness, globally asymptotically stable, and permanence of system (8) in next sections.


## 3. Main Theoretical Results

### 3.1. Boundedness of System (8)

Theorem 1.
Letxt,yt,zt,qt be the solution of system (8), and it always can find a positive constant M, which makes xt≤M,yt≤M,zt≤M, and qt≤M hold.Proof.
For the solutionXt=xt,yt,zt,qt of model (8) with initial value x0,y0,z0,q0, let Vt,Xt=e1xt+e2yt+zt+qt, where V∈V0. Since(19)dxtdt≤rxtk0−xtk1−xt−bx2t≤xtr−bxt,qt=∫−∞tFt−sxsds,we can obtain xt,qt≤r/b. Choose 0<λ<minm1−rd/b,m2,a, we can get the following inequality relations:(20)D+Vt+λVt≤−e1bx2t+λe1xt+axt+e1rxt+λe2+re2db−e2m1yt+λ−m2zt+λ−aqt≤M0,t≠n+l−1T,t≠nT,where(21)M0=λe1+a+e1r24e1b.
In addition, clearly,(22)Vn+l−1T+≤Vn+l−1T,t=n+l−1T;VnT+≤VnT+p,t=nT.
By using Lemma 3, we can obtain(23)D+Vt+λVt≤−e1bx2t+λe1xt+axt+e1rxt+λe2+re2db−e2m1yt+λ−m2zt+λ−aqt≤M0.
To sum up the above discussion, for sufficiently larget, Vt,Xt is bounded; i.e., xt≤M,yt≤M,zt≤M, and qt≤M hold for sufficiently large t. □


### 3.2. Globally Asymptotically Stable Solution0,0,z∗t,0

Theorem 2.
Let0,0,z∗t,0 be the eradication solution of prey and intermediate predator of system (8), which is globally asymptotically stable if and only if the following two inequalities hold:(24)rk0Tk1+a1c1mm2A1+A2<ln11−δ1,−m1T+a2c2mm2B1+B2<ln11−δ2,where(25)A=1−1−δ3exp−m2Tm,A1=ln1−c1pm1−exp−m2mlTb1A+c1pm,A2=ln1−c1pm1−δ3mexp−lm2mT−exp−m2mTb1A+c1pm1−δ3mexp−lm2mT,B1=ln1−c2pm1−exp−lm2mTb2A+c2pm,B2=ln1−c2pm1−δ3mexp−lm2mT−exp−m2mTb2A+c2pm1−δ3mexp−lm2mT.Proof.
LetXt=xt,yt,zt,qt be the solution of system (8), by taking the technique of small amplitude perturbations for the solution, and we can derive the local stability of periodic solutions 0,0,z∗t,0. Let xt=u1t,yt=u2t,zt=z∗t+u3t,andqt=u4t, where uiti=1,2,3,4 is a small perturbation. When t≠nT, n+l−1T, by applying Taylor series expansion and neglecting higher order terms, system (8) can be read as follows:(26)u1′t=rk0k1−a1z∗tmb1+c1z∗tu1t,u2′t=−a2z∗tmb2+c2z∗t−m1u2t,u3′t=e1a1z∗tmu1tb1+c1z∗t+e2a2z∗tmu2tb2+c2z∗t−m2u3t,u4′t=au1t−au4t,t≠nT,t≠n+l−1T,u1t+=1−δ1u1t,u2t+=1−δ2u2t,u3t+=1−δ3u3t,u4t+=u4t,t=n+l−1T,u1t+=u1t,u2t+=u2t,u3t+=u3t,u4t+=u4t,t=nT,so we can get(27)u1tu2tu3tu4t=Φtu10u20u30u40,0≤t<T,where Φt is the fundamental matrix of (26). At the same time, Φt satisfies(28)dΦtdt=rk0k1−a1z∗tmb1+c1z∗t0000−a2z∗tmb2+c2z∗t−m100e1a1z∗tmb1+c1z∗te2a2z∗tmb2+c2z∗t−m20a00−aΦt.Let t=0, then Φ0=I is a identity matrix, so we can get the linearized structure of system (8) as follows:(29)u1n+l−1T+u2n+l−1T+u3n+l−1T+u4n+l−1T+=1−δ100001−δ200001−δ300001u1n+l−1Tu2n+l−1Tu3n+l−1Tu4n+l−1T,u1nT+u2nT+u3nT+u4nT+=1000010000100001u1nTu2nTu3nTu4nT.
By using the eigenvalues of following matrixM, we can derive the stability of the solutions 0,0,z∗t,0:(30)M=1−δ100001−δ200001−δ3000011000010000100001ΦT.
Ifλi<1 i=1,2,3,4, then the periodic solution 0,0,z∗t,0 is locally stable, and their expressions are as follows:(31)λ1=1−δ1exp∫0Trk0k1−a1z∗tmb1+c1z∗tdt,λ2=1−δ2exp∫0T−a2z∗tmb2+c2z∗t−m1dt,λ3=1−δ3exp−m2T,λ4=exp−aT.
The analysis is similar to literature [35], i.e., by using Floquet theory, if λi<1i=1,2,3,4, then 0,0,z∗t,0 is locally asymptotically stable. So we need to calculate λi<1i=1,2; herein, λ3 and λ4 are already less than one. Therefore, we have(32)λ1<1⟺T<k1rk0ln11−δ1−a1c1mm2A1+A2,λ2<1⟺T>1−m1ln11−δ2−a2c2mm2B1+B2,where(33)∫0Tz∗tmdt=∫n−1Tn+l−1Tpexp−m2t−n−1T1−1−δ3exp−m2Tmdt+∫n+l−1TnTp1−δ3exp−m2t−n−1T1−1−δ3exp−m2Tmdt.
Inequality (32) holds if and only if(34)rk0Tk1+a1c1mm2A1+A2<ln11−δ1,−m1T+a2c2mm2B1+B2<ln11−δ2.
Hence, when the inequality (34) is established, the system (8) is locally asymptotically stable.
After the above discussion, we can obtain the conditions of local asymptotic stability of periodic solutions0,0,z∗t,0. Now, we explore global attraction of the periodic solutions; for ε>0, we have(35)ξ1≜1−δ1exp∫n+l−1Tn+lTrk0k1−a1z∗t−εmb1+c1z∗t−εdt<1,ξ2≜1−δ2exp∫0TdM−m1−a2z∗t−εmb2+c2z∗t−εdt<1.
Noticing thatdzt/dt≥−m2zt, according to Lemma 2 and comparison theorem of the impulsive equation, for all sufficiently large t, we can obtain(36)zt>z∗t−ε.
For allt>0, assuming that inequality (36) is established, so we consider the system (8) and inequality (36), we can get(37)dxtdt≤rk0k1−a1z∗t−εmb1+c1z∗t−εxt,t≠n+l−1T,xt+=1−δ1x1t,t=n+l−1T,integrating (37) on n+l−1T,n−1T, and from (8), we have(38)xn+lT≤xn+l−1T+exp∫n+l−1Tn+lTrk0k1−a1z∗t−εmb1+c1z∗t−εdt,=xn+l−1T1−δ1exp∫n+l−1Tn+lTrk0k1−a1z∗t−εmb1+c1z∗t−εdt,=xn+l−1Tξ1.
After some algebraic calculations of (38), which produces xn+lT≤xlTξ1n, when n⟶∞, it follows that xn+lT⟶0. When n+l−1T<t≤n+lT, we get 0<xt<xn+l−1T1−δ1exprT, so xt⟶0 as n⟶∞. In addition, if xt⟶0, which produces qt⟶0. Then, for yt, we have(39)dytdt≤dM−m1−a2z∗t−εmb2+c2z∗t−εyt,t≠n+l−1T,yt+=1−δ2y1t,t=n+l−1T.
The analysis is similar toxt, and it can also obtain yt⟶0 as n⟶∞.
Next, our goal is to explorezt to z∗t t⟶∞. Let 0<ε<r3/γγ=e1a1/b1+e2a2/b2, we can find a T′>0, which produces 0<xt,yt,qt<ε as t>T′. Actually, assume that 0<xt,yt,qt<ε for all t≥0. Then, we have the following inequality by system (8):(40)−m2zt≤dztdt≤−m2zt+γεztm,γ=e1a1b1+e2a2b2.
Because0<m≤1, the inequality (40) can be changed to −m2zt≤dzt/dt≤−m2+γεzt. We now consider comparative systems (41) and (42):(41)dv1tdt=−m2v1t,t≠nT,t≠n+l−1T,v1t+=1−δ3v1t,t=n+l−1T,v1t+=v1t+p,t=nT,v10+=z0+,(42)dv2tdt=−m2+γεv2t,t≠nT,t≠n+l−1T,v2t+=1−δ3v2t,t=n+l−1T,v2t+=v2t+p,t=nT,v20+=z0+,let v1t and v2t be any solution of system (41) and (42), respectively, where(43)v2∗t=pexp−m2+γεt−n−1T1−1−δ3exp−m2+γεT,n−1T<t≤n+l−1T,p1−δ3exp−m2+γεt−n−1T1−1−δ3exp−m2+γεT,n+l−1T<t≤nT.
By using Lemmas1 and 2, which can derive v1t≤zt≤v2t. Furthermore, when t⟶∞, v1t⟶z∗t and v2t⟶v2∗t. Then, ∃T1>0 for any ε1>0, and it yields the following inequalities:(44)z∗t−ε1<zt<z∗t+ε,t≥T1.
Obviously,v2∗t⟶z∗t as ε⟶0. Therefore, for t large enough, it yields z∗t−ε1<zt<z∗t+ε1, which implies zt⟶z∗t as t⟶∞. □


### 3.3. Permanence of the System

Theorem 3.
If the following two inequalities hold, then system (8) is permanent:(45)rk0Tk1+a1c1mm2A1+A2>ln11−δ1,−m1T+a2c2mm2B1+B2>ln11−δ2,where A,A1,A2,B1, and B2 are given in Theorem 2.Proof.
From Theorem1, we can know that there is a constant M which makes xt,yt,zt,qt<M for t≥0. Since z′t≥−m2zt, we give the following differential equation (46):(46)v3′t=−m2v3t,t≠nT,t≠n+l−1T,v3t+=1−δ3v3t,t=n+l−1T,v3t+=v3t+p,t=nT,v30+=z0+>0.
By Lemmas1 and 2, for sufficiently large t≥0, we can get zt≥v3t>z∗t−ε and zt≥p1−δ3exp−m2T/1−1−δ3exp−m2T−ε2=m3ε2>0. So our next goal is to find positive m1 and m2 such that xt≥m1,yt≥m2, and qt≥m1 for t large enough. Now we are exploring how to find m1.
The following two steps prove the existence of such am1. Firstly, let 0<m1¯<b1m2/e1a1Mm−1,b2m2/e2a2Mm−1 and sufficiently small ε1>0, then we have(47)γ1≜1−δ1exp∫n+l−1Tn+lTrk0−m1¯k1−b+a3M−a1b1v4∗t+η1mdt>1.
It is easy to prove thatxt<m1¯ can not hold for all t≥0. Otherwise, we have(48)z′t=zt−m2+e1a1xtztm−1b1+e2a2ytztm−1b2,≤zt−m2+e1a1m1¯ztm−1b1+e2a2Mztm−1b2≜ztB,t≠nT,t≠n+l−1T,zt+=1−δ3zt,t=n+l−1T,zt+=zt+p,t=nT,z0+=z0.
So we can obtainv4t⟶v4∗t as t⟶∞ (v4t is the solution of (49)) and zt≤v4t:(49)v4′t=Bv4t,t≠nT,t≠n+l−1T,v4t+=1−δ3v4t,t=n+l−1T,v4t+=v4t+p,t=nT,v40+=z0+>0,and its periodic solution is as follows:(50)v∗t=pexpBt−n−1T1−1−δ3expBT,n−1T<t≤n+l−1T,p1−δ3expBt−n−1T1−1−δ3expBT,n+l−1T<t≤nT.
We can easily getzt≤v4t as t⟶+∞ by using the comparison theorem. So there is a T1>0, when t>T1, zt≤v4t<v4∗t+η1. Therefore, we have(51)x′t≥xtrk0−m1¯k1−b+a3M−a1b1v4∗t+η1m,t≠nT,t≠n+l−1T,xt+=1−δ1xt,t=n+l−1T.
N
1
∈
N and N1+l−1T≥T1. Integrating (51) on n+l−1T,n+lT,n≥N1, we have(52)xn+lT≥xn+l−1T+exp∫n+l−1Tn+lTrk0−m1¯k1−b+a3M−a1b1v4∗t+η1mdt=xn+l−1T1−δ1exp∫n+l−1Tn+lTrk0−m1¯k1−b+a3M−a1b1v4∗t+η1mdt=xn+l−1Tγ1.
Some algebraic calculations of (52), which produces xN1+kT≥xN1Tγ1k⟶∞ as k⟶∞. Obviously, this is contradictory to the boundedness of xt. Therefore, there is a t1>0 such that xt1≥m1¯.
Secondly, ifxt,yt,qt≥m1¯ for all t≥t1, then the proof of this theorem is completed. Otherwise, xt>m1¯ can not hold for all t≥t1. So we need to find out which solutions leave the set Xt∈R+4:xt≥m1¯, and let these solutions return to the set again. Let t∗=inft≥t1xt<m1¯, which has the following cases:Case 1: ift∗=n1+l−1T, then xt≥m1¯ and yt≥m1¯ for t∈t1,t∗ and 1−δ1m1¯≤xt∗+=1−δ1xt∗<m1¯ and 1−δ2m1¯≤yt∗+=1−δ2yt∗<m1¯. Choose n2,n3∈N such that n2−1T>lnη1/M+p/B1,1−δ1n2γ1n3expγ2n2Tγ1n3>1−δ1n2expγ2n2+1Tγ1n3>1, where γ2=rk0−m1¯/k1−b+a3M−a1/b1Mm<0. Let T¯=n2T+n3T, and it is deduced that there must be a t′∈t∗,t∗+T¯ such that xt′>m1¯; considering (49) with v4n1T+=x3n1T+, it is easy to get(53)v4t=1−δ3n−n1+lzn1T+−p1−1−δ3expBTexpBt−n1T+v4∗t,n−1T<t≤n+l−1T,1−δ3n−n1+l−1zn1T+−p1−1−δ3expB1TexpBt−n1T+v4∗t,n+l−1T<t≤nT,andn1+1≤n≤n1+n2+n3+1. So we can get v4t−v4∗t<M+pexpBt−n1T<η1 and zt≤v4t<v4∗t+η1 for t∗+n2T≤t≤t∗+T¯, which indicate that (51) holds for t∗+n2T≤t≤t∗+T¯. We can obtain xt∗+T¯≥xt∗+n2Tγ1n3 through a discussion similar to the first step. From the system (8), we have(54)x′t≥xtrk0−m1¯k1−b+a3M−a1b1Mm=γ2xt,t≠nT,t≠n+l−1T,xt+=1−δ1xt,t=n+l−1T.Integrating (54) on t∈t∗,t∗+n2T, then we get xt∗+n2T≥m1¯1−δ1n2expγ2n2T and xt∗+T¯≥m1¯1−δ1n2γ1n3expγ2n2T>m1¯; it is a contradiction with xt<m1¯. So xt>m1¯ holds for all t∈t1,t∗. Let t¯=inft≥t∗xt≥m1¯, then we can obtain xt≤m1¯,xt¯=m1¯ and xt≥m1¯1−δ1n2+n3expn2+n3γ2T≜m1 for t∗<t<t¯. Then, we derive that xt≥m1 holds for t∗<t<t¯. In addition, for t>t¯, since xt¯≥m1¯, the same discussion can be continued.Case 2: the case oft∗=n1+l−1T is discussed in Case 1. Now, assume that t∗≠n1+l−1T, which produces xt≥m1¯ for t1≤t<t∗ and x1t∗=m1¯. If t∗∈n1′+l−1T,n1′+lT,n1′∈N, which have two possible cases for t∈t∗,n1′+lT.Case 2(a): for allt∈t∗,n1′+lT, xt≤m1¯ is valid. By the same discussion with case 1, which exists t1′∈n1′+lT,n1′+lT+T¯ such that xt1′>m1¯; in order to simplify the proof, we will not elaborate on it. Let t˜=inft≥t∗xt>m1¯, then we have xt˜≤m1¯ for t∗≤t<t˜. We get xt≥m1¯1−δ1n2+n3+1expn2+n3γ2T≜m1 and xt≥m1 for t∈t∗,t˜. For t>t˜, the same arguments can be followed since x1t˜≥m1¯.Case 2(b): there existst∗<t<n1′+lT such that xt>m1¯. Let t˜=inft≥t∗xt>m1¯, then xt˜≤m1¯ for t∈t∗,t˜ and x1t˜=m1¯. When t∈t∗,t˜, (54) holds. Integrating (54) on t∗,t˜, we obtain xt≥xt∗expγ2t−t∗≥m1¯expγ2T>m1. Since xt≥m1¯ for t>t∗, the same arguments can be continued. Hence, xt≥m1 for all t≥t˜.
To sum up all the above discussions, we derive thatxt≥m1 holds for all t≥t˜. Using similar methods for yt, we can get yt≥m2 for all t≥t2. Obviously, qt≥m1 for all t>t˜1. Let m=maxm1,m2,m3; by Theorem 1 and Definition 2, we obtain that m<xt,yt,zt,qt<M. So, the system (8) proved to be permanent. □


## 3.1. Boundedness of System (8)

Theorem 1.
Letxt,yt,zt,qt be the solution of system (8), and it always can find a positive constant M, which makes xt≤M,yt≤M,zt≤M, and qt≤M hold.Proof.
For the solutionXt=xt,yt,zt,qt of model (8) with initial value x0,y0,z0,q0, let Vt,Xt=e1xt+e2yt+zt+qt, where V∈V0. Since(19)dxtdt≤rxtk0−xtk1−xt−bx2t≤xtr−bxt,qt=∫−∞tFt−sxsds,we can obtain xt,qt≤r/b. Choose 0<λ<minm1−rd/b,m2,a, we can get the following inequality relations:(20)D+Vt+λVt≤−e1bx2t+λe1xt+axt+e1rxt+λe2+re2db−e2m1yt+λ−m2zt+λ−aqt≤M0,t≠n+l−1T,t≠nT,where(21)M0=λe1+a+e1r24e1b.
In addition, clearly,(22)Vn+l−1T+≤Vn+l−1T,t=n+l−1T;VnT+≤VnT+p,t=nT.
By using Lemma 3, we can obtain(23)D+Vt+λVt≤−e1bx2t+λe1xt+axt+e1rxt+λe2+re2db−e2m1yt+λ−m2zt+λ−aqt≤M0.
To sum up the above discussion, for sufficiently larget, Vt,Xt is bounded; i.e., xt≤M,yt≤M,zt≤M, and qt≤M hold for sufficiently large t. □


## 3.2. Globally Asymptotically Stable Solution0,0,z∗t,0

Theorem 2.
Let0,0,z∗t,0 be the eradication solution of prey and intermediate predator of system (8), which is globally asymptotically stable if and only if the following two inequalities hold:(24)rk0Tk1+a1c1mm2A1+A2<ln11−δ1,−m1T+a2c2mm2B1+B2<ln11−δ2,where(25)A=1−1−δ3exp−m2Tm,A1=ln1−c1pm1−exp−m2mlTb1A+c1pm,A2=ln1−c1pm1−δ3mexp−lm2mT−exp−m2mTb1A+c1pm1−δ3mexp−lm2mT,B1=ln1−c2pm1−exp−lm2mTb2A+c2pm,B2=ln1−c2pm1−δ3mexp−lm2mT−exp−m2mTb2A+c2pm1−δ3mexp−lm2mT.Proof.
LetXt=xt,yt,zt,qt be the solution of system (8), by taking the technique of small amplitude perturbations for the solution, and we can derive the local stability of periodic solutions 0,0,z∗t,0. Let xt=u1t,yt=u2t,zt=z∗t+u3t,andqt=u4t, where uiti=1,2,3,4 is a small perturbation. When t≠nT, n+l−1T, by applying Taylor series expansion and neglecting higher order terms, system (8) can be read as follows:(26)u1′t=rk0k1−a1z∗tmb1+c1z∗tu1t,u2′t=−a2z∗tmb2+c2z∗t−m1u2t,u3′t=e1a1z∗tmu1tb1+c1z∗t+e2a2z∗tmu2tb2+c2z∗t−m2u3t,u4′t=au1t−au4t,t≠nT,t≠n+l−1T,u1t+=1−δ1u1t,u2t+=1−δ2u2t,u3t+=1−δ3u3t,u4t+=u4t,t=n+l−1T,u1t+=u1t,u2t+=u2t,u3t+=u3t,u4t+=u4t,t=nT,so we can get(27)u1tu2tu3tu4t=Φtu10u20u30u40,0≤t<T,where Φt is the fundamental matrix of (26). At the same time, Φt satisfies(28)dΦtdt=rk0k1−a1z∗tmb1+c1z∗t0000−a2z∗tmb2+c2z∗t−m100e1a1z∗tmb1+c1z∗te2a2z∗tmb2+c2z∗t−m20a00−aΦt.Let t=0, then Φ0=I is a identity matrix, so we can get the linearized structure of system (8) as follows:(29)u1n+l−1T+u2n+l−1T+u3n+l−1T+u4n+l−1T+=1−δ100001−δ200001−δ300001u1n+l−1Tu2n+l−1Tu3n+l−1Tu4n+l−1T,u1nT+u2nT+u3nT+u4nT+=1000010000100001u1nTu2nTu3nTu4nT.
By using the eigenvalues of following matrixM, we can derive the stability of the solutions 0,0,z∗t,0:(30)M=1−δ100001−δ200001−δ3000011000010000100001ΦT.
Ifλi<1 i=1,2,3,4, then the periodic solution 0,0,z∗t,0 is locally stable, and their expressions are as follows:(31)λ1=1−δ1exp∫0Trk0k1−a1z∗tmb1+c1z∗tdt,λ2=1−δ2exp∫0T−a2z∗tmb2+c2z∗t−m1dt,λ3=1−δ3exp−m2T,λ4=exp−aT.
The analysis is similar to literature [35], i.e., by using Floquet theory, if λi<1i=1,2,3,4, then 0,0,z∗t,0 is locally asymptotically stable. So we need to calculate λi<1i=1,2; herein, λ3 and λ4 are already less than one. Therefore, we have(32)λ1<1⟺T<k1rk0ln11−δ1−a1c1mm2A1+A2,λ2<1⟺T>1−m1ln11−δ2−a2c2mm2B1+B2,where(33)∫0Tz∗tmdt=∫n−1Tn+l−1Tpexp−m2t−n−1T1−1−δ3exp−m2Tmdt+∫n+l−1TnTp1−δ3exp−m2t−n−1T1−1−δ3exp−m2Tmdt.
Inequality (32) holds if and only if(34)rk0Tk1+a1c1mm2A1+A2<ln11−δ1,−m1T+a2c2mm2B1+B2<ln11−δ2.
Hence, when the inequality (34) is established, the system (8) is locally asymptotically stable.
After the above discussion, we can obtain the conditions of local asymptotic stability of periodic solutions0,0,z∗t,0. Now, we explore global attraction of the periodic solutions; for ε>0, we have(35)ξ1≜1−δ1exp∫n+l−1Tn+lTrk0k1−a1z∗t−εmb1+c1z∗t−εdt<1,ξ2≜1−δ2exp∫0TdM−m1−a2z∗t−εmb2+c2z∗t−εdt<1.
Noticing thatdzt/dt≥−m2zt, according to Lemma 2 and comparison theorem of the impulsive equation, for all sufficiently large t, we can obtain(36)zt>z∗t−ε.
For allt>0, assuming that inequality (36) is established, so we consider the system (8) and inequality (36), we can get(37)dxtdt≤rk0k1−a1z∗t−εmb1+c1z∗t−εxt,t≠n+l−1T,xt+=1−δ1x1t,t=n+l−1T,integrating (37) on n+l−1T,n−1T, and from (8), we have(38)xn+lT≤xn+l−1T+exp∫n+l−1Tn+lTrk0k1−a1z∗t−εmb1+c1z∗t−εdt,=xn+l−1T1−δ1exp∫n+l−1Tn+lTrk0k1−a1z∗t−εmb1+c1z∗t−εdt,=xn+l−1Tξ1.
After some algebraic calculations of (38), which produces xn+lT≤xlTξ1n, when n⟶∞, it follows that xn+lT⟶0. When n+l−1T<t≤n+lT, we get 0<xt<xn+l−1T1−δ1exprT, so xt⟶0 as n⟶∞. In addition, if xt⟶0, which produces qt⟶0. Then, for yt, we have(39)dytdt≤dM−m1−a2z∗t−εmb2+c2z∗t−εyt,t≠n+l−1T,yt+=1−δ2y1t,t=n+l−1T.
The analysis is similar toxt, and it can also obtain yt⟶0 as n⟶∞.
Next, our goal is to explorezt to z∗t t⟶∞. Let 0<ε<r3/γγ=e1a1/b1+e2a2/b2, we can find a T′>0, which produces 0<xt,yt,qt<ε as t>T′. Actually, assume that 0<xt,yt,qt<ε for all t≥0. Then, we have the following inequality by system (8):(40)−m2zt≤dztdt≤−m2zt+γεztm,γ=e1a1b1+e2a2b2.
Because0<m≤1, the inequality (40) can be changed to −m2zt≤dzt/dt≤−m2+γεzt. We now consider comparative systems (41) and (42):(41)dv1tdt=−m2v1t,t≠nT,t≠n+l−1T,v1t+=1−δ3v1t,t=n+l−1T,v1t+=v1t+p,t=nT,v10+=z0+,(42)dv2tdt=−m2+γεv2t,t≠nT,t≠n+l−1T,v2t+=1−δ3v2t,t=n+l−1T,v2t+=v2t+p,t=nT,v20+=z0+,let v1t and v2t be any solution of system (41) and (42), respectively, where(43)v2∗t=pexp−m2+γεt−n−1T1−1−δ3exp−m2+γεT,n−1T<t≤n+l−1T,p1−δ3exp−m2+γεt−n−1T1−1−δ3exp−m2+γεT,n+l−1T<t≤nT.
By using Lemmas1 and 2, which can derive v1t≤zt≤v2t. Furthermore, when t⟶∞, v1t⟶z∗t and v2t⟶v2∗t. Then, ∃T1>0 for any ε1>0, and it yields the following inequalities:(44)z∗t−ε1<zt<z∗t+ε,t≥T1.
Obviously,v2∗t⟶z∗t as ε⟶0. Therefore, for t large enough, it yields z∗t−ε1<zt<z∗t+ε1, which implies zt⟶z∗t as t⟶∞. □


## 3.3. Permanence of the System

Theorem 3.
If the following two inequalities hold, then system (8) is permanent:(45)rk0Tk1+a1c1mm2A1+A2>ln11−δ1,−m1T+a2c2mm2B1+B2>ln11−δ2,where A,A1,A2,B1, and B2 are given in Theorem 2.Proof.
From Theorem1, we can know that there is a constant M which makes xt,yt,zt,qt<M for t≥0. Since z′t≥−m2zt, we give the following differential equation (46):(46)v3′t=−m2v3t,t≠nT,t≠n+l−1T,v3t+=1−δ3v3t,t=n+l−1T,v3t+=v3t+p,t=nT,v30+=z0+>0.
By Lemmas1 and 2, for sufficiently large t≥0, we can get zt≥v3t>z∗t−ε and zt≥p1−δ3exp−m2T/1−1−δ3exp−m2T−ε2=m3ε2>0. So our next goal is to find positive m1 and m2 such that xt≥m1,yt≥m2, and qt≥m1 for t large enough. Now we are exploring how to find m1.
The following two steps prove the existence of such am1. Firstly, let 0<m1¯<b1m2/e1a1Mm−1,b2m2/e2a2Mm−1 and sufficiently small ε1>0, then we have(47)γ1≜1−δ1exp∫n+l−1Tn+lTrk0−m1¯k1−b+a3M−a1b1v4∗t+η1mdt>1.
It is easy to prove thatxt<m1¯ can not hold for all t≥0. Otherwise, we have(48)z′t=zt−m2+e1a1xtztm−1b1+e2a2ytztm−1b2,≤zt−m2+e1a1m1¯ztm−1b1+e2a2Mztm−1b2≜ztB,t≠nT,t≠n+l−1T,zt+=1−δ3zt,t=n+l−1T,zt+=zt+p,t=nT,z0+=z0.
So we can obtainv4t⟶v4∗t as t⟶∞ (v4t is the solution of (49)) and zt≤v4t:(49)v4′t=Bv4t,t≠nT,t≠n+l−1T,v4t+=1−δ3v4t,t=n+l−1T,v4t+=v4t+p,t=nT,v40+=z0+>0,and its periodic solution is as follows:(50)v∗t=pexpBt−n−1T1−1−δ3expBT,n−1T<t≤n+l−1T,p1−δ3expBt−n−1T1−1−δ3expBT,n+l−1T<t≤nT.
We can easily getzt≤v4t as t⟶+∞ by using the comparison theorem. So there is a T1>0, when t>T1, zt≤v4t<v4∗t+η1. Therefore, we have(51)x′t≥xtrk0−m1¯k1−b+a3M−a1b1v4∗t+η1m,t≠nT,t≠n+l−1T,xt+=1−δ1xt,t=n+l−1T.
N
1
∈
N and N1+l−1T≥T1. Integrating (51) on n+l−1T,n+lT,n≥N1, we have(52)xn+lT≥xn+l−1T+exp∫n+l−1Tn+lTrk0−m1¯k1−b+a3M−a1b1v4∗t+η1mdt=xn+l−1T1−δ1exp∫n+l−1Tn+lTrk0−m1¯k1−b+a3M−a1b1v4∗t+η1mdt=xn+l−1Tγ1.
Some algebraic calculations of (52), which produces xN1+kT≥xN1Tγ1k⟶∞ as k⟶∞. Obviously, this is contradictory to the boundedness of xt. Therefore, there is a t1>0 such that xt1≥m1¯.
Secondly, ifxt,yt,qt≥m1¯ for all t≥t1, then the proof of this theorem is completed. Otherwise, xt>m1¯ can not hold for all t≥t1. So we need to find out which solutions leave the set Xt∈R+4:xt≥m1¯, and let these solutions return to the set again. Let t∗=inft≥t1xt<m1¯, which has the following cases:Case 1: ift∗=n1+l−1T, then xt≥m1¯ and yt≥m1¯ for t∈t1,t∗ and 1−δ1m1¯≤xt∗+=1−δ1xt∗<m1¯ and 1−δ2m1¯≤yt∗+=1−δ2yt∗<m1¯. Choose n2,n3∈N such that n2−1T>lnη1/M+p/B1,1−δ1n2γ1n3expγ2n2Tγ1n3>1−δ1n2expγ2n2+1Tγ1n3>1, where γ2=rk0−m1¯/k1−b+a3M−a1/b1Mm<0. Let T¯=n2T+n3T, and it is deduced that there must be a t′∈t∗,t∗+T¯ such that xt′>m1¯; considering (49) with v4n1T+=x3n1T+, it is easy to get(53)v4t=1−δ3n−n1+lzn1T+−p1−1−δ3expBTexpBt−n1T+v4∗t,n−1T<t≤n+l−1T,1−δ3n−n1+l−1zn1T+−p1−1−δ3expB1TexpBt−n1T+v4∗t,n+l−1T<t≤nT,andn1+1≤n≤n1+n2+n3+1. So we can get v4t−v4∗t<M+pexpBt−n1T<η1 and zt≤v4t<v4∗t+η1 for t∗+n2T≤t≤t∗+T¯, which indicate that (51) holds for t∗+n2T≤t≤t∗+T¯. We can obtain xt∗+T¯≥xt∗+n2Tγ1n3 through a discussion similar to the first step. From the system (8), we have(54)x′t≥xtrk0−m1¯k1−b+a3M−a1b1Mm=γ2xt,t≠nT,t≠n+l−1T,xt+=1−δ1xt,t=n+l−1T.Integrating (54) on t∈t∗,t∗+n2T, then we get xt∗+n2T≥m1¯1−δ1n2expγ2n2T and xt∗+T¯≥m1¯1−δ1n2γ1n3expγ2n2T>m1¯; it is a contradiction with xt<m1¯. So xt>m1¯ holds for all t∈t1,t∗. Let t¯=inft≥t∗xt≥m1¯, then we can obtain xt≤m1¯,xt¯=m1¯ and xt≥m1¯1−δ1n2+n3expn2+n3γ2T≜m1 for t∗<t<t¯. Then, we derive that xt≥m1 holds for t∗<t<t¯. In addition, for t>t¯, since xt¯≥m1¯, the same discussion can be continued.Case 2: the case oft∗=n1+l−1T is discussed in Case 1. Now, assume that t∗≠n1+l−1T, which produces xt≥m1¯ for t1≤t<t∗ and x1t∗=m1¯. If t∗∈n1′+l−1T,n1′+lT,n1′∈N, which have two possible cases for t∈t∗,n1′+lT.Case 2(a): for allt∈t∗,n1′+lT, xt≤m1¯ is valid. By the same discussion with case 1, which exists t1′∈n1′+lT,n1′+lT+T¯ such that xt1′>m1¯; in order to simplify the proof, we will not elaborate on it. Let t˜=inft≥t∗xt>m1¯, then we have xt˜≤m1¯ for t∗≤t<t˜. We get xt≥m1¯1−δ1n2+n3+1expn2+n3γ2T≜m1 and xt≥m1 for t∈t∗,t˜. For t>t˜, the same arguments can be followed since x1t˜≥m1¯.Case 2(b): there existst∗<t<n1′+lT such that xt>m1¯. Let t˜=inft≥t∗xt>m1¯, then xt˜≤m1¯ for t∈t∗,t˜ and x1t˜=m1¯. When t∈t∗,t˜, (54) holds. Integrating (54) on t∗,t˜, we obtain xt≥xt∗expγ2t−t∗≥m1¯expγ2T>m1. Since xt≥m1¯ for t>t∗, the same arguments can be continued. Hence, xt≥m1 for all t≥t˜.
To sum up all the above discussions, we derive thatxt≥m1 holds for all t≥t˜. Using similar methods for yt, we can get yt≥m2 for all t≥t2. Obviously, qt≥m1 for all t>t˜1. Let m=maxm1,m2,m3; by Theorem 1 and Definition 2, we obtain that m<xt,yt,zt,qt<M. So, the system (8) proved to be permanent. □


## 4. Extension of Theoretical Results

In this section, we can see that Section3 of this paper can be extended to the results of the literature [26].From Theorem 2, we know if(55)rk0Tk1+a1c1mm2A1+A2<ln11−δ1,(56)−m1T+a2c2mm2B1+B2<ln11−δ2,hold, then the eradication periodic solution 0,0,z∗t,0 of species of xt and yt is globally asymptotically stable, where A,A1,A2,B1, and B2 are given in Theorem 2. If m=1,δii=1,2,3=0,l=1, inequality (55) is equivalent to(57)rk0Tk1+a1m2c1ln1−c1p1−exp−m2Tb11−exp−m2T+c1p<0.In addition, inequality (56) is equivalent to(58)−m1T+a2m2c2ln1−c2p1−exp−m2Tb21−exp−m2T+c2p<0.Remark 2.
Theorem 2.1 of the literature [26] can be obtained by combining inequalities (57) and inequalities (58).
From the above discussion, the following corollary can be obtained.Corollary 1.
Ifm=1,δii=1,2,3=0,l=1, then the eradication periodic solution 0,0,z∗t,0 of species of xt and yt of system (8) is locally asymptotically stable if(59)rk0Tk1+a1m2c1ln1−c1p1−exp−m2Tb11−exp−m2T+c1p<0,−m1T+a2m2c2ln1−c2p1−exp−m2Tb21−exp−m2T+c2p<0.From Theorem3, we know that if(60)rk0Tk1+a1c1mm2A1+A2>ln11−δ1,−m1T+a2c2mm2B1+B2>ln11−δ2,hold, then system (8) is permanent. Similar to the discussion of inequality (55)–(58), we can obtain the following corollary.Corollary 2.
Ifm=1,δii=1,2,3=0,andl=1, the system (8) is permanent provided that(61)rk0Tk1+a1m2c1ln1−c1p1−exp−m2Tb11−exp−m2T+c1p>0,−m1T+a2m2c2ln1−c2p1−exp−m2Tb21−exp−m2T+c2p>0.Remark 3.
Corollary1 and Corollary 2 are the same as Theorem 2.1 and Theorem 2.4 of the literature [26], respectively. Our theoretical results generalize the results of the literature [26].Remark 4.
After calculation and comparison, we can know that there are some mistakes in the form of Theorem 2.1, Theorem 2.3, and Theorem 2.4 of the literature [26].


## 5. Examples and Simulations

The effect of numerical experiments on this system is to help us visually observe the dynamic behavior of species changing with parameters, so some dynamical properties of the system (8) are showed by numerical simulation in this section, which is based on the parameters of the system in accordance with actual biological significance.Firstly, letr=0.6,k1=20,k0=15,a1=0.5,a2=0.5,a3=0.5,e1=0.8,e2=0.85,c1=0.6, c2=0.5,b1=1.5,b2=0.75,m1=0.1,m2=0.15, d=0.9,b=0.15,a=1.8,l=0.7,δ1=0.65,δ2=0.45,δ3=0.2,m=0.92,andp=10. By verification, these parameters satisfy biological significance and the conditions of Theorem 2, so the eradication periodic solution 0,0,z∗t,0 of species xt and yt is globally asymptotically stable under this numerical scenario. By the numerical experiment, Figures 1(a)–1(e) and Figure 2 show that the system is in asymptotic stability. Figures 1(a)–1(e) illustrate the existence of the solution 0,0,z∗t,0, and the dynamic behavior of species xt and yt with different initial values is described in Figure 2. In addition, let m=0.61 while keeping other parameters unchanged, we will be surprised to find that the system becomes persistent (see Figure 1(f)) through numerical experiments, so Figure 1(f) shows that the parameter m has a certain influence on the dynamic properties of the system (8).Figure 1
The eradication periodic solution0,0,z∗t,0 of species xt and yt with x0=0.3,y0=0.2,z0=0.3,andq0=0.3. (a)–(c) Time series of species xt, yt, and zt, respectively; (d) time series of three species xt,yt,andzt, which are in the same quadrant; (e) phase portrait of xt,yt,andzt; (f) a special case of dynamic behaviors as m=0.61.
(a)
(b)
(c)
(d)
(e)
(f)Figure 2
Dynamic behavior of0,0,z∗t,0 with different initial values, where the values of parameters are consistent with the settings in Figures 1(a) and 1(b) phase portrait of species xt and zt and yt and zt for different initial values, respectively.
(a)
(b)Secondly, letr=0.6,k1=19,k0=15,a1=0.8,a2=0.9,a3=0.2,e1=0.33,e2=0.1,c1=0.6, c2=0.33,b1=1.5,b2=0.75,m1=0.25,m2=0.25,d=0.9,b=0.119,a=1.8,l=0.6,δ1=0.1,δ2=0.35,δ3=0.25,m=0.9,andp=0.5.By verification, these parameters satisfy biological significance and the conditions of Theorem3, so the system (8) is permanent under this numerical scenario. By numerical simulation, the same conclusion is verified (see Figures 3(a)–3(d)). Similarly, let m=0.3 while keeping other parameters unchanged, we will be surprised to find that species xt and yt tend to extinct (see Figure 3(f)); i.e., the system has changed from a persistent situation to an extinct one, so we can see that the influence of parameter m on the system is effective.Figure 3
(a–c) Time series of persistence ofxt,yt,andzt, respectively, which initial values of x0=1, y0=0.1, z0=1, and q0=1; (d) time series of three species xt,yt,andzt in the same quadrant; (e) phase portrait of xt,yt,andzt; (f) special case of dynamic behaviors as m=0.34.
(a)
(b)
(c)
(d)
(e)
(f)Thirdly, by numerical analysis of bifurcation diagrams of impulsive periodT, the dynamic properties of the system (8) can be showed. Let x0=2,y0=2,z0=3, and q0=1. r=0.6,k1=20,k0=15,a1=0.8,a2=0.9,a3=0.2,e1=0.8, e2=0.1,c1=0.6,c2=0.5,b1=0.25,b2=0.75,m1=0.1,m2=0.15,d=0.25,b=0.15,a=1.8,l=0.7,δ1=0.65,δ2=0.45,δ3=0.25,m=0.9,andp=0.5. The bifurcation diagrams of impulsive period T is presented in Figure 4, whose range of values is [17, 23]. From Figures 4(a)–4(c), we can observe that impulsive period T has a great influence on the system. Since their bifurcation diagrams show many of the same phenomena, we will only explain the bifurcation diagram of xt in detail here. Figures 5 and 6 are detailed demonstrations of Figure 4(a). When T is small, the bifurcation diagram exhibits periodic oscillations of the system; for example, it has a periodic oscillation as T=16.5 (see Figure 6(a)). When T=17.85,19.68,20.46, and 20.62, bifurcation occurs, respectively. More concretely, if T>17.85, we can see from Figure 5(a) that the phenomenon of bifurcation diagram changes from period-doubling bifurcations’ cascade to cascade; if T>19.68, 4T period solution occurs (see Figure 5(a)), and a example is showed in Figure 6(b) T=19.7; if T>20.46, 8T period solution appears (see Figure 5(a)), and a example is showed in Figure 6(c) (T=20.51); if T>20.62, 16T period solution occurs(see Figure 5(a)). As T increased by more than 20.66, the system was in “crisis.” Period bifurcation reappears when T>21.22 (see Figure 5(b)). In the neighborhood of T=22.1, symmetry-breaking bifurcations occur after period bifurcation (see Figure 5(c)), which is a special simple bifurcation that produces multiple steady states. When T=22.61, a typical chaotic oscillation appears (see Figure 5(c)); after this situation, chaos appears such as Figure 6(d) T=22.8. That is to say, when T is less than 17.85, the harvesting and release of the system is frequent; when 17.85<T<20.66, period bifurcation appears; if 20.66<T<23, the dynamical properties of the system (8) become very complex; that is, the development of the system is unpredictable. All these results indicate that impulsive period T has a serious effect on the dynamic behavior of the system.Figure 4
Bifurcation diagrams of parameter (T). (a–c) Bifurcation diagram of xt,yt,andzt, respectively.
(a)
(b)
(c)Figure 5
The detailed demonstrations of Figure4(a). (a, b) The phenomenon of period bifurcation. (c) The phenomenon of crisis, symmetry-breaking bifurcations, and chaotic oscillation.
(a)
(b)
(c)Figure 6
(a) Periodic oscillation whenT=16.5; (b) 4T period solution when p=19.7; (c) 8T period solution when p=20.51; (d) chaos when T=22.8.
(a)
(b)
(c)
(d)Fourthly, we explore the bifurcation diagram of impulsive releasep. Let T=0.5, while other parameters are consistent with Figure 4, whose range of values is [0, 2.5] (see Figure 7). As we can see from Figure 5, if the amount of release p is small, the system (8) has periodic oscillations such as p=0.23 (see Figure 8(a)). If the releasing amount p of predator zt is comparatively large, then species xt and yt will die out. From Figure 7(b), we can see that predator yt will die out when T=1.77, the number of xt drops rapidly to zero when p>27.25, and this phenomenon shows the role of zt as an advanced predator, predator zt, control predator yt, and prey xt. If the release amount p is appropriate, there will be abundant phenomena in the system (8); for example, it produces chaos (see Figure 8(c)) when p=0.65, and it has 1T-period solution (see Figure 8(b)) when p=1.12. Furthermore, if the value of parameter m in Figure 8(b) is changed, different dynamic properties can be obtained such as a 4T-period solution when m=0.35 (Figure 9(a)), a 6T-period solution when m=0.31 (see Figure 9(b)), and a phenomenon of periodic oscillation when m=0.95 (see Figure 9(c)). All these results show that the amount of impulsive release p and mutual interference of predator m produces complex dynamic behavior to the system (8).Figure 7
(a–c) The bifurcation diagrams of impulsive releasep of xt,yt,andzt, respectively, with initial values (1, 1, 1, 1) and on [0, 2.5].
(a)
(b)
(c)Figure 8
(a) Periodic oscillations whenp=0.23; (b) 1T-period solution with p=1.12; (c) chaos with p=0.65.
(a)
(b)
(c)Figure 9
The value of parameterm in Figure 8(b) is changed: (a) 4T-period solution when m=0.35; (b) 6T-period solution when m=0.31; (c) periodic oscillation when m=0.95.
(a)
(b)
(c)Finally, the effect of distributed time delaya on the system (8) is explored. Several recent experiments have become deterministic evidence for chaos, and the results of these literatures show that the study of chaos is very meaningful [37–39]. Among them, the largest Lyapunov exponents provide a very effective diagnosis for the chaotic system, and the largest Lyapunov exponent is the average change rate of the two trajectories which are close in phase space as time goes on by exponential separation or polymerization. In any system, as long as there is a largest Lyapunov exponent greater than zero, it considers the system to be a chaotic system [38]. Rosenstein et al. [39] introduces the general method and process of calculating largest Lyapunov exponents. The corresponding largest Lyapunov exponents of a of system (8) is calculated and output in Figure 10(a), and an example of a strange attractor is given when a=1.72 (Figure 10(b)). All these results show that distributed time delay a has an important effect on the system.Figure 10
(a) The largest Lyapunov exponents of the system (8) with varying a in [6, 13]; (b) a strange attractor when a=1.72.
(a)
(b)


## 6. Conclusion

The dynamical analysis of a complex predator-prey system with mutual interference and the effect of time delay is presented in this paper, and numerical experiments of corresponding theoretical analysis are also presented. The periodic solution0,0,z∗t,0 is global asymptotically stable, which has been proved in Theorem 2. The persistence of system (8) was proved in Theorem 3. As special cases, Corollary 1 and Corollary 2 generalize the theorems in the literature [26], and we have obtained more abundant theoretical results and numerical experiments.After theoretical analysis and numerical simulation, we can know that impulsive periodT, release amount p, mutual interference parameters m, and distributed time delay a can produce different dynamical properties for system (8). For example, the solution 0,0,z∗t,0 and the persistence of system (8) are given in Figures 1 and 3, respectively. Specially, the dynamical effect of mutual interference parameters m is investigated. From Figures 1(d), 1(f), 3(d), 3(f), 8(b), and 9, we can know that different values of m can produce different dynamical behaviors with the system (8), so it is necessary for us to study the influence of mutual interference in such a system.In addition, the figure of largest Lyapunov exponents is obtained by computer simulation, which can prove the existence of chaotic phenomena in the system (8). Based on the above theoretical results and numerical experiments, we can know that the development of the system can be effectively controlled by choosing the appropriate impulsive period T, the advanced predator release amount p, and the reasonable harvest rate δi i=1,2,3. The comprehensive theoretical analysis and numerical simulation can give some guidance to the development of a real ecosystem and the balance of nature. Controlling the number of insect populations is critical because sudden outbreaks of insects can disrupt ecological balance and affect economic growth. Chemical control and biological control are very important methods for controlling agricultural pests. The usual measure for chemical control is spraying pesticides. Chemical control can kill pests quickly and effectively and reduce economic losses, but it also causes certain pollution to the environment. In order to reduce the environmental pollution caused by chemical control, biological control can be achieved by stocking or releasing natural enemies, but its effect is not obvious and rapid. Therefore, integrated control is a very effective means of balancing. The results obtained in this work provide a useful platform for pest control, fishery production, species breeding, and other fields.



---

*Source: 7629763-2020-01-13.xml*