# Global Stability Analysis of a Nonautonomous Stage-Structured Competitive System with Toxic Effect and Double Maturation Delays

Authors: Chao Liu; Yuanke Li
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We investigate a nonautonomous two-species competitive system with stage structure and double time delays due to maturation for two species, where toxic effect of toxin liberating species on nontoxic species is considered and the inhibiting effect is zero in absence of either species. Positivity and boundedness of solutions are analytically studied. By utilizing some comparison arguments, an iterative technique is proposed to discuss permanence of the species within competitive system. Furthermore, existence of positive periodic solutions is investigated based on continuation theorem of coincidence degree theory. By constructing an appropriate Lyapunov functional, sufficient conditions for global stability of the unique positive periodic solution are analyzed. Numerical simulations are carried out to show consistency with theoretical analysis.


---

## Body

## 1. Introduction

In recent years, many research efforts have been made on competitive Lotka-Volterra system with stage structure and time delay. By incorporating a constant time delay into single species model, a stage-structured model is proposed in the pioneering work [1], where time delay reflects a delayed birth of immature population and a reduced survival of immature population to their maturity. The model system takes the following form:
(1)x˙i(t)=αxm(t)-γxi(t)-αe-γτxm(t-τ),x˙m(t)=αe-γτxm(t-τ)-βxm2(t),
where xi(t) and xm(t) represent the immature population and mature population density at time t, respectively. α>0 denotes the birth rate of immature population; γ>0 stands for the death rate of immature population. β>0 is the death and overcrowding rate of mature population. τ denotes time of immature population to maturity. The term αe-γτxm(t-τ) represents the immature species which are born at time t-τ and survive at time t with immature death rate γ and therefore represents transformation of immature species to mature species. It is found that all ecologically relevant solutions tend to the positive equilibrium solution as time t→∞, and various aspects of the above proposed system including positivity and boundedness of solutions are discussed in [1].Zeng et al. propose a nonautonomous competitive two-species model with stage structure in one species in [2], where conditions of permanence are obtained. Furthermore, existence and asymptotic stability of periodic solution are proved under some assumptions if the proposed model turns out to be a periodic system. A two-species Lotka-Volterra type competition model with stage structure for both species is proposed and investigated in [3], where the individuals of each species are classified as immature and mature. By constructing a suitable Lyapunov function, sufficient conditions are derived for the global stability of nonnegative equilibria of the proposed model in the case of constant coefficients. Furthermore, a set of easily verifiable sufficient conditions are obtained for the existence of positive periodic solution when coefficients are assumed to be positively continuous periodic functions. In [4], there is a time delayed periodic system which describes the competition among mature populations. The evolutionary behavior of model system is analyzed and some sufficient conditions for competitive coexistence and exclusion are obtained.A nonautonomous competitive Lotka-Volterra system is studied in [5]; it reveals a computable necessary and sufficient condition for the system to be totally permanent when the growth rates have averages and the interaction coefficients are nonnegative constants. Along with this research, permanence for a class of competitive Lotka-Volterra systems is discussed in [6] which extends the work done in [5], and a computable necessary and sufficient condition is found for the permanence of all subsystems of the system and its small perturbation on the interaction matrix. In [7], a two-species competitive model with stage structure is discussed, and the dynamics of coupled system of semilinear parabolic equations with time delays are investigated, which show that the introduction of diffusion does not affect the permanence and extinction of the species, though the introduction of stage structure brings negative effect on it. In [8], sufficient conditions are obtained for the existence of a unique, globally attractive, strictly positive (componentwise), almost periodic solution of a nonautonomous, almost periodic competitive two-species model with a stage structure in one species. An example together with its numeric simulations shows the feasibility of our main results, which generalize the main results of Zeng et al. [2]. According to two types of well-known periodic single species population growth models with time delay, two corresponding periodic competitive systems with multiple delays are proposed in [9], and the same criteria for the existence and globally asymptotic stability of positive periodic solutions of the above two competitive systems are derived. In [10], a discrete periodic competitive model with stage structure is established, and some sufficient and realistic conditions are obtained for existence of a positive periodic solution of the proposed system. In [11], a periodic nonautonomous competitive stage-structured system with infinite delay is considered, where the adult members of n-species are in competition. For each of the n-species the model incorporates a time delay which represents the time from birth to maturity of that species. Infinite delay is introduced which denotes the influential effect of the entire past history of the system on the current competition interactions. By using the comparison principle, if the growth rates are sufficiently large, then the solutions are uniformly permanent. Then, by using Horn’s fixed point theorem, the existence of positive periodic solution of the system with finite delay is discussed. Finally, it is proved that even the system with infinite delay admits a positive periodic solution.In [12], a nonautonomous predator-prey system with discrete time delay is studied, where there is epidemic disease in the predator. By using some techniques of the differential inequalities and delay differential inequalities, the permanence of system is discussed under some appropriate conditions. When all the coefficients of the system are periodic, the existence and global attractivity of the positive periodic solution are studied by Mawhin’s continuation theorem and constructing a suitable Lyapunov functional. Furthermore, when the coefficients of the system are not absolutely periodic but almost periodic, sufficient conditions are also derived for the existence and asymptotic stability of the almost periodic solution. In [13], general n-species nonautonomous Lotka-Volterra competitive systems with pure-delays and feedback controls are discussed. New sufficient conditions, for which a part of the n-species remains permanent, are established by applying the method of multiple Lyapunov functionals and introducing a new analysis technique.By utilizing Brouwer fixed point theorem and constructing a suitable Lyapunov function, the periodic solution and global stability for a nonautonomous competitive Lotka-Volterra diffusion system are investigated in [14]; it can be found that the system has a unique periodic solution which is globally stable under some appropriate conditions. In [15], a delay differential equation model for the interaction between two species is investigated. The maturation delay for each species is modelled as a distribution, to allow for the possibility that individuals may take different amount of time to maturity. Positivity and boundedness of the solutions are studied, and global stability is analyzed for each equilibrium. A Lotka-Volterra competitive system with infinite delay and feedback controls is proposed in [16]. By using Mawhin’s continuation theorem of coincidence degree theory, an impulsive nonautonomous Lotka-Volterra predator-prey system with harvesting terms is investigated in [17]. Some sufficient conditions for the existence of multiple positive almost periodic solutions for the system under consideration are discussed. Furthermore, existence of multiple positive almost periodic solutions to other types of population systems can be studied by using the same method obtained in this paper. By using the method of multiple Lyapunov functionals and by developing a new analysis technique, some sufficient conditions are obtained that guarantee that some of the n species are driven to extinction. A three-dimensional nonautonomous competitive Lotka-Volterra system is considered in [18]; it is shown that if the growth rates are positive, bounded, and continuous functions, and the averages of the growth rates satisfy certain inequalities, then any positive solution has the property that one of its components vanishes. In [19], an almost periodic multispecies Lotka-Volterra mutualism system with time delays and impulsive effects is investigated. By using the theory of comparison theorem and constructing a suitable Lyapunov functional, sufficient conditions which guarantee the existence and uniqueness and global asymptotical stability of almost periodic solution of this system are obtained.It is well known that the effect of toxins on ecological systems is an important issue from mathematical and experimental points of view [20, 21]. The first mathematical model to represent the toxic liberating interaction between two competing species is introduced by Maynard Smith [22]. The model is based upon a two-species Lotka-Volterra competition model with an additional term to take into account the effect of toxic substances released by one species to another, which takes the following form:
(2)N˙1(t)=N1(t)[α1-β1N1(t)-c1N2(t)-ρ1N1(t)N2(t)],N˙2(t)=N2(t)[α2-β2N2(t)-c2N1(t)-ρ2N1(t)N2(t)],
where N1(t) and N2(t) represent the density of two competing species at time t, respectively. α1 and α2 denote the birth rate of N1(t) species and N2(t) species, respectively. β1 and β2 are the rate of intraspecific competition term for the first and second species, respectively. c1 and c2 stand for the rate of interspecific competition, respectively. ρ1 and ρ2 represent the toxic inhibition rate for the first species by the second species and vice versa. By considering that ρ denotes the rate of toxic inhibition for the nontoxic species N1(t) released by the toxin liberating species N2(t) and all other parameters share the same biological interpretations mentioned in model system (2), work done in [22] is extended in [23] and the generalized model system is as follows:
(3)N˙1(t)=N1(t)[α1-β1N1(t)-c1N2(t)-ρN1(t)N22(t)],N˙2(t)=N2(t)[α2-β2N2(t)-c2N1(t)],
where the toxic substance producing action follows the mathematical term ρN12(t)N22(t) [23].It should be noted that models of the persistence and extinction of a population or community in a polluted environment have been investigated in [23]. But all of those papers have a basic assumption that the capacity of the environment is so large that the change of toxicant in the environment that comes from uptake and egestation by the organisms can be neglected. This assumption is not made in [24, 25], some sufficient conditions on persistence or extinction of a population have been obtained, and the threshold between the two has also been obtained for most situations. In [26, 27], there are modified delay differential equation models of the growth of two species of plankton having competitive and allelopathic effects on each other. By using the continuation theorem of coincidence degree theory, a set of easily verifiable sufficient conditions are obtained for the existence of positive periodic solutions for this model. Recently, some discussions and investigations of the nonautonomous competitive model with toxic effects are made. A periodic competitive stage-structured Lotka-Volterra model with the effects of toxic substances is investigated in [28]. It is shown that toxic substances play an important role in the extinction of species. A set of sufficient conditions guarantee that one of the components is driven to extinction while the other is globally attractive. The dynamical behavior of a two-species competitive system affected by toxic substances is investigated in [21], where each species produces a substance toxic to the other species. Boundedness and local and global stabilities are also addressed. It should be noted that toxic interaction follows the mathematical term suggested in model system (2) and each mature individual produces a substance toxic to the other mature individuals only when the other mature individual is present, and the immature individual is not affected by the toxicant [21]. However, to the author’s best knowledge, dynamical behavior and stability analysis of nonautonomous stage-structured competitive system with toxin liberating species and nontoxic species have not been investigated. Generally speaking, it takes some time for a species to reach maturity to produce the toxicant; then toxin liberating mature individual produces a substance toxic to the nontoxic mature individuals only. The inhibiting effect is zero in absence of either species, and the immature individual of each species is not affected by the toxicant. Furthermore, the species compete each other for the limited life resource within closed environment, but this competition only happens among the mature individuals and does not involve the immature individuals. Consequently, it is necessary to investigate the dynamic effect of stage structure and toxic effect on the population dynamics of two-species competitive system with toxin liberating species and nontoxic species.The rest section of this paper is organized as follows: a nonautonomous two-species competitive model is established in the second section. Stage structure and maturation delay for each species are introduced, and toxic effect of toxin liberating species on nontoxic species is considered. In the third section, qualitative analyses are performed to investigate the effect of stage structure and toxic substances on the dynamical behavior of two-species competitive model system. The positivity and boundedness of solutions are analytically studied. By utilizing some comparison arguments, an iterative technique is proposed to discuss permanence of the species within competitive system. Furthermore, existence of positive periodic solutions is considered based on continuation theorem of coincidence degree theory. By constructing an appropriate Lyapunov functional, sufficient conditions for global stability of the unique positive periodic solution are analyzed. Numerical simulations are provided to support the theoretical findings obtained in this paper. Finally, this paper ends with a conclusion.


