# On the Dynamical Behavior of Toxic-Phytoplankton-Zooplankton Model with Delay

Authors: Mehbuba Rehim; Weixin Wu; Ahmadjan Muhammadhaji
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

A toxin producing phytoplankton-zooplankton model with inhibitory exponential substrate and time delay has been formulated and analyzed. Since the liberation of toxic substances by phytoplankton species is not an instantaneous process but is mediated by some time lag required for maturity of the species and the zooplankton mortality due to the toxic phytoplankton bloom occurs after some time laps of the bloom of toxic phytoplankton, we induced a discrete time delay to both of the consume response function and distribution of toxic substance term. Furthermore, based on the fact that the predation rate decreases at large toxic-phytoplankton density, the system is modelled via a Tissiet type functional response. We study the dynamical behaviour and investigate the conditions to guarantee the coexistence of two species. Analytical methods and numerical simulations are used to obtain information about the qualitative behaviour of the models.


---

## Body

## 1. Introduction

Phytoplankton are one of the most important components of the marine ecosystem. They not only form a basis for all aquatic food chains but also perform a very useful service by producing a huge amount of oxygen for human and other living animals after absorbing carbon-dioxide from surrounding environments [1]. Zooplankton are microscopic animals that eat other plankton and serve as a most favorable food source for fish and other aquatic animals. During the recent years, many authors have studied the system between zooplankton and phytoplankton. Authors in [2] have dealt with a nutrient-plankton model in an aquatic environment in the context of phytoplankton bloom. In [3] the effect of seasonality and periodicity on plankton dynamics is investigated. In [4], two plankton ecosystem models with explicit representation of viruses and virally infected phytoplankton are presented.The most common features of the phytoplankton population is rapid increase of biomass due to rapid cell proliferation and almost equally rapid decrease in populations, separated by some fixed time period. This type of rapid change in phytoplankton population density is known as “bloom” [5]. Due to the accumulation of high biomass or to the presence of toxicity, some of these blooms, more adequately called “harmful algal blooms” [6], are noxious to marine ecosystems or to human health and can produce great socioeconomic damage. There has been a global increase in harmful plankton blooms in last two decades [7–9].Because of the difficulty in measuring plankton biomass, mathematical modeling of plankton population is an important alternative method of improving our knowledge of the physical and biological processes relating to plankton ecology [10–15]. In [16], nutrient-plankton-zooplankton interaction models with a toxic substance which inhibits either the growth of phytoplankton, zooplankton, or both trophic levels are proposed and studied. In [17], the authors have constructed a mathematical model for describing the interaction between a nontoxic and toxic phytoplankton with a single nutrient.Based on the fact that(1) the release of toxin from phytoplankton species is not an instantaneous process but is mediated by some time lag required for maturity of the species and (2) the zooplankton may die after some time lapse of the bloom of toxic phytoplankton (see http://www.mote.org/, http://www.mdsg.umd.edu/), models incorporating time delay in diverse biological models are extensively reviewed by Beretta and Kuang [18], Gopalsamy [19], Cooke and Grossman [20], and Cushing [21]. The discrete time delay has potential to change the qualitative behavior of the dynamical systems [22–29]. Chattopadhyay et al. proposed a delay model incorporating time lag in toxin liberation by phytoplankton to avoid predation by zooplankton [28]. In [28], the authors introduced distribution delay and discrete to toxin liberation term. Due to discrete time delay in toxin liberation, the local existence of periodic solution through Hopf bifurcation has been obtained in [28].In [5, 10, 28–31], the following plausible toxic-phytoplankton-zooplankton system has been studied:
(1)dPdt=rP(1-P(t)K)-μf(P(t))Z(t),dZdt=βf(P(t))Z(t)-DZ(t)-ρgPtZt,
where P=P(t) is the density of the toxin-producing phytoplankton (TPP) population and Z=Z(t) is the density of zooplankton population at any instant of time t. In model (1), r is the intrinsic growth rate and K is the environmental carrying capacity of TPP population. The term f(P) describes the functional response for the grazing of phytoplankton by zooplankton, and μ is the maximum uptake rate for zooplankton. β denotes the ratio of biomass conversion and D is the natural death rate of zooplankton. The function g(P) represents the distribution of toxic substance which ultimately contributes to the death of zooplankton populations, and the parameter ρ denotes the rate of toxic substances produced by per unit biomass of phytoplankton.Model (1) has been studied for the following cases: (1) when f(P) is linear but g(P) is Holling type II [10, 28, 29] or Holling type III [10]; (2) when f(P) and g(P) are both Holling type II [5, 10] or Holling type III [10]; (3) when g(P) is Holling type II while f(P) is Holling type III [10].In above cases,f(P) and g(P) are both increasing functions of P over the entire interval 0<P<∞. However, in some cases, very high substrate concentrations in the lakes actually inhabit the growth of phytoplankton cells. Moreover, with the substrate concentrations increasing unlimitedly, some kind of microorganism will die eventually [32]. To describe the above phenomenon accurately, we consider f(P) from a different point of view. Adopting the idea used in [32], we assume that there exists a constant 0<P*<∞ such that f(P) is increasing over the interval 0≤P<P* and is decreasing on the interval P*<P<∞. More precisely, we use the so-called Tissiet functional response of the form of f(P)=μPe-P/K/(α+P) (see, e.g., [32]). This type of functional response takes care of the fact that the predation rate decreases at large toxic-phytoplankton density.The remainder of this paper is organized as follows. The model is described in Section2. In Section 3, we state and prove the positivity and boundedness of the solutions. Then, in Section 4, equilibria, their existence, and local asymptotic stability are considered. Under the aids of numerical simulation method, we further analyse the model and determine if there is a parameter range for the delay parameter τ, where more complicated dynamics occur. In Section 5, the permanence of the system is discussed by some analytic techniques on limit sets of differential dynamical systems. Finally, a brief discussion is presented in Section 6.


## 2. State of the Model

In a real ecological context, the interaction between phytoplankton and zooplankton will not be essentially instantaneous. Instead, the response of zooplankton to contacts with phytoplankton is likely to be delayed due to a gestation period. Another fact, during the interaction between phytoplankton and zooplankton, is that the liberation of toxic substances by phytoplankton must be mediated by some time lag which is required for the maturity of toxic-phytoplankton. Letτ1 stand for the time delay in conversion of food to viable biomass for the species, and τ2 is the discrete time period required for the maturity of phytoplankton cells to liberate toxic substances. Based on model 3 in [5], we intend to study a model system with the assumption that f(P) and g(P) are described by same type of function, namely, Tissiet functional response. Then we arrive at the following model:
(2)dPdt=rP(t)(1-P(t)K)-μP(t)Z(t)α+P(t)e-P(t)/K,dZdt=βP(t-τ1)Z(t-τ1)α+P(t-τ1)e-P(t-τ1)/K-DZ(t)-ρP(t-τ2)Z(t-τ2)α+P(t-τ2)e-P(t-τ2)/K,
where the parameter r is the intrinsic growth rate and K is the environmental carrying capacity of TPP population. The constant μ is the maximum per capita grazing rate, β denotes the ratio of biomass conversion, ρ denotes the rate of toxic substances produced by per unit biomass, D is the natural death rate of the zooplankton. All the parameters in system (2) are positive constants with their usual ecological meanings.Here we observe that if there is no delay (i.e.,τ1=τ2=0) and β-D-ρ<0, then Z′(t)<0. Hence, throughout our analysis, we assume that β-D-ρ>0.For the sake of simplicity of mathematical analysis, in this paper, we consider model (2) for the special case: the time delay in conversion of food to viable biomass for the species equal to the discrete time period required for the maturity of phytoplankton cells to liberate toxic substances; that is, τ1=τ2=τ.By performing the following scaling for model (2):
(3)t*=rt,P*(t*)=P(t)K,Z*(t*)=μZ(t)rK,α*=αK,β*=βμ,ρ*=ρμ,τ*=rτ,D*=Dμ,
we obtain a dimensionless system in the state variables P*, Z*, which can be written, by removing the stars, in the form
(4)P′(t)=P(t)(1-P(t))-P(t)Z(t)α+P(t)e-P(t),Z′t=βP(t-τ)Z(t-τ)α+P(t-τ)e-P(t-τ)-DZ(t)-ρP(t-τ)Z(t-τ)α+P(t-τ)e-P(t-τ),
with initial conditions:
(5)P(θ)=ϕ(θ)≥0,Z(θ)=φ(θ)≥0,θ∈[-τ,0],
where (ϕ(θ),φ(θ))∈C([-τ,0],R+2).


## 3. Positive and Boundedness

In this section, we consider the positive and boundedness of the solutions of model (4) with initial condition (5). One has the following theorem.Theorem 1.
Letϕ(θ)≥0, φ(θ)≥0 on -τ≤θ<0 and P(0)>0, Z(0)>0. Then(a)
all the solutions of (4) with initial condition (5) exist on [0,σ) for some constant σ>0 and are unique and positive for 0<t<σ,(b)
limsup⁡t→+∞P(t)≤1, limsup⁡t→+∞Z(t)≤M, where M=(β-ρ)(1+D)2/4D,(c)
ifβ<ρ+4αD/(1+D)2, then liminf⁡t→+∞P(t)≥m0, where m0=1-M/α. Further, the subset
(6)Γ={(ϕ(θ),φ(θ))∈C∣m0≤ϕ(θ)≤1,φ(θ)≥0}is positively invariant with respect to (4).Proof.
By Theorems  2.1 and 2.3 in Hale and Lunel [33], solutions of (4) with initial data ∈IntC+2 exist on 0<t<σ for some σ>0 and are unique. Suppose (P(t),Z(t)) is a solution of (4) for t∈[0,σ). Without loss of generality, assume that t∈[0,σ) is the maximum internal of the solution and σ=∞ if the solution exists for any t>0. Integrating the first equation of (4) gives
(7)P(t)=ϕ(0)exp⁡[∫0t(1-P(u)-Z(u)α+P(u)e-P(u))du]>0,t∈0,σ.
To prove the Z(t)>0 for any t∈[0,σ), use the method of contradiction. Suppose there exists a t*∈[0,σ) such that
(8)Zt*=0,Z′t*≤0,Zt>0foranyt∈[0,t*).
From the second equation of system (4), we have
(9)Z′(t*)=(β-ρ)P(t*-τ)Z(t*-τ)e-P(t*-τ)α+P(t*-τ)>0.
This contradicts Z′(t*)≤0. Hence, Z(t)>0 for all t∈[0,σ). This completes the proof of conclusion (a) in Theorem 1.
It follows from the first equation of (4) that P′(t)≤P(t)(1-P(t)), which implies that limsup⁡t→∞P(t)≤1. Define W(t)=(β-ρ)P(t)+Z(t+τ), t≥0. Then from (4) we obtain
(10)dWdt=-DW(t)+(β-ρ)P(t)(1+D-P(t))≤-DW(t)+(β-ρ)1+D24.
Applying the theorem of differential inequality we obtain that
(11)0<WP,Z≤(β-ρ)1+D24D1-e-Dt<WP,Z≤+W(p(θ),Z(0))e-Dt.
Therefore, limsup⁡t→+∞Z(t)≤M. Thus, there is a constant T>0, such that Z(t)≤M for all t≥T. This completes the proof of conclusion (b) in Theorem 1.
From the first equation of system (4) we get
(12)P′t≥Pt1-Pt-MαPt=Pt1-Mα-Pt,t≥T
which implies that liminf⁡t→+∞P(t)≥1-M/α=m0.
For anyψ=(ϕ,φ)∈Γ, let (P(t),Z(t)) be the solution of (4) with the initial function ψ. If there is a t1>0 such that P(t1)>1, then P˙(t0)>0 for some t0∈(0,t1) and P(t0)=1. Hence, it follows from the first equation of (4) that
(13)P˙(t0)=P(t0)(1-P(t0))-P(t0)Z(t0)(α+P(t0))eP(t0)<0,
which is a contradiction to P˙(t0)>0. So, P(t)≤1 for all t≥0.
It is easy to prove that if there is at2>0 such that P(t2)=m0, then for all t>t2 one has P(t)>m0. Which implies that P(t)≥m0 for all t≥0. This completes the proof of conclusion (c) in Theorem 1.Remark. The conclusions (b) and (c) in Theorem 1 indicate that if the ratio of biomass conversion is less than certain value, then phytoplankton population will be persistent.


## 4. Equilibrium, Stability, and Hopf Bifurcation

### 4.1. Existence of Equilibria

It is easy to see that model (4) has two boundary equilibria E0(0,0) and E1(1,0). To discuss the existence of the positive equilibria, we work on the equation
(14)Z=(1-P)(α+p)eP,(β-ρ)Pα+Pe-P-D=0.
Denote the left-hand side of the second equation in (14):
(15)F(P)=(β-ρ)Pα+Pe-P-D,P∈[0,1].
Then F(0)=-D<0, limP→+∞F(P)=-D,
(16)F′P=(β-ρ)[(α+P)-P(1+(α+P))]α+P2eP=β-ρα+P2eP-P2-αP+α,P∈[0,1].
Let G(P)=-P2-αP+α, we have G(p1)=G(p2)=0, where
(17)p1=-α2+4α-α2<0,0<p2=α2+4α-α2<1.
From this we know that the function F(P) is monotone increasing in the interval [0,p2] and monotone decreasing in the interval [p2,+∞). Hence, F(P) reaches its maximum Fm=((β-ρ)/(α+p2))e-p2-D on the interval (0,+∞). Therefore, we have the following results.Theorem 2.
The following statements hold.(1)
IfF(p2)<0, then equation F(P)=0 has no roots in the interval [0,+∞) (Figure 1(a)). In this case, for system (4), there is no positive equilibrium.(2)
IfF(p2)>0 and F(1)>0, then equation F(P)=0 has only one root, P=p1*, in the interval [0,1], where 0<p1*<1 (Figure 1(c)). In this case, for system (4), there exists a single positive equilibrium E*(p1*,z1*).(3)
IfF(p2)>0 and F(1)<0, then equation F(P)=0 has two distinct roots, P=p12* and P=p21*, in the interval [0,1], where 0<p12*<p2<p21*<1 (Figure 1(d)). In this case, for system (4), there are two distinct positive equilibria E*(p12*,z12*) and E*(p21*,z21*).(4)
IfF(p2)=0, then equation F(P)=0 has a unique root, p2*=p2, in the interval [0,+∞), 0<p2*=p2<1 (Figure 1(b)). In this case, for system (4), there is a unique positive equilibrium E*(p2*,z2*).(5)
IfF(p2)>0 and F(1)=0, then equation F(P)=0 has two distinct roots, P=p3* and P=1, in the interval (0,+∞), where 0<p3*<p2<1 (Figure 1(e)). In this case, for system (4), there is also unique positive equilibrium E*(p3*,z3*),where zi*, i=1,2,3,12,21, are determined by the first equation in (14).The graph of functionF(P)=0.
(a)
(b)
(c)
(d)
(e)


### 4.2. Local Stability and Hopf Bifurcation

LetE*(p*,z*) denote any one of the equilibrium points of system (4). Linearizing (4) about E* we obtain
(18)P′t=[1-2p*-z*α+p*(αα+p*-p*)e-p*]P(t)-p*α+p*e-p*Z(t),z′t=(β-ρ)(1-p*)(αα+p*-p*)P(t-τ)-DZ(t)+(β-ρ)p*α+p*e-p*Z(t-τ).
In what follows, existence of the interior equilibria, dynamical properties of the zooplankton-free equilibrium, and the interior equilibrium are investigated.Theorem 3.
For any time delayτ≥0, E1(1,0) is locally asymptotically stable if F(1)<0 and is unstable if F(1)>0; the trivial solution of the linearized system of (4) about E1(1,0) is stable for F(1)=0.Proof.
From (18) we know that, at any equilibrium point E(p*,z*), the characteristic equation is
(19)λ2+a0λ+a1+(a2-Dλ)e-λτ=0,
where
(20)a0=D-P*α+P*H(P*),a1=-Dp*α+p*H(p*),a2=D[p*α+p*H(p*)+(1-p*)(αα+p*-p*)],
in which H(p*)=[-p*2-(α+1)p*+1].
We first consider local stability ofE1(1,0). For E1(1,0), (19) becomes
(21)(λ+1)(λ+D-(β-ρ)(α+1)ee-λτ)=0.
One of the roots of (21) is λ1=-1<0, and other roots are given by solution of the following equality:
(22)(D+λ)eλτ=(β-ρ)(α+1)e.
It follows from [34] that (i)
ifF(1)<0, then all roots of (22) have negative parts for any time delay τ≥0. Hence, E1(1,0) is locally asymptotically stable for any time delay τ≥0;(ii)
ifF(1)>0, then (22) has roots which have positive real parts for any time delay τ≥0. Hence, E1(1,0) is unstable for any time delay τ≥0;(iii)
ifF(1)=0, this is a critical case and (22) is equivalent to
(23)λ+D(1-e-λτ)=0.From [34], it has that except λ=0, any root of (23) has negative real part for any time delay τ≥0. Hence, the trivial solution of the linearized system of (4) about E1(1,0) is stable for any time delay τ≥0.Theorem 4.
Suppose that the positive equilibriaE*(pi*,zi*), i=1,2.12,21, of system (4) exist. For system (4), the following statements hold. (i)
Ifpi*>Δα and 2p*H(pi*)+(1-pi*)G(pi*)>0, then E*(pi*,zi*) are locally asymptotically stable for τ<τ0(i), E*(pi*,zi*) are unstable for τ>τ0(i), and there is a periodic solution around E*(pi*,zi*) for τ=τ0(i). If 0<pi*<Δα, then E*(pi*,zi*) is unstable for any τ≥0, where τ0(i) is determined by (30), i=1,12.(ii)
E*(p21*,z21*) is unstable for any τ≥0.(iii)
Ifp2*>Δα, then the trivial solution of the linearized system of (4) about E*(p2*,z2*) is stable for any τ≥0. If 0<p2*<Δα, then the trivial solution of the linearized system of (4) about E*(p2*,z2*) is unstable for any τ≥0,where Δα=((α+1)2+4-α+1)/2, 0<Δα<1.Proof.
Whenτ=0, the transcendental equation (19) becomes
(24)λ2+(a0-D)λ+a1+a2=0,
where
(25)a0-D=-pi*α+pi*Hpi*,a1+a2=D1-pi*α+pi*Gpi*,i=1,2,12,21.H(pi*)>0 if 0<pi*<Δα, and H(pi*)<0 if Δα<pi*≤1, i=1,2,12,21.
By the argument in Section4.1, as long as E*(pi*,zi*) exists, it must be G(pi*)>0, i=1,12; thus, a1+a2>0; as long as E*(p21*,z21*) exists, it must be G(p21*)<0; thus, a1+a2<0; as long as E*(p2*,z2*) exists, it must be G(p2*)=0; thus, a1+a2=0; hence, from Routh-Hurwitz theorem, we have that if pi*>Δα, i=1.12, then E*(p1*,z1*) and E*(p12*,z12*) are locally asymptotically stable. If 0<pi*<Δα, i=1,12, then E*(p1*,z1*) and E*(p12*,z12*) are unstable; if pi*=Δα, i=1,12, then E*(p1*,z1*) and E*(p12*,z12*) are nonhyperbolic equilibria; E*(p21*,z21*) is unstable; E*(p2*,z2*) is a critical case.
Suppose thatλ=iω, ω>0, is a root of (19) for some τ>0. Substituting λ=iω in the characteristic equation (19), and separating real and imaginary parts, then we get
(26)-ω2+a1+a2cos⁡(ωτ)-ωDsin(ωτ)=0,a0ω-ωDcos⁡(ωτ)-a2sin(ωτ)=0.
Eliminating τ from (26), we get a biquadratic equation in τ as
(27)ω4+(a02-2a1-D2)ω2+a12-a22=0.
Its roots are
(28)ω±2=12{D2+2a1-a02±a02-2a1-D22-4(a12-a22)}.
By (20) we obtain
(29)a02-2a1-D2=D-pi*α+pi*Hpi*2+2Dpi*α+pi*H(pi*)-D2=pi*α+pi*Hpi*2>0,a12-a22=-D21-pi*α+pi*2Gpi*·2p*Hpi*+1-pi*Gpi*,mmmmmclmmmmcci=1,2,12,21.(1) We consider the stability of E*(pi*,zi*),  i=1,12.
SinceG(pi*)>0, if 2p*H(pi*)+(1-Pi*)G(pi*)>0,  i=1,12, then there is only one root, ω=ω+,  ω+>0, such that
(30)ω±2=12{D2+2a1-a02+a02-2a1-D22-4(a12-a22)}.
That is, (19) has one imaginary solution, λ=iω+,  ω+>0. From (26), we obtain the following set of τ for which there are imaginary root:
(31)τn(i)=θω++2nπω+,i=1,12,n=0,1,2…,
where 0≤θ≤2π and
(32)sinθ=-Dω+2+a0a2ω++a1Da22+D2,cos⁡θ=ω+2a2-Da0ω+-a1a2a22+D2.
If pi*>Δα, then H(pi*)<0, when τ=0, E*(pi*,zi*), i=1,12 are locally asymptotically stable. Hence, if pi*>Δα and 2p*H(pi*)+(1-pi*)G(pi*)>0, then E*(pi*,zi*) are locally asymptotically stable for τ<τ0(i) (n=0,  i=1,12).
Differentiating (19) with respect to τ, we obtain
(33)[2λ+a0-De-λτ-(a2-Dλ)τe-λτ]dλdτ-(a2-Dλ)λe-λτ=0.
It follows that
(34)dλdτ=(a2-Dλ)λe-λτ2λ+a0-De-λτ-(a2-Dλ)τe-λτ,(35)dλdτ-1=(a0+2λ)eλτ-D(a2-Dλ)λ-τλ.
From (35), we have
(36)Redλdτ-1λ=iω+=Rea0+2λeλτ-Da2-Dλλλ=iω+=1D2ω+4+a22ω+2a2ω+2ω+cos⁡(ω+τ)+a0sin(ω+τ)mmmmmmmmmmmm+Dω+2a0cos⁡⁡ω+τmmmmmmmmmmmmmmmcm-2ω+sin(ω+τ)-D=1D2ω+4+a22ω+2(2ω+4-2a1ω+2+a02ω+2-D2ω+2)=ω+2D2ω+4+a22ω+2(2ω+2+a02-2a1-D2).
Thus, we have
(37)sign[Redλdτλ=iω+]=sign{Redλdτ-1λ=iω+}>0.
Hence, there is a Hopf bifurcation at ω=ω+, τ=τ0(i). Therefore, if pi*>Δα and 2p*H(pi*)+(1-pi*)G(pi*)>0, then E*(pi*,zi*) is locally asymptotically stable for τ<τ0(i), E*(pi*,zi*) is unstable for τ>τ0(i), and there is a periodic solution around E*(pi*,zi*) for τ=τ0(i). If pi*>Δα and 2P*H(pi*)+(1-pi*)G(pi*)<0 or 0<pi*<Δα, then E*(pi*,zi*) is unstable for any τ≥0,  i=1,12.
(
2
) We consider the stability of E*(p21*,z21*).
SinceG(p21*)<0, if p21*>Δα, then (19) has only one imaginary solution. From proof of E*(p1*,z1*), sign[Re(dλ/dτ)λ=iω+]>0. Hence, E*(p21*,z21*) is unstable for any τ≥0. If 0<p21*<Δα and 2p*H(p21*)+(1-p21*)G(p21*)<0, then (19) has only one imaginary solution. If 0<p21*<Δα and 2P*H(p21*)+(1-p21*)G(p21*)>0, then (19) does not have imaginary solution. Wnen τ=0,  E*(p21*,z21*) is unstable. Therefore, E*(p21*,z21*) is unstable for any τ≥0.
(
3
) We consider the stability of E*(p2*,z2*).
Consider the characteristic equation (19). λ(τ)=0 is a root of (19) for all τ≥0. Suppose λ=u+iv is a root of (19), then substituting λ=u+iv in the characteristic equation, we have
(38)u+iv2+a0(u+iv)+a1+(a2-D(u+iv))e-(u+iv)τ=0
and separating real and imaginary parts, then we get
(39)u2-v2+a0u+a1+e-uτ[(a2-Du)cos⁡(vτ)-Dvsinvτ]=0,2uv+a0v-e-uτ[Dvcos⁡(vτ)+(a2-Du)sin(vτ)]=0.
From (39), we obtain
(40)u2-v2+a0u+a12+v22u+a02=e-2uτ[D2v2+a1+Du2].
(i) Assume p2*≥Δα. Suppose (19) has a root λ=u+iv,  u>0 for all τ≥0. Then, from (40), we have
(41)u2-v2+a0u+a12+v22u+a02-[D2v2+a1+Du2]=(e-2uτ-1)[D2v2+a1+Du2].
The left-hand size of (41) is equal to
(42)u2-v22+4u2v2+v2(a02-2a1-D2)+2a1u(a0-D)+u2(a02+2a1-D2).
Obviously,
(43)a02-2a1-D2=p2*α+p2*Hp2*2>0,2a1u(a0-D)>0,
and by assumption,
(44)a02+2a1-D2=p2*α+p2*H(p2*)[p2*α+p2*H(p2*)-4D]>0.
Therefore,
(45)u2-v2+a0u+a12+v22u+a02-[D2v2+a1+Du2]<0.
According to (41), this is impossible, since we are assuming u>0. Hence, all roots of (19) have nonpositive real parts, and this implies that the trivial solutions of the linearized system of (4) about E*(p2*,z2*) are stable.
(ii) Assume0<p2*<Δα. Consider the following real function:
(46)f(λ,τ)=λ2+a0λ+a1+(a2-Dλ)e-λτ.
It is clear that f(0,τ)=0 and limλ→+∞f(λ,τ)=+∞. There exists a λ0>0 such that if λ≥λ0, f(λ,τ)≥0, we also have
(47)∂f(λ,τ)∂λ=2λ+a0-[τ(a2-Dλ)+D]e-λτ,∂f(0,τ)∂λ=a0-τa2-D=-P2*α+p2*(H(p2*)+Dτ).
Since 0<p2*<Δα yields H(p2*)>0, which gives ∂f(0,τ)/∂λ<0(τ≥0), hence, there exists a δ(τ)>0 such that when 0<λ≤δ(τ), f(λ,τ)<0. Therefore, there exists at least one λ*, δ<λ*≤λ0, such that f(λ*,τ)=0, that is, (19) has at least a positive root. This indicates that trivial solution of the linearized system of (18) about E*(p2*,z2*) is unstable.


### 4.3. Numerical Simulation

For numerical simulation, let us first consider case (i) in Theorem4. We choose parameter values as α=3, β=1, ρ=0.3, and D=0.05. For this choice, F(p2)≈0,0630>0, F(1)≈0,0144>0. The stability of the interior equilibrium point E1*(0.3315,3.1025) depends upon the magnitude of delay τ. From (30) and (31) we obtain ω+2=0.3486, τ0=3.5. Consider Δα≈0.2361<0.3315, 2P*H(p1*)+(1-p1*)G(p1*)≈0.97932>0. When τ=0.1<τ0, Figures 2(a) and 2(b) show that the system is approaching the equilibrium point E1*. When τ increases to τ=3.57, E1* becomes unstable (Figure 2(c)). Stable oscillations appear when τ=3.5 and the Hopf bifurcation periodic solutions are depicted in Figure 3.((a)-(b)) Behavior and phase portrait of system (4) with τ=0.1<τ0. The positive equilibrium E1*(0.3315,3.1025) is stable. (c) The oscillations of P and Z versus time t when τ=3.57>τ0 show that the positive equilibrium E1* is unstable. Here τ0=3.5 and the initial value is (ϕ(0),φ(0))=(0.7,3).
(a)
(b)
(c)Figure 3
Behavior and phase portrait of system (4) when τ=τ0=3.5. Hopf bifurcation occurs from the interior equilibrium E1*. The initial value is (ϕ(0),φ(0))=(0.33,3.1).For case (ii) in Theorem4, consider the following set of parametric values: α=1, β=1, ρ=0.1, and D=0.179. Then we have F(p2)=0.0063>0, F(1)=-0.0135<0. By Theorem 4, the system (4) has two equilibria E12*(p12*,z12*) and E21*(p21*,z21*). From direct calculation, we get that p12*=0.4599, z12*=1.2489, p21*=0.8128, z21*=0.7650, τ0=5.11, Δα=0.414, and 2p*H(p12*)+(1-p12*)G(p12*)≈0.0567>0. For τ=1<τ0, the numerical simulation shows that equilibrium point E12*(0.4599,1.2489) is stable (Figure 4). When τ increases to τ=8>τ0, the equilibrium E12*(0.4599,1.2489) becomes unstable (Figure 6). There is a periodic solution around the equilibrium E12*=(0.4599,1.2489) when τ=5.11 (Figure 5).Figure 4
Behavior and phase portrait of system (4) with τ=1<τ0 shows the equilibrium E12*(0.4599,1.2489) is stable. Here τ0=5.11 and the initial value is (ϕ(0),φ(0))=(0.41,1.25).Figure 5
Behavior and phase portrait of system (4) with τ=5.11=τ0. Hopf bifurcation occurs from the interior equilibrium E12*. The initial value is (ϕ(0),φ(0))=(0.45,1.23).Figure 6
Behavior and phase portrait of system (4) with τ=8>τ0. The interior equilibrium E12* is unstable. Here τ0=5.11 and the initial value is (ϕ(0),φ(0))=(0.46,1.25).For case (iii) in Theorem4, let us take the parametric values α=1, β=1, ρ=0.1, and D=0.1853. For these set of values, by direct calculation, we get F(p2)=0, p2*=p2=0.618, and z2*=1.1467p2*>Δα=0.414. By Theorem 1, system (4) has only one equilibrium E2*(p2*,z2*), and from Theorem 4 we know that the equilibrium E2*(p2*,z2*) is stable (Figure 7(a)). Take parametric values as α=0.4, β=1, ρ=0.1, and D=0.2724. Then we have F(p2)=0, P2*=p2=0.4633, Z2*=0.7374, and P2*<Δα=0.5207. By Theorem 2, system (4) has only one equilibrium E2*(0.4633,0.7374). The numerical simulation shows E2*(0.4633,0.7374) is unstable (Figure 7(b)).Behavior and phase portrait of system (4). (a) The equilibrium E2*(0.618,1.1467) is stable for any τ>0 when p2*>Δα=0.414. (b) The equilibrium E2*(0.4633,0.5207) is unstable for any τ>0 when p2*<Δα=0.414. The initial value is (ϕ(0),φ(0))=(0.41,1.25).
(a)
(b)


## 4.1. Existence of Equilibria

It is easy to see that model (4) has two boundary equilibria E0(0,0) and E1(1,0). To discuss the existence of the positive equilibria, we work on the equation
(14)Z=(1-P)(α+p)eP,(β-ρ)Pα+Pe-P-D=0.
Denote the left-hand side of the second equation in (14):
(15)F(P)=(β-ρ)Pα+Pe-P-D,P∈[0,1].
Then F(0)=-D<0, limP→+∞F(P)=-D,
(16)F′P=(β-ρ)[(α+P)-P(1+(α+P))]α+P2eP=β-ρα+P2eP-P2-αP+α,P∈[0,1].
Let G(P)=-P2-αP+α, we have G(p1)=G(p2)=0, where
(17)p1=-α2+4α-α2<0,0<p2=α2+4α-α2<1.
From this we know that the function F(P) is monotone increasing in the interval [0,p2] and monotone decreasing in the interval [p2,+∞). Hence, F(P) reaches its maximum Fm=((β-ρ)/(α+p2))e-p2-D on the interval (0,+∞). Therefore, we have the following results.Theorem 2.
The following statements hold.(1)
IfF(p2)<0, then equation F(P)=0 has no roots in the interval [0,+∞) (Figure 1(a)). In this case, for system (4), there is no positive equilibrium.(2)
IfF(p2)>0 and F(1)>0, then equation F(P)=0 has only one root, P=p1*, in the interval [0,1], where 0<p1*<1 (Figure 1(c)). In this case, for system (4), there exists a single positive equilibrium E*(p1*,z1*).(3)
IfF(p2)>0 and F(1)<0, then equation F(P)=0 has two distinct roots, P=p12* and P=p21*, in the interval [0,1], where 0<p12*<p2<p21*<1 (Figure 1(d)). In this case, for system (4), there are two distinct positive equilibria E*(p12*,z12*) and E*(p21*,z21*).(4)
IfF(p2)=0, then equation F(P)=0 has a unique root, p2*=p2, in the interval [0,+∞), 0<p2*=p2<1 (Figure 1(b)). In this case, for system (4), there is a unique positive equilibrium E*(p2*,z2*).(5)
IfF(p2)>0 and F(1)=0, then equation F(P)=0 has two distinct roots, P=p3* and P=1, in the interval (0,+∞), where 0<p3*<p2<1 (Figure 1(e)). In this case, for system (4), there is also unique positive equilibrium E*(p3*,z3*),where zi*, i=1,2,3,12,21, are determined by the first equation in (14).The graph of functionF(P)=0.
(a)
(b)
(c)
(d)
(e)


## 4.2. Local Stability and Hopf Bifurcation

LetE*(p*,z*) denote any one of the equilibrium points of system (4). Linearizing (4) about E* we obtain
(18)P′t=[1-2p*-z*α+p*(αα+p*-p*)e-p*]P(t)-p*α+p*e-p*Z(t),z′t=(β-ρ)(1-p*)(αα+p*-p*)P(t-τ)-DZ(t)+(β-ρ)p*α+p*e-p*Z(t-τ).
In what follows, existence of the interior equilibria, dynamical properties of the zooplankton-free equilibrium, and the interior equilibrium are investigated.Theorem 3.
For any time delayτ≥0, E1(1,0) is locally asymptotically stable if F(1)<0 and is unstable if F(1)>0; the trivial solution of the linearized system of (4) about E1(1,0) is stable for F(1)=0.Proof.
From (18) we know that, at any equilibrium point E(p*,z*), the characteristic equation is
(19)λ2+a0λ+a1+(a2-Dλ)e-λτ=0,
where
(20)a0=D-P*α+P*H(P*),a1=-Dp*α+p*H(p*),a2=D[p*α+p*H(p*)+(1-p*)(αα+p*-p*)],
in which H(p*)=[-p*2-(α+1)p*+1].
We first consider local stability ofE1(1,0). For E1(1,0), (19) becomes
(21)(λ+1)(λ+D-(β-ρ)(α+1)ee-λτ)=0.
One of the roots of (21) is λ1=-1<0, and other roots are given by solution of the following equality:
(22)(D+λ)eλτ=(β-ρ)(α+1)e.
It follows from [34] that (i)
ifF(1)<0, then all roots of (22) have negative parts for any time delay τ≥0. Hence, E1(1,0) is locally asymptotically stable for any time delay τ≥0;(ii)
ifF(1)>0, then (22) has roots which have positive real parts for any time delay τ≥0. Hence, E1(1,0) is unstable for any time delay τ≥0;(iii)
ifF(1)=0, this is a critical case and (22) is equivalent to
(23)λ+D(1-e-λτ)=0.From [34], it has that except λ=0, any root of (23) has negative real part for any time delay τ≥0. Hence, the trivial solution of the linearized system of (4) about E1(1,0) is stable for any time delay τ≥0.Theorem 4.
Suppose that the positive equilibriaE*(pi*,zi*), i=1,2.12,21, of system (4) exist. For system (4), the following statements hold. (i)
Ifpi*>Δα and 2p*H(pi*)+(1-pi*)G(pi*)>0, then E*(pi*,zi*) are locally asymptotically stable for τ<τ0(i), E*(pi*,zi*) are unstable for τ>τ0(i), and there is a periodic solution around E*(pi*,zi*) for τ=τ0(i). If 0<pi*<Δα, then E*(pi*,zi*) is unstable for any τ≥0, where τ0(i) is determined by (30), i=1,12.(ii)
E*(p21*,z21*) is unstable for any τ≥0.(iii)
Ifp2*>Δα, then the trivial solution of the linearized system of (4) about E*(p2*,z2*) is stable for any τ≥0. If 0<p2*<Δα, then the trivial solution of the linearized system of (4) about E*(p2*,z2*) is unstable for any τ≥0,where Δα=((α+1)2+4-α+1)/2, 0<Δα<1.Proof.
Whenτ=0, the transcendental equation (19) becomes
(24)λ2+(a0-D)λ+a1+a2=0,
where
(25)a0-D=-pi*α+pi*Hpi*,a1+a2=D1-pi*α+pi*Gpi*,i=1,2,12,21.H(pi*)>0 if 0<pi*<Δα, and H(pi*)<0 if Δα<pi*≤1, i=1,2,12,21.
By the argument in Section4.1, as long as E*(pi*,zi*) exists, it must be G(pi*)>0, i=1,12; thus, a1+a2>0; as long as E*(p21*,z21*) exists, it must be G(p21*)<0; thus, a1+a2<0; as long as E*(p2*,z2*) exists, it must be G(p2*)=0; thus, a1+a2=0; hence, from Routh-Hurwitz theorem, we have that if pi*>Δα, i=1.12, then E*(p1*,z1*) and E*(p12*,z12*) are locally asymptotically stable. If 0<pi*<Δα, i=1,12, then E*(p1*,z1*) and E*(p12*,z12*) are unstable; if pi*=Δα, i=1,12, then E*(p1*,z1*) and E*(p12*,z12*) are nonhyperbolic equilibria; E*(p21*,z21*) is unstable; E*(p2*,z2*) is a critical case.
Suppose thatλ=iω, ω>0, is a root of (19) for some τ>0. Substituting λ=iω in the characteristic equation (19), and separating real and imaginary parts, then we get
(26)-ω2+a1+a2cos⁡(ωτ)-ωDsin(ωτ)=0,a0ω-ωDcos⁡(ωτ)-a2sin(ωτ)=0.
Eliminating τ from (26), we get a biquadratic equation in τ as
(27)ω4+(a02-2a1-D2)ω2+a12-a22=0.
Its roots are
(28)ω±2=12{D2+2a1-a02±a02-2a1-D22-4(a12-a22)}.
By (20) we obtain
(29)a02-2a1-D2=D-pi*α+pi*Hpi*2+2Dpi*α+pi*H(pi*)-D2=pi*α+pi*Hpi*2>0,a12-a22=-D21-pi*α+pi*2Gpi*·2p*Hpi*+1-pi*Gpi*,mmmmmclmmmmcci=1,2,12,21.(1) We consider the stability of E*(pi*,zi*),  i=1,12.
SinceG(pi*)>0, if 2p*H(pi*)+(1-Pi*)G(pi*)>0,  i=1,12, then there is only one root, ω=ω+,  ω+>0, such that
(30)ω±2=12{D2+2a1-a02+a02-2a1-D22-4(a12-a22)}.
That is, (19) has one imaginary solution, λ=iω+,  ω+>0. From (26), we obtain the following set of τ for which there are imaginary root:
(31)τn(i)=θω++2nπω+,i=1,12,n=0,1,2…,
where 0≤θ≤2π and
(32)sinθ=-Dω+2+a0a2ω++a1Da22+D2,cos⁡θ=ω+2a2-Da0ω+-a1a2a22+D2.
If pi*>Δα, then H(pi*)<0, when τ=0, E*(pi*,zi*), i=1,12 are locally asymptotically stable. Hence, if pi*>Δα and 2p*H(pi*)+(1-pi*)G(pi*)>0, then E*(pi*,zi*) are locally asymptotically stable for τ<τ0(i) (n=0,  i=1,12).
Differentiating (19) with respect to τ, we obtain
(33)[2λ+a0-De-λτ-(a2-Dλ)τe-λτ]dλdτ-(a2-Dλ)λe-λτ=0.
It follows that
(34)dλdτ=(a2-Dλ)λe-λτ2λ+a0-De-λτ-(a2-Dλ)τe-λτ,(35)dλdτ-1=(a0+2λ)eλτ-D(a2-Dλ)λ-τλ.
From (35), we have
(36)Redλdτ-1λ=iω+=Rea0+2λeλτ-Da2-Dλλλ=iω+=1D2ω+4+a22ω+2a2ω+2ω+cos⁡(ω+τ)+a0sin(ω+τ)mmmmmmmmmmmm+Dω+2a0cos⁡⁡ω+τmmmmmmmmmmmmmmmcm-2ω+sin(ω+τ)-D=1D2ω+4+a22ω+2(2ω+4-2a1ω+2+a02ω+2-D2ω+2)=ω+2D2ω+4+a22ω+2(2ω+2+a02-2a1-D2).
Thus, we have
(37)sign[Redλdτλ=iω+]=sign{Redλdτ-1λ=iω+}>0.
Hence, there is a Hopf bifurcation at ω=ω+, τ=τ0(i). Therefore, if pi*>Δα and 2p*H(pi*)+(1-pi*)G(pi*)>0, then E*(pi*,zi*) is locally asymptotically stable for τ<τ0(i), E*(pi*,zi*) is unstable for τ>τ0(i), and there is a periodic solution around E*(pi*,zi*) for τ=τ0(i). If pi*>Δα and 2P*H(pi*)+(1-pi*)G(pi*)<0 or 0<pi*<Δα, then E*(pi*,zi*) is unstable for any τ≥0,  i=1,12.
(
2
) We consider the stability of E*(p21*,z21*).
SinceG(p21*)<0, if p21*>Δα, then (19) has only one imaginary solution. From proof of E*(p1*,z1*), sign[Re(dλ/dτ)λ=iω+]>0. Hence, E*(p21*,z21*) is unstable for any τ≥0. If 0<p21*<Δα and 2p*H(p21*)+(1-p21*)G(p21*)<0, then (19) has only one imaginary solution. If 0<p21*<Δα and 2P*H(p21*)+(1-p21*)G(p21*)>0, then (19) does not have imaginary solution. Wnen τ=0,  E*(p21*,z21*) is unstable. Therefore, E*(p21*,z21*) is unstable for any τ≥0.
(
3
) We consider the stability of E*(p2*,z2*).
Consider the characteristic equation (19). λ(τ)=0 is a root of (19) for all τ≥0. Suppose λ=u+iv is a root of (19), then substituting λ=u+iv in the characteristic equation, we have
(38)u+iv2+a0(u+iv)+a1+(a2-D(u+iv))e-(u+iv)τ=0
and separating real and imaginary parts, then we get
(39)u2-v2+a0u+a1+e-uτ[(a2-Du)cos⁡(vτ)-Dvsinvτ]=0,2uv+a0v-e-uτ[Dvcos⁡(vτ)+(a2-Du)sin(vτ)]=0.
From (39), we obtain
(40)u2-v2+a0u+a12+v22u+a02=e-2uτ[D2v2+a1+Du2].
(i) Assume p2*≥Δα. Suppose (19) has a root λ=u+iv,  u>0 for all τ≥0. Then, from (40), we have
(41)u2-v2+a0u+a12+v22u+a02-[D2v2+a1+Du2]=(e-2uτ-1)[D2v2+a1+Du2].
The left-hand size of (41) is equal to
(42)u2-v22+4u2v2+v2(a02-2a1-D2)+2a1u(a0-D)+u2(a02+2a1-D2).
Obviously,
(43)a02-2a1-D2=p2*α+p2*Hp2*2>0,2a1u(a0-D)>0,
and by assumption,
(44)a02+2a1-D2=p2*α+p2*H(p2*)[p2*α+p2*H(p2*)-4D]>0.
Therefore,
(45)u2-v2+a0u+a12+v22u+a02-[D2v2+a1+Du2]<0.
According to (41), this is impossible, since we are assuming u>0. Hence, all roots of (19) have nonpositive real parts, and this implies that the trivial solutions of the linearized system of (4) about E*(p2*,z2*) are stable.
(ii) Assume0<p2*<Δα. Consider the following real function:
(46)f(λ,τ)=λ2+a0λ+a1+(a2-Dλ)e-λτ.
It is clear that f(0,τ)=0 and limλ→+∞f(λ,τ)=+∞. There exists a λ0>0 such that if λ≥λ0, f(λ,τ)≥0, we also have
(47)∂f(λ,τ)∂λ=2λ+a0-[τ(a2-Dλ)+D]e-λτ,∂f(0,τ)∂λ=a0-τa2-D=-P2*α+p2*(H(p2*)+Dτ).
Since 0<p2*<Δα yields H(p2*)>0, which gives ∂f(0,τ)/∂λ<0(τ≥0), hence, there exists a δ(τ)>0 such that when 0<λ≤δ(τ), f(λ,τ)<0. Therefore, there exists at least one λ*, δ<λ*≤λ0, such that f(λ*,τ)=0, that is, (19) has at least a positive root. This indicates that trivial solution of the linearized system of (18) about E*(p2*,z2*) is unstable.


## 4.3. Numerical Simulation

For numerical simulation, let us first consider case (i) in Theorem4. We choose parameter values as α=3, β=1, ρ=0.3, and D=0.05. For this choice, F(p2)≈0,0630>0, F(1)≈0,0144>0. The stability of the interior equilibrium point E1*(0.3315,3.1025) depends upon the magnitude of delay τ. From (30) and (31) we obtain ω+2=0.3486, τ0=3.5. Consider Δα≈0.2361<0.3315, 2P*H(p1*)+(1-p1*)G(p1*)≈0.97932>0. When τ=0.1<τ0, Figures 2(a) and 2(b) show that the system is approaching the equilibrium point E1*. When τ increases to τ=3.57, E1* becomes unstable (Figure 2(c)). Stable oscillations appear when τ=3.5 and the Hopf bifurcation periodic solutions are depicted in Figure 3.((a)-(b)) Behavior and phase portrait of system (4) with τ=0.1<τ0. The positive equilibrium E1*(0.3315,3.1025) is stable. (c) The oscillations of P and Z versus time t when τ=3.57>τ0 show that the positive equilibrium E1* is unstable. Here τ0=3.5 and the initial value is (ϕ(0),φ(0))=(0.7,3).
(a)
(b)
(c)Figure 3
Behavior and phase portrait of system (4) when τ=τ0=3.5. Hopf bifurcation occurs from the interior equilibrium E1*. The initial value is (ϕ(0),φ(0))=(0.33,3.1).For case (ii) in Theorem4, consider the following set of parametric values: α=1, β=1, ρ=0.1, and D=0.179. Then we have F(p2)=0.0063>0, F(1)=-0.0135<0. By Theorem 4, the system (4) has two equilibria E12*(p12*,z12*) and E21*(p21*,z21*). From direct calculation, we get that p12*=0.4599, z12*=1.2489, p21*=0.8128, z21*=0.7650, τ0=5.11, Δα=0.414, and 2p*H(p12*)+(1-p12*)G(p12*)≈0.0567>0. For τ=1<τ0, the numerical simulation shows that equilibrium point E12*(0.4599,1.2489) is stable (Figure 4). When τ increases to τ=8>τ0, the equilibrium E12*(0.4599,1.2489) becomes unstable (Figure 6). There is a periodic solution around the equilibrium E12*=(0.4599,1.2489) when τ=5.11 (Figure 5).Figure 4
Behavior and phase portrait of system (4) with τ=1<τ0 shows the equilibrium E12*(0.4599,1.2489) is stable. Here τ0=5.11 and the initial value is (ϕ(0),φ(0))=(0.41,1.25).Figure 5
Behavior and phase portrait of system (4) with τ=5.11=τ0. Hopf bifurcation occurs from the interior equilibrium E12*. The initial value is (ϕ(0),φ(0))=(0.45,1.23).Figure 6
Behavior and phase portrait of system (4) with τ=8>τ0. The interior equilibrium E12* is unstable. Here τ0=5.11 and the initial value is (ϕ(0),φ(0))=(0.46,1.25).For case (iii) in Theorem4, let us take the parametric values α=1, β=1, ρ=0.1, and D=0.1853. For these set of values, by direct calculation, we get F(p2)=0, p2*=p2=0.618, and z2*=1.1467p2*>Δα=0.414. By Theorem 1, system (4) has only one equilibrium E2*(p2*,z2*), and from Theorem 4 we know that the equilibrium E2*(p2*,z2*) is stable (Figure 7(a)). Take parametric values as α=0.4, β=1, ρ=0.1, and D=0.2724. Then we have F(p2)=0, P2*=p2=0.4633, Z2*=0.7374, and P2*<Δα=0.5207. By Theorem 2, system (4) has only one equilibrium E2*(0.4633,0.7374). The numerical simulation shows E2*(0.4633,0.7374) is unstable (Figure 7(b)).Behavior and phase portrait of system (4). (a) The equilibrium E2*(0.618,1.1467) is stable for any τ>0 when p2*>Δα=0.414. (b) The equilibrium E2*(0.4633,0.5207) is unstable for any τ>0 when p2*<Δα=0.414. The initial value is (ϕ(0),φ(0))=(0.41,1.25).
(a)
(b)


## 5. Permanence

In this section, we will use the same methods as [35] to prove the permanence of system (4).Theorem 5.
Suppose thatβ<ρ+4αD/(1+D)2, D>1, and F(1)>0. Then, for any time delay τ≥0, system (4) is permanent.Proof.
From Theorem1, we see that it is enough to consider the solution (P(t),Z(t)) with initial condition ψ=(ϕ,φ)∈Γ. Theorem 1 implies that the omega limit set ω(ψ) of (P(t),Z(t))(t≥0) is nonempty, compact, and invariant and ω(ψ)⊂Γ. It follows from definition of permanence and Theorem 1, that we only need to show
(48)liminf⁡t→∞Z(t)≥υ.
Here υ is some positive constant which dose not depend on the initial function ψ.
Let us first show(49)liminf⁡t→∞Z(t)>0.
In fact, if (49) is not true, then, from liminf⁡t→∞Z(t)=0, we see that there exists a positive time sequence {tn}:tn→+∞(n→+∞) such that
(50)limtn→+∞Z(tn)=0,Z˙(tn)≤0,Z(t)≥Z(tn)tn-τ≤t≤tn.
Note that the solution (P(t),Z(t)) is bounded on [0,+∞) by Theorem 1. It follows from (4) that (P(t),Z(t)) is uniformly continuous on [0,+∞). Hence, from Ascoli’s theorem there is a subsequence of {tn}, still denoted by {tn}, such that
(51)limtn→+∞(P(tn),Z(tn))=(P~(t),Z~(t))
holds uniformly on R in the wide sense. From Theorem 1, we have that (P~(t),Z~(t))∈Γ for any t∈R and that, for any τ∈R, the function (P~(t+τ),Z~(t+τ)) of t is the solution of (4) with the initial function (P~(τ),Z~(τ)). Here we note that Z~(0)=0 and m0≤P~(t)≤1 for any t∈R.
We claim that(P~(t),Z~(t))=(1,0) for any t∈R. From Theorem 1, we know that if ϕ(0)>0, then the solution (P(t),Z(t)) of (4) exists and P(t)>0 and Z(t)>0(t≥0). Thus, from Z~(0)=0, we have that Z~(t)=0 for any t<0. Then, it follows from (4) that Z(t)≡0 for any t∈R, and P′~(t)=P~(t)(1-P~(t)) for any t≥τ. Hence,
(52)P~(t)=11+(1/P(0)-1)e-t.
Since P~(t) is bounded for t∈R, we must have P~(0)=1, which implies that P~(t)=1 for t∈R. It follows from (4) and the invariant of Γ that (P~(t),Z~(t))=(1,0) for any t∈R. This shows that above claim holds. In particular, we have that
(53)limn→+∞Z(tn-τ)=Z~(-τ)=0,limn→+∞P(tn-τ)=P~(-τ)=1,limn→+∞(β-ρ)P(tn-τ)(α+P(tn-τ))eP(tn-τ)=β-ρ(α+1)e>D.
For sufficiently small ε>0 and sufficiently large N, n>N, we have
(54)(β-ρ)P(tn-τ)(α+P(tn-τ))eP(tn-τ)>β-ρ(α+1)e-ε>D.
Hence, for large n, from (50) and above inequality, we obtain
(55)Z˙nt=(β-ρ)P(tn-τ)(α+P(tn-τ))eP(tn-τ)Z(tn-τ)-Z(tn)≥((β-ρ)P(tn-τ)(α+P(tn-τ))eP(tn-τ)-1)Z(tn)>((β-ρ)(α+1)e-ε-1)Z(tn)>0
which is a contradiction to Z˙n(t)≤0. This completes the proof of liminf⁡t→∞Z(t)>0.
Next, let us show(56)liminf⁡t→∞Z(t)≥υ>0.
For any initial functions sequence ψn={(ϕ(n),φ(n))}⊂Γ, let (P(n)(t),Z(n)(t)) be the solution of (4) with the initial function ψn. Let ωn(ψn) be the omega limit set of (P(n)(t),Z(n)(t)). We have that there exists some compact and invariant set ω*⊂Γ such that dist⁡(ωn(ψn),ω*)→0 as n→+∞. Here, dist⁡(ωn(ψn),ω*) means Hausdor ff distance [35, 36].
If (48) does not hold, for some initial function sequence ψn={(ϕ(n),φ(n))}⊂Γ such that φ(n)(0)>0, we have that there is some ψ~=(ϕ~,φ~)∈ω* such that φ~(θ0)=0 for some θ0∈[-τ,0]. Now, let (P~(t),Z~(t)) be the solution of (4) with the initial function ψ~. Then, by the invariant of ω*, we have that (P~(t+θ),Z~(t+θ))∈ω* for all t∈R and θ∈[-τ,0]. Note φ~(θ0)=0 and the positivity of all solutions, we easily have that Z~(t)=0 for all t≤θ0. Hence, it follows from (4) that Z~(θ)=0, θ∈[-τ,0], and Z~(t)=0 for all t∈R. This implies that P~(t)=[1+(P-1(0)-1)e-t]-1, Z~(t)=0, for all t∈R. If P~(0)<1, we see that the negative semiorbit (P~(t+θ),Z~(t+θ))∈ω*(t≤0) is unbounded. This is contradiction.
IfP~(0)=1, we have that P~(t)=1, Z~(t)=0 for all t∈R. This shows that ψ~=(1,0)=E0∈ω*. Let us show E0 is factually isolated [35, 36]. that is, there exists some neighborhood U of E0 in Γ such that E0 is the largest invariant set in U. In fact, let us choose
(57)U={ψ=(ϕ,φ)∈Γ,ψ-E0<ε}
for some sufficiently small positive constant ε and ε<1-D(α+1)e/(β-ρ). We will show that E0 is the largest invariant set in U for some ε.
If not, for any sufficiently smallε, there exists some invariant set Ω(Ω⊂U) such that Ω∖E0 is not empty. Let ψ=(ϕ,φ)∈Ω∖E0 and (P(t+θ),Z(t+θ)) be the solution of (4) with the initial function ψ. Then, (P(t+θ),Z(t+θ))∈Ω for all t∈R.
Ifφ(0)=0, by the invariance of Ω and Theorem 1, we also have the contradiction that ψ=E0 or that the negative semiorbit (P(t+θ),Z(t+θ))(t<0) of (4) through ψ is unbounded.
Ifφ(0)>0, from the Theorem 1, we see that Z(t)>0 for all t≥0. Now, let us consider the continuous function:
(58)Qt=Zt+η∫t-τtZθdθ,
for some constant η>1. Because of (P(t+θ),Z(t+θ))∈U(t∈R), we have 1-ε≤P(t)≤1(t∈R). The time derivative of Q(t) along the solution (P(t),Z(t)) satisfies
(59)Q˙t=Z˙(t)+η(Z(t)-Z(t-τ))=(η-D)Z(t)+[(β-ρ)P(t-τ)(α+P(t-τ))eP(t-τ)-η]Z(t-τ)≥(η-D)Z(t)+[(β-ρ)(1-ε)(α+1)e-η]Z(t-τ).
Since ε<1-D(α+1)e/(β-ρ) and F(1)>0, we have that (β-ρ)(1-ε)/(α+1)e>D. We can choose η>D, such that D<η<(β-ρ)(1-ε)/(α+1)e. From (49), there exists a constant t1>0 such that Z(t)≥δ>0 for some constant δ and t≥t1. Hence, it follows from (59)(60)Q˙(t)≥(η-D)δ>0,
for all t≥t1. Thus, Q(t)→+∞ as t→+∞. This contradicts Theorem 1 and shows that E0 is isolated.
We easily see that the semigroup defined by the solution of (4) satisfies the conditions of Lemma  4.3 in [36] with Ω=E0. Thus, by Lemma  4.3 in [36], we have that there is some ξ=(ξ1,ξ2), such that ξ∈ω*⋂‍(Ωs(E0)∖E0). Here, Ωs(E0) denotes the stable set of E0.
Ifξ2(0)=0, again by the invariant of Ω and Theorem 1, we also have the contradiction that ξ=E0 or that the negative semiorbit (P~t+θ,Z~t+θt<0,θ∈-τ,0) of (4) through ξ is unbounded.
Ifξ2(0)>0, from Theorem 1, we see that P~(t)>0, Z~(t)>0 for all t>0. It follows from ξ∈ω*⋂‍(Ωs(E0)∖E0) that limt→+∞P~(t)=1, limt→+∞Z~(t)=0, which contradicts (49). This shows that (48) holds. Thus, (4) is permanent. This completes the proof of Theorem 5.


## 6. Discussion

In a real ecological context, the interaction between phytoplankton and zooplankton will not be essentially instantaneous. Instead, the response of zooplankton to contact with phytoplankton is likely to be delayed due to a gestation period. Another fact, during the interaction between phytoplankton and zooplankton, is that the liberation of toxic substances by phytoplankton must be mediated by some time lag which is required for the maturity of toxic-phytoplankton. And in some cases, very high substrate concentrations in the lakes actually inhabit the growth of phytoplankton cells. Moreover, with the substrate concentrations increasing unlimitedly, some kind of microorganism will die eventually [32]. Based on the above fact, in this paper, we introduce time delay to a phytoplankton-zooplankton interaction model with exponential substrate uptake and exponential distribution of toxic substance term. The model (4) accounts for some natural phenomenon. By using comparison principle for functional differential equations and traditional analysis technique for transcendental equations [34], we give a detailed analysis on boundedness of solutions of system (4) and local asymptotic stability of the equilibria of system (4). Our results show that time delay is factually harmless for the local asymptotic stability of the zooplankton free equilibrium of (4), but it is not always harmless for the positive equilibrium; that is to say, because of the time delay the positive equilibrium becomes unstable (Theorem 4). Based on some known techniques on limit sets of differential dynamical systems, we show that, for any time delay, the phytoplankton-zooplankton interaction model is permanent if and only if one positive equilibrium exists.Under the aids of the numerical simulation we further investigate the delayed model system (4). Figures 1–6 illustrate that the delay τ plays crucial role in determining the asymptotic behavior of solutions of model (4). The stability or oscillatory coexistence depends upon not only the parametric restriction but also upon the gestation delay (liberation delay). Our system exhibits stable behavior when τ<τ0, where τ0 is the threshold value of the parameter τ. This threshold value of τ0 is determined by our numerical simulations. This value of τ0 is 3.5 for case (i) of Theorem 4 and is 5.1 for case (ii). These stable solutions are shown in Figures 2(a), 2(b), and 4. When τ crosses τ0, there is a delay induced instability demonstrated in Figures 2(c) and 5 for τ>τ0. Stable oscillations appear when τ=τ0 and the Hopf bifurcation periodic solution occurs (see Figures 3 and 5). Thus, there is a range of gestation delay (liberation delay) which initially imparts stability, then induces instability, and ultimately leads to periodic behavior.



---

*Source: 756315-2015-02-17.xml*