# Dynamical Behaviour of a Nutrient-Plankton Model with Holling Type IV, Delay, and Harvesting

Authors: Xin-You Meng; Jiao-Guo Wang; Hai-Feng Huo
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

In this paper, a Holling type IV nutrient-plankton model with time delay and linear plankton harvesting is investigated. The existence and local stability of all equilibria of model without time delay are given. Regarding time delay as bifurcation parameter, such system around the interior equilibrium loses its local stability, and Hopf bifurcation occurs when time delay crosses its critical value. In addition, the properties of the bifurcating periodic solutions are investigated based on normal form theory and center manifold theorem. What is more, the global continuation of the local Hopf bifurcation is discussed by using a global Hopf bifurcation result. Furthermore, the optimal harvesting is obtained by the Pontryagin’s Maximum Principle. Finally, some numerical simulations are given to confirm our theoretical analysis.


---

## Body

## 1. Introduction

In the aquatic ecosystem, all aquatic food chains are based on plankton, which are composed of phytoplankton and zooplankton. As pointed out in [1], phytoplankton in particular have great contribution for our earth; for example, they provide food for marine life and oxygen for human life, and they also absorb more than half the carbon dioxide which may be contributing to global warming.In the past decade, the issues of phytoplankton-zooplankton models have been investigated by many experts [2–13]. Chattopadhy et al. [2] showed the effect of the delay forming the period of toxic substances on blooms. Wang et al. [4] considered a toxic phytoplankton-zooplankton model with delay-dependent coefficient and investigated influence of delay and harvesting for the system. Yang et al. [8] discussed the impact of reaction-diffusion on the dynamic behaviours of plankton. Sharma et al. [11] studied the bifurcation behaviour analysis of a plankton model with multiple delays. Nevertheless, these models mainly included interactions of phytoplankton-zooplankton.However, nutrients play an significant role in the growth of plankton. The concentration of nutrients can better explain their mechanism in the phytoplankton-zooplankton system. As far as we know, the nutrient-plankton systems have been considered by a number of scholars [14–21]. Chakrabortyour et al. [14] demonstrated that the toxin producing phytoplankton (TPP) species control the bloom of other nontoxic phytoplankton species when nutrient-deficient conditions are favorable for TPP species to release toxic chemicals. Zhang and Wang [15] discussed Hopf bifurcation and bistability of nutrient-phytoplankton-zooplankton model. Fan et al. [16] proposed a nutrient model for a water ecosystem. Chakraborty et al. [17] investigated spatial dynamics of a nutrient-phytoplankton system with toxic effect on phytoplankton. Wang et al. [18] analyzed a nutrient-plankton model with delayed nutrient cycling. Jang and Allen [19] studied deterministic and stochastic nutrient-phytoplankton-zooplankton models with periodic toxin producing phytoplankton. Das and Ray [20] debated the effect of delay on nutrient cycling in phytoplankton-zooplankton interactions in estuarine system. Sharma et al. [21] discussed the effect of the discrete time delay on the nutrient-plankton interaction and proposed the model as follows:(1)dNdt=N0-αN-bNx+k1b1x,dxdt=α1Nx-b1x-βxyγ+x,dydt=β1xyγ+x-α2y-ρxt-τyγ+xt-τ,where N(t), x(t), and y(t) denote the concentration of nutrient, phytoplankton, and zooplankton population at time t, respectively. τ is the time delay which is incorporated with the assumption that the liberation of toxin is not instantaneous; rather it is mediated by some time lag. Li and Xiao [22] showed a predator-prey model with Holling type IV: p(u)=cu/(m+u2), where the predator is specifically not taking nutrient, and the function form of biomass conversion of herbivore is Holling type IV function. They also assumed that the extra mortality in zooplankton due to toxic substance term is also expressed by Holling type IV function.Generally speaking, there are some economic benefits for some plankton in aquatic food chain. Some phytoplankton such as phylum Diatom, phylum Pyrrophyta, and phylum Chrysophyta can be harvested as bait for aquaculture, and some zooplankton such as jellyfish, krill, andAcetes can be harvested for human food. Thus, these plankton play a crucial role in marine conservation and fishery management. As we all know, harvesting is an essential element in our life. On the one hand, the harvesting is one of the economic sources of human economy. On the other hand, it is a method to control the balance of ecosystem. Thus, zooplankton-phytoplankton system with harvesting has attracted the attention of scholars [23–29]. Chakraborty and Das [23] proposed a two-zooplankton one-phytoplankton system with harvesting and analyzed the impact of harvesting on the coexistence and competitive exclusion of competitive predators. Lv et al. [24] proposed a phytoplankton-zooplankton system with harvesting and pointed out that overharvesting could lead to population extinction and appropriate harvesting should guarantee the sustainability of the population types in an aquatic ecosystem.Some works have already been made to understand the nutrient-plankton system with time delay and toxic phytoplankton, but the study of nutrient-plankton system which contains Holling type IV function, time delay, toxic phytoplankton, and linear harvesting is rarely done so far. Taking account of the effect of the time delay in maturation of zooplankton into nutrient-plankton model can make the model more realistic. The main aim of the present study is to see whether the time delay in maturation of zooplankton and harvesting of plankton can exert an impact on the nutrient-plankton system.Based on the works [21, 22] and the above discussions, the following model is proposed by introducing harvesting and time delay(2)N˙t=N0-αN-bNP+k1P,P˙t=b1NP-kP-βPZr+P2-q1EP,Z˙t=β1Pt-τZr+P2t-τ-μZ-θPZr+P2-q2EZ,where N(t), P(t), and Z(t) denote the number of nutrient, phytoplankton, and zooplankton population at time t, respectively. N0 is the constant input of nutrient concentration and α is washout rate. b and b1(b1≤b) are the nutrient uptake rate and conversion rate of nutrient for growth of the phytoplankton, respectively. k and μ are the mortality rates of the phytoplankton and zooplankton population, respectively. β and β1(β1≤β) are the maximal ingestion rate and conversion rate of zooplankton. θ denotes the rate of toxin liberation. k1(k1≤k) is the nutrient recycle rate after death of phytoplankton population. The parameter r is the half saturation constant for a Holling type IV function. τ denotes the time delay in maturation of zooplankton. The constants q1 and q2 are the catchability coefficients of the phytoplankton and zooplankton population, and E is the harvesting effort.The initial conditions of system (2) take the following form(3)Nϱ=ϕ1ϱ≥0,Pϱ=ϕ2ϱ≥0,Zϱ=ϕ3ϱ≥0,ϱ∈-τ,0,ϕ10≥0,ϕ20≥0,ϕ30≥0,where ϕ1(ϱ),ϕ2(ϱ),ϕ3(ϱ)∈C([-τ,0],R+3), the Banach space of continuous functions mapping the interval [-τ,0] into R+3, where R+3={(x1,x2,x3):xi≥0,i=1,2,3}.The organisation of this paper is as follows: in Section2, preliminary dynamical results of system without delay are analyzed, such as the existence of equilibria, boundness, and the local stability of equilibria. In Section 3, with the presence of time delay, the local stability and the existence of Hopf bifurcation of system (2) are investigated. Based on normal form theory and center manifold theorem, direction of Hopf bifurcation and stability of the bifurcating periodic solutions are also discussed. What is more, the global existence of bifurcating periodic solution is given by using a global Hopf bifurcation result due to Wu [30]. In Section 4, the optimal control problem is formulated by Pontryagin’s Maximum Principle. In Section 5, numerical simulations are given to illustrate the results.


## 2. Preliminary Results of System without Delay

### 2.1. Positivity and Boundedness

Whenτ=0, system (2) is(4)N˙t=N0-αN-bNP+k1P,P˙t=b1NP-kP-βPZr+P2-q1EP,Z˙t=β1PZr+P2-μZ-θPZr+P2-q2EZ.The initial conditions for system (4) take the following form(5)N0=ϕ10≥0,P0=ϕ20≥0,Z0=ϕ30≥0.Lemma 1.
Ifr≥1, solutions of system (4) with initial conditions are defined on [0,+∞) and remain positive for all t≥0.Proof.
From the first equation of system (4), we obtain(6)N˙t≥-αN-bNP.Hence, we obtain(7)Nt≥ϕ10exp⁡∫0t-α-bPsds>0.From the second equation of system (4), we have(8)Pt≥ϕ20exp⁡∫0tb1Ns-k-βZsr-q1Eds>0,forr≥1.Since ϕ3(0)≥0, the third equation of system (4) shows that by a standard comparison argument(9)Zt≥ϕ30e-μ+q2Et>0.Theorem 2.
All solutions of system (4) with initial conditions are bounded for all t≥0.Proof.
In order to check the boundedness of solutions of system (4), we define the function (10)Wt=Nt+bb1Pt+βbβ1-θb1Zt.As derivative of W(t) with respect to system (4), we obtain(11)W˙t=N0-αN+k1P-bb1kP-bb1q1EP-βbβ1-θb1μ+q2EZ=N0-αN-bb1k+q1EP-βbβ1-θb1μ+q2EZ+k1P≤N0-αN-bb1k+q1E-b1kbP-βbβ1-θb1μ+q2EZ≤-mWt+N0,where m=min⁡{α,k+q1E-b1k/b,μ+q2E}. Therefore, by the comparison theory [31], we have(12)Wt≤W0-N0me-mt+N0m.Thus, N(t),P(t), and Z(t) are bounded.


### 2.2. Equilibria and Their Local Stability