## 2. Model Formulation

In this paper, the effect of stage structure and toxic substances on the dynamical behavior of two-species competitive model system is investigated under the following five hypotheses, which are given as follows.(H1)
Two competing species, that is, nontoxic species and toxin liberating species, are considered in this paper. It is assumed that each species is divided into two-stage groups, and the immature and mature individuals are divided by a fixed period.x1(t) and y1(t) represent immature population density of nontoxic species and toxin liberating species at time t, respectively; x2(t) and y2(t) denote mature population density of nontoxic species and toxin liberating species at time t, respectively.(H2)
ω-periodic continuous functions α1(t)>0 and α2(t)>0 denote the birth rate of immature population of nontoxic species and toxin liberating species at time t, respectively. ω-periodic continuous functions γ1(t)>0 and γ2(t)>0 stand for the death rate of immature population of nontoxic species and toxin liberating species at time t, respectively. ω-periodic continuous functions β1(t)>0 and β2(t)>0 are the death and overcrowding rate of mature population of nontoxic species and toxin liberating species at time t, respectively.(H3)
τ1 denotes time of immature nontoxic species to maturity. The term α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1) represents the immature nontoxic species which are born at time t-τ1 and survive at time t with immature death rate. τ2 denotes time of immature toxin liberating species to maturity. The term α2(t-τ2)e-∫t-τ2t‍γ2(s)dsy2(t-τ2) represents the immature toxin liberating species which are born at time t-τ2 and survive at time t with immature death rate.(H4)
For toxin liberating species, it takes some time to attain its level of maturity to produce the toxic substances, and toxin liberating mature individual produces a substance toxic to the nontoxic mature individuals only. The inhibiting effect is zero in absence of either species, and the immature individual of each species is not affected by the toxicant. Based on model system (3), the toxic effect released by toxin liberating species on nontoxic species is described by the mathematical term ρ(t)x22(t)y22(t), where toxic inhibition rate is represented by an ω-periodic continuous function ρ(t)>0.(H5)
Nontoxic species and toxin liberating species compete each other for the common resource within closed environment, but this competition only happens among the mature individuals and does not involve the immature individuals.ω-periodic continuous function c1(t)>0 represents interspecific competition rate for the mature nontoxic species by the mature toxin liberating species, and ω-periodic continuous function c2(t)>0 represents interspecific competition rate for the mature toxin liberating species by the mature nontoxic species.Based on hypotheses (H1)–(H5), a nonautonomous stage-structured competitive model with toxic effect and double maturation delays is established as follows:(4)x˙1(t)=α1(t)x2(t)-γ1(t)x1(t)-α1(t-τ1)e-∫t-τ1tγ1(s)dsx2(t-τ1),x˙2(t)=α1(t-τ1)e-∫t-τ1tγ1(s)dsx2(t-τ1)-β1(t)x22(t)-c1(t)x2(t)y2(t)-ρ(t)x22(t)y22(t),y˙1(t)=α2(t)y2(t)-γ2(t)y1(t)-α2(t-τ2)e-∫t-τ2tγ2(s)dsy2(t-τ2),y˙2(t)=α2(t-τ2)e-∫t-τ2tγ2(s)dsy2(t-τ2)-β2(t)y22(t)-c2(t)x2(t)y2(t).In this paper, model system (4) is investigated with the following initial conditions:
(5)xi(θ)=ϕi(θ)>0,-τ1≤θ≤0,i=1,2,yi(θ)=ψi(θ)>0,-τ2≤θ≤0,i=1,2.For the continuity of the initial conditions, it is required that(6)x1(0)=∫-τ10‍α1(θ)ϕ2(θ)e∫0θ‍γ1(s)dsdθ,y1(0)=∫-τ20‍α2(θ)ψ2(θ)e∫0θ‍γ2(s)dsdθ.


## 3. Qualitative Analysis of Model System

In this section, qualitative analysis of the nonautonomous model system (4) is performed, which is utilized to discuss dynamic effect of toxic effect and maturation delay on population dynamics. The positivity and boundedness of solutions are analytically studied. By utilizing some comparison arguments, an iterative technique is proposed to discuss permanence of the species within competitive system. Furthermore, existence of positive periodic solutions is investigated based on continuation theorem of coincidence degree theory. By constructing an appropriate Lyapunov functional, sufficient conditions for global stability of the unique positive periodic solution are analyzed.Some mathematical notations are adopted for convenience of the following statement:(7)fL=min⁡t∈[0,ω]|f(t)|,fM=max⁡t∈[0,ω]|f(t)|,
where f(t) is a ω-periodic continuous function.


### 3.1. Permanence of Solutions