In this subsection, we will discuss the existence of equilibria of system (4) and their local asymptotical stability.The given system (4) has three nonnegative equilibria, defined as the plankton extinction equilibrium E0N0/α,0,0, the zooplankton free equilibrium E1(N1,P1,0), and the coexistence equilibrium E∗(N∗,P∗,Z∗).The zooplankton free equilibrium isE1(N1,P1,0), where(13)N1=k+q1Eb1,P1=αk+q1E-b1N0k1b1-bk-bq1E.Since k1b1-bk-bq1E<0, if b1N0>α(k+q1E), then there is only one possible equilibrium involving phytoplankton but not zooplankton.The coexistence equilibrium isE∗(N∗,P∗,Z∗), where N∗, P∗, and Z∗ satisfy the following equations(14)N0-αN-bNP+k1P=0,b1N-k+q1E-βZr+P2=0,β1-θPr+P2-μ+q2E=0.From the third equation of (14), we have(15)μ+q2EP2-β1-θP+rμ+q2E=0.If β1-θ>0, and △=(β1-θ)2-4r(μ+q2E)2>0, then (15) has exactly two positive real roots(16)P1∗=β1-θ-△2μ+q2E,P2∗=β1-θ+△2μ+q2E.And if β1-θ>0, and △=0, then (15) has two equal positive real roots(17)P3∗=β1-θ2μ+q2E.Again from the first and second equations of (14), we have(18)N1∗=2N0μ+q2E+k1β1-θ-△2αμ+q2E+bβ1-θ-△,Z1∗=β1-θβ1-θ-△2βμ+q2E22b1N0μ+q2E+b1k1β1-θ-△2αμ+q2E+bβ1-θ-△-k-q1E,N2∗=2N0μ+q2E+k1β1-θ+△2αμ+q2E+bβ1-θ+△,Z2∗=β1-θβ1-θ+△2βμ+q2E22b1N0μ+q2E+b1k1β1-θ+△2αμ+q2E+bβ1-θ+△-k-q1E,N3∗=2N0μ+q2E+k1β1-θ2αμ+q2E+bβ1-θ,Z3∗=β1-θ22βμ+q2E22b1N0μ+q2E+b1k1β1-θ2αμ+q2E+bβ1-θ-k-q1E.If β1-θ>0, k<2b1N0μ+q2E+b1k1β1-θ-△/2αμ+q2E+bβ1-θ-△-q1E, then Z1∗>0. If β1-θ>0, k<2b1N0(μ+q2E)+b1k1(β1-θ+△)/2α(μ+q2E)+b(β1-θ+△)-q1E, then Z2∗>0. And if β1-θ>0, k<2b1N0(μ+q2E)+b1k1(β1-θ)/2α(μ+q2E)+b(β1-θ)-q1E, then Z3∗>0.In order to discuss the existence and the local stability of equilibria of system (4), some assumptions are given.( H 1 ): Δ>0, β1-θ>0, Z1∗>Z2∗(orZ1∗<Z2∗) and 2b1N0(μ+q2E)+b1k1(β1-θ+△)/2α(μ+q2E)+b(β1-θ+△) − q1E<k<2b1N0(μ+q2E)+b1k1(β1-θ-△)/2α(μ+q2E)+b(β1-θ-△)-q1E (or 2b1N0(μ+q2E)+b1k1(β1-θ-△)/2α(μ+q2E)+b(β1-θ+△) − q1E<k<2b1N0(μ+q2E)+b1k1(β1-θ+△)/2α(μ+q2E)+b(β1-θ+△)-q1E).( H 2 ): Δ>0, β1-θ>0 and k<min{2b1N0μ+q2E+b1k1β1-θ-△/2αμ+q2E+bβ1-θ-△ − q1E, 2b1N0μ+q2E+b1k1β1-θ+△/2αμ+q2E+bβ1-θ+△ − q1E}.( H 3 ): Δ=0, β1-θ>0 and k<2b1N0(μ+q2E)+b1k1(β1-θ)/2α(μ+q2E)+b(β1-θ)-q1E.( H 1 ′ ): min{q1E2αμ+q2E+bβ1-θ-Δ/2N0μ+q2E+k1β1-θ-Δ, q1E2αμ+q2E+bβ1-θ+Δ/2N0μ+q2E+k1β1-θ+Δ, bq1E/k1} < b1 < max{q1E2αμ+q2E+bβ1-θ-Δ/2N0μ+q2E+k1β1-θ-Δ, q1E2αμ+q2E+bβ1-θ+Δ/2N0μ+q2E+k1β1-θ+Δ, bq1E/k1} and r>(β1-θ)(β1-θ+Δ)/2(μ+q2E)2.( H 2 ′ ): b1<min{q1E2αμ+q2E+bβ1-θ-Δ/2N0μ+q2E+k1β1-θ-Δ, q1E2αμ+q2E+bβ1-θ+Δ/2N0μ+q2E+k1β1-θ+Δ, bq1E/k1} and r>(β1-θ)(β1-θ+Δ)/2(μ+q2E)2.( H 3 ′ ) : b 1 < min ⁡ q 1 E 2 α μ + q 2 E + b β 1 - θ / 2 N 0 μ + q 2 E + k 1 β 1 - θ , b q 1 E / k 1and  r>(β1-θ)2/2(μ+q2E)2.From the above analysis, we obtain the following results.Theorem 3.
(i) The plankton extinction equilibriumE0N0/α,0,0 always exists.
(ii) The zooplankton free equilibriumE1 exists if b1N0>α(k+q1E).
(iii) The coexistence equilibriumE1∗(N1∗,P1∗,Z1∗) (or E2∗(N2∗,P2∗,Z2∗)) exists if (H1) holds; the coexistence equilibria E1∗ and E2∗ exist if (H2) holds; the coexistence equilibrium E3∗(N3∗,P3∗,Z3∗) exists if (H3) holds.According to the results of Theorem3 and the above assumptions, we have the results as follows.Theorem 4.
(1) Ifb1N0<α(k+q1E), then the plankton extinction equilibrium E0 is locally asymptotically stable.
(2) Ifb1N0>α(k+q1E) and μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E, then the zooplankton free equilibrium E1 is locally asymptotically stable.
(3) The coexistence equilibriumE1∗ (or E2∗) is locally asymptotically stable if conditions (H1) and (H1′) hold. The coexistence equilibria E1∗ and E2∗ are locally asymptotically stable if conditions (H2) and (H2′) hold. The coexistence equilibrium E3∗ is locally asymptotically stable if conditions (H3) and (H3′) hold.Proof.
(1) The characteristic equation onE0(N0/α,0,0) is given(19)λ+αλ-b1N0α+k+q1Eλ+μ+Eq2=0.It is clear that (19) has negative roots λ1=-α, and λ2=-(μ+Eq2). So if b1N0<α(k+q1E), then(20)λ3=b1N0α-k-q1E<0.Thus, the plankton extinction equilibrium E0=(N0/α,0,0) is locally asymptotically stable. (2) The characteristic equation about E1(N1,P1,0) is given by(21)λ-β1-θP1r+P12+μ+q2Eλ2+α+bP1λ+b1P1bN1-k1=0.If μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E, then(22)λ1=β1-θP1r+P12-μ-q2E<0.Further, the other roots λ2 and λ3 satisfy that λ2+λ3=-(α+b1P1)<0. If λ2λ3=b1P1(-bN1+k1)>0, then λ2<0,λ3<0. Therefore, if b1N0>α(k+q1E) and μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E, the zooplankton free equilibrium E1=(N1,P1,0) is locally asymptotically stable.
(3) The characteristic equation aboutE1∗(N1∗,P1∗,Z1∗) is given(23)λ3+M1λ2+M2λ+M3=0,where(24)m1=α+bP1∗,m2=βZ1∗r-P1∗2r+P1∗22+k,m3=q1E-b1N1∗,m4=βP1∗r+P1∗2,m5=β1-θZ1∗r-P1∗2r+P1∗22,m6=bq1E-k1b1,M1=m1+m2+m3,M2=m1m2+m4m5+αm3+P1∗m6,M3=m1m4m5.Since M1M2-M3=m12m2+αm1m3+P1∗m1m6+(m2+m3)M2, if (H1) and (H1′) hold, then M1>0, M1M2-M3>0. Therefore all roots of (23) have negative real parts. Based on Routh-Hurwitz criterion we obtain that the coexistence equilibrium E1∗(N1∗,P1∗,Z1∗) (or E2∗(N2∗,P2∗,Z2∗)) is locally asymptotically stable.
Similarly, if(H2) and (H2′) hold, the coexistence equilibria E1∗(N1∗,P1∗,Z1∗) and E2∗(N2∗,P2∗,Z2∗) are locally asymptotically stable; if (H3) and (H3′) hold, the coexistence equilibrium E3∗(N3∗,P3∗,Z3∗) is locally asymptotically stable.Remark 5.
From biological point of view, phytoplankton and zooplankton cannot survive, and the state is stable if the assumptionb1N0<α(k+q1E) holds; phytoplankton survive but zooplankton cannot survive, and the state is stable if the assumptions b1N0>α(k+q1E) and μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E hold; both phytoplankton and zooplankton survive, and the state is stable if the assumptions (H1)-(H3′) hold.


## 2.1. Positivity and Boundedness

Whenτ=0, system (2) is(4)N˙t=N0-αN-bNP+k1P,P˙t=b1NP-kP-βPZr+P2-q1EP,Z˙t=β1PZr+P2-μZ-θPZr+P2-q2EZ.The initial conditions for system (4) take the following form(5)N0=ϕ10≥0,P0=ϕ20≥0,Z0=ϕ30≥0.Lemma 1.
Ifr≥1, solutions of system (4) with initial conditions are defined on [0,+∞) and remain positive for all t≥0.Proof.
From the first equation of system (4), we obtain(6)N˙t≥-αN-bNP.Hence, we obtain(7)Nt≥ϕ10exp⁡∫0t-α-bPsds>0.From the second equation of system (4), we have(8)Pt≥ϕ20exp⁡∫0tb1Ns-k-βZsr-q1Eds>0,forr≥1.Since ϕ3(0)≥0, the third equation of system (4) shows that by a standard comparison argument(9)Zt≥ϕ30e-μ+q2Et>0.Theorem 2.
All solutions of system (4) with initial conditions are bounded for all t≥0.Proof.
In order to check the boundedness of solutions of system (4), we define the function (10)Wt=Nt+bb1Pt+βbβ1-θb1Zt.As derivative of W(t) with respect to system (4), we obtain(11)W˙t=N0-αN+k1P-bb1kP-bb1q1EP-βbβ1-θb1μ+q2EZ=N0-αN-bb1k+q1EP-βbβ1-θb1μ+q2EZ+k1P≤N0-αN-bb1k+q1E-b1kbP-βbβ1-θb1μ+q2EZ≤-mWt+N0,where m=min⁡{α,k+q1E-b1k/b,μ+q2E}. Therefore, by the comparison theory [31], we have(12)Wt≤W0-N0me-mt+N0m.Thus, N(t),P(t), and Z(t) are bounded.


## 2.2. Equilibria and Their Local Stability