Theorem 1.
Solutions of model system (4) with initial conditions (5) and (6) are positive for all t>0.Proof.
Firstly, we show thatx2(t)>0 for all t>0. Otherwise, if it is false, since x2(t)>0 for all t∈[-τ1,0], then it can be derived that there exists a t1>0 such that x2(t1)=0.
Definet0=inf⁡{t>0∣x2(t)=0}. According to the definition of t0, it can be obtained that
(8)x˙2(t0)≤0.
It follows from the second equation of model system (4) that
(9)x˙2(t0)={α1(t0-τ1)e-∫t0-τ1t0‍γ1(s)dsϕ2(t0-τ1)>0,0≤t0≤τ1,α1(t0-τ1)e-∫t0-τ1t0‍γ1(s)dsx2(t0-τ1)>0,t>τ1,
and it is easy to show that x˙2(t0)>0, which is a contradiction to (8). Hence, x2(t)>0 for all t>0.
By a direct computation, it follows from the first equation of model system (4) that
(10)x1(t)=∫t-τ1t‍α1(s)e∫ts‍γ1(m)dsx2(s)ds,
since x2(t)>0 for all t>0; it is easy to show that x1(t)>0 for all t>0 based on (10).
By utilizing the similar proof, it can be obtained thaty1(t)>0 and y2(t)>0 for all t>0. Consequently, solutions of model system (4) with initial conditions (5) and (6) are positive for all t>0.Theorem 2.
Solutions of model system (4) with initial conditions (5) and (6) are ultimately bounded.Proof.
Letw(t)=x1(t)+x2(t)+y1(t)+y2(t), where (x1(t),x2(t),y1(t),y2(t)) is an arbitrary positive solution of model system (4) with the initial conditions (5) and (6).
Calculating the derivative ofw(t) along the solution of model system (4) gives that
(11)w˙(t)≤(α1(t)+γ1(t))x2(t)-γ1(t)(x1(t)+x2(t))-β1(t)x22(t)+(α2(t)+γ2(t))y2(t)-γ2(t)(y1(t)+y2(t))-β2(t)y22(t)≤(α1M+γ1M)x2(t)-β1Lx22(t)-r1L(x1(t)+x2(t))+(α2M+γ2M)y2(t)-β2Ly22(t)-r2L(y1(t)+y2(t))≤-γLw(t)+(α1M+γ1M)x2(t)-β1Lx22(t)+(α2M+γ2M)y2(t)-β2Ly22(t)≤-γLw(t)+(α1M+γ1M)24β1L+(α2M+γ2M)24β2L,
where γL=min⁡{γ1L,γ2L}.
By using the standard comparison principle [20], it follows from (11) that
(12)w(t)≤β2L(α1M+γ1M)2+β1L(α2M+γ2M)24γLβ1Lβ2L,
which implies that any solution of model system (4) with initial conditions (5) and (6) is ultimately bounded.Lemma 3 (see [29]).
Consider the following differential equation:(13)x˙(t)=ax(t-τ)-bx(t)-cx2(t),
where a,b,c,τ>0 and x(t)>0 for -τ≤t≤0; we have that(i)
ifa>b, then lim⁡t→+∞x(t)=(a-b)/c;(ii)
ifa<b, then lim⁡t→+∞x(t)=0.Lemma 4 (see [29]).
Consider the following differential equation:(14)x˙(t)=dx(t-σ)-ex2(t),
where d,e,σ>0 and x(t)>0 for -σ≤t≤0; we have
(15)lim⁡t→+∞x(t)=de.Definition 5 (see [30]).
Model system(16)X˙(t)=f(t,Xt(θ)),
where t≥0, θ∈[-τ,0], X∈Rn. Model system (16) is said to be permanent if, for any solution X(t,ϕ), there exists a constant m>0 and T=T(ϕ) such that X(t)>m for all t>T.Definition 6 (see [30]).
The domainD∈Cn is said to be an ultimately bounded domain, if D is a closed, bounded subset of Cn, and there exists a constant T=T(ϕ) such that Xt(θ)∈D for all t>T.Theorem 7.
Ifα1Lβ2L>c1Mα2M and α2Lβ1L>c2Mα1M, then model system (4) is permanent with initial conditions (5) and (6).Proof.
According to the second equation of model system (4) and Theorem 1, we get that
(17)x˙2(t)≤α1Me-γ1Lτ1x2(t-τ1)-β1Lx22(t).
By virtue of Lemma4 and (17), there exists a positive time T1 such that, for sufficiently small ϵ>0 and t≥T1, it yields
(18)x2(t)≤α1Me-γ1Lτ1β1L+ϵ:=M2(1).
By rearranging (10), it can be obtained that
(19)x1(t)=e-∫0t‍γ1(s)ds∫t-τ1t‍α1(s)e∫0s‍γ1(m)dmx2(s)ds.
For anyt≥T1, it follows from (18) and (19) that
(20)x1(t)≤a1MM2(1)(1-e-γ1Mτ1)γ1L≔M1(1).
Based on the fourth equation of model system (4) and Theorem 1, it can be obtained that
(21)y˙2(t)≤α2Me-γ2Lτ2y2(t-τ2)-β2Ly22(t)
holds for t≥T1.
By virtue of Lemma4 and (21), there exists T2>T1 such that, for sufficiently small ϵ>0 and t≥T2, it yields
(22)y2(t)≤α2Me-γ2Lτ2β2L+ϵ≔M4(1).
By direct computing, it follows from the third equation of model system (4) that
(23)y1(t)=e-∫0t‍γ2(s)ds∫t-τ2t‍α2(s)e∫0s‍γ2(m)dmy2(s)ds.
For anyt≥T2, it follows from (22) and (23) that
(24)y1(t)≤α2MM4(1)(1-e-γ2Mτ2)γ2L≔M3(1).
Furthermore, it follows from the second equation of model system (4) that
(25)x˙2(t)≥α1Le-γ1Mτ1x2(t-τ1)-(β1M+ρ1M(M4(1))2)x22(t)-c1MM4(1)x2(t).
Based on Lemma3 and (25), there exists T3>T2 and for any t≥T3 and sufficiently small ϵ>0,
(26)x2(t)≥α1Le-γ1Mτ1-c1MM4(1)β1M+ρ1M(M4(1))2-ϵ≔m2(1)
holds provided that
(27)α1Le-γ1Mτ1>c1MM4(1).
According to (19), for any t≥T3, we get that
(28)x1(t)≥α1L(1-e-γ1Lτ1)m2(1)γ1M≔m1(1).
For anyt≥T3, it follows from the fourth equation of model system (4) that
(29)y˙2(t)≥α2Le-γ2Mτ2y2(t-τ2)-β2My22(t)-c2MM2(1)y2(t).
Based on Lemma3 and (29), there exists T4>T3 and for any t≥T4 and sufficiently small ϵ>0,
(30)y2(t)≥α2Le-γ2Mτ2-c2MM2(1)β2M-ϵ≔m4(1)
holds provided that
(31)α2Le-γ2Mτ2>c2MM2(1).
According to (23) and (30), it can be obtained that
(32)y1(t)≥α2Lm4(1)(1-e-γ2Lτ2)γ2M≔m3(1).
According to the second equation of model system (4), we get that
(33)x˙2(t)≤α1Me-γ1Lτ1x2(t-τ1)-(β1L+ρL(m4(1))2)x22(t)-c1Lm4(1)x2(t).
By virtue of Lemma3 and (33), there exists T5>T4 such that, for sufficiently small ϵ>0 and t≥T5, it yields
(34)x2(t)≤α1Me-γ1Lτ1-c1Lm4(1)β1L+ρL(m4(1))2+ϵ≔M2(2),
which holds provided that
(35)α1Me-γ1Lτ1>c1Lm4(1).
For anyt≥T5, it follows from (19) and (34) that
(36)x1(t)≤α1MM2(2)(1-e-γ1Mτ1)γ1L≔M1(2).
Based on the fourth equation of model system (4), it can be obtained that
(37)y˙2(t)≤α2Me-γ2Lτ2y2(t-τ2)-β2Ly22(t)-c2Lm2(1)y2(t)
holds for t≥T5.
By virtue of Lemma3 and (37), there exists T6>T5 such that, for sufficiently small ϵ>0 and t≥T6, it yields
(38)y2(t)≤α2Me-γ2Lτ2-c2Lm2(1)β2L+ϵ≔M4(2),
which holds provided that
(39)α2Me-γ2Lτ2>c2Lm2(1).
For anyt≥T6, it follows from (23) and (38) that
(40)y1(t)≤α2MM4(2)(1-e-γ2Mτ2)γ2L≔M3(2).
Furthermore, it follows from the second equation of model system (4) that
(41)x˙2(t)≥α1Le-γ1Mτ2x2(t-τ1)-(β1M+ρM(M4(2))2)x22(t)-c1MM4(2)x2(t).
Based on Lemma3 and (41), there exists T7>T6 and, for any t≥T7 and sufficiently small ϵ>0,
(42)x2(t)≥α1Le-γ1Mτ1-c1MM4(2)β1M+ρM(M4(2))2-ϵ≔m2(2)
holds provided that
(43)α1Le-γ1Mτ1>c1MM4(2).
According to (19) and (42), for any t≥T7, we get that
(44)x1(t)≥α1L(1-e-γ1Lτ1)m2(2)γ1M≔m1(2).
For anyt≥T7, it follows from the fourth equation of model system (4) that
(45)y˙2(t)≥α2Le-γ2Mτ2y2(t-τ2)-β2My22(t)-c2MM2(2)y2(t).
Based on Lemma3 and (45), there exists T8>T7 and, for any t≥T8 and sufficiently small ϵ>0,
(46)y2(t)≥α2Le-γ2Mτ2-c2MM2(2)β2M-ϵ≔m4(2)
holds provided that
(47)α2Le-γ2Mτ2>c2MM2(2).
According to (23) and (46), it can be obtained that
(48)y1(t)≥α2Lm4(2)(1-e-γ2Lτ2)γ2M≔m3(2).
By using simple computation, it is easy to show that six inequalities (27), (31), (35), (39), (43), and (47) hold if the following two inequalities α1Lβ2L>c1Mα2M and α2Lβ1L>c2Mα1M hold.
Furthermore, eight sequences will be obtained by repeating the discussion in this manner, which are given as follows:(49)M1(n+1)=α1MM2(n+1)(1-e-γ1Mτ1)γ1L,M2(n+1)=α1Me-γ1Lτ1-c1Lm4(n)β1L+ρL(m4(n))2+ϵ,M3(n+1)=α2MM4(n+1)(1-e-γ2Mτ2)γ2L,M4(n+1)=α2Me-γ2Lτ2-c2Lm2(n)β2L+ϵ,m1(n+1)=α1L(1-e-γ1Lτ1)m2(n+1)γ1M,m2(n+1)=α1Le-γ1Mτ1-c1MM4(n+1)β1M+ρM(M4(n+1))2-ϵ,m3(n+1)=α2Lm4(n+1)(1-e-γ2Lτ2)γ2M,m4(n+1)=α2Le-γ2Mτ2-c2MM2(n+1)β2M-ϵ.
It is easy to show thatMi(n)>0 and the sequences {Mi(n)} (i=1,2,3,4) are decreasing as n increases, which implies that lim⁡n→∞Mi(n)=Mi* exists; furthermore, it is easy to show that mi(n)<Mi(n) and the sequences {mi(n)} (i=1,2,3,4) are increasing as n increases, which implies that lim⁡n→∞mi(n)=mi* exists. Consequently, it follows from (49) that
(50)M1*=α1MM2*(1-e-γ1Mτ1)γ1L,M2*=α1Me-γ1Lτ1-c1Lm4*β1L+ρLm4*2,M3*=α2MM4*(1-e-γ2Mτ2)γ2L,M4*=α2Me-γ2Lτ2-c2Lm2*β2L,m1*=α1L(1-e-γ1Lτ1)m2*γ1M,m2*=α1Le-γ1Mτ1-c1MM4*β1M+ρMM4*2,m3*=α2Lm4*(1-e-γ2Lτ2)γ2M,m4*=α2Le-γ2Mτ2-c2MM2*β2M.
Based on Definition5 and (50), it can be concluded that model system (4) is persistent if α1Lβ2L>c1Mα2M and α2Lβ1L>c2Mα1M hold.


### 3.2. Existence of Positive Periodic Solutions

Definition 8 (see [31]).
LetL:Dom⁡L⊂X→Y be a linear mapping and let N:X→Y be a continuous mapping, where X and Y are real Banach spaces. If dimKerL=codimImL<+∞ and Im⁡L is closed in Y, then L is called a Fredholm mapping of index zero.IfL is Fredholm mapping of index zero and there exist continuous projectors P:X→X and Q:Y→Y such that
(51)Im⁡P=Ker⁡L,Im⁡L=Ker⁡Q=Im⁡(I-Q),
then restriction Lp of L to Dom⁡L∩Ker⁡P:(I-P)X→Im⁡L is invertible.Definition 9 (see [31]).
Denote the inverse ofLp by Kp. Supposing that Ω is an open bounded subset of X, if QN(Ω-) is bounded and Kp(I-Q)N:Ω→X is compact, then the mapping N is called L-compact on Ω-. Since Im⁡Q is isomorphic to Ker⁡L, there exists an isomorphism J:Im⁡Q→Ker⁡L.Lemma 10 (see [31]).
LetΩ⊂X be an open bounded set, let L be a Fredholm mapping of index zero, and let N be L-compact on Ω-. If the following three conditions hold: (i)
Lx≠λNx for any λ∈(0,1) and x∈∂Ω∩Dom⁡L,(ii)
QNx≠0 for any x∈∂Ω∩Ker⁡L,(iii)
deg⁡{JQN,Ω∩Ker⁡L,0}≠0,then Lx=Nx has at least one solution in Ω-∩Dom⁡L.Theorem 11.
Ifα1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2, α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ2, then model system (4) with initial conditions (5) and (6) has at least one positive ω-periodic solution.Proof.
Consider the subsystem of model system (4):
(52)x˙2(t)=α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1)-β1(t)x22(t)-c1(t)x2(t)y2(t)-ρ(t)x22(t)y22(t),y˙2(t)=α2(t-τ2)e-∫t-τ2t‍γ2(s)dsy2(t-τ2)-β2(t)y22(t)-c2(t)x2(t)y2(t).
Letu1(t)=ln⁡[x2(t)], u2(t)=ln⁡[y2(t)].
By substitutingu1(t) and u2(t) into (52), it can be obtained that
(53)u˙1(t)=α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)-u1(t)-β1(t)eu1(t)-c1(t)eu2(t)-ρ(t)eu1(t)+2u2(t),u˙2(t)=α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)-u2(t)-β2(t)eu2(t)-c2(t)eu1(t).
It should be noted that if model system (53) has one ω-periodic solution (u1*(t),u2*(t))T, then (x2*(t),y2*(t))T=(eu1*(t),eu2*(t))T is a positive ω-periodic solution of model system (52).
In order to utilize Lemma10 in a straightforward manner, we define
(54)X=Y={(u1(t),u2(t))Ti∈C(R,R2):ui(t+ω)=ui(t),i=1,2(u1(t),u2(t))T},∥(u1(t),u2(t))T∥=max⁡t∈[0,ω]|u1(t)|+max⁡t∈[0,ω]|u2(t)|,
where |·| denotes the Euclidean norm; it is easy to show that both X and Y are Banach spaces with the norm ∥·∥; then define
(55)Dom⁡L∩X⟶X,L(u1(t),u2(t))T=(du1(t)dt,du2(t)dt)T,
where Dom⁡L={(u1(t),u2(t))T∈C(R,R2)}, N[u1u2]=[f1(t)f2(t)], and
(56)f1(t)=α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)-u1(t)-β1(t)eu1(t)-c1(t)eu2(t)-ρ(t)eu1(t)+2u2(t),f2(t)=α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)-u2(t)-β2(t)eu2(t)-c2(t)eu1(t).
Furthermore, we define(57)P[u1u2]=Q[u1u2]=[1ω∫0ω‍u1(t)dt1ω∫0ω‍u2(t)dt],[u1u2]∈X=Y.
According to the above definitions, it is not difficult to verify thatKer⁡L={x∣x∈X,x=h,h∈R2}, Im⁡L={y∈Y∣∫0ω‍y(t)dt=0} are closed in Y, dimKerL=codimImL=2, and both P and Q are continuous projectors such that Im⁡P=Ker⁡L and Ker⁡Q=Im⁡L=Im⁡(I-Q).
Based on the above analysis, it can be obtained thatL is a Fredholm mapping of index zero.
Furthermore, the inverseKp:Im⁡L→Dom⁡L∩Ker⁡P of Lp exists and takes the following form:
(58)Kp(y)=∫0t‍y(s)ds-1ω∫0ω‍∫0t‍y(s)dsdt.
Hence,QN:X→Y and Kp(I-Q)N:X→X can be defined as follows, respectively,
(59)QNx=[1ω∫0ω‍f1(t)dt1ω∫0ω‍f2(t)dt],Kp(I-Q)Nx=∫0t‍Nx(s)ds-1ω∫0ω‍∫0t‍Nx(s)dsdt-(tω-12)∫0ω‍Nx(s)ds.
It is easy to show thatQN and Kp(I-Q)N are continuous. In order to facilitate the proof based on Lemma 10, we also need to find an appropriate open and bounded subset Ω, which can be found by the following two steps.
Step  1. According to the operator equation Lx=λNx for λ∈(0,1), the upper and lower bound of u1(t) and u2(t) will be estimated as follows:
(60)du1(t)dt=λf1(t),du2(t)dt=λf2(t),
where f1(t) and f2(t) have been defined in (56).
Suppose that(u1(t),u2(t))T∈X is a solution of model system (60) for some λ∈(0,1). By integrating (60) over the interval [0,ω], it can be obtained that
(61)∫0ω‍α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)-u1(t)dt=∫0ω‍β1(t)eu1(t)+c1(t)eu2(t)+ρ(t)eu1(t)+2u2(t)dt,(62)∫0ω‍α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)-u2(t)dt=∫0ω‍β2(t)eu2(t)+c2(t)eu1(t)dt.
Based on definition(u1(t),u2(t))T∈X, there exist ξi,ηi∈[0,ω] such that
(63)ui(ξi)=min⁡t∈[0,ω]ui(t),ui(ηi)=max⁡t∈[0,ω]ui(t),i=1,2.
Multiplying the first equation of (60) by eu1(t) and integrating it over [0,ω] give that
(64)∫0ω‍α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)dt=∫0ω‍β1(t)e2u1(t)+c1(t)eu1(t)+u2(t)hhhhhh+ρ(t)e2(u1(t)+u2(t))dt.
It follows from (64) that
(65)β1L∫0ω‍e2u1(t)dt+c1L∫0ω‍eu1(t)+u2(t)dt+ρL∫0ω‍e2(u1(t)+u2(t))dt≤α1Me-γ1Lτ1∫0ω‍eu1(t)dt.
On the other hand, by using the inequality,(66)(∫0ω‍eu1(t)dt)2≤ω∫0ω‍e2u1(t)dt.
Based on (65) and (66), it can be obtained that
(67)β1L(∫0ω‍eu1(t)dt)2≤ωβ1L∫0ω‍e2u1(t)dt≤ωα1Me-γ1Lτ1∫0ω‍eu1(t)dt,
which derives that
(68)β1L∫0ω‍eu1(t)dt≤ωα1Me-γ1Lτ1,u1(ξ1)≤ln⁡α1Me-γ1Lτ1β1L.
It follows from (60) and (68) that
(69)∫0ω‍|u1′(t)|dt<2∫0ω‍β1(t)eu1(t)≤2ωα1Mβ1Me-γ1Lτ1β1L.
According to (68) and (69), it can be obtained that
(70)u1(t)≤u1(ξ1)+∫0ω|u1′(t)|dt≤ln⁡α1Me-γ1Lτ1β1L+2ωα1Mβ1Me-γ1Lτ1β1L.
Multiplying the second equation of (60) by eu2(t) and integrating it over [0,ω] give that
(71)∫0ω‍α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)dt=∫0ω‍β2(t)e2u2(t)+c2(t)eu1(t)+u2(t)dt.
It follows from (71) that
(72)β2L∫0ω‍e2u2(t)dt+c2L∫0ω‍eu1(t)+u2(t)dt≤α2Me-γ2Lτ2∫0ω‍eu2(t)dt.
Based on (66) and (72), it can be obtained that
(73)β2L(∫0ω‍eu1(t)dt)2≤ωα2Me-γ2Lτ2∫0ω‍eu2(t)dt,
which derives that
(74)∫0ω‍eu2(t)dt≤ωα2Me-γ2Lτ2β2L,u2(ξ2)≤ln⁡α2Me-γ2Lτ2β2L.
It follows from (60) and (74) that
(75)∫0ω‍|u2′(t)|dt<2∫0ω‍β2(t)eu2(t)dt≤2ωα2Mβ2Me-γ2Lτ2β2L.
According to (74) and (75), it can be obtained that
(76)u2(t)≤u2(ξ2)+∫0ω‍|u2′(t)|dt≤ln⁡α2Me-γ2Lτ2β2L+2ωα2Mβ2Me-γ2Lτ2β2L.
It should be noted that(77)∫0ω‍α1(t)e-∫tt+τ1‍γ1(s)dseu1(t)dt=∫0ω‍α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)dt.
Based on (64), it can be obtained that
(78)(α1Le-γ2Mτ1-c1Mα2Me-γ2Lτ2β2L)∫0ω‍eu1(t)dt≤(β1M+2ρMα2Me-γ2Lτ2β2L)∫0ω‍e2u1(t)dt,
which derives that
(79)∫0ω‍eu1(t)dt≥ω(α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2)β1Mβ2L+2α2MρMe-γ2Lτ2,u1(η1)≥ln⁡α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2β1Mβ2L+2α2MρMe-γ2Lτ2
hold provided that α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2.
According to (69) and (79), it can be obtained that
(80)u1(t)≥u1(η1)-∫0ω‍|u1′(t)|dt≥ln⁡α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2β1Mβ2L+2α2MρMe-γ2Lτ2-2ωα1Mβ1Me-γ2Lτ2β1L.
By virtue of (70) and (80), if α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2, then
(81)max⁡t∈[0,ω]|u1(t)|<max⁡{|ln⁡α1Me-γ1Lτ1β1L|+2ωα1Mβ1Me-γ1Lτ1β1L,|ln⁡α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2β1Mβ2L+2α2MρMe-γ2Lτ2|+2ωα1Mβ1Me-γ1Lτ1β1L}≔H1.
Similarly, it is easy to show that(82)∫0ω‍α2(t)e-∫tt+τ2‍γ2(s)dseu2(t)dt=∫0ω‍α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)dt,
which derives that
(83)u2(η2)≥ln⁡α2Lβ1Le-γ2Mτ2-α1Mc2Me-γ1Lτ1β1Lβ2M
holds provided that α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ1.
According to (75) and (83), it is derived that
(84)u2(t)≥u2(η2)-∫0ω‍|u2′(t)|dt≥ln⁡α2Lβ1Le-γ2Mτ2-α1Mc2Me-γ1Lτ1β1Lβ2M-2ωα2Mβ2Me-γ2Lτ2β2L.
By virtue of (76) and (84), if α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ2, then
(85)max⁡t∈[0,ω]|u2(t)|<max⁡{|ln⁡α2Me-γ2Lτ2β2L|+2ωα2Mβ2Me-γ2Lτ2β2L,|ln⁡α2Lβ1Le-γ2Mτ2-α1Mc2Me-γ1Lτ1β1Lβ2M|+2ωα2Mβ2Me-γ2Lτ2β2L}≔H2.
It is obvious thatH1 and H2 in (81) and (85) are independent of λ.
Step  2. In order to construct an appropriate open and bounded subset Ω, denote H=H1+H2+H0, where H0 is sufficiently large such that the unique solution (u*,v*)T of the following algebraic equations:
(86)1ω∫0ω‍f1(t)dt=0,1ω∫0ω‍f2(t)dt=0,
satisfies ∥(u*,v*)T∥=|u*|+|v*|<H.
SelectΩ={(u1(t),u2(t))T∈X:∥(u1,u2)T∥<H}, which implies that condition (i) of Lemma 10 holds.
When(u1(t),u2(t))T∈∂Ω∩Ker⁡L=∂Ω∩R2, (u1,u2)T is a constant vector in R2 with |u1|+|u2|=H. Consequently, it can be concluded that
(87)QN[u1u2]=[1ω∫0ω‍f1(t)dt1ω∫0ω‍f2(t)dt]≠[00],
which implies that condition (ii) of Lemma 10 is satisfied.
TakeJ=I:Im⁡Q→Ker⁡L, (u1,u2)T→(u1,u2)T. It follows from straightforward computation that
(88)deg⁡(JQN(u1,u2)T,Ω∩Ker⁡L,(0,0)T)=1,
where (u1*,u2*) is the unique solution of (86). Hence, the condition (iii) of Lemma 10 holds.
Furthermore, it is easy to see that the set{Kp(I-Q)Nx∣x∈Ω-} is equicontinuous and uniformly bounded. By using the Arzela-Ascoli theorem [31], it can be shown that Kp(I-Q)N:Ω-→X is compact and N is L-compact.
Consequently, all conditions (i)–(iii) of Lemma10 hold for Ω. It follows from Lemma 10 that model system (53) has at least one ω-periodic solution (u1*(t),u2*(t))T, and model system (52) has at least one ω-periodic solution (x2*(t),y2*(t))T=(eu1*(t),eu2*(t))T.
Let(x2*(t),y2*(t))T be a positive ω-periodic solution of model system (52); it follows from (19) and (23) that
(89)x1*(t)=e-∫0t‍γ1(s)ds∫t-τ1t‍α1(s)e∫0s‍γ1(m)dmx2*(s)ds,y1*(t)=e-∫0t‍γ2(s)ds∫t-τ2t‍α2(s)e∫0s‍γ2(m)dmy2*(s)ds
are ω-periodic continuous function.
Based on the above analysis, if the following two inequalities hold:(90)α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2,α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ2,
then model system (4) with initial conditions (5) and (6) has at least one positive ω-periodic solution (x1*(t),x2*(t),y1*(t),y2*(t))T.


### 3.3. Global Stability Analysis

Theorem 12.
Ifliminf⁡t→+∞Gk(t)>0, k=1,2, then model system (4) with initial conditions (5) and (6) has a unique positive ω-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T, where
(91)G1=-α1(t)e-∫tt+τ1‍γ1(m)dm-2c1(t)M4*+2β1(t)m2*-8M2*M4*2ρ(t),G2(t)=-qα2(t)e-∫tt+τ2‍γ2(m)dm+2qβ2(t)m4*+2qc1(t)m2*,
and q is a positive constant and mi, Mi, i=1,2,3,4, have been defined in (50).Proof.
Suppose that(x1*(t),x2*(t),y1*(t),y2*(t))T is a positive ω-periodic solution of model system (4) with initial conditions (5) and (6).
Construct a Lyapunov functional as follows:(92)V1(t)=|x2(t)-x2*(t)|+∫t-τ1t‍α1(s)e-∫ss+τ1‍γ1(s)dm×|x2(s)-x2*(s)|ds.
By calculating the upper right derivative ofV1(t) along the positive ω-periodic solution of model system (4), it can be obtained that
(93)D+V1(t)=sgn⁡[x2(t)-x2*(t)]×{α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1)-α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2*(t-τ1)-β1(t)x22(t)+β1(t)x2*2(t)-c1(t)x2(t)y2(t)+c1(t)x2*(t)y2*(t)-ρ(t)x22(t)y22(t)+ρ(t)x2*2(t)y2*2(t)α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1)}+α1(t)e-∫tt+τ1‍γ1(m)dm|x2(t)-x2*(t)|-α1(t-τ1)e-∫t-τ1t‍γ1(m)dm×|x2(t-τ1)-x2*(t-τ1)|≤-β1(t)[x2(t)+x2*(t)]|x2(t)-x2*(t)|+α1(t)e-∫tt+τ1‍γ1(m)dm|x2(t)-x2*(t)|-c1(t)sgn⁡[x2(t)-x2*(t)]×[x2(t)y2(t)-x2*(t)y2*(t)]-ρ(t)sgn⁡[x2(t)-x2*(t)]×[x2(t)y2(t)+x2*(t)y2*(t)]×[x2(t)y2(t)-x2*(t)y2*(t)]≤-{-α1(t)e-∫tt+τ1‍γ1(m)dmβ1(t)[(x2(t)+x2*(t))]+c1(t)[y2(t)-y2*(t)]+ρ(t)(y2(t)-y2*(t))×[(x2(t)-x2*(t))y2(t)+(y2(t)-y2*(t))x2*(t)]-α1(t)e-∫tt+τ1‍γ1(m)dm}|x2(t)-x2*(t)|.
Similarly, construct another Lyapunov functional as follows:(94)V2(t)=|y2(t)-y2*(t)|+∫t-τ2t‍α2(s)e-∫ss+τ2‍γ2(s)dm|y2(s)-y2*(s)|ds.
By calculating the upper right derivative ofV2(t) along the positive ω-periodic solution of model system (4), it can be obtained that
(95)D+V2(t)=sgn⁡[y2(t)-y2*(t)]×{α2(t-τ2)e-∫t-τ2t‍γ2(s)dshhhhh×[y2(t-τ2)-y2*(t-τ2)]hhhhh-β2(t)[y2(t)+y2*(t)][y2(t)-y2*(t)]hhhhh-c2(t)[x2(t)y2(t)-x2*(t)y2*(t)]α2(t-τ2)e-∫t-τ2t‍γ2(s)ds}+α2(t)e-∫tt+τ2‍γ2(m)dm|y2(t)-y2*(t)|-α2(t-τ2)e-∫t-τ2t‍γ2(m)dm×|y2(t-τ2)-y2*(t-τ2)|≤-{β2(t)[y2(t)+y2*(t)]+c1(t)[x2(t)-x2*(t)]}×|y2(t)-y2*(t)|+α2(t)e-∫tt+τ2‍γ2(m)dm|y2(t)-y2*(t)|.
Let(96)V(t)=V1(t)+qV2(t),
where q is a positive constant.
By calculating the upper right derivative ofV(t) along the positive ω-periodic solution of model system (4) based on (93) and (95), it can be obtained as follows:
(97)D+V(t)≤|x2(t)-x2*(t)|×{α1(t)e-∫tt+τ1‍γ1(m)dm-β1(t)[x2(t)+x2*(t)]hhhh+[y2(t)-y2*(t)][+(y2(t)-y2*(t))x2*(t))c1(t)+ρ(t)hhhhhhhhhhhhhhhhhhhh×((x2(t)-x2*(t))y2(t)hhhhhhhhhhhhhihhhhiihhh+(y2(t)-y2*(t))hhhhhhhhhhhhhhhhhhhhh×x2*(t))]α1(t)e-∫tt+τ1‍γ1(m)dm}+|y2(t)-y2*(t)|{qα2(t)e-∫tt+τ2‍γ2(m)dm-qβ2(t)hhhhhhhhhhiihhhhh×[y2(t)+y2*(t)]hhhhhhhhhhhhiihhh-qc1(t)[x2(t)-x2*(t)]qα2(t)e-∫tt+τ2‍γ2(m)dm-qβ2(t)}.
According to Theorem7, there exists a positive value T>0, when t≥T; we get that
(98)m2*-ϵ<x2(t)<M2*+ϵ,m2*-ϵ<x2*(t)<M2*+ϵ,m4*-ϵ<y2(t)<M4*+ϵ,m4*-ϵ<y2*(t)<M4*+ϵ
hold for sufficiently small ϵ>0.
Based on (98), when t>T+max⁡{τ1,τ2}, it is derived that
(99)D+V(t)≤-(G1(t)-ϵ)|x2(t)-x2*(t)|-(G2(t)-ϵ)|y2(t)-y2*(t)|,
where G1(t) and G2(t) have been defined in Theorem 12.
Ifliminf⁡t→+∞Gk(t)>0 for k=1,2, then there exist two constants δ1>0 and δ2>0 such that for t≥T*≔T+2max⁡{τ1,τ2}(100)G1(t)≥δ1>0,G2(t)≥δ2>0.
Consequently, fort≥T*, we have
(101)D+V(t)≤-δ12|x2(t)-x2*(t)|-δ22|y2(t)-y2*(t)|.
By integrating both sides of (101) on the interval [T*,t], it can be obtained that, for t≥T*,
(102)V(t)+δ12∫T*t‍|x2(s)-x2*(s)|ds+δ22∫T*t‍|y2(s)-y2*(s)|ds≤V(T*).
Hence,V(t) is bounded on the interval [T*,+∞) and
(103)∫T*t‍|x2(s)-x2*(s)|ds<+∞,∫T*t‍|y2(s)-y2*(s)|ds<+∞.
According to Barbalat’s Lemma [31], it can be concluded that
(104)lim⁡t→∞|x2(t)-x2*(t)|=0,lim⁡t→∞|y2(t)-y2*(t)|=0.
It follows from (19) and (23) that
(105)|x1(t)-x1*(t)|≤∫t-τ1t‍α1(s)e∫ts‍γ1(m)dm|x2(s)-x2*(s)|ds≤∫t-τ1t‍α1M|x2(s)-x2*(s)|ds,|y1(t)-y1*(t)|≤∫t-τ2t‍α2(s)e∫ts‍γ2(m)dm|y2(s)-y2*(s)|ds≤∫t-τ2t‍α2M|y2(s)-y2*(s)|ds.
Based on (104) and (105), it can be concluded that
(106)lim⁡t→∞|x1(t)-x1*(t)|=0,lim⁡t→∞|y1(t)-y1*(t)|=0.
Therefore, it follows from (104) and (106) that model system (4) with initial conditions (5) and (6) has a unique positive ω-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T.