In this subsection, we will discuss the existence of equilibria of system (4) and their local asymptotical stability.The given system (4) has three nonnegative equilibria, defined as the plankton extinction equilibrium E0N0/α,0,0, the zooplankton free equilibrium E1(N1,P1,0), and the coexistence equilibrium E∗(N∗,P∗,Z∗).The zooplankton free equilibrium isE1(N1,P1,0), where(13)N1=k+q1Eb1,P1=αk+q1E-b1N0k1b1-bk-bq1E.Since k1b1-bk-bq1E<0, if b1N0>α(k+q1E), then there is only one possible equilibrium involving phytoplankton but not zooplankton.The coexistence equilibrium isE∗(N∗,P∗,Z∗), where N∗, P∗, and Z∗ satisfy the following equations(14)N0-αN-bNP+k1P=0,b1N-k+q1E-βZr+P2=0,β1-θPr+P2-μ+q2E=0.From the third equation of (14), we have(15)μ+q2EP2-β1-θP+rμ+q2E=0.If β1-θ>0, and △=(β1-θ)2-4r(μ+q2E)2>0, then (15) has exactly two positive real roots(16)P1∗=β1-θ-△2μ+q2E,P2∗=β1-θ+△2μ+q2E.And if β1-θ>0, and △=0, then (15) has two equal positive real roots(17)P3∗=β1-θ2μ+q2E.Again from the first and second equations of (14), we have(18)N1∗=2N0μ+q2E+k1β1-θ-△2αμ+q2E+bβ1-θ-△,Z1∗=β1-θβ1-θ-△2βμ+q2E22b1N0μ+q2E+b1k1β1-θ-△2αμ+q2E+bβ1-θ-△-k-q1E,N2∗=2N0μ+q2E+k1β1-θ+△2αμ+q2E+bβ1-θ+△,Z2∗=β1-θβ1-θ+△2βμ+q2E22b1N0μ+q2E+b1k1β1-θ+△2αμ+q2E+bβ1-θ+△-k-q1E,N3∗=2N0μ+q2E+k1β1-θ2αμ+q2E+bβ1-θ,Z3∗=β1-θ22βμ+q2E22b1N0μ+q2E+b1k1β1-θ2αμ+q2E+bβ1-θ-k-q1E.If β1-θ>0, k<2b1N0μ+q2E+b1k1β1-θ-△/2αμ+q2E+bβ1-θ-△-q1E, then Z1∗>0. If β1-θ>0, k<2b1N0(μ+q2E)+b1k1(β1-θ+△)/2α(μ+q2E)+b(β1-θ+△)-q1E, then Z2∗>0. And if β1-θ>0, k<2b1N0(μ+q2E)+b1k1(β1-θ)/2α(μ+q2E)+b(β1-θ)-q1E, then Z3∗>0.In order to discuss the existence and the local stability of equilibria of system (4), some assumptions are given.( H 1 ): Δ>0, β1-θ>0, Z1∗>Z2∗(orZ1∗<Z2∗) and 2b1N0(μ+q2E)+b1k1(β1-θ+△)/2α(μ+q2E)+b(β1-θ+△) − q1E<k<2b1N0(μ+q2E)+b1k1(β1-θ-△)/2α(μ+q2E)+b(β1-θ-△)-q1E (or 2b1N0(μ+q2E)+b1k1(β1-θ-△)/2α(μ+q2E)+b(β1-θ+△) − q1E<k<2b1N0(μ+q2E)+b1k1(β1-θ+△)/2α(μ+q2E)+b(β1-θ+△)-q1E).( H 2 ): Δ>0, β1-θ>0 and k<min{2b1N0μ+q2E+b1k1β1-θ-△/2αμ+q2E+bβ1-θ-△ − q1E, 2b1N0μ+q2E+b1k1β1-θ+△/2αμ+q2E+bβ1-θ+△ − q1E}.( H 3 ): Δ=0, β1-θ>0 and k<2b1N0(μ+q2E)+b1k1(β1-θ)/2α(μ+q2E)+b(β1-θ)-q1E.( H 1 ′ ): min{q1E2αμ+q2E+bβ1-θ-Δ/2N0μ+q2E+k1β1-θ-Δ, q1E2αμ+q2E+bβ1-θ+Δ/2N0μ+q2E+k1β1-θ+Δ, bq1E/k1} < b1 < max{q1E2αμ+q2E+bβ1-θ-Δ/2N0μ+q2E+k1β1-θ-Δ, q1E2αμ+q2E+bβ1-θ+Δ/2N0μ+q2E+k1β1-θ+Δ, bq1E/k1} and r>(β1-θ)(β1-θ+Δ)/2(μ+q2E)2.( H 2 ′ ): b1<min{q1E2αμ+q2E+bβ1-θ-Δ/2N0μ+q2E+k1β1-θ-Δ, q1E2αμ+q2E+bβ1-θ+Δ/2N0μ+q2E+k1β1-θ+Δ, bq1E/k1} and r>(β1-θ)(β1-θ+Δ)/2(μ+q2E)2.( H 3 ′ ) : b 1 < min ⁡ q 1 E 2 α μ + q 2 E + b β 1 - θ / 2 N 0 μ + q 2 E + k 1 β 1 - θ , b q 1 E / k 1and  r>(β1-θ)2/2(μ+q2E)2.From the above analysis, we obtain the following results.Theorem 3.
(i) The plankton extinction equilibriumE0N0/α,0,0 always exists.
(ii) The zooplankton free equilibriumE1 exists if b1N0>α(k+q1E).
(iii) The coexistence equilibriumE1∗(N1∗,P1∗,Z1∗) (or E2∗(N2∗,P2∗,Z2∗)) exists if (H1) holds; the coexistence equilibria E1∗ and E2∗ exist if (H2) holds; the coexistence equilibrium E3∗(N3∗,P3∗,Z3∗) exists if (H3) holds.According to the results of Theorem3 and the above assumptions, we have the results as follows.Theorem 4.
(1) Ifb1N0<α(k+q1E), then the plankton extinction equilibrium E0 is locally asymptotically stable.
(2) Ifb1N0>α(k+q1E) and μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E, then the zooplankton free equilibrium E1 is locally asymptotically stable.
(3) The coexistence equilibriumE1∗ (or E2∗) is locally asymptotically stable if conditions (H1) and (H1′) hold. The coexistence equilibria E1∗ and E2∗ are locally asymptotically stable if conditions (H2) and (H2′) hold. The coexistence equilibrium E3∗ is locally asymptotically stable if conditions (H3) and (H3′) hold.Proof.
(1) The characteristic equation onE0(N0/α,0,0) is given(19)λ+αλ-b1N0α+k+q1Eλ+μ+Eq2=0.It is clear that (19) has negative roots λ1=-α, and λ2=-(μ+Eq2). So if b1N0<α(k+q1E), then(20)λ3=b1N0α-k-q1E<0.Thus, the plankton extinction equilibrium E0=(N0/α,0,0) is locally asymptotically stable. (2) The characteristic equation about E1(N1,P1,0) is given by(21)λ-β1-θP1r+P12+μ+q2Eλ2+α+bP1λ+b1P1bN1-k1=0.If μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E, then(22)λ1=β1-θP1r+P12-μ-q2E<0.Further, the other roots λ2 and λ3 satisfy that λ2+λ3=-(α+b1P1)<0. If λ2λ3=b1P1(-bN1+k1)>0, then λ2<0,λ3<0. Therefore, if b1N0>α(k+q1E) and μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E, the zooplankton free equilibrium E1=(N1,P1,0) is locally asymptotically stable.
(3) The characteristic equation aboutE1∗(N1∗,P1∗,Z1∗) is given(23)λ3+M1λ2+M2λ+M3=0,where(24)m1=α+bP1∗,m2=βZ1∗r-P1∗2r+P1∗22+k,m3=q1E-b1N1∗,m4=βP1∗r+P1∗2,m5=β1-θZ1∗r-P1∗2r+P1∗22,m6=bq1E-k1b1,M1=m1+m2+m3,M2=m1m2+m4m5+αm3+P1∗m6,M3=m1m4m5.Since M1M2-M3=m12m2+αm1m3+P1∗m1m6+(m2+m3)M2, if (H1) and (H1′) hold, then M1>0, M1M2-M3>0. Therefore all roots of (23) have negative real parts. Based on Routh-Hurwitz criterion we obtain that the coexistence equilibrium E1∗(N1∗,P1∗,Z1∗) (or E2∗(N2∗,P2∗,Z2∗)) is locally asymptotically stable.
Similarly, if(H2) and (H2′) hold, the coexistence equilibria E1∗(N1∗,P1∗,Z1∗) and E2∗(N2∗,P2∗,Z2∗) are locally asymptotically stable; if (H3) and (H3′) hold, the coexistence equilibrium E3∗(N3∗,P3∗,Z3∗) is locally asymptotically stable.Remark 5.
From biological point of view, phytoplankton and zooplankton cannot survive, and the state is stable if the assumptionb1N0<α(k+q1E) holds; phytoplankton survive but zooplankton cannot survive, and the state is stable if the assumptions b1N0>α(k+q1E) and μ>(β1-θ)(αk+αq1E-b1N0)(k1b1-bk-bq1E)/r(k1b1-bk-bq1E)2+(αk+αq1E-b1N0)2-q1E hold; both phytoplankton and zooplankton survive, and the state is stable if the assumptions (H1)-(H3′) hold.


## 3. Hopf Bifurcation and Analysis of System with Delay

### 3.1. Hopf Bifurcation

In this part, we will only investigate the effect of time delay on the dynamics of system (2) about the coexistence equilibrium E1∗. For convenience, we denote the coexistence equilibrium E1∗= E∗. We first translate the coexistence equilibrium E∗ to the origin by setting N¯=N-N∗,P¯=P-P∗, and Z¯=Z-Z∗ and drop the bar for convenience of notation. Then, the linearized system of system (2) can be obtained(25)N˙t=a11N+a12P,P˙t=a21N+a22P+a23Z,Z˙t=a31Pt-τ+a32P,where(26)a11=-α+bP∗,a12=-bN∗+k1,a21=b1P∗,a22=b1N∗-βZ∗r-P∗2r+P∗22-k-q1E,a23=-βP∗r+P∗2,a31=β1Z∗r-P∗2r+P∗22,a32=-θZ∗r-P∗2r+P∗22.Obviously, the characteristic equation of system (25) around E∗ takes the following form(27)λ3+r1λ2+r2λ+r5+r3+r4λe-λτ=0,where(28)r1=-a11-a22,r2=a11a22-a23a32-a12a21,r3=a11a23a31,r4=-a23a31,r5=a11a23a32.Forτ≠0, according to Corollary 2.4 [32], the sum of orders of the zeros of (27) in the open right half-plane can change only when a zero appears on or crosses the imaginary axis. It is obvious that λ=0 is not a root of (27) when β1≠θ, and then the imaginary axis cannot be crossed by λ=0 for any τ>0. And due to Corollary 2.4 [32], a stability switch (or a cross of the imaginary axis) occurs with λ=±iω for some real ω>0.Now, we assume thatλ=iω(ω>0) is a root of (27). Separating real and imaginary parts, we have(29)r3cos⁡ωτ+r4ωsin⁡ωτ=r1ω2-r5,r4ωcos⁡ωτ-r3sin⁡ωτ=ω3-r2ω.Adding up the squares of both equations, we obtain(30)ω6+r12-2r2ω4+r22-2r1r5-r42ω2+r52-r32=0.By denoting z=ω2, (30) becomes(31)z3+p1z2+p2z+p3=0,where p1=r12-2r2,p2=r22-2r1r5-r42,p3=r52-r32.Leth(z)=z3+p1z2+p2z+p3. According to the distribution of the roots of (31) and the results of [33], we can easily get the following lemmas.Lemma 6.
(i) Ifp3<0, (31) has at least one positive root.
(ii) Ifp3≥0 and p12-3p2≤0, (31) does not have any positive root.
(iii) Ifp3≥0 and q12-3q2>0, (31) has positive roots if and only if z1∗=1/3(-p1+p12-3p2)>0 and h(z1∗)≤0.Thus, if the condition(H1) holds, then h(0)=p3<0. Equation (31) has at least one positive root. Without loss of generality, we assume that it has three real positive roots, which are defined by z1,z2,z3, respectively. For convenience, we assume z1<z2<z3. Then (30) has three positive roots ωk=zk,k=1,2,3. From (29), we have(32)τkj=1ωkarccos⁡r4ωk4+r1r3-r2r4ωk2-r3r5r42ωk2+r32+2jπ,k=1,2,3;j=0,1,2,….Then±iωk is a pair of purely imaginary roots of (27) with τkj. Define(33)τ0=τk0=mink∈1,2,3⁡τk0,ω0=ωk0.And we have the following transversality condition.Lemma 7.
Letλ(τ)=α(τ)+iω(τ) be the root of (27) satisfying λ(τkj)=0, ω(τkj)=ωk, k=1,2,3, j=0,1,2,…. Then(34)sign⁡Re⁡dλτdτλ=iωk=sign⁡h′zk=>0,k=1,3,<0,k=2.Proof.
Differentiating both sides of (27) with respect to τ, we obtain(35)dλdτ=λe-λτr4λ+r33λ2+2r1λ+r2-τr4λ+r3e-λτ+r4e-λτ.Then(36)dλdτ-1=3λ2+2r1λ+r2eλτ+r4λr4λ+r3-τλ.From (30), we have(37)Redλdτ-1τ=τkj=Redλdτ-1λ=iωk=Re-3ω2+2r1iωk+r2eiωkτkj+r4r4iωk+r3iωk-τkjiωk=ωk2r42ωk4+r32ωk23ωk4+2r12-4r2ωk2+r22-2r1r5-r42=1r42ωk2+r32h′zk.
Thus we have(38)sign⁡Redλdτλ=iωk=sign⁡Redλdτ-1λ=iωk=sign⁡1r4ωk2+r32h′zk.Since h(0)<0, h(∞)=∞, we obtain h′(z1)>0, h′(z2)<0, h′(z3)>0. This ends the proof.To summarize the above discussions, we have the following result.Theorem 8.
Suppose that(H1) and (H1′) hold. System (2) at E∗ is locally asymptotically stable whenever τ∈[0,τ0] and unstable when τ>τ0. System (2) undergoes a Hopf bifurcation for τ=τ0.


### 3.2. Properties of Hopf Bifurcation