### 3.4. Numerical Simulation

In this subsection, numerical simulations are carried out to substantiate the analytical findings obtained this paper. In order to facilitate the numerical simulations,ω-periodic continuous functions introduced in model system (4) are selected as follows: α1(t)=2.1+sin(t)/10, γ1(t)=0.2+sin(t)/200, β1(t)=1+sin(t)/300, c1(t)=0.2+sin(t)/300, ρ(t)=0.05+sin(t)/400, α2(t)=4.1+sin(t)/18, γ2(t)=0.1+sin(t)/580, β2(t)=0.3+sin(t)/30, and c2(t)=0.15+sin(t)/270. The maturation delay for nontoxic species and toxin liberating species is given as follows: τ1=0.1 and τ2=0.2, respectively. By using straightforward computation, it can be found that α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2, α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ1; then model system (4) has at least one positive ω-periodic solution based on Theorem 11. Further computations show that G1(t)≥0.4281 and G2(t)≥1.4006. Consequently, it follows from Theorem 12 that model system (4) has a unique positive 2π-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T, whose dynamical responses are plotted in Figure 1. Furthermore, the unique positive 2π-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T is plotted in the x1-x2 plane and y1-y2 plane, which can be found in Figures 2 and 3, respectively. It should be noted that two different initial solutions are included to show the attractivity of different solutions.Figure 1
Dynamical responses of the unique positive2π-periodic globally stable solution of model system (4).Figure 2
The unique positive2π-periodic globally stable solution, which is plotted in the x1-x2 plane. Two different initial solutions are included to show the attractivity of different solutions, which are plotted in red and blue color, respectively.Figure 3
The unique positive2π-periodic globally stable solution, which is plotted in the y1-y2 plane. Two different initial solutions are included to show the attractivity of different solutions, which are plotted in red and blue color, respectively.


## 3.1. Permanence of Solutions

Theorem 1.
Solutions of model system (4) with initial conditions (5) and (6) are positive for all t>0.Proof.
Firstly, we show thatx2(t)>0 for all t>0. Otherwise, if it is false, since x2(t)>0 for all t∈[-τ1,0], then it can be derived that there exists a t1>0 such that x2(t1)=0.
Definet0=inf⁡{t>0∣x2(t)=0}. According to the definition of t0, it can be obtained that
(8)x˙2(t0)≤0.
It follows from the second equation of model system (4) that
(9)x˙2(t0)={α1(t0-τ1)e-∫t0-τ1t0‍γ1(s)dsϕ2(t0-τ1)>0,0≤t0≤τ1,α1(t0-τ1)e-∫t0-τ1t0‍γ1(s)dsx2(t0-τ1)>0,t>τ1,
and it is easy to show that x˙2(t0)>0, which is a contradiction to (8). Hence, x2(t)>0 for all t>0.
By a direct computation, it follows from the first equation of model system (4) that
(10)x1(t)=∫t-τ1t‍α1(s)e∫ts‍γ1(m)dsx2(s)ds,
since x2(t)>0 for all t>0; it is easy to show that x1(t)>0 for all t>0 based on (10).
By utilizing the similar proof, it can be obtained thaty1(t)>0 and y2(t)>0 for all t>0. Consequently, solutions of model system (4) with initial conditions (5) and (6) are positive for all t>0.Theorem 2.
Solutions of model system (4) with initial conditions (5) and (6) are ultimately bounded.Proof.
Letw(t)=x1(t)+x2(t)+y1(t)+y2(t), where (x1(t),x2(t),y1(t),y2(t)) is an arbitrary positive solution of model system (4) with the initial conditions (5) and (6).
Calculating the derivative ofw(t) along the solution of model system (4) gives that
(11)w˙(t)≤(α1(t)+γ1(t))x2(t)-γ1(t)(x1(t)+x2(t))-β1(t)x22(t)+(α2(t)+γ2(t))y2(t)-γ2(t)(y1(t)+y2(t))-β2(t)y22(t)≤(α1M+γ1M)x2(t)-β1Lx22(t)-r1L(x1(t)+x2(t))+(α2M+γ2M)y2(t)-β2Ly22(t)-r2L(y1(t)+y2(t))≤-γLw(t)+(α1M+γ1M)x2(t)-β1Lx22(t)+(α2M+γ2M)y2(t)-β2Ly22(t)≤-γLw(t)+(α1M+γ1M)24β1L+(α2M+γ2M)24β2L,
where γL=min⁡{γ1L,γ2L}.
By using the standard comparison principle [20], it follows from (11) that
(12)w(t)≤β2L(α1M+γ1M)2+β1L(α2M+γ2M)24γLβ1Lβ2L,
which implies that any solution of model system (4) with initial conditions (5) and (6) is ultimately bounded.Lemma 3 (see [29]).
Consider the following differential equation:(13)x˙(t)=ax(t-τ)-bx(t)-cx2(t),
where a,b,c,τ>0 and x(t)>0 for -τ≤t≤0; we have that(i)
ifa>b, then lim⁡t→+∞x(t)=(a-b)/c;(ii)
ifa<b, then lim⁡t→+∞x(t)=0.Lemma 4 (see [29]).
Consider the following differential equation:(14)x˙(t)=dx(t-σ)-ex2(t),
where d,e,σ>0 and x(t)>0 for -σ≤t≤0; we have
(15)lim⁡t→+∞x(t)=de.Definition 5 (see [30]).
Model system(16)X˙(t)=f(t,Xt(θ)),
where t≥0, θ∈[-τ,0], X∈Rn. Model system (16) is said to be permanent if, for any solution X(t,ϕ), there exists a constant m>0 and T=T(ϕ) such that X(t)>m for all t>T.Definition 6 (see [30]).
The domainD∈Cn is said to be an ultimately bounded domain, if D is a closed, bounded subset of Cn, and there exists a constant T=T(ϕ) such that Xt(θ)∈D for all t>T.Theorem 7.
Ifα1Lβ2L>c1Mα2M and α2Lβ1L>c2Mα1M, then model system (4) is permanent with initial conditions (5) and (6).Proof.
According to the second equation of model system (4) and Theorem 1, we get that
(17)x˙2(t)≤α1Me-γ1Lτ1x2(t-τ1)-β1Lx22(t).
By virtue of Lemma4 and (17), there exists a positive time T1 such that, for sufficiently small ϵ>0 and t≥T1, it yields
(18)x2(t)≤α1Me-γ1Lτ1β1L+ϵ:=M2(1).
By rearranging (10), it can be obtained that
(19)x1(t)=e-∫0t‍γ1(s)ds∫t-τ1t‍α1(s)e∫0s‍γ1(m)dmx2(s)ds.
For anyt≥T1, it follows from (18) and (19) that
(20)x1(t)≤a1MM2(1)(1-e-γ1Mτ1)γ1L≔M1(1).
Based on the fourth equation of model system (4) and Theorem 1, it can be obtained that
(21)y˙2(t)≤α2Me-γ2Lτ2y2(t-τ2)-β2Ly22(t)
holds for t≥T1.
By virtue of Lemma4 and (21), there exists T2>T1 such that, for sufficiently small ϵ>0 and t≥T2, it yields
(22)y2(t)≤α2Me-γ2Lτ2β2L+ϵ≔M4(1).
By direct computing, it follows from the third equation of model system (4) that
(23)y1(t)=e-∫0t‍γ2(s)ds∫t-τ2t‍α2(s)e∫0s‍γ2(m)dmy2(s)ds.
For anyt≥T2, it follows from (22) and (23) that
(24)y1(t)≤α2MM4(1)(1-e-γ2Mτ2)γ2L≔M3(1).
Furthermore, it follows from the second equation of model system (4) that
(25)x˙2(t)≥α1Le-γ1Mτ1x2(t-τ1)-(β1M+ρ1M(M4(1))2)x22(t)-c1MM4(1)x2(t).
Based on Lemma3 and (25), there exists T3>T2 and for any t≥T3 and sufficiently small ϵ>0,
(26)x2(t)≥α1Le-γ1Mτ1-c1MM4(1)β1M+ρ1M(M4(1))2-ϵ≔m2(1)
holds provided that
(27)α1Le-γ1Mτ1>c1MM4(1).
According to (19), for any t≥T3, we get that
(28)x1(t)≥α1L(1-e-γ1Lτ1)m2(1)γ1M≔m1(1).
For anyt≥T3, it follows from the fourth equation of model system (4) that
(29)y˙2(t)≥α2Le-γ2Mτ2y2(t-τ2)-β2My22(t)-c2MM2(1)y2(t).
Based on Lemma3 and (29), there exists T4>T3 and for any t≥T4 and sufficiently small ϵ>0,
(30)y2(t)≥α2Le-γ2Mτ2-c2MM2(1)β2M-ϵ≔m4(1)
holds provided that
(31)α2Le-γ2Mτ2>c2MM2(1).
According to (23) and (30), it can be obtained that
(32)y1(t)≥α2Lm4(1)(1-e-γ2Lτ2)γ2M≔m3(1).
According to the second equation of model system (4), we get that
(33)x˙2(t)≤α1Me-γ1Lτ1x2(t-τ1)-(β1L+ρL(m4(1))2)x22(t)-c1Lm4(1)x2(t).
By virtue of Lemma3 and (33), there exists T5>T4 such that, for sufficiently small ϵ>0 and t≥T5, it yields
(34)x2(t)≤α1Me-γ1Lτ1-c1Lm4(1)β1L+ρL(m4(1))2+ϵ≔M2(2),
which holds provided that
(35)α1Me-γ1Lτ1>c1Lm4(1).
For anyt≥T5, it follows from (19) and (34) that
(36)x1(t)≤α1MM2(2)(1-e-γ1Mτ1)γ1L≔M1(2).
Based on the fourth equation of model system (4), it can be obtained that
(37)y˙2(t)≤α2Me-γ2Lτ2y2(t-τ2)-β2Ly22(t)-c2Lm2(1)y2(t)
holds for t≥T5.
By virtue of Lemma3 and (37), there exists T6>T5 such that, for sufficiently small ϵ>0 and t≥T6, it yields
(38)y2(t)≤α2Me-γ2Lτ2-c2Lm2(1)β2L+ϵ≔M4(2),
which holds provided that
(39)α2Me-γ2Lτ2>c2Lm2(1).
For anyt≥T6, it follows from (23) and (38) that
(40)y1(t)≤α2MM4(2)(1-e-γ2Mτ2)γ2L≔M3(2).
Furthermore, it follows from the second equation of model system (4) that
(41)x˙2(t)≥α1Le-γ1Mτ2x2(t-τ1)-(β1M+ρM(M4(2))2)x22(t)-c1MM4(2)x2(t).
Based on Lemma3 and (41), there exists T7>T6 and, for any t≥T7 and sufficiently small ϵ>0,
(42)x2(t)≥α1Le-γ1Mτ1-c1MM4(2)β1M+ρM(M4(2))2-ϵ≔m2(2)
holds provided that
(43)α1Le-γ1Mτ1>c1MM4(2).
According to (19) and (42), for any t≥T7, we get that
(44)x1(t)≥α1L(1-e-γ1Lτ1)m2(2)γ1M≔m1(2).
For anyt≥T7, it follows from the fourth equation of model system (4) that
(45)y˙2(t)≥α2Le-γ2Mτ2y2(t-τ2)-β2My22(t)-c2MM2(2)y2(t).
Based on Lemma3 and (45), there exists T8>T7 and, for any t≥T8 and sufficiently small ϵ>0,
(46)y2(t)≥α2Le-γ2Mτ2-c2MM2(2)β2M-ϵ≔m4(2)
holds provided that
(47)α2Le-γ2Mτ2>c2MM2(2).
According to (23) and (46), it can be obtained that
(48)y1(t)≥α2Lm4(2)(1-e-γ2Lτ2)γ2M≔m3(2).
By using simple computation, it is easy to show that six inequalities (27), (31), (35), (39), (43), and (47) hold if the following two inequalities α1Lβ2L>c1Mα2M and α2Lβ1L>c2Mα1M hold.
Furthermore, eight sequences will be obtained by repeating the discussion in this manner, which are given as follows:(49)M1(n+1)=α1MM2(n+1)(1-e-γ1Mτ1)γ1L,M2(n+1)=α1Me-γ1Lτ1-c1Lm4(n)β1L+ρL(m4(n))2+ϵ,M3(n+1)=α2MM4(n+1)(1-e-γ2Mτ2)γ2L,M4(n+1)=α2Me-γ2Lτ2-c2Lm2(n)β2L+ϵ,m1(n+1)=α1L(1-e-γ1Lτ1)m2(n+1)γ1M,m2(n+1)=α1Le-γ1Mτ1-c1MM4(n+1)β1M+ρM(M4(n+1))2-ϵ,m3(n+1)=α2Lm4(n+1)(1-e-γ2Lτ2)γ2M,m4(n+1)=α2Le-γ2Mτ2-c2MM2(n+1)β2M-ϵ.
It is easy to show thatMi(n)>0 and the sequences {Mi(n)} (i=1,2,3,4) are decreasing as n increases, which implies that lim⁡n→∞Mi(n)=Mi* exists; furthermore, it is easy to show that mi(n)<Mi(n) and the sequences {mi(n)} (i=1,2,3,4) are increasing as n increases, which implies that lim⁡n→∞mi(n)=mi* exists. Consequently, it follows from (49) that
(50)M1*=α1MM2*(1-e-γ1Mτ1)γ1L,M2*=α1Me-γ1Lτ1-c1Lm4*β1L+ρLm4*2,M3*=α2MM4*(1-e-γ2Mτ2)γ2L,M4*=α2Me-γ2Lτ2-c2Lm2*β2L,m1*=α1L(1-e-γ1Lτ1)m2*γ1M,m2*=α1Le-γ1Mτ1-c1MM4*β1M+ρMM4*2,m3*=α2Lm4*(1-e-γ2Lτ2)γ2M,m4*=α2Le-γ2Mτ2-c2MM2*β2M.
Based on Definition5 and (50), it can be concluded that model system (4) is persistent if α1Lβ2L>c1Mα2M and α2Lβ1L>c2Mα1M hold.


## 3.2. Existence of Positive Periodic Solutions

Definition 8 (see [31]).
LetL:Dom⁡L⊂X→Y be a linear mapping and let N:X→Y be a continuous mapping, where X and Y are real Banach spaces. If dimKerL=codimImL<+∞ and Im⁡L is closed in Y, then L is called a Fredholm mapping of index zero.IfL is Fredholm mapping of index zero and there exist continuous projectors P:X→X and Q:Y→Y such that
(51)Im⁡P=Ker⁡L,Im⁡L=Ker⁡Q=Im⁡(I-Q),
then restriction Lp of L to Dom⁡L∩Ker⁡P:(I-P)X→Im⁡L is invertible.Definition 9 (see [31]).
Denote the inverse ofLp by Kp. Supposing that Ω is an open bounded subset of X, if QN(Ω-) is bounded and Kp(I-Q)N:Ω→X is compact, then the mapping N is called L-compact on Ω-. Since Im⁡Q is isomorphic to Ker⁡L, there exists an isomorphism J:Im⁡Q→Ker⁡L.Lemma 10 (see [31]).
LetΩ⊂X be an open bounded set, let L be a Fredholm mapping of index zero, and let N be L-compact on Ω-. If the following three conditions hold: (i)
Lx≠λNx for any λ∈(0,1) and x∈∂Ω∩Dom⁡L,(ii)
QNx≠0 for any x∈∂Ω∩Ker⁡L,(iii)
deg⁡{JQN,Ω∩Ker⁡L,0}≠0,then Lx=Nx has at least one solution in Ω-∩Dom⁡L.Theorem 11.
Ifα1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2, α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ2, then model system (4) with initial conditions (5) and (6) has at least one positive ω-periodic solution.Proof.
Consider the subsystem of model system (4):
(52)x˙2(t)=α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1)-β1(t)x22(t)-c1(t)x2(t)y2(t)-ρ(t)x22(t)y22(t),y˙2(t)=α2(t-τ2)e-∫t-τ2t‍γ2(s)dsy2(t-τ2)-β2(t)y22(t)-c2(t)x2(t)y2(t).
Letu1(t)=ln⁡[x2(t)], u2(t)=ln⁡[y2(t)].
By substitutingu1(t) and u2(t) into (52), it can be obtained that
(53)u˙1(t)=α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)-u1(t)-β1(t)eu1(t)-c1(t)eu2(t)-ρ(t)eu1(t)+2u2(t),u˙2(t)=α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)-u2(t)-β2(t)eu2(t)-c2(t)eu1(t).
It should be noted that if model system (53) has one ω-periodic solution (u1*(t),u2*(t))T, then (x2*(t),y2*(t))T=(eu1*(t),eu2*(t))T is a positive ω-periodic solution of model system (52).
In order to utilize Lemma10 in a straightforward manner, we define
(54)X=Y={(u1(t),u2(t))Ti∈C(R,R2):ui(t+ω)=ui(t),i=1,2(u1(t),u2(t))T},∥(u1(t),u2(t))T∥=max⁡t∈[0,ω]|u1(t)|+max⁡t∈[0,ω]|u2(t)|,
where |·| denotes the Euclidean norm; it is easy to show that both X and Y are Banach spaces with the norm ∥·∥; then define
(55)Dom⁡L∩X⟶X,L(u1(t),u2(t))T=(du1(t)dt,du2(t)dt)T,
where Dom⁡L={(u1(t),u2(t))T∈C(R,R2)}, N[u1u2]=[f1(t)f2(t)], and
(56)f1(t)=α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)-u1(t)-β1(t)eu1(t)-c1(t)eu2(t)-ρ(t)eu1(t)+2u2(t),f2(t)=α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)-u2(t)-β2(t)eu2(t)-c2(t)eu1(t).
Furthermore, we define(57)P[u1u2]=Q[u1u2]=[1ω∫0ω‍u1(t)dt1ω∫0ω‍u2(t)dt],[u1u2]∈X=Y.
According to the above definitions, it is not difficult to verify thatKer⁡L={x∣x∈X,x=h,h∈R2}, Im⁡L={y∈Y∣∫0ω‍y(t)dt=0} are closed in Y, dimKerL=codimImL=2, and both P and Q are continuous projectors such that Im⁡P=Ker⁡L and Ker⁡Q=Im⁡L=Im⁡(I-Q).
Based on the above analysis, it can be obtained thatL is a Fredholm mapping of index zero.
Furthermore, the inverseKp:Im⁡L→Dom⁡L∩Ker⁡P of Lp exists and takes the following form:
(58)Kp(y)=∫0t‍y(s)ds-1ω∫0ω‍∫0t‍y(s)dsdt.
Hence,QN:X→Y and Kp(I-Q)N:X→X can be defined as follows, respectively,
(59)QNx=[1ω∫0ω‍f1(t)dt1ω∫0ω‍f2(t)dt],Kp(I-Q)Nx=∫0t‍Nx(s)ds-1ω∫0ω‍∫0t‍Nx(s)dsdt-(tω-12)∫0ω‍Nx(s)ds.
It is easy to show thatQN and Kp(I-Q)N are continuous. In order to facilitate the proof based on Lemma 10, we also need to find an appropriate open and bounded subset Ω, which can be found by the following two steps.
Step  1. According to the operator equation Lx=λNx for λ∈(0,1), the upper and lower bound of u1(t) and u2(t) will be estimated as follows:
(60)du1(t)dt=λf1(t),du2(t)dt=λf2(t),
where f1(t) and f2(t) have been defined in (56).
Suppose that(u1(t),u2(t))T∈X is a solution of model system (60) for some λ∈(0,1). By integrating (60) over the interval [0,ω], it can be obtained that
(61)∫0ω‍α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)-u1(t)dt=∫0ω‍β1(t)eu1(t)+c1(t)eu2(t)+ρ(t)eu1(t)+2u2(t)dt,(62)∫0ω‍α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)-u2(t)dt=∫0ω‍β2(t)eu2(t)+c2(t)eu1(t)dt.
Based on definition(u1(t),u2(t))T∈X, there exist ξi,ηi∈[0,ω] such that
(63)ui(ξi)=min⁡t∈[0,ω]ui(t),ui(ηi)=max⁡t∈[0,ω]ui(t),i=1,2.
Multiplying the first equation of (60) by eu1(t) and integrating it over [0,ω] give that
(64)∫0ω‍α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)dt=∫0ω‍β1(t)e2u1(t)+c1(t)eu1(t)+u2(t)hhhhhh+ρ(t)e2(u1(t)+u2(t))dt.
It follows from (64) that
(65)β1L∫0ω‍e2u1(t)dt+c1L∫0ω‍eu1(t)+u2(t)dt+ρL∫0ω‍e2(u1(t)+u2(t))dt≤α1Me-γ1Lτ1∫0ω‍eu1(t)dt.
On the other hand, by using the inequality,(66)(∫0ω‍eu1(t)dt)2≤ω∫0ω‍e2u1(t)dt.
Based on (65) and (66), it can be obtained that
(67)β1L(∫0ω‍eu1(t)dt)2≤ωβ1L∫0ω‍e2u1(t)dt≤ωα1Me-γ1Lτ1∫0ω‍eu1(t)dt,
which derives that
(68)β1L∫0ω‍eu1(t)dt≤ωα1Me-γ1Lτ1,u1(ξ1)≤ln⁡α1Me-γ1Lτ1β1L.
It follows from (60) and (68) that
(69)∫0ω‍|u1′(t)|dt<2∫0ω‍β1(t)eu1(t)≤2ωα1Mβ1Me-γ1Lτ1β1L.
According to (68) and (69), it can be obtained that
(70)u1(t)≤u1(ξ1)+∫0ω|u1′(t)|dt≤ln⁡α1Me-γ1Lτ1β1L+2ωα1Mβ1Me-γ1Lτ1β1L.
Multiplying the second equation of (60) by eu2(t) and integrating it over [0,ω] give that
(71)∫0ω‍α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)dt=∫0ω‍β2(t)e2u2(t)+c2(t)eu1(t)+u2(t)dt.
It follows from (71) that
(72)β2L∫0ω‍e2u2(t)dt+c2L∫0ω‍eu1(t)+u2(t)dt≤α2Me-γ2Lτ2∫0ω‍eu2(t)dt.
Based on (66) and (72), it can be obtained that
(73)β2L(∫0ω‍eu1(t)dt)2≤ωα2Me-γ2Lτ2∫0ω‍eu2(t)dt,
which derives that
(74)∫0ω‍eu2(t)dt≤ωα2Me-γ2Lτ2β2L,u2(ξ2)≤ln⁡α2Me-γ2Lτ2β2L.
It follows from (60) and (74) that
(75)∫0ω‍|u2′(t)|dt<2∫0ω‍β2(t)eu2(t)dt≤2ωα2Mβ2Me-γ2Lτ2β2L.
According to (74) and (75), it can be obtained that
(76)u2(t)≤u2(ξ2)+∫0ω‍|u2′(t)|dt≤ln⁡α2Me-γ2Lτ2β2L+2ωα2Mβ2Me-γ2Lτ2β2L.
It should be noted that(77)∫0ω‍α1(t)e-∫tt+τ1‍γ1(s)dseu1(t)dt=∫0ω‍α1(t-τ1)e-∫t-τ1t‍γ1(s)dseu1(t-τ1)dt.
Based on (64), it can be obtained that
(78)(α1Le-γ2Mτ1-c1Mα2Me-γ2Lτ2β2L)∫0ω‍eu1(t)dt≤(β1M+2ρMα2Me-γ2Lτ2β2L)∫0ω‍e2u1(t)dt,
which derives that
(79)∫0ω‍eu1(t)dt≥ω(α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2)β1Mβ2L+2α2MρMe-γ2Lτ2,u1(η1)≥ln⁡α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2β1Mβ2L+2α2MρMe-γ2Lτ2
hold provided that α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2.
According to (69) and (79), it can be obtained that
(80)u1(t)≥u1(η1)-∫0ω‍|u1′(t)|dt≥ln⁡α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2β1Mβ2L+2α2MρMe-γ2Lτ2-2ωα1Mβ1Me-γ2Lτ2β1L.
By virtue of (70) and (80), if α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2, then
(81)max⁡t∈[0,ω]|u1(t)|<max⁡{|ln⁡α1Me-γ1Lτ1β1L|+2ωα1Mβ1Me-γ1Lτ1β1L,|ln⁡α1Lβ2Le-γ2Mτ1-α2Mc1Me-γ2Lτ2β1Mβ2L+2α2MρMe-γ2Lτ2|+2ωα1Mβ1Me-γ1Lτ1β1L}≔H1.
Similarly, it is easy to show that(82)∫0ω‍α2(t)e-∫tt+τ2‍γ2(s)dseu2(t)dt=∫0ω‍α2(t-τ2)e-∫t-τ2t‍γ2(s)dseu2(t-τ2)dt,
which derives that
(83)u2(η2)≥ln⁡α2Lβ1Le-γ2Mτ2-α1Mc2Me-γ1Lτ1β1Lβ2M
holds provided that α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ1.
According to (75) and (83), it is derived that
(84)u2(t)≥u2(η2)-∫0ω‍|u2′(t)|dt≥ln⁡α2Lβ1Le-γ2Mτ2-α1Mc2Me-γ1Lτ1β1Lβ2M-2ωα2Mβ2Me-γ2Lτ2β2L.
By virtue of (76) and (84), if α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ2, then
(85)max⁡t∈[0,ω]|u2(t)|<max⁡{|ln⁡α2Me-γ2Lτ2β2L|+2ωα2Mβ2Me-γ2Lτ2β2L,|ln⁡α2Lβ1Le-γ2Mτ2-α1Mc2Me-γ1Lτ1β1Lβ2M|+2ωα2Mβ2Me-γ2Lτ2β2L}≔H2.
It is obvious thatH1 and H2 in (81) and (85) are independent of λ.
Step  2. In order to construct an appropriate open and bounded subset Ω, denote H=H1+H2+H0, where H0 is sufficiently large such that the unique solution (u*,v*)T of the following algebraic equations:
(86)1ω∫0ω‍f1(t)dt=0,1ω∫0ω‍f2(t)dt=0,
satisfies ∥(u*,v*)T∥=|u*|+|v*|<H.
SelectΩ={(u1(t),u2(t))T∈X:∥(u1,u2)T∥<H}, which implies that condition (i) of Lemma 10 holds.
When(u1(t),u2(t))T∈∂Ω∩Ker⁡L=∂Ω∩R2, (u1,u2)T is a constant vector in R2 with |u1|+|u2|=H. Consequently, it can be concluded that
(87)QN[u1u2]=[1ω∫0ω‍f1(t)dt1ω∫0ω‍f2(t)dt]≠[00],
which implies that condition (ii) of Lemma 10 is satisfied.
TakeJ=I:Im⁡Q→Ker⁡L, (u1,u2)T→(u1,u2)T. It follows from straightforward computation that
(88)deg⁡(JQN(u1,u2)T,Ω∩Ker⁡L,(0,0)T)=1,
where (u1*,u2*) is the unique solution of (86). Hence, the condition (iii) of Lemma 10 holds.
Furthermore, it is easy to see that the set{Kp(I-Q)Nx∣x∈Ω-} is equicontinuous and uniformly bounded. By using the Arzela-Ascoli theorem [31], it can be shown that Kp(I-Q)N:Ω-→X is compact and N is L-compact.
Consequently, all conditions (i)–(iii) of Lemma10 hold for Ω. It follows from Lemma 10 that model system (53) has at least one ω-periodic solution (u1*(t),u2*(t))T, and model system (52) has at least one ω-periodic solution (x2*(t),y2*(t))T=(eu1*(t),eu2*(t))T.
Let(x2*(t),y2*(t))T be a positive ω-periodic solution of model system (52); it follows from (19) and (23) that
(89)x1*(t)=e-∫0t‍γ1(s)ds∫t-τ1t‍α1(s)e∫0s‍γ1(m)dmx2*(s)ds,y1*(t)=e-∫0t‍γ2(s)ds∫t-τ2t‍α2(s)e∫0s‍γ2(m)dmy2*(s)ds
are ω-periodic continuous function.
Based on the above analysis, if the following two inequalities hold:(90)α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2,α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ2,
then model system (4) with initial conditions (5) and (6) has at least one positive ω-periodic solution (x1*(t),x2*(t),y1*(t),y2*(t))T.


## 3.3. Global Stability Analysis

Theorem 12.
Ifliminf⁡t→+∞Gk(t)>0, k=1,2, then model system (4) with initial conditions (5) and (6) has a unique positive ω-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T, where
(91)G1=-α1(t)e-∫tt+τ1‍γ1(m)dm-2c1(t)M4*+2β1(t)m2*-8M2*M4*2ρ(t),G2(t)=-qα2(t)e-∫tt+τ2‍γ2(m)dm+2qβ2(t)m4*+2qc1(t)m2*,
and q is a positive constant and mi, Mi, i=1,2,3,4, have been defined in (50).Proof.
Suppose that(x1*(t),x2*(t),y1*(t),y2*(t))T is a positive ω-periodic solution of model system (4) with initial conditions (5) and (6).
Construct a Lyapunov functional as follows:(92)V1(t)=|x2(t)-x2*(t)|+∫t-τ1t‍α1(s)e-∫ss+τ1‍γ1(s)dm×|x2(s)-x2*(s)|ds.
By calculating the upper right derivative ofV1(t) along the positive ω-periodic solution of model system (4), it can be obtained that
(93)D+V1(t)=sgn⁡[x2(t)-x2*(t)]×{α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1)-α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2*(t-τ1)-β1(t)x22(t)+β1(t)x2*2(t)-c1(t)x2(t)y2(t)+c1(t)x2*(t)y2*(t)-ρ(t)x22(t)y22(t)+ρ(t)x2*2(t)y2*2(t)α1(t-τ1)e-∫t-τ1t‍γ1(s)dsx2(t-τ1)}+α1(t)e-∫tt+τ1‍γ1(m)dm|x2(t)-x2*(t)|-α1(t-τ1)e-∫t-τ1t‍γ1(m)dm×|x2(t-τ1)-x2*(t-τ1)|≤-β1(t)[x2(t)+x2*(t)]|x2(t)-x2*(t)|+α1(t)e-∫tt+τ1‍γ1(m)dm|x2(t)-x2*(t)|-c1(t)sgn⁡[x2(t)-x2*(t)]×[x2(t)y2(t)-x2*(t)y2*(t)]-ρ(t)sgn⁡[x2(t)-x2*(t)]×[x2(t)y2(t)+x2*(t)y2*(t)]×[x2(t)y2(t)-x2*(t)y2*(t)]≤-{-α1(t)e-∫tt+τ1‍γ1(m)dmβ1(t)[(x2(t)+x2*(t))]+c1(t)[y2(t)-y2*(t)]+ρ(t)(y2(t)-y2*(t))×[(x2(t)-x2*(t))y2(t)+(y2(t)-y2*(t))x2*(t)]-α1(t)e-∫tt+τ1‍γ1(m)dm}|x2(t)-x2*(t)|.
Similarly, construct another Lyapunov functional as follows:(94)V2(t)=|y2(t)-y2*(t)|+∫t-τ2t‍α2(s)e-∫ss+τ2‍γ2(s)dm|y2(s)-y2*(s)|ds.
By calculating the upper right derivative ofV2(t) along the positive ω-periodic solution of model system (4), it can be obtained that
(95)D+V2(t)=sgn⁡[y2(t)-y2*(t)]×{α2(t-τ2)e-∫t-τ2t‍γ2(s)dshhhhh×[y2(t-τ2)-y2*(t-τ2)]hhhhh-β2(t)[y2(t)+y2*(t)][y2(t)-y2*(t)]hhhhh-c2(t)[x2(t)y2(t)-x2*(t)y2*(t)]α2(t-τ2)e-∫t-τ2t‍γ2(s)ds}+α2(t)e-∫tt+τ2‍γ2(m)dm|y2(t)-y2*(t)|-α2(t-τ2)e-∫t-τ2t‍γ2(m)dm×|y2(t-τ2)-y2*(t-τ2)|≤-{β2(t)[y2(t)+y2*(t)]+c1(t)[x2(t)-x2*(t)]}×|y2(t)-y2*(t)|+α2(t)e-∫tt+τ2‍γ2(m)dm|y2(t)-y2*(t)|.
Let(96)V(t)=V1(t)+qV2(t),
where q is a positive constant.
By calculating the upper right derivative ofV(t) along the positive ω-periodic solution of model system (4) based on (93) and (95), it can be obtained as follows:
(97)D+V(t)≤|x2(t)-x2*(t)|×{α1(t)e-∫tt+τ1‍γ1(m)dm-β1(t)[x2(t)+x2*(t)]hhhh+[y2(t)-y2*(t)][+(y2(t)-y2*(t))x2*(t))c1(t)+ρ(t)hhhhhhhhhhhhhhhhhhhh×((x2(t)-x2*(t))y2(t)hhhhhhhhhhhhhihhhhiihhh+(y2(t)-y2*(t))hhhhhhhhhhhhhhhhhhhhh×x2*(t))]α1(t)e-∫tt+τ1‍γ1(m)dm}+|y2(t)-y2*(t)|{qα2(t)e-∫tt+τ2‍γ2(m)dm-qβ2(t)hhhhhhhhhhiihhhhh×[y2(t)+y2*(t)]hhhhhhhhhhhhiihhh-qc1(t)[x2(t)-x2*(t)]qα2(t)e-∫tt+τ2‍γ2(m)dm-qβ2(t)}.
According to Theorem7, there exists a positive value T>0, when t≥T; we get that
(98)m2*-ϵ<x2(t)<M2*+ϵ,m2*-ϵ<x2*(t)<M2*+ϵ,m4*-ϵ<y2(t)<M4*+ϵ,m4*-ϵ<y2*(t)<M4*+ϵ
hold for sufficiently small ϵ>0.
Based on (98), when t>T+max⁡{τ1,τ2}, it is derived that
(99)D+V(t)≤-(G1(t)-ϵ)|x2(t)-x2*(t)|-(G2(t)-ϵ)|y2(t)-y2*(t)|,
where G1(t) and G2(t) have been defined in Theorem 12.
Ifliminf⁡t→+∞Gk(t)>0 for k=1,2, then there exist two constants δ1>0 and δ2>0 such that for t≥T*≔T+2max⁡{τ1,τ2}(100)G1(t)≥δ1>0,G2(t)≥δ2>0.
Consequently, fort≥T*, we have
(101)D+V(t)≤-δ12|x2(t)-x2*(t)|-δ22|y2(t)-y2*(t)|.
By integrating both sides of (101) on the interval [T*,t], it can be obtained that, for t≥T*,
(102)V(t)+δ12∫T*t‍|x2(s)-x2*(s)|ds+δ22∫T*t‍|y2(s)-y2*(s)|ds≤V(T*).
Hence,V(t) is bounded on the interval [T*,+∞) and
(103)∫T*t‍|x2(s)-x2*(s)|ds<+∞,∫T*t‍|y2(s)-y2*(s)|ds<+∞.
According to Barbalat’s Lemma [31], it can be concluded that
(104)lim⁡t→∞|x2(t)-x2*(t)|=0,lim⁡t→∞|y2(t)-y2*(t)|=0.
It follows from (19) and (23) that
(105)|x1(t)-x1*(t)|≤∫t-τ1t‍α1(s)e∫ts‍γ1(m)dm|x2(s)-x2*(s)|ds≤∫t-τ1t‍α1M|x2(s)-x2*(s)|ds,|y1(t)-y1*(t)|≤∫t-τ2t‍α2(s)e∫ts‍γ2(m)dm|y2(s)-y2*(s)|ds≤∫t-τ2t‍α2M|y2(s)-y2*(s)|ds.
Based on (104) and (105), it can be concluded that
(106)lim⁡t→∞|x1(t)-x1*(t)|=0,lim⁡t→∞|y1(t)-y1*(t)|=0.
Therefore, it follows from (104) and (106) that model system (4) with initial conditions (5) and (6) has a unique positive ω-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T.