From Theorem8, we have obtained the conditions for the existence of Hopf bifurcation when τ=τ0. In this subsection, we shall consider the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions of system (2) by using the center manifold and normal form theory presented by Hassard [34].Letx(t)=N(τt)-N∗,y(t)=P(τt)-P∗,z(t)=Z(τt)-Z∗, τ=τ0+σ,ϕ(t)=(x(t),y(t),z(t))T, and ϕt(θ~)=ϕ(t+θ~) for θ~∈[-1,0]. Using Taylor series expansion about E∗=(N∗,P∗,Z∗), system (2) can be rewritten as(39)ϕ˙t=τG1ϕt+G2ϕt-1+f,where(40)G1=a11a120a21a22a230a320,G2=0000000a310,f(σ,ϕ)=(f1(σ,ϕ),f2(σ,ϕ),f3(μ,ϕ))T and(41)f1σ,ϕ=a110xtyt,f2σ,ϕ=b110xtyt+b011ytzt+b020y2t,f3σ,ϕ=c110ytzt+c011ztyt-1+c200y2t+c002y2t-1,where(42)a110=-b,b110=b1,b011=-βr-P∗2r+P∗22,b020=6βZ∗P∗r-P∗2r+P∗23,c110=-θr-P∗2r+P∗22,c011=β1r-P∗2r+P∗22,c200=6θZ∗P∗r-P∗2r+P∗23,c002=-6β1Z∗P∗r-P∗2r+P∗23.Linearization of the delayed system (41) around the origin is given(43)ϕ˙t=τG1ϕt+G2ϕt-1.Define(44)Lσϕ=τ0+σG1ϕt+G2ϕt-1,where ϕ∈C1([-1,0],R3).By the Riesz representation theorem [35], there exists 3×3 matrix-valued function ζ(θ~,σ) of bounded variation for θ~∈[-1,0] such that(45)Lσϕ=∫-10dζθ~,σϕθ~,forϕ∈C-1,0,R3.In fact, we can choose(46)ζθ~,σ=τ0+σG1δθ~+G2δθ~+1,where δ(θ~) is a Dirac delta function.Next, forϕ∈C1([-1,0],R3), we define(47)Aσϕθ~=dϕθ~dθ~,θ~∈-1,0,∫-10dζθ~,σϕθ~,θ~=0,and(48)Rσϕθ~=0,θ~∈-1,0,fσϕ,θ~=0,and then system (43) can be rewritten as(49)u˙t=Aσut+Rσut,where u=(x,y,z)T and ut=u(t+θ~) for θ~∈[-1,0].Forψ∈C1([-1,0],R3), the adjoint operator A∗ of A is defined as(50)A∗σϕs=-dψsds,s∈0,1,∫-10ψ-ξdζξ,0,s=0,and a bilinear inner product(51)ψs,ϕθ~=ψ¯0ϕ0-∫θ~=-10∫ξ=0θ~ψ¯ξ-θ~dζθ~ϕξdξ,where ζ(θ~)=ζ(θ~,0). From the above discussion, we can know that ±iω0 are the eigenvalues of A(0), and thus ±iω0 are also the eigenvalues of A∗(0).Letq(θ~)=(1,q2,q3)Teiω0τ0θ~ be the eigenvector of A(0) corresponding to the eigenvalue iω0τ0, and q∗(θ~)=D1(1,q2∗,q3∗)Te-iω0τ0θ~ be the eigenvector of A∗(0) corresponding to the eigenvalue -iω0τ0. With the conditions A(0)q(0)=iω0q(0) and A∗(0)q∗(0)=-iω0q∗(0), we obtain(52)iω0I-G1+G2e-iω0τ0q0=0,and(53)-iω0I-G1T+G2Teiω0τ0q∗0=0.where I is a third order identity matrix. Hence, we can obtain q2=iω0-a11/a12,q3=(iω0-a22)q2-a21/a23, q2∗=-iω0+a11/a21,q3∗=-a23q3/iω0.In order to satisfy thatq∗,q=1, we need to choose the value of D1 in the following part. For convenience, we denote Γ=(1,q2,q3)T,Γ∗=(1,q2∗,q3∗)T. By virtue of (52) and (53), it derives that(54)q∗,q=q∗¯0q0-∫θ~=-10∫ξ=0θ~q¯∗ξ-θ~dζθ~qξdξ,=D1¯Γ∗¯ΓT-D1¯∫θ~=-10∫ξ=0θ~Γ∗¯e-iω0τ0ξ-θ~dζθ~ΓTeiω0τ0ξdξ=D1¯Γ∗¯ΓT-D1¯∫-10Γ∗¯θ~eiω0τ0θ~dζθ~ΓT=D1¯Γ∗¯ΓT+D1¯τ0e-iω0τ0Γ∗¯G2ΓT=D1¯1+q2q2∗¯+q3q3∗¯+τ0q3∗¯q2a31e-iω0τ0.Hence, we can choose D1¯=[1+q2q2∗¯+q3q3∗¯+τ0q3∗¯q2a31e-iω0τ0]-1.Nextly, we will compute the coordinates to detail center manifoldC0 at σ≠0. Let μt be solution of (49) when σ=0. We define(55)zt=q∗,ut,Wt,θ~=ut-zq-z¯q¯=ut-2Re⁡ztqθ~.On the center manifoldC0, it derives that(56)Wt,θ~=Wzt,z¯t,θ~=W20θ~z22+W11θ~zz¯+W02θ~z¯22+⋯,where z and z¯ are local coordinates for C0 in the direction of q∗ and q∗¯, respectively.We only consider real solutions asW is real. If μt is real, then we get(57)q∗,W=q∗,ut-zq-z¯q¯=q∗,ut-q∗,q-q∗,q¯.For the solution μt of (49), since σ=0, from (49) and (51), we have(58)z˙t=q∗,u˙t=q∗,A0ut+R0ut=A∗0q∗,ut+q∗¯0fut,0≔iω0τ0zt+gz,z¯,where(59)gz,z¯=q∗¯0f0Wz,z¯,0+2Re⁡ztqθ~=g20z22+g11zz¯+g02z¯22+g21z2z¯2+···.It follows from (55) and (56) that(60)utθ~=W20θ~z22+W11θ~zz¯+W02θ~z¯22+ΓTeiω0τ0θ~z+Γ∗Te-iω0τ0θ~z¯+···.From (41) and (59), it derives that(61)gz,z¯=τ0D1¯1,q2∗¯,q3∗¯a110x0y0b110x0y0+b011y0z0+b020y20c110y0z0+c011zty-1+c200y20+c002y2-1=τ0D1¯A1W201z22+W111zz¯+W021z¯22+z+z¯W202z22+W112zz¯+W022z¯22+q2z+q2¯z¯+A2W203z22+W113zz¯+W023z¯22+q3z+q3¯z¯W202z22+W112zz¯+W022z¯22+q2z+q2¯z¯+A3W202z22+W112zz¯+W022z¯22+q2z+q2¯z¯2+A4W2030z22+W1130zz¯+W0230z¯22+q3z+q3¯z¯W202-1z22+W112-1zz¯+W022-1z¯22+q2ze-iω0τ0+q2¯z¯eiω0τ0+A5W202-1z22+W112-1zz¯+W022-1z¯22+q2ze-iω0τ0+q2¯z¯eiω0τ02.Here A1=a110+b110q2∗¯, A2=b011q2∗¯+c110q3∗¯, A3=b020q2∗¯+c200q3∗¯, A4=c011q3∗¯, A5=c002q3∗¯.According to comparing its coefficients with (59), it derives that(62)g20=2τ0D1¯A1q2+A2q2q3+A3q22+A4q2q3e-iω0τ0+A5q22e-2iω0τ0,g11=τ0D1¯A1q2+A2q2q3¯+A3q2q2¯+A4q2q3¯e-iω0τ0+A5q1q2¯,g02=2τ0D1¯A1q2¯+A2q2¯q3¯+A3q2¯2+A4q2¯q3¯eiω0τ0+A5q22e2iω0τ0,g21=τ0D1¯A1W201q2¯+2q2W111+2W112+W202+τ0D1¯A2W202q3¯+2q3W112+2q2W113+q2¯W203+τ0D1¯A3W202q2¯+2q2W112+2q2W112+W202q2¯+τ0D1¯A4W203q2¯eiω0τ0+2q2W113e-iω0τ0+2q3W112-1+W202-1q3¯+τ0D1¯A52W202-1q2¯eiω0τ0+4q2W112-1e-iω0τ0.Since W11 and W20 are still unknown, we still compute them next. By virtue of (49) and (55), we have(63)W˙=u˙t-z˙q-z¯˙q¯=A0W-2Re⁡q¯∗0f0qθ~,θ~∈-1,0A0W-2Re⁡q¯∗0f0q0+f0z,z¯,θ~=0≜A0W0+Hz,z¯,θ~.Here,(64)Hz,z¯,θ~=H20θ~z22+H11θ~zz¯+H02θ~z¯22+···.By computing and comparing coefficients, we can obtain(65)A0-2iω0τ0W20θ~=-H20θ~,A0W11θ~=-H11θ~.From (63), we know that(66)Hz,z¯,θ~=-q∗0f0qθ~-q∗0f0¯q¯θ~=-gz,z¯qθ~-g¯z,z¯q¯θ~.Using (59) in (66) and comparing coefficient with (64), it derives that(67)H20θ~=-g20qθ~-g¯02q¯θ~,(68)H11θ~=-g11qθ~-g¯11q¯θ~.From the definition of A(0), based on (65) and (67), we can obtain (69)W˙20θ~=2iω0τ0W20θ~+g20qθ~+g¯02q¯θ~.Because of q(θ~)=(1,q2,q3)Teiω0τ0θ, we have(70)W20θ~=ig20q0ω0τ0eiω0τ0θ~+ig¯20q¯03ω0τ0e-iω0τ0θ~+E1e2iω0τ0θ~.Similarly, it derives that from (65) and (66)(71)W11θ~=-ig11q0ω0eiω0θ~+ig¯11q¯0ω0e-iω0θ~+E2.Here E1=(E1(1),E1(2),E1(3))T, E2=(E2(1),E2(2),E2(3))T are three dimensional constant vectors and can be determined by setting θ~=0 in H(z,z¯,θ~).By above definitions and conditions, we can compute(72)E11=2Υ11Υ1,E12=2Υ12Υ1,E13=2Υ13Υ1,E21=Υ21Υ2,E22=Υ22Υ2,E23=Υ23Υ2,where(73)Υ1=2iω0-a110-c110+c001b001-ω02-2iω0a010b100,Υ11=a1012iω0ς12-b001ς13+ς11-4ω02-2ω0b010-b001c100+c001,Υ12=2iω0-a1002iω0ς12+b001ς13+2iω0b100ς11,Υ13=ς132iω0-a1002iω0-b010-a010b100+c100+c001ς122iω0-a100+b100ς11,Υ2=-a100b001c100+c001,Υ21=-ς21b001c100+c001+ς23a010b001,Υ22=-ς23a100b001,Υ23=ς23a100b010-a010b100+c100+c001-ς22a100+-ς21b100,where(74)ς11=a110q2,ς12=b110q2+b011q2q3+b020q22,ς13=2c110q2q3+c011q2q3eiω0τ0+c200q22+c002q22e-2iω0τ0ς21=2a110q2,ς22=2b110q2+2b011q2q3¯+2b020q2q2¯,ς23=2c100q2q3¯+2c011q2q3¯e-2iω0τ0+2c200q2q2¯+2c002q2q2¯.Based on above analysis, the following values can be given(75)C10=i2ω0τ0g20g11-2g112-g0223+g212,μ2=-Re⁡C0Re⁡λ′τ0,β2=2Re⁡C0,T2=-Im⁡C0+μ2Im⁡λ′τ0ω0τ0.By computing (75), we have the following theorem.Theorem 9.
Properties of bifurcating periodic solution in the center manifold at critical valueτ0 are determined as follows:
(i) Ifμ2>0(μ2<0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τ0(τ<τ0).
(ii) Bifurcating periodic solutions are stable (unstable) ifβ2<0(β2>0).
(iii) Period of the periodic solution increases (decreases) ifT2>0(T2<0).


### 3.3. Global Hopf Bifurcation

In the above parts, we have obtained that system (2) undergoes Hopf bifurcation at E∗(N∗,P∗,Z∗) when τ=τkj. We will study the global continuation of these local bifurcating periodic solutions by using the global Hopf bifurcation theorem for delay differential equations [30], and state that system (2) admits periodic solutions globally for delay τ in a finite interval [τkj,+∞],k=1,2,3;j=1,2,…. This includes the values that are not necessarily near the local Hopf bifurcation values.For convenience, we give some notations,R+3={(N,P,Z):N>0,P>0,Z>0}, X=C([-τ,0],R+3), and x(t)=(N(t),P(t),Z(t))T. We can rewrite system (2) as the following function(76)dxtθ^dt=Fxtθ^,τ,T,where xt(θ^)=(N(t+θ^),P(t+θ^),Z(t+θ^))∈X.According to the work of Wu [30], we need to introduce the following notions(77)F^=FR+3×R+×R+:R+3×R+×R+→R+3,ΣF=clx,τ,T∈X×R+×R+:xis  a  T-periodic  solution  of76,NF=x¯,τ,T:Fx¯,τ,T=0.Let C(E∗,τkj,2π/ωk) be the connected component of C(E∗,τkj,2π/ωk) in Σ, and denote Projτ(E∗,τkj,2π/ωk) as its projection on τ component. ±iωk is a pair purely imaginary roots of (27) with τ=τkj. By Theorem 8, we know that C(E∗,τkj,2π/ωk) is nonempty.Lemma 10.
Whenever the equilibriumE∗ exists, all the nontrivial solutions of system (2) are uniformly bounded.Proof.
We assume that(N(t),P(t),Z(t)) is an arbitrary positive periodic solution of system (2) subject to initial condition (3). Let t1 and t2 be real quantities defined by N(t1)=min{N(t)}, N(t2)=max{N(t)}, P(t1)=min{P(t)}, P(t2)=max{P(t)}, Z(t1)=min{Z(t)}, and Z(t2)=max{Z(t)}. From the first and second equations of system (2), we have(78)0=N0-αNt2-bNt2Pt2+k1Pt2<N0-bNt2Pt2+k1Pt20=b1Nt2Pt2-kPt2-βPt2Zt2r+P2t2-q1EPt2<b1Nt2Pt2-kPt2-q1EPt2,which clearly implies that(79)Nt2Pt2<-N0k+q1Ek1b1-kb-q1bE≡θ1.Since k1b1<kb, θ1>0.
Similarly, from the second equation of system (2), we have(80)0=b1Nt2Pt2-kPt2-βPt2Zt2r+P2t2-q1EPt2<b1Nt2Pt2-βPt2Zt2r+P2t2.Thus, we obtain(81)βPt2Zt2r+Pt22<b1Nt2Pt2<b1θ1.
We defined the functionW1(t)=N(t)+b/b1P(t)+bβ/b1β1Z(t), for all t>0. Then it follows from system (2) that(82)W˙1=N0-αN+k1p-bb1kP-bb1βPZr+P2-bb1q1EP+bβb1β1β1Pt-τZr+P2t-τ-bβb1β1θPZr+P2-bβb1β1μ+q2EZ≤N0-αN-bb1k+q1EP+bβb1β1β1Pt-τZr+P2t-τ-βbb1β1μ+q2EZ+k1P≤N0-αN-bb1k+q1EP-bβb1β1μ+q2EZ+bθ1≤-m~W1t+N0+bθ1,where m~=min⁡{α,k+q1E,μ+q2E}. Therefore, by comparison with theory [31], we have(83)W1t≤W10-N0+bθ1m~e-m~t+N0+bθ1m~.Thus, all the nontrivial solutions of system (2) are uniformly bounded for t≥0.Lemma 11.
System (2) has no nonconstant τ-periodic solution if the condition (H1) and(84)H4:sup⁡a11+a22+a23,a31+a32+a11+a12,a21+a22<0are satisfied.Proof.
In order to explain that system (2) has no nonconstant periodic solution of periodic τ, we only prove that system (4) has no nonconstant periodic solution (see [36], Lemma 5.5). We rewrite system (4) as u˙(t)=f(u(t)). Suppose that u(t) is nonconstant periodic solution for system (4); then u(t)=(u1,u2,u3) is nonconstant periodic solution of the following equations(85)u˙1=N0-αu1-bu1u2+k1u2,u˙2=b1u1u2-ku2-βu2u3r+u22-q1Eu2,u˙3=β1u2u3r+u22-μu3-θu2u3r+u22-q2Eu3.From the Lemma 10, we know that the periodic orbits of system (85) corresponding to u(t) are included in the following area(86)G=u∈R3:m¯<uj<M¯,j=1,2,3,where m¯ and M¯ are the boundedness of the periodic solutions of system (2) in the G.
The Jacobian matrix of system (85) is(87)Ju=a11a120a21a22a230a31+a320.By using the method of [37], its second additive compound matrix is(88)J2u=a11+a22a230a31+a32a11a120a21a22.The second compound system is(89)x˙1=a11+a22x1+a23x2,x˙2=a31+a32x1+a11x2+a12x3,x˙3=a21x2+a22x3.We select vector module for R3 as follows(90)x1,x2,x3=max⁡x1,x2,x3.Then, the Lozinski uˇ measure μ(J[2](u)) of second additive compound matrix corresponding to a vector module above is(91)μJ2u=max⁡a11+a22+a23,a31+a32+a11+a12,a21+a22.Thus, if (H4) holds, there exists an ε such that μ(t)≤-ε<0 for all u∈G.This shows that the conclusion of Lemma11 follows from Corollary 3.5 in [38].Theorem 12.
Assume(H1) and (H4) are fulfilled; then there exists τ∈I, such that system (2) has at least one positive periodic solution when τ varies among τ1j, τ2j, and τ3j, j≥1, where(92)I=τkj>0:τkjare  the  Hopf  bifurcation  points  of  system2.Proof.
The characteristic matrix of system (2) at an equilibrium E¯=(N¯,P¯,Z¯) is given(93)ΔE¯,τ,Tλ=a¯11-λa¯120a¯21a¯22-λa¯230a¯31e-λτ+a¯32-λ,where(94)a¯11=-α+bP¯,a¯12=-bN¯+k1,a¯21=b1P¯,a¯22=b1N¯-βZ¯r-P¯2r+P¯22-k-q1E,a¯23=-βP¯r+P¯2,a¯31=β1Z¯r-P¯2r+P¯22,a¯32=-θZ¯r-P¯2r+P¯22.(E¯,τ¯,T¯) is called a center if F^(E¯,τ¯,T¯)=0 and det⁡ΔE¯,τ¯,T¯t2πi/T=0. A center is said to be isolated if it is the only center in some neighbourhood of (E¯,τ¯,T¯). It is obvious that (76) has an equilibrium x∗=E∗ in R+3 under the condition (H1). It follows that E∗,τkj,2π/ωk is an isolated center from the discussion about the local Hopf bifurcation theorem in Section 3.1. By Lemma 7 and Theorem 8, there exist constants ϵ>0 and δ>0, as well as a smooth curve λ:τkj-δ,τkj+δ→C such that(95)det⁡ΔE¯,τ¯,T¯λτ=0,λ-iωk<ϵwhen τ∈[τkj-δ,τkj+δ] holds, and(96)λτkj=iωk,dRe⁡λdττ=τkj≠0.Denotes Tk=2π/ωk, and(97)Ωϵ,T=u,T:0<u<ϵ,T-Tk<ϵ.Apparently, if τ-τkj≤δ, (u,T)∈∂Ωϵ,T such that det⁡ΔE∗,τkj,Tλu+i2π/T=0 if and only if u=0,τ=τkj. Hence all the assumptions (A1)-(A4) of Theorem 3.2 in [30] are satisfied for m=1. Moreover, if we define(98)H±E∗,τkj,2πωku,T=det⁡ΔE∗,τkj+δ,Tu+i2πT,then we have the crossing number(99)γ1E∗,τkj,2πωk=degBH-E∗,τkj,2πωk,Ωϵ,2π/ωk-degBH+E∗,τkj,2πωk,Ωϵ,2π/ωk=-1,dRe⁡λτ1,3jdτ1,3j>0,1,dRe⁡λτ2jdτ2j<0.By Theorem 3.2 of Wu [30], we conclude that the connected component C(E∗,τkj,2π/ωk) through (E∗,τkj,2π/ωk) in ∑‍ is nonempty. Meanwhile, we have that
(i)C(E∗,τkj,2π/ωk) is unbounded or
(ii)C(E∗,τkj,2π/ωk) is bounded, CE∗,τkj,2π/ωk⋂N(F) is finite, and ∑(E¯,τ,T)∈C(E∗,τkj,2π/ωk)⋂N(F)‍γ1(E¯,τ,T)=0.
By expression (32), we all know that there exists an integer n>0, such that 2π<τkjωk<2(n+1)π, τkj∈J, which implies that τkj/n+1<2π/ωk<τkj.
Therefore, we have thatτ/n+1<T<τ if (E¯,τ,T)∈C(E∗,τkj,2π/ωk). This fact and Lemma 11 show that the projection of C(E∗,τkj,2π/ωk) onto the T-space is bounded. On the flip side, Lemma 10 implies that the projection of C(E∗,τkj,2π/ωk) onto the x-space is uniformly bounded for τ∈Projτ(E∗,τkj,2π/ωk). Otherwise, C(E∗,τkj,2π/ωk) is bounded and Projτ(E∗,τkj,2π/ωk)⊂[τk,∞], which contradicts the fact that nonconstant periodic solution and T are both uniformly bounded for τ∈Projτ(E∗,τkj,2π/ωk). The proof is completed.


## 3.1. Hopf Bifurcation

In this part, we will only investigate the effect of time delay on the dynamics of system (2) about the coexistence equilibrium E1∗. For convenience, we denote the coexistence equilibrium E1∗= E∗. We first translate the coexistence equilibrium E∗ to the origin by setting N¯=N-N∗,P¯=P-P∗, and Z¯=Z-Z∗ and drop the bar for convenience of notation. Then, the linearized system of system (2) can be obtained(25)N˙t=a11N+a12P,P˙t=a21N+a22P+a23Z,Z˙t=a31Pt-τ+a32P,where(26)a11=-α+bP∗,a12=-bN∗+k1,a21=b1P∗,a22=b1N∗-βZ∗r-P∗2r+P∗22-k-q1E,a23=-βP∗r+P∗2,a31=β1Z∗r-P∗2r+P∗22,a32=-θZ∗r-P∗2r+P∗22.Obviously, the characteristic equation of system (25) around E∗ takes the following form(27)λ3+r1λ2+r2λ+r5+r3+r4λe-λτ=0,where(28)r1=-a11-a22,r2=a11a22-a23a32-a12a21,r3=a11a23a31,r4=-a23a31,r5=a11a23a32.Forτ≠0, according to Corollary 2.4 [32], the sum of orders of the zeros of (27) in the open right half-plane can change only when a zero appears on or crosses the imaginary axis. It is obvious that λ=0 is not a root of (27) when β1≠θ, and then the imaginary axis cannot be crossed by λ=0 for any τ>0. And due to Corollary 2.4 [32], a stability switch (or a cross of the imaginary axis) occurs with λ=±iω for some real ω>0.Now, we assume thatλ=iω(ω>0) is a root of (27). Separating real and imaginary parts, we have(29)r3cos⁡ωτ+r4ωsin⁡ωτ=r1ω2-r5,r4ωcos⁡ωτ-r3sin⁡ωτ=ω3-r2ω.Adding up the squares of both equations, we obtain(30)ω6+r12-2r2ω4+r22-2r1r5-r42ω2+r52-r32=0.By denoting z=ω2, (30) becomes(31)z3+p1z2+p2z+p3=0,where p1=r12-2r2,p2=r22-2r1r5-r42,p3=r52-r32.Leth(z)=z3+p1z2+p2z+p3. According to the distribution of the roots of (31) and the results of [33], we can easily get the following lemmas.Lemma 6.
(i) Ifp3<0, (31) has at least one positive root.
(ii) Ifp3≥0 and p12-3p2≤0, (31) does not have any positive root.
(iii) Ifp3≥0 and q12-3q2>0, (31) has positive roots if and only if z1∗=1/3(-p1+p12-3p2)>0 and h(z1∗)≤0.Thus, if the condition(H1) holds, then h(0)=p3<0. Equation (31) has at least one positive root. Without loss of generality, we assume that it has three real positive roots, which are defined by z1,z2,z3, respectively. For convenience, we assume z1<z2<z3. Then (30) has three positive roots ωk=zk,k=1,2,3. From (29), we have(32)τkj=1ωkarccos⁡r4ωk4+r1r3-r2r4ωk2-r3r5r42ωk2+r32+2jπ,k=1,2,3;j=0,1,2,….Then±iωk is a pair of purely imaginary roots of (27) with τkj. Define(33)τ0=τk0=mink∈1,2,3⁡τk0,ω0=ωk0.And we have the following transversality condition.Lemma 7.
Letλ(τ)=α(τ)+iω(τ) be the root of (27) satisfying λ(τkj)=0, ω(τkj)=ωk, k=1,2,3, j=0,1,2,…. Then(34)sign⁡Re⁡dλτdτλ=iωk=sign⁡h′zk=>0,k=1,3,<0,k=2.Proof.
Differentiating both sides of (27) with respect to τ, we obtain(35)dλdτ=λe-λτr4λ+r33λ2+2r1λ+r2-τr4λ+r3e-λτ+r4e-λτ.Then(36)dλdτ-1=3λ2+2r1λ+r2eλτ+r4λr4λ+r3-τλ.From (30), we have(37)Redλdτ-1τ=τkj=Redλdτ-1λ=iωk=Re-3ω2+2r1iωk+r2eiωkτkj+r4r4iωk+r3iωk-τkjiωk=ωk2r42ωk4+r32ωk23ωk4+2r12-4r2ωk2+r22-2r1r5-r42=1r42ωk2+r32h′zk.
Thus we have(38)sign⁡Redλdτλ=iωk=sign⁡Redλdτ-1λ=iωk=sign⁡1r4ωk2+r32h′zk.Since h(0)<0, h(∞)=∞, we obtain h′(z1)>0, h′(z2)<0, h′(z3)>0. This ends the proof.To summarize the above discussions, we have the following result.Theorem 8.
Suppose that(H1) and (H1′) hold. System (2) at E∗ is locally asymptotically stable whenever τ∈[0,τ0] and unstable when τ>τ0. System (2) undergoes a Hopf bifurcation for τ=τ0.


## 3.2. Properties of Hopf Bifurcation

From Theorem8, we have obtained the conditions for the existence of Hopf bifurcation when τ=τ0. In this subsection, we shall consider the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions of system (2) by using the center manifold and normal form theory presented by Hassard [34].Letx(t)=N(τt)-N∗,y(t)=P(τt)-P∗,z(t)=Z(τt)-Z∗, τ=τ0+σ,ϕ(t)=(x(t),y(t),z(t))T, and ϕt(θ~)=ϕ(t+θ~) for θ~∈[-1,0]. Using Taylor series expansion about E∗=(N∗,P∗,Z∗), system (2) can be rewritten as(39)ϕ˙t=τG1ϕt+G2ϕt-1+f,where(40)G1=a11a120a21a22a230a320,G2=0000000a310,f(σ,ϕ)=(f1(σ,ϕ),f2(σ,ϕ),f3(μ,ϕ))T and(41)f1σ,ϕ=a110xtyt,f2σ,ϕ=b110xtyt+b011ytzt+b020y2t,f3σ,ϕ=c110ytzt+c011ztyt-1+c200y2t+c002y2t-1,where(42)a110=-b,b110=b1,b011=-βr-P∗2r+P∗22,b020=6βZ∗P∗r-P∗2r+P∗23,c110=-θr-P∗2r+P∗22,c011=β1r-P∗2r+P∗22,c200=6θZ∗P∗r-P∗2r+P∗23,c002=-6β1Z∗P∗r-P∗2r+P∗23.Linearization of the delayed system (41) around the origin is given(43)ϕ˙t=τG1ϕt+G2ϕt-1.Define(44)Lσϕ=τ0+σG1ϕt+G2ϕt-1,where ϕ∈C1([-1,0],R3).By the Riesz representation theorem [35], there exists 3×3 matrix-valued function ζ(θ~,σ) of bounded variation for θ~∈[-1,0] such that(45)Lσϕ=∫-10dζθ~,σϕθ~,forϕ∈C-1,0,R3.In fact, we can choose(46)ζθ~,σ=τ0+σG1δθ~+G2δθ~+1,where δ(θ~) is a Dirac delta function.Next, forϕ∈C1([-1,0],R3), we define(47)Aσϕθ~=dϕθ~dθ~,θ~∈-1,0,∫-10dζθ~,σϕθ~,θ~=0,and(48)Rσϕθ~=0,θ~∈-1,0,fσϕ,θ~=0,and then system (43) can be rewritten as(49)u˙t=Aσut+Rσut,where u=(x,y,z)T and ut=u(t+θ~) for θ~∈[-1,0].Forψ∈C1([-1,0],R3), the adjoint operator A∗ of A is defined as(50)A∗σϕs=-dψsds,s∈0,1,∫-10ψ-ξdζξ,0,s=0,and a bilinear inner product(51)ψs,ϕθ~=ψ¯0ϕ0-∫θ~=-10∫ξ=0θ~ψ¯ξ-θ~dζθ~ϕξdξ,where ζ(θ~)=ζ(θ~,0). From the above discussion, we can know that ±iω0 are the eigenvalues of A(0), and thus ±iω0 are also the eigenvalues of A∗(0).Letq(θ~)=(1,q2,q3)Teiω0τ0θ~ be the eigenvector of A(0) corresponding to the eigenvalue iω0τ0, and q∗(θ~)=D1(1,q2∗,q3∗)Te-iω0τ0θ~ be the eigenvector of A∗(0) corresponding to the eigenvalue -iω0τ0. With the conditions A(0)q(0)=iω0q(0) and A∗(0)q∗(0)=-iω0q∗(0), we obtain(52)iω0I-G1+G2e-iω0τ0q0=0,and(53)-iω0I-G1T+G2Teiω0τ0q∗0=0.where I is a third order identity matrix. Hence, we can obtain q2=iω0-a11/a12,q3=(iω0-a22)q2-a21/a23, q2∗=-iω0+a11/a21,q3∗=-a23q3/iω0.In order to satisfy thatq∗,q=1, we need to choose the value of D1 in the following part. For convenience, we denote Γ=(1,q2,q3)T,Γ∗=(1,q2∗,q3∗)T. By virtue of (52) and (53), it derives that(54)q∗,q=q∗¯0q0-∫θ~=-10∫ξ=0θ~q¯∗ξ-θ~dζθ~qξdξ,=D1¯Γ∗¯ΓT-D1¯∫θ~=-10∫ξ=0θ~Γ∗¯e-iω0τ0ξ-θ~dζθ~ΓTeiω0τ0ξdξ=D1¯Γ∗¯ΓT-D1¯∫-10Γ∗¯θ~eiω0τ0θ~dζθ~ΓT=D1¯Γ∗¯ΓT+D1¯τ0e-iω0τ0Γ∗¯G2ΓT=D1¯1+q2q2∗¯+q3q3∗¯+τ0q3∗¯q2a31e-iω0τ0.Hence, we can choose D1¯=[1+q2q2∗¯+q3q3∗¯+τ0q3∗¯q2a31e-iω0τ0]-1.Nextly, we will compute the coordinates to detail center manifoldC0 at σ≠0. Let μt be solution of (49) when σ=0. We define(55)zt=q∗,ut,Wt,θ~=ut-zq-z¯q¯=ut-2Re⁡ztqθ~.On the center manifoldC0, it derives that(56)Wt,θ~=Wzt,z¯t,θ~=W20θ~z22+W11θ~zz¯+W02θ~z¯22+⋯,where z and z¯ are local coordinates for C0 in the direction of q∗ and q∗¯, respectively.We only consider real solutions asW is real. If μt is real, then we get(57)q∗,W=q∗,ut-zq-z¯q¯=q∗,ut-q∗,q-q∗,q¯.For the solution μt of (49), since σ=0, from (49) and (51), we have(58)z˙t=q∗,u˙t=q∗,A0ut+R0ut=A∗0q∗,ut+q∗¯0fut,0≔iω0τ0zt+gz,z¯,where(59)gz,z¯=q∗¯0f0Wz,z¯,0+2Re⁡ztqθ~=g20z22+g11zz¯+g02z¯22+g21z2z¯2+···.It follows from (55) and (56) that(60)utθ~=W20θ~z22+W11θ~zz¯+W02θ~z¯22+ΓTeiω0τ0θ~z+Γ∗Te-iω0τ0θ~z¯+···.From (41) and (59), it derives that(61)gz,z¯=τ0D1¯1,q2∗¯,q3∗¯a110x0y0b110x0y0+b011y0z0+b020y20c110y0z0+c011zty-1+c200y20+c002y2-1=τ0D1¯A1W201z22+W111zz¯+W021z¯22+z+z¯W202z22+W112zz¯+W022z¯22+q2z+q2¯z¯+A2W203z22+W113zz¯+W023z¯22+q3z+q3¯z¯W202z22+W112zz¯+W022z¯22+q2z+q2¯z¯+A3W202z22+W112zz¯+W022z¯22+q2z+q2¯z¯2+A4W2030z22+W1130zz¯+W0230z¯22+q3z+q3¯z¯W202-1z22+W112-1zz¯+W022-1z¯22+q2ze-iω0τ0+q2¯z¯eiω0τ0+A5W202-1z22+W112-1zz¯+W022-1z¯22+q2ze-iω0τ0+q2¯z¯eiω0τ02.Here A1=a110+b110q2∗¯, A2=b011q2∗¯+c110q3∗¯, A3=b020q2∗¯+c200q3∗¯, A4=c011q3∗¯, A5=c002q3∗¯.According to comparing its coefficients with (59), it derives that(62)g20=2τ0D1¯A1q2+A2q2q3+A3q22+A4q2q3e-iω0τ0+A5q22e-2iω0τ0,g11=τ0D1¯A1q2+A2q2q3¯+A3q2q2¯+A4q2q3¯e-iω0τ0+A5q1q2¯,g02=2τ0D1¯A1q2¯+A2q2¯q3¯+A3q2¯2+A4q2¯q3¯eiω0τ0+A5q22e2iω0τ0,g21=τ0D1¯A1W201q2¯+2q2W111+2W112+W202+τ0D1¯A2W202q3¯+2q3W112+2q2W113+q2¯W203+τ0D1¯A3W202q2¯+2q2W112+2q2W112+W202q2¯+τ0D1¯A4W203q2¯eiω0τ0+2q2W113e-iω0τ0+2q3W112-1+W202-1q3¯+τ0D1¯A52W202-1q2¯eiω0τ0+4q2W112-1e-iω0τ0.Since W11 and W20 are still unknown, we still compute them next. By virtue of (49) and (55), we have(63)W˙=u˙t-z˙q-z¯˙q¯=A0W-2Re⁡q¯∗0f0qθ~,θ~∈-1,0A0W-2Re⁡q¯∗0f0q0+f0z,z¯,θ~=0≜A0W0+Hz,z¯,θ~.Here,(64)Hz,z¯,θ~=H20θ~z22+H11θ~zz¯+H02θ~z¯22+···.By computing and comparing coefficients, we can obtain(65)A0-2iω0τ0W20θ~=-H20θ~,A0W11θ~=-H11θ~.From (63), we know that(66)Hz,z¯,θ~=-q∗0f0qθ~-q∗0f0¯q¯θ~=-gz,z¯qθ~-g¯z,z¯q¯θ~.Using (59) in (66) and comparing coefficient with (64), it derives that(67)H20θ~=-g20qθ~-g¯02q¯θ~,(68)H11θ~=-g11qθ~-g¯11q¯θ~.From the definition of A(0), based on (65) and (67), we can obtain (69)W˙20θ~=2iω0τ0W20θ~+g20qθ~+g¯02q¯θ~.Because of q(θ~)=(1,q2,q3)Teiω0τ0θ, we have(70)W20θ~=ig20q0ω0τ0eiω0τ0θ~+ig¯20q¯03ω0τ0e-iω0τ0θ~+E1e2iω0τ0θ~.Similarly, it derives that from (65) and (66)(71)W11θ~=-ig11q0ω0eiω0θ~+ig¯11q¯0ω0e-iω0θ~+E2.Here E1=(E1(1),E1(2),E1(3))T, E2=(E2(1),E2(2),E2(3))T are three dimensional constant vectors and can be determined by setting θ~=0 in H(z,z¯,θ~).By above definitions and conditions, we can compute(72)E11=2Υ11Υ1,E12=2Υ12Υ1,E13=2Υ13Υ1,E21=Υ21Υ2,E22=Υ22Υ2,E23=Υ23Υ2,where(73)Υ1=2iω0-a110-c110+c001b001-ω02-2iω0a010b100,Υ11=a1012iω0ς12-b001ς13+ς11-4ω02-2ω0b010-b001c100+c001,Υ12=2iω0-a1002iω0ς12+b001ς13+2iω0b100ς11,Υ13=ς132iω0-a1002iω0-b010-a010b100+c100+c001ς122iω0-a100+b100ς11,Υ2=-a100b001c100+c001,Υ21=-ς21b001c100+c001+ς23a010b001,Υ22=-ς23a100b001,Υ23=ς23a100b010-a010b100+c100+c001-ς22a100+-ς21b100,where(74)ς11=a110q2,ς12=b110q2+b011q2q3+b020q22,ς13=2c110q2q3+c011q2q3eiω0τ0+c200q22+c002q22e-2iω0τ0ς21=2a110q2,ς22=2b110q2+2b011q2q3¯+2b020q2q2¯,ς23=2c100q2q3¯+2c011q2q3¯e-2iω0τ0+2c200q2q2¯+2c002q2q2¯.Based on above analysis, the following values can be given(75)C10=i2ω0τ0g20g11-2g112-g0223+g212,μ2=-Re⁡C0Re⁡λ′τ0,β2=2Re⁡C0,T2=-Im⁡C0+μ2Im⁡λ′τ0ω0τ0.By computing (75), we have the following theorem.Theorem 9.
Properties of bifurcating periodic solution in the center manifold at critical valueτ0 are determined as follows:
(i) Ifμ2>0(μ2<0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τ0(τ<τ0).
(ii) Bifurcating periodic solutions are stable (unstable) ifβ2<0(β2>0).
(iii) Period of the periodic solution increases (decreases) ifT2>0(T2<0).


## 3.3. Global Hopf Bifurcation

In the above parts, we have obtained that system (2) undergoes Hopf bifurcation at E∗(N∗,P∗,Z∗) when τ=τkj. We will study the global continuation of these local bifurcating periodic solutions by using the global Hopf bifurcation theorem for delay differential equations [30], and state that system (2) admits periodic solutions globally for delay τ in a finite interval [τkj,+∞],k=1,2,3;j=1,2,…. This includes the values that are not necessarily near the local Hopf bifurcation values.For convenience, we give some notations,R+3={(N,P,Z):N>0,P>0,Z>0}, X=C([-τ,0],R+3), and x(t)=(N(t),P(t),Z(t))T. We can rewrite system (2) as the following function(76)dxtθ^dt=Fxtθ^,τ,T,where xt(θ^)=(N(t+θ^),P(t+θ^),Z(t+θ^))∈X.According to the work of Wu [30], we need to introduce the following notions(77)F^=FR+3×R+×R+:R+3×R+×R+→R+3,ΣF=clx,τ,T∈X×R+×R+:xis  a  T-periodic  solution  of76,NF=x¯,τ,T:Fx¯,τ,T=0.Let C(E∗,τkj,2π/ωk) be the connected component of C(E∗,τkj,2π/ωk) in Σ, and denote Projτ(E∗,τkj,2π/ωk) as its projection on τ component. ±iωk is a pair purely imaginary roots of (27) with τ=τkj. By Theorem 8, we know that C(E∗,τkj,2π/ωk) is nonempty.Lemma 10.
Whenever the equilibriumE∗ exists, all the nontrivial solutions of system (2) are uniformly bounded.Proof.
We assume that(N(t),P(t),Z(t)) is an arbitrary positive periodic solution of system (2) subject to initial condition (3). Let t1 and t2 be real quantities defined by N(t1)=min{N(t)}, N(t2)=max{N(t)}, P(t1)=min{P(t)}, P(t2)=max{P(t)}, Z(t1)=min{Z(t)}, and Z(t2)=max{Z(t)}. From the first and second equations of system (2), we have(78)0=N0-αNt2-bNt2Pt2+k1Pt2<N0-bNt2Pt2+k1Pt20=b1Nt2Pt2-kPt2-βPt2Zt2r+P2t2-q1EPt2<b1Nt2Pt2-kPt2-q1EPt2,which clearly implies that(79)Nt2Pt2<-N0k+q1Ek1b1-kb-q1bE≡θ1.Since k1b1<kb, θ1>0.
Similarly, from the second equation of system (2), we have(80)0=b1Nt2Pt2-kPt2-βPt2Zt2r+P2t2-q1EPt2<b1Nt2Pt2-βPt2Zt2r+P2t2.Thus, we obtain(81)βPt2Zt2r+Pt22<b1Nt2Pt2<b1θ1.
We defined the functionW1(t)=N(t)+b/b1P(t)+bβ/b1β1Z(t), for all t>0. Then it follows from system (2) that(82)W˙1=N0-αN+k1p-bb1kP-bb1βPZr+P2-bb1q1EP+bβb1β1β1Pt-τZr+P2t-τ-bβb1β1θPZr+P2-bβb1β1μ+q2EZ≤N0-αN-bb1k+q1EP+bβb1β1β1Pt-τZr+P2t-τ-βbb1β1μ+q2EZ+k1P≤N0-αN-bb1k+q1EP-bβb1β1μ+q2EZ+bθ1≤-m~W1t+N0+bθ1,where m~=min⁡{α,k+q1E,μ+q2E}. Therefore, by comparison with theory [31], we have(83)W1t≤W10-N0+bθ1m~e-m~t+N0+bθ1m~.Thus, all the nontrivial solutions of system (2) are uniformly bounded for t≥0.Lemma 11.
System (2) has no nonconstant τ-periodic solution if the condition (H1) and(84)H4:sup⁡a11+a22+a23,a31+a32+a11+a12,a21+a22<0are satisfied.Proof.
In order to explain that system (2) has no nonconstant periodic solution of periodic τ, we only prove that system (4) has no nonconstant periodic solution (see [36], Lemma 5.5). We rewrite system (4) as u˙(t)=f(u(t)). Suppose that u(t) is nonconstant periodic solution for system (4); then u(t)=(u1,u2,u3) is nonconstant periodic solution of the following equations(85)u˙1=N0-αu1-bu1u2+k1u2,u˙2=b1u1u2-ku2-βu2u3r+u22-q1Eu2,u˙3=β1u2u3r+u22-μu3-θu2u3r+u22-q2Eu3.From the Lemma 10, we know that the periodic orbits of system (85) corresponding to u(t) are included in the following area(86)G=u∈R3:m¯<uj<M¯,j=1,2,3,where m¯ and M¯ are the boundedness of the periodic solutions of system (2) in the G.
The Jacobian matrix of system (85) is(87)Ju=a11a120a21a22a230a31+a320.By using the method of [37], its second additive compound matrix is(88)J2u=a11+a22a230a31+a32a11a120a21a22.The second compound system is(89)x˙1=a11+a22x1+a23x2,x˙2=a31+a32x1+a11x2+a12x3,x˙3=a21x2+a22x3.We select vector module for R3 as follows(90)x1,x2,x3=max⁡x1,x2,x3.Then, the Lozinski uˇ measure μ(J[2](u)) of second additive compound matrix corresponding to a vector module above is(91)μJ2u=max⁡a11+a22+a23,a31+a32+a11+a12,a21+a22.Thus, if (H4) holds, there exists an ε such that μ(t)≤-ε<0 for all u∈G.This shows that the conclusion of Lemma11 follows from Corollary 3.5 in [38].Theorem 12.
Assume(H1) and (H4) are fulfilled; then there exists τ∈I, such that system (2) has at least one positive periodic solution when τ varies among τ1j, τ2j, and τ3j, j≥1, where(92)I=τkj>0:τkjare  the  Hopf  bifurcation  points  of  system2.Proof.
The characteristic matrix of system (2) at an equilibrium E¯=(N¯,P¯,Z¯) is given(93)ΔE¯,τ,Tλ=a¯11-λa¯120a¯21a¯22-λa¯230a¯31e-λτ+a¯32-λ,where(94)a¯11=-α+bP¯,a¯12=-bN¯+k1,a¯21=b1P¯,a¯22=b1N¯-βZ¯r-P¯2r+P¯22-k-q1E,a¯23=-βP¯r+P¯2,a¯31=β1Z¯r-P¯2r+P¯22,a¯32=-θZ¯r-P¯2r+P¯22.(E¯,τ¯,T¯) is called a center if F^(E¯,τ¯,T¯)=0 and det⁡ΔE¯,τ¯,T¯t2πi/T=0. A center is said to be isolated if it is the only center in some neighbourhood of (E¯,τ¯,T¯). It is obvious that (76) has an equilibrium x∗=E∗ in R+3 under the condition (H1). It follows that E∗,τkj,2π/ωk is an isolated center from the discussion about the local Hopf bifurcation theorem in Section 3.1. By Lemma 7 and Theorem 8, there exist constants ϵ>0 and δ>0, as well as a smooth curve λ:τkj-δ,τkj+δ→C such that(95)det⁡ΔE¯,τ¯,T¯λτ=0,λ-iωk<ϵwhen τ∈[τkj-δ,τkj+δ] holds, and(96)λτkj=iωk,dRe⁡λdττ=τkj≠0.Denotes Tk=2π/ωk, and(97)Ωϵ,T=u,T:0<u<ϵ,T-Tk<ϵ.Apparently, if τ-τkj≤δ, (u,T)∈∂Ωϵ,T such that det⁡ΔE∗,τkj,Tλu+i2π/T=0 if and only if u=0,τ=τkj. Hence all the assumptions (A1)-(A4) of Theorem 3.2 in [30] are satisfied for m=1. Moreover, if we define(98)H±E∗,τkj,2πωku,T=det⁡ΔE∗,τkj+δ,Tu+i2πT,then we have the crossing number(99)γ1E∗,τkj,2πωk=degBH-E∗,τkj,2πωk,Ωϵ,2π/ωk-degBH+E∗,τkj,2πωk,Ωϵ,2π/ωk=-1,dRe⁡λτ1,3jdτ1,3j>0,1,dRe⁡λτ2jdτ2j<0.By Theorem 3.2 of Wu [30], we conclude that the connected component C(E∗,τkj,2π/ωk) through (E∗,τkj,2π/ωk) in ∑‍ is nonempty. Meanwhile, we have that
(i)C(E∗,τkj,2π/ωk) is unbounded or
(ii)C(E∗,τkj,2π/ωk) is bounded, CE∗,τkj,2π/ωk⋂N(F) is finite, and ∑(E¯,τ,T)∈C(E∗,τkj,2π/ωk)⋂N(F)‍γ1(E¯,τ,T)=0.
By expression (32), we all know that there exists an integer n>0, such that 2π<τkjωk<2(n+1)π, τkj∈J, which implies that τkj/n+1<2π/ωk<τkj.
Therefore, we have thatτ/n+1<T<τ if (E¯,τ,T)∈C(E∗,τkj,2π/ωk). This fact and Lemma 11 show that the projection of C(E∗,τkj,2π/ωk) onto the T-space is bounded. On the flip side, Lemma 10 implies that the projection of C(E∗,τkj,2π/ωk) onto the x-space is uniformly bounded for τ∈Projτ(E∗,τkj,2π/ωk). Otherwise, C(E∗,τkj,2π/ωk) is bounded and Projτ(E∗,τkj,2π/ωk)⊂[τk,∞], which contradicts the fact that nonconstant periodic solution and T are both uniformly bounded for τ∈Projτ(E∗,τkj,2π/ωk). The proof is completed.


## 4. Optimal Harvesting

The emphasis of this section is on the profit making aspect of plankton; optimal control addresses the problem of finding control law for a given system such that a certain optimality criterion is satisfied [25]. In the economic perspective, the commercial development of recyclable resources is a critical issue and it is necessary to determine the optimal tradeoff between present and future harvests.The harvesting problem covers cost of harvesting comprising a set of differential equations that describe the paths of the control variables that minimize the cost functional. It is assumed that price is a function that decreases with increasing biomass. Thus, to maximize the total discounted net revenues from the plankton, the optimal control problem can be formulated as(100)JE=∫0tfp1-v1q1EPq1EP+p2-v2q2EZq2EZ-cEdt,where c is the constant harvesting cost per unit effort, p1 and p2 are the constant price per unit biomass of the phytoplankton and zooplankton, and v1, v2 are the economic constants. The convexity of the objective function with respect to E and the compactness of range values of the state variables can be combined to give the existence of the optimal control.This problem is subject to system (2) and the control constrains 0<E(t)<Emax; here, Emax is the maximum effect capacity of the harvesting industry. We want to find the optimal harvesting effect E∗ such that(101)JE∗=max⁡JE∣E∈U,where U is the control set defined by(102)U=E∣Eis  measurable  and0≤E≤Emaxfor  allt.To find the maximum value ofJ(E), we set X(t)=(N(t),P(t),Z(t)) and Xτ(t)=(Nτ(t),Pτ(t),Zτ(t)), where Nτ(t)=N(t-τ),Pτ=P(t-τ),Zτ(t)=Z(t-τ). Now, the Hamiltonian function of this optimal control problem is given(103)Ht,Xt,Xτt,Et,λt=p1-v1q1EtPtq1EtPt+p2-v2q2EtZtq2EtZt-cEt+λ1tN0-αNt-bNtPt+k1Pt+λ2tb1NtPt-kPt-βPtZtr+P2t-q1EtPt+λ3tβ1Pt-τZtr+P2t-τ-μZt-θPtZtr+P2t-q2EtZt,where λi=λi(t),i=1,2,3, are adjoint variables. Transversality conditions are λi(tf)=0,i=1,2,3.Pontryagin’s Maximum Principle with delay given in [28] is used to derive the necessary conditions for the optimal control. If we consider X(t) and Xτ(t) defined above, then there exists a continuous function λ(t) on [0,tf] satisfying the following three functions. The state equation is(104)X′t=HλXt,Xτt,Et,λt,the optimality condition is(105)0=HEXt,Xτt,Et,λt,and the adjoint equation is(106)λ′t=-HXt,Xt,Xτt,Et,λt-HXτt+τ,Xt+τ,Xτt+τ,Et+τ,λt+τ,where Hλ, HE, HX, and HXτ denote the derivative with respect to λ, E, X, Xτ, respectively. Now we use the necessary conditions to Hamiltonian H.Then there exist adjoint equations(107)λ1′t=-∂H∂Nt-∂H∂Nτt+τ=λ1α+bPt-λ2b1Pt,λ2′t=-∂H∂Pt-∂H∂Pτt+τ=2v1q12E2tPt-p1q1Et+λ1bNt-k1+λ2-b1Nt+k+βZtr-P2tr+P2t2+q1Et+λ3θZtr-P2tr+P2t2+λ3t+τ-β1Zt+τr-P2tr+P2t2,λ3′t=-∂H∂Zt-∂H∂Zτt+τ=2v2q22E2Z-p2q2Et+λ2βPtr+P2t+λ3t-β1Pt-τr+P2t-τ+μ+θPtr+P2t+q2Et.The condition for the optimal harvesting can be obtained from (105), which is (108)Et=p1q1Pt+p2q2Zt-c-λ2q1Pt-λ3q2Zt2v1q12P2t+v2q22Z2t.Based on the above analysis, we have the following theorem.Theorem 13.
There exist an optimal controlE(t) and corresponding solutions to system (2) that maximin J(E) over U. Furthermore, there exist adjoint variables λ1(t),λ2(t), and λ3(t) which satisfy (107) with the transversality conditions λi(tf)=0,i=1,2,3. Moreover the optimal control is given by expression (108).


## 5. Numerical Simulation

In this section, numerical simulations are carried out to verify the theoretical findings by the use of MATLAB. Most values of parameters in system (2) are taken from [21]. N0=1, α=0.485, b=0.388, k1=0.0555, b1=0.35, k=0.07, β=0.76, q1=0.3, E=0.5, β1=0.66, μ=0.03, θ=0.085, q2=0.1, r=5.Whenτ=0, we observe that (H1) holds and the coexisting equilibrium is E∗=(1.37,0.78,1.90). Condition (H1′) is also valid. Therefore, the equilibrium E∗ is locally asymptotically stable (see Figure 1). Then, when τ≠0, we obtained h(0)=p3=-0.0002<0. Hence, (30) has at least one positive root. We obtain ω0=0.13, and the critical value τ0=8.43. According to Theorem 8, system (2) remains stable when 0<τ<τ0, but system (2) is unstable when τ>τ0. The stability behaviour of system (2) around E∗ for τ=8.3 is shown in Figure 2. Hopf bifurcation occurs when τ=9>τ0, shown in Figure 3.Figure 1
Waveform plot and phase of system (4) at the equilibrium (1.37,0.78,1.90). (a) Waveform plot, (b) phase diagram.
(a)
Waveform plot (b)
Phase diagramFigure 2
The equilibriumE∗ of system (2) is locally asymptotically stable with initial value (1.2,1,1.7) when τ=8.3<τ0. (a) Waveform plot, (b) phase diagram.
(a)
Waveform plot (b)
Phase diagramFigure 3
The form of limit cycle of system (2) around E∗ with initial value (1.2,1,1.6) when τ=9>τ0. (a) Waveform plot, (b) phase diagram.
(a)
Waveform plot (b)
Phase diagramAgain, through the above parameter values, the properties of bifurcating periodic solutions atτ0 are given(109)C10=-307.56-114.2,μ2=53958,β2=-615.12,T2=19131.Therefore, Hopf bifurcation is supercritical, bifurcating periodic solutions are stable (see Figure4), and period of the periodic solution increases according to Theorem 9.Figure 4
A stable periodic orbits atτ=13. (a) Regular oscillations of system (2), (b) phase diagram.
(a)
Regular oscillations (b)
Phase diagramNextly, we consider a set of parametric values to illustrate optimal control problem as follows:N0=1, α=0.485,b=0.388, k1=0.0555,b1=0.35, k=0.07,β=0.76, q1=0.3,β1=0.66, μ=0.03,θ=0.085, q2=0.1,r=5, c=0.1,a=0.2, δ=0.2,v1=0.4, v2=0.7,p1=1,p2=5. We first use forward forth order Runge-Kutta method for the state variables of (4) with nonnegative initial conditions. Then we solve the adjoint variables from (107) by using backward forth order Runge-Kutta method. It is observed that as the constant input of nutrient concentration or the half saturation constant increases, the optimal harvesting effort increases (see (a,b) of Figure 5). However, as the rate of toxin liberation increases, the optimal harvesting effort decreases (see (c) of Figure 5). We observe that the optimal harvesting effort tends to be stable. We see that the nutrient increases firstly and then tends to be stable as time varies in the presence of harvesting (see (a) of Figure 6), and the phytoplankton and zooplankton population decrease and even tend to be stable as time changes in the presence of harvesting (see (b,c) of Figure 6). All the final values of adjoint variables are zero (see Figure 7). Finally, using the same method in subject (100), we research optimal harvesting policy for model (2). Optimal harvesting effort is mainly dependent on gestation delay τ of zooplankton displayed in Figure 8. We can observe that the optimal harvesting effort increases monotonically.Figure 5
Optimal plankton harvesting with respect to different parameters: (a) the constant input of nutrient concentrationN0, (b) the half saturation constant r, (c) the rate of toxin liberation θ.
(a)
N 0 (b)
r (c)
θFigure 6
Diagrams for the state variables with control for system (4): (a) the nutrient, (b) the phytoplankton population, (c) the zooplankton population.
(a)
The nutrient (b)
The phytoplankton (c)
The zooplanktonFigure 7
Diagrams for the adjoint variables of system (4): (a) λ1, (b) λ2, (c) λ3.
(a)
λ 1 (b)
λ 2 (c)
λ 3Figure 8
Variation of the optimal effortE with respect to time with different τ in system (2).


## 6. Conclusion and Discussion

Sharma et al. [21] discussed the effect of the time delay about the liberation of toxin on nutrient-plankton interaction. In this paper, the dynamic of nutrient-plankton interaction was considered with the gestation delay and linear harvesting. It is assumed that the mortality caused by toxin to zooplankton population is described as Holling type IV function. We obtained that the system lost local stability around E∗ and Hopf bifurcation occurred when time delay crossed its corresponding critical value. By using normal form theory and center manifold theorem, we obtained the properties of Hopf bifurcation. Next, as time delay varies in some regions, the global existence result of the periodic solution of Hopf bifurcation was given by a global Hopf bifurcation result due to Wu [30]. Pontryagin’s Maximum Principle was used to obtain the optimal harvesting policy. We analyzed the effect of different parameter variables on the harvesting effort (see of Figure 5) and obtained the impact of time delay τ on the harvesting effort (see of Figure 8). In the sense of biology, optimal harvesting effort increased with increasing of the constant input of nutrient concentration; optimal harvesting effort stayed the same as the half saturation constant changes; optimal harvesting effort decreased with the increasing toxin liberation of plankton. And optimal harvesting effort increased with the longer time delay in maturation of zooplankton.In the future research, in order to make the model more practical, harvesting might be delayed because zooplankton avoid harvesting, so this factor can be considered into our nutrient-plankton system, which is(110)N˙t=N0-αN-bNP+k1P,P˙t=b1NP-kP-βPZr+P2-q1EP,Z˙t=β1Pt-τ1Zr+P2t-τ1-μZ-θPZr+P2-q2EZt-τ2,where τ1 and τ2 are the time delays in maturation of zooplankton and harvesting because zooplankton avoid harvesting, respectively.



---

*Source: 9232590-2018-12-09.xml*