## 3.4. Numerical Simulation

In this subsection, numerical simulations are carried out to substantiate the analytical findings obtained this paper. In order to facilitate the numerical simulations,ω-periodic continuous functions introduced in model system (4) are selected as follows: α1(t)=2.1+sin(t)/10, γ1(t)=0.2+sin(t)/200, β1(t)=1+sin(t)/300, c1(t)=0.2+sin(t)/300, ρ(t)=0.05+sin(t)/400, α2(t)=4.1+sin(t)/18, γ2(t)=0.1+sin(t)/580, β2(t)=0.3+sin(t)/30, and c2(t)=0.15+sin(t)/270. The maturation delay for nontoxic species and toxin liberating species is given as follows: τ1=0.1 and τ2=0.2, respectively. By using straightforward computation, it can be found that α1Lβ2Le-γ2Mτ1>α2Mc1Me-γ2Lτ2, α2Lβ1Le-γ2Mτ2>α1Mc2Me-γ1Lτ1; then model system (4) has at least one positive ω-periodic solution based on Theorem 11. Further computations show that G1(t)≥0.4281 and G2(t)≥1.4006. Consequently, it follows from Theorem 12 that model system (4) has a unique positive 2π-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T, whose dynamical responses are plotted in Figure 1. Furthermore, the unique positive 2π-periodic globally stable solution (x1*(t),x2*(t),y1*(t),y2*(t))T is plotted in the x1-x2 plane and y1-y2 plane, which can be found in Figures 2 and 3, respectively. It should be noted that two different initial solutions are included to show the attractivity of different solutions.Figure 1
Dynamical responses of the unique positive2π-periodic globally stable solution of model system (4).Figure 2
The unique positive2π-periodic globally stable solution, which is plotted in the x1-x2 plane. Two different initial solutions are included to show the attractivity of different solutions, which are plotted in red and blue color, respectively.Figure 3
The unique positive2π-periodic globally stable solution, which is plotted in the y1-y2 plane. Two different initial solutions are included to show the attractivity of different solutions, which are plotted in red and blue color, respectively.


## 4. Conclusion

In this paper, a nonautonomous dynamical model is proposed to investigate population dynamics of competitive system with toxin liberating species and nontoxic species, where stage structure and maturation delay for two species are considered. It is well known that the effect of toxin ecological systems is an important issue from mathematical and experimental points of view. Generally speaking, it takes some time for a species to reach maturity to produce the toxicant, the toxin liberating mature individual produces a substance toxic to the nontoxic mature individuals only, and the inhibiting effect is zero in absence of either species. Furthermore, the species compete each other for the limited life resource within closed environment, but this competition only happens among the mature individuals and does not involve the immature individuals. Consequently, it is necessary to investigate the dynamic effect of stage structure and toxic substances on population dynamics of two-species competitive system.Qualitative analysis of the proposed model system is discussed in the third section of this paper. It follows from Theorems1 and 2 that solutions of model system (4) with initial conditions are positive and ultimately bounded. By utilizing some comparison arguments, an iterative technique is proposed to discuss permanence of the species; model system (4) is persistent, which can be found in Theorem 7. Furthermore, existence of positive periodic solutions is considered in Theorem 11 based on continuation theorem of coincidence degree theory, which shows that model system (4) has at least one positive ω-solution. By constructing an appropriate Lyapunov functional, sufficient conditions for global stability of the unique positive periodic solution are analyzed; that is, liminf⁡t→+∞Gk(t)>0,k=1,2, which can be found in Theorem 12. Finally, numerical simulations are provided to show dynamical responses of the unique positive 2π-periodic globally stable solution, which are plotted in Figure 1. Furthermore, the unique positive 2π-periodic globally stable solution is plotted in the x1-x2 plane and y1-y2 plane, which can be found in Figures 2 and 3, respectively. Since biological phenomenon associated with stage structure and toxin substances extensively exists within competitive system, theoretical results obtained in this paper are theoretically beneficial to discuss dynamic effect of maturation delay and toxic effect on population dynamics; it makes this work done in this paper has some positive and new features.



---

*Source: 689573-2014-10-19.xml*