# Dynamics Induced by Delay in a Nutrient-Phytoplankton Model with Multiple Delays

Authors: Chuanjun Dai; Hengguo Yu; Qing Guo; He Liu; Qi Wang; Zengling Ma; Min Zhao
Journal: Complexity
 (2019)

---

## Abstract

A nutrient-phytoplankton model with multiple delays is studied analytically and numerically. The aim of this paper is to study how the delay factors influence dynamics of interaction between nutrient and phytoplankton. The analytical analysis indicates that the positive equilibrium is always globally asymptotically stable when the delay does not exist. On the contrary, the positive equilibrium loses its stability via Hopf instability induced by delay and then the corresponding periodic solutions emerge. Especially, the stability switches for positive equilibrium occur as the delay is increased. Furthermore, the numerical simulations show that periodic-2 and periodic-3 solutions can appear due to the existence of delays. Numerical results are consistent with the analytical results. Our results demonstrate that the delay has a great impact on the nutrient-phytoplankton dynamics.


---

## Body

## 1. Introduction

Some phytoplankton, for example, Cyanobacteria, can form dense and sometimes toxic blooms in freshwater and marine environments, which threaten ecological balance, drinking water, fisheries, and even human health [1]. However, the mechanism, by which phytoplankton blooms occur, is currently not very clear, which contribute to the difficulty to prevent or mitigate the proliferation of phytoplankton blooms. These have stimulated lots of researches aiming to understand the growth mechanisms of phytoplankton.In recent years, dynamics in phytoplankton growth have drawn increasing attention from experimental ecologists, as well as mathematical ecologists. Some results from experiments and field observations imply that many factors affecting the dynamics of phytoplankton growth are bound to exist, such as nutrient [2], light [3], temperature [4], iron supply [5], zooplankton [6]. Especially, due to the effects of limiting factors including temperature, light, and day length, it has been indicated by Rhee and Gotham [7] that the population dynamics of phytoplankton in aquatic environments can change with season, latitude, and depth. Among factors affecting phytoplankton growth, nutrient has been an essential element [8–10], mainly including nitrogen and phosphate. Results reported by Ryther [11] indicated that phytoplankton indeed consumes lots of nitrogen and phosphate in their growth process, but reducing the nitrogen content in aquatic cannot slow the eutrophication. Using data from 17 lakes, Smith [8] analysed the influence of ratio of total nitrogen to phosphorus on the growth of blue-green algae (Cyanophya) and showed that controlling the ratio can help us improve the quality of aquatic environment very well. Obviously, the production process of phytoplankton is more complex.However, due to the complexity and nonlinearity of aquatic ecosystem, there are some difficulties in understanding nutrient-phytoplankton dynamics only depending on experiment or field observation, which makes it necessary to use models to provide quantitative insights into dynamic mechanism of phytoplankton growth. For different aquatic environments, we can use various modifications of the classical prey-predator models by introducing functional responses to model nutrient-phytoplankton dynamics [12–14]. For example, Huppert et al. [15] describe the dynamics of nutrient-driven phytoplankton blooms by a simple model and identify, using the model analysis, an important threshold effect that a bloom will only be triggered when nutrients exceed a certain defined level. Additionally, most nutrient-phytoplankton models reveal that phytoplankton population and nutrient population can coexist at equilibrium globally under some conditions [16, 17]. However, Sherratt and Smith [18] have reported that a constant population density may not exist in reality because of the existence of some factors, such as noise and physical factors. Actually, experiments and field observations show that the changes of phytoplankton population density usually possess oscillatory behaviour [19, 20].For the single cell phytoplankton species, in most studies of nutrient-phytoplankton models, it is usually assumed that the processes, such as conversion process of nutrient, in the dynamics of phytoplankton growth are instantaneous [14–17, 19–23]. It may be doubtful whether there exists the delay in the growth of phytoplankton over the large area or not. Yet, J. Caperon [24] studied time lag in population growth response ofIsochrysis galbana, a phytoplankton species, to a variable nitrate environment by both experiments and model, and demonstrated the existence of delay in the growth ofIsochrysis galbana. Hence, the delay may indeed exist in the phytoplankton growth, which means that it is necessary to consider delay in nutrient-phytoplankton models. An approach that has been attempted by researchers to model the dynamics of phytoplankton is the role of delay since delay appears as an important component in biosystems and ecosystems [25–30].Actually, growing evidence shows that there exists time lag in some conversion processes from one state to another in some systems, and delay is an important factor because it can affect the dynamics of these systems. Volterra [31] considered time delay in a prey-predator model first and found oscillatory behaviour for the spatial distribution. For a long time, it has been recognized that delays can give rise to destabilizing effect of the dynamics of systems, where periodic solutions, as well as chaos, may emerge [32–35]. Models incorporating delays in diverse biological and ecological models are extensively studied [36–42]. Especially, the characteristic equation with respect to the linearized system of delay differential equations plays a key role in dynamic analysis, by which we can obtain some information on the stability of equilibrium. In addition, using the normal form theory, one can carry out the bifurcation analysis, such as the direction and stability of periodic solutions arising through Hopf bifurcation [43, 44].The main purpose of this paper is to consider the effects of multiple delays on the nutrient-phytoplankton dynamics. In [15], Huppert et al. presented a simple model to investigate effect of nutrient on phytoplankton blooms, and much better results are obtained. Here, this model is extended into a “two preys-one predator” type to describe nitrogen- phosphorus-phytoplankton dynamics, as follows:(1)dNdt=I1-q1N-α1NAdPdt=I2-q2P-α2PAdAdt=β1α1Nt-τ1A+β2α2Pt-τ2A-mAwhere N, P, and A represent nitrogen, phosphorus, and phytoplankton population density at time t, respectively; I1 is the nitrogen nutrients input flowing into the system and I2 is the phosphorus nutrients input flowing into the system; q1 is the loss rate of the nitrogen nutrients, and q2 is the loss rate of the phosphorus nutrients; α1 is nitrogen nutrient uptake rate of phytoplankton, and α2 is phosphorus nutrient uptake rate of phytoplankton; β1 and β2 denote the efficiency of nutrient utilization; τ1 and τ2 are time delay parameters; m is the mortality rate of phytoplankton. Although the function, which describes nutrient uptake dynamics, is not a Michaelis-Menten function, but Lotka-Volterra type, Huppert et al. [15] have indicated that the Lotka-Volterra term is a good first approximation to the Michaelis-Menten type. From biological viewpoint, all parameters are nonnegative. Nt, Pt, and A(t)≥0 are continuous on -τ≤t<0, where τ=max⁡(τ1,τ2) and N0, P0, and A0>0.The paper is organized as follows. In Section2, we analyze the existence and stability of positive equilibrium in model (1) without delays. In Section 3, we discuss stability of positive equilibrium and Hopf bifurcation under five different cases for delay effect. Subsequently, the direction of bifurcation and the stability of periodic solutions arising through Hopf bifurcation are given in Section 4. In order to analyze further how delay effects influence nutrient-phytoplankton dynamics, a series of numerical simulations are carried out in Section 5. Finally, the paper ends with conclusion in Section 6.


## 2. Existence and Stability of Positive Equilibrium in Model (1) without Delays

In this section, it is presented first that the first octant is positive invariant in model (1) without delays and the following lemma holds.Lemma 1.
All the solutions of model (1) with initial conditions that initiate in {R+3} are positive invariant in the absence of delays.Proof.
From the first equation of model (1), we have(2)dNdt=I1-q1N-α1NA≥-q1+α1AN.Hence, Nt≥N0exp⁡-∫0t(q1+α1A)ds>0 under N0>0.
Likewise, from the second equation of model (1), we have Pt≥P0exp⁡-∫0t(q2+α2A)ds>0 under P0>0.
In the absence of delays in model (1), from the third equation of model (1), if A0>0, it can be obtained that(3)At=A0exp⁡∫0tβ1α1N+β2α2P-mds>0
Obviously, all the solutions of model (1) without delays are positive invariant if the initial conditions initiate in {R+3}.
Then, we complete the proof.For model (1), it is obvious that the extinction equilibrium, (I1/q1,I2/q2,0), exists. Moreover, in order to discuss the existence of positive equilibrium, the following function is defined:(4)fx=α1α2mx2+mq1-α1β1I1α2+mq2-α2β2I2α1x+mq1q2-β1α1I1q2-β2α2I2q1,and then we can obtain(5)N∗=I1q1+α1A∗,P∗=I2q2+α2A∗,where A∗ is the positive root of (4).For the functionf(x), we have f0=0 when the condition, m=β1α1I1/q1+β2α2I2/q2, holds, and then there is no positive equilibrium in model (1). Obviously, f0<0 holds if m<β1α1I1/q1+β2α2I2/q2, and then there exists a unique positive root with respect to fx=0, which means that there exists a unique positive equilibrium in model (1) under this condition. However, when m>β1α1I1/q1+β2α2I2/q2, it can be verified directly that f0>0 and mq1-α1β1I1α2+mq2-α2β2I2α1>0, which implies that there is no positive equilibrium in model (1). Thus, summarizing these results, the following theorem can be obtained.Theorem 2.
If0<m<β1α1I1/q1+β2α2I2/q2 holds, then there exists a unique positive equilibrium in model (1); otherwise, there is no positive equilibrium in model (1) if m≥β1α1I1/q1+β2α2I2/q2 holds.Letting the unique positive equilibrium beE∗=(N∗,P∗,A∗), then the following theorem holds for model (1) without delays.Theorem 3.
If the unique positive equilibrium exists in model (1) without delays, then it is globally asymptotically stable.Proof.
We construct a Lyapunov function, as follows:(6)V=β1∫N∗Ns-N∗sds+β2∫P∗Ps-P∗sds+∫A∗As-A∗sdsIn the model (1) without delays,(7)dVdt=β1N-N∗NdNdt+β2P-P∗PdPdt+A-A∗AdAdt=β1N-N∗NI1-q1N-α1NA+β2P-P∗PI2-q2P-α2PA+A-A∗Aβ1α1NA+β2α2PA-mA=β1N-N∗Nq1N∗+α1N∗A∗-q1N+α1NA+β2P-P∗Pq2P∗+α2P∗A∗-q2P+α2PA+A-A∗Aβ1α1NA+β2α2PA-β1α1N∗+β2α2P∗A=-β1q1+α1A∗N-N∗2N-β2q2+α2A∗P-P∗2P
Obviously,dV/dt≤0 holds under existence of positive equilibrium and dV/dt=0 holds if and only if N=N∗ and P=P∗. The largest invariant subset of the set of the point where dV/dt=0 is E∗=(N∗,P∗,A∗). Therefore, according to LaSalle’s theorem, E∗=(N∗,P∗,A∗) is globally asymptotically stable.
Then, we complete the proof.Letting the extinction equilibrium beE0=N0,P0,0=(I1/q1,I2/q2,0), then we can obtain the following theorem in model (1) in the absence of delay.Theorem 4.
In the absence of delays, letm∗=β1α1I1/q1+β2α2I2/q2, so that(i)
ifm>m∗, then the extinction equilibrium E0 is locally asymptotically stable;(ii)
ifm<m∗, then the extinction equilibrium E0 is unstable;(iii)
ifm=m∗, then the model (1) undergoes transcritical bifurcation at the extinction equilibrium E0.Proof.
For simplicity, let(8)fwX,m=I1-q1N-α1NAI2-q2P-α2PAβ1α1NA+β2α2PA-mAandX=N,P,AT.
The Jacobian matrix atE0 is (9)
J
E
0
=
-
q
1
0
-
α
1
I
1
q
1
0
-
q
2
-
α
2
I
2
q
2
0
0
-
m
-
m
∗
.The eigenvalues are -q1,-q2,-m-m∗.
Obviously, ifm>m∗, then the extinction equilibrium E0 is locally asymptotically stable.
Ifm<m∗, then the extinction equilibrium E0 is unstable.
Whenm=m∗, the Jacobian matrix at E0 is (10)
J
E
0
=
-
q
1
0
-
α
1
I
1
q
1
0
-
q
2
-
α
2
I
2
q
2
0
0
0
.Then J(E0) has a geometrically simple zero eigenvalue with right eigenvector Φ=α1I1q22,α2I2q12,-q12q22T and left eigenvector Ψ=(0,0,1).
Now(11)Dmfw=00-Aand(12)ΨDXDmfwΦE0=q12q22≠0,ΨDXXfwΦ,Φ=Ψ∑i=13eiΦTDXDXfiTΦE0=-2q12q22β1I1α12q22+β2I2α22q12≠0
According to [45], the model (1) undergoes transcritical bifurcation at the extinction equilibrium E0 in the absence of delays.
Then, we complete the proof.Actually, whenm>m∗ holds, the positive equilibrium does not exist, and the extinction equilibrium E0 is globally asymptotically stable. Then, the following theorem holds.Theorem 5.
Ifm>β1α1I1/q1+β2α2I2/q2 holds, then the extinction equilibrium E0=N0,P0,0=(I1/q1,I2/q2,0) is globally asymptotically stable.Proof.
We construct a Lyapunov function, as follows:(13)V=β1∫N0Ns-N0sds+β2∫P0Ps-P0sds+A
In the model (1) without delays,(14)dVdt=β1N-N0NdNdt+β2P-P0PdPdt+dAdt=β1N-N0NI1-q1N-α1NA+β2P-P0PI2-q2P-α2PA+β1α1NA+β2α2PA-mA=β1N-N0Nq1N0-q1N-α1NA+β2P-P0Pq2P0-q2P-α2PA+β1α1NA+β2α2PA-mA=-β1q1N-N02N-β2q2P-P02P+β1α1N0+β2α2P0-mA=-β1q1N-N02N-β2q2P-P02P+β1α1I1q1+β2α2I2q2-mA
Obviously,dV/dt<0 holds if m>β1α1I1/q1+β2α2I2/q2. Therefore, the extinction equilibrium E0=(N0,P0,0)=(I1/q1,I2/q2,0) is globally asymptotically stable when m>β1α1I1/q1+β2α2I2/q2.
Then, we complete the proof.


## 3. Local Stability Analysis and the Hopf Bifurcation

In this section, we first state the following positive invariant theorem.Lemma 6.
All the solutions of model (1) with initial conditions that initiate in {R+3} are positive invariant.Proof.
We consider(N,P,A) a noncontinuable solution of model (1); see [46], defined on [-τ,B), where B∈(0,+∞]. Then we can use the method from [47] to prove that, for all t∈[0,B), Nt>0, Pt>0, and At>0. Suppose that is not true. Then, there exists 0<T<B such that, for all t∈[0,T), Nt>0, Pt>0, and At>0 and either NT=0, PT=0, or AT=0. According to Lemma 1, for all t∈[0,T), we have(15)Nt>N0exp⁡-∫0tq1+α1Ads,Pt>P0exp⁡-∫0tq2+α2Ads,and(16)At=A0exp⁡∫0tβ1α1Ns-τ1+β2α2Ps-τ2-mds.
As(N,P,A) is defined and continuous on [-τ,T], there is a M≥0 such that, for all t∈[-τ,T),(17)Nt>N0exp⁡-∫0tq1+α1Ads≥N0exp⁡-TM,Pt>P0exp⁡-∫0tq2+α2Ads≥P0exp⁡-TMand(18)At=A0exp⁡∫0tβ1α1Ns-τ1+β2α2Ps-τ2-mds≥A0exp⁡-TM.
Taking the limit, ast→T, we can get(19)NT≥N0exp⁡-TM>0,PT≥P0exp⁡-TM>0and(20)AT≥A0exp⁡-TM>0,which contradicts the fact that either NT=0, PT=0, or AT=0. Thus, for all t∈0,B, Nt>0, Pt>0, and At>0.
Therefore, all the solutions of model (1) are positive invariant if the initial conditions initiate in {R+3}.
Then, we complete the proof.Next, we will discuss the stability of the unique positive equilibrium and existence of Hopf bifurcation in model (1) for five different cases: τ1>0, τ2=0; τ1=0, τ2>0; τ1=τ2=τ; τ1∈0,τ10, τ2>0; τ1>0, τ2∈0,τ20.According to Theorem2, let u1=N(t)-N∗, u2=P(t)-P∗, u3=A(t)-A∗; the linearized form of model (1) can be obtained as follows:(21)u˙1=a11u1t+a13u3tu˙2=a22u2t+a23u3tu˙3=a31u1t-τ1+a32u2t-τ2where (22)a11=-q1-α1A∗;a13=-α1N∗;a22=-q2-α2A∗;a23=-α2P∗;a31=β1α1A∗;a32=β2α2A∗.Then, we obtain the associated characteristic equation of model (21) as follows:(23)λ3+Bλ2+Cλ+Dλe-λτ1+Ee-λτ1+Fλe-λτ2+Ge-λτ2=0,where(24)B=-a11+a22;C=a11a22;D=-a13a31;E=a13a31a22;F=-a23a32;G=a11a23a32.Case 1.
τ 1 > 0, τ2=0.
Due toτ1>0, τ2=0, (23) becomes (25)λ3+Bλ2+C+Fλ+Dλ+Ee-λτ1+G=0.
Assumingλ=iω1(ω1>0) is the pure imaginary root of (25), then the following can be obtained:(26)-ω13+C+Fω1=Esin⁡ω1τ1-Dω1cos⁡ω1τ1,-Bω12+G=-Ecos⁡ω1τ1-Dω1sin⁡ω1τ1,Then(27)ω16+B2-2C+Fω14+C+F2-2BG-D2ω12+G2-E2=0.
Now, we define a function as follows:(28)f1v1=v13+B2-2C+Fv12+C+F2-2BG-D2v1+G2-E2.(i)
If (H21):G2-E2<0 holds, then, (28) has at least one positive root. Without loss of generality, we denote v11, v12, and v13 as the roots of (28); hence ω1k=v1k, k=1, 2, 3, if v1k>0.(ii)
IfH22:G2-E2>0 holds, let M1=B2-2(C+F) and M2=C+F2-2BG-D2. When Δ=M12-3M2≤0, then (28) has no positive roots. However, when Δ>0, f1′v1=0 has two real roots, denoted as(29)x1∗=-M1+Δ3,x2∗=-M1-Δ3.Obviously, limx→+∞⁡f1(x)=+∞. If H22 and Δ=M12-3M2>0 holds, then (28) has two positive real roots if and only if x1∗>0 and f1(x1∗)<0. In addition, we denote two positive roots of (28) as χ1 and χ2; then, (27) has two positive roots, namely, ω1a=χ1 and ω1b=χ2. Furthermore, we can have the following results.Proposition 7.
(i)
If(H21) holds, then (28) has at least one positive root. (ii)
If(H22) and Δ=M12-3M2≤0 holds, then, (28) has no positive root. (iii)
If(H22) and Δ=M12-3M2>0 holds, then, (28) has two positive roots if and only if x1∗>0, and f1(x1∗)<0.
Then, according to (26), the critical delay can be obtained as follows:(30)τ1pj=1ω1parccos⁡Dω1p4+BE-CD+DFω1p2-EGD2ω1p2+E2+2jπ,j=0,1,2,…,p=1,2,3,a,b
Lettingτ10=minp=1,2,3,a,b⁡τ1p0, from [37] we know that Redλ/dτ1≠0 also needs to be proved. Differentiating left side of (26) with respect to τ1, we have(31)dλdτ1-1=3λ2+2Bλ+C+F+De-λτ1λDλ+Ee-λτ1-τ1λ,so that the following can be obtained:(32)Redλdτ1λ=iω10-1=ω102Δf1′ω102,where(33)Δ=-Dω102cos⁡ω10τ1+Eω10sin⁡ω10τ12+Eω10cosω10τ1+Dω102sinω10τ12.If (i) in Proposition7 and (H23):f1′(ω102)≠0 hold, then we have Redλ/dτ1λ=iω10-1≠0. However, if (iii) in Proposition 7 holds, assuming τ1aj<τ1bj, then we obtain f1′(χ1)>0 and f1′(χ2)<0. Hence, we have dRe⁡λτ/dτ|τ=τ1aj>0, (dRe⁡λ(τ)/dτ)|τ=τ1bj<0, and j=0,1,2,…. Therefore, we have the following results.Theorem 8.
For model (1) with τ1>0, τ2=0,(i)
If(H21) and(H23) both hold, then the positive equilibrium E∗ is locally asymptotically stable for τ1∈0,τ10 and Hopf bifurcation occurs at τ1=τ10.(ii)
If (ii) in Proposition7 holds, then the positive equilibrium E∗ is locally asymptotically stable for all τ1≥0.(iii)
If (iii) in Proposition7 holds, there exists a nonnegative integer n, such that the positive equilibrium E∗ is locally asymptotically stable whenever τ1∈[0,τ1a0)∪(τ1b0,τ1a1)∪⋯∪(τ1bn-1,τ1an) and is unstable whenever τ1∈[τ1a0,τ1b0)∪(τ1a1,τ1b1)∪⋯∪(τ1an-1,τ1bn-1)∪(τ1an,+∞). Then, model (1) undergoes Hopf bifurcation around E∗ at every τ1=τ1ajandτ1bj, j=0,1,2,….Case 2.
τ 1 = 0, τ2>0.
Sinceτ1=0, τ2>0, (23) becomes (34)λ3+Bλ2+C+Dλ+Fλ+Ge-λτ2+E=0.
Similar to Case1, let λ=iω2(ω2>0) be the pure imaginary root of (34); then we obtain(35)-ω23+C+Dω2=Gsin⁡ω2τ2-Fω2cos⁡ω2τ2,-Bω22+E=Gcos⁡ω2τ2-Fω2sin⁡ω2τ2.That is,(36)ω26+B2-2C+Dω24+C+D2-2BE-F2ω22+E2-G2=0.
Lettingv2=ω22, we define the following function:(37)f2v2=v23+B2-2C+Dv22+C+D2-2BE-F2v2+E2-G2.(i)
If(H31):E2-G2<0 holds, then (37) has at least one positive root. Without loss of generality, we denote v21, v22, v23 as the roots of (37); then ω2k=v2k, k=1,2,3, if v2k>0.(ii)
IfH32:E2-G2>0 holds, let M1=B2-2(C+D) and M2=C+D2-2BE-F2. When Δ=M12-3M2≤0, then (37) has no positive roots. However, when Δ>0, f2′v2=0 has two real roots, denoted as(38)x1∗∗=-M1+Δ3,x2∗∗=-M1-Δ3.
Obviously,limx→+∞⁡f2(x)=+∞. If H32 and Δ=M12-3M2>0 holds, then (37) has two positive real roots if and only if x1∗∗>0 and f2(x1∗∗)<0. In addition, we denote two positive roots of (37) as χ1∗ and χ2∗; then (36) has two positive roots, namely, ω2a=χ1∗ and ω2b=χ2∗. Furthermore, we can obtain the following results.Proposition 9.
(i)
If(H31) holds, then (36) has at least one positive root. (ii)
If(H32) and Δ=M12-3M2≤0 holds, then, (36) has no positive root. (iii)
If(H32) and Δ=M12-3M2>0 holds, then, (36) has two positive roots if and only if x1∗∗>0 and f2(x1∗∗)<0.
Then the critical delay can be derived by (35):(39)τ2pj=1ω2parccos⁡Dω2p4+BE-CD+DFω2p2-EGD2ω2p2+E2+2jπ,j=0,1,2,…,p=1,2,3,a,b
Letτ20=minp=1,2,3,a,b⁡τ2p0. Differentiating left side of (34) with respect to τ2, we obtain(40)dλdτ2-1=3λ2+2Bλ+C+D+Fe-λτ2λFλ+Ge-λτ2-τ2λ,Hence, we obtain the following:(41)Redλdτ2λ=iω20-1=ω202Δf2′ω202,where(42)Δ=-Fω202cos⁡ω20τ2+Gω20sin⁡ω20τ22+Gω20cos⁡ω20τ2+Fω202sin⁡ω20τ22.If (i) in Proposition9 and (H33):f2′(ω202)≠0 both hold, then Redλ/dτ2λ=iω20-1≠0 is obtained. However, if (iii) in Proposition 9 holds, assuming τ2aj<τ2bj, we obtain f2′(χ1∗)>0 and f2′(χ2∗)<0. Then, we have dRe⁡λτ/dτ|τ=τ2aj>0, dRe⁡λτ/dτ|τ=τ2bj<0, and j=0,1,2,…. Therefore, we have the following theorem.Theorem 10.
For model (1) with τ1=0, τ2>0, (i)
IfH31 and H33 hold, then, the positive equilibrium E∗ is locally asymptotically stable for τ2∈0,τ20 and Hopf bifurcation occur at τ2=τ20.(ii)
If (ii) in Proposition9 holds, then the positive equilibrium E∗ is locally asymptotically stable for all τ2≥0.(iii)
If (iii) in Proposition9 holds, then there exists a nonnegative integer n, such that the positive equilibrium E∗ is locally asymptotically stable whenever τ2∈[0,τ2a0)∪(τ2b0,τ2a1)∪⋯∪(τ2bn-1,τ2an) and is unstable whenever τ2∈[τ2a0,τ2b0)∪(τ2a1,τ2b1)∪⋯∪(τ2an-1,τ2bn-1)∪(τ2an,+∞). Then, model (1) undergoes Hopf bifurcation around E∗ for every τ2=τ2ajandτ2bj, j=0,1,2,….Case 3.
τ 1 = τ 2 = τ.
Whenτ1=τ2=τ, (23) becomes (43)λ3+Bλ2+Cλ+D+Fλe-λτ+E+Ge-λτ=0.
Letingλ=iω3(ω3>0) be the pure imaginary root of (43), then(44)-ω33+Cω3=E+Gsin⁡ω3τ-D+Fω3cos⁡ω3τ,-Bω32=-E+Gcos⁡ω3τ-D+Fω3sin⁡ω3τ,that is,(45)ω36+B2-2Cω34+C2-D+F2ω32-E+G2=0.Let v3=ω32 and define the following function:(46)f3v3=v33+B2-2Cv32+C2-D+F2v3-E+G2,From (46), we can clearly see that f3(0)=-E+G2<0; hence (46) has at least one positive root. Without loss of generality, we denote v31, v32, and v33 as the roots of (46); then we have ω3k=v3k, k=1,2,3, if v3k>0 holds. Hence, the critical delay can be derived by (44):(47)τ3kj=1ω3karccos⁡D+Fω3k4+BE+Gω3k2-CD+Fω3k2D+F2ω3k2+E+G2+2πj,j=0,1,2…,k=1,2,3.
Letτ30=mink=1,2,3⁡τ3k0. Differentiating left side of (43) with respect to τ, we obtain(48)dλdτ-1=3λ2+2Bλ+C+D+Fe-λτλD+Fλ+E+Ge-λτ-τλ.Hence, we obtain the following:(49)Redλdτλ=iω30-1=ω302Δf3′ω302,where(50)Δ=-D+Fω302cos⁡ω30τ+E+Gω30sin⁡ω30τ2+E+Gω30cos⁡ω30τ+D+F·ω302sin⁡ω30τ2.If (H41): f3′(ω302)≠0 holds, then Re⁡dλ/dτλ=iω30-1≠0 is obtained; hence, we have the following result.Theorem 11.
For model (1), when τ1=τ2=τ, if (H41) holds, then the positive equilibrium E∗ is locally asymptotically stable for τ∈0,τ30 and Hopf bifurcation occurs at τ=τ30.Case 4.
τ 1 ∈ 0 , τ 10, τ2>0.
Under this case,τ2 is considered as a parameter. The same as Case 1, let λ=iω2∗ be the root of (23); then we have the following:(51)R51cos⁡ω2∗τ2-R52sin⁡ω2∗τ2=R53,R51sin⁡ω2∗τ2+R52cos⁡ω2∗τ2=R54,where(52)R51=Fω2∗;R53=ω2∗3-Cω2∗-Dω2∗cos⁡ω2∗τ1+Esin⁡ω2∗τ1;R52=G;R54=Bω2∗2-Dω2∗sin⁡ω2∗τ1-Ecos⁡ω2∗τ1.
According to (51), the following holds:(53)F1ω2∗+F2ω2∗sin⁡ω2∗τ1+F3ω2∗cos⁡ω2∗τ1=0,where(54)F1ω2∗=ω2∗6+B2-2Cω2∗4+C2+D2-F2ω2∗2+E2-G2,F2ω2∗=2E-BDω2∗3-2CEω2∗,F3ω2∗=-2Dω2∗4+2CD-BEω2∗2.
Assuming (H51): (53) has finite positive root and denoting as ω2k∗, k=1,2,…,l1, then the critical value can be represented as follows:(55)τ2k∗j=1ω2k∗arccos⁡R51·R53+R52·R54R512+R522+2πj,j=0,1,2…;k=1,2,…l1.
Letτ20∗=min⁡τ2k∗(0), (k=1,2,…l1). Differentiating left side of (23) with respect to τ2, the following is obtained:(56)dλdτ2-1=3λ2+2Bλ+C+D-Dλτ1-τ1Ee-λτ1+Fe-λτ2λFλ+Ge-λτ2-τ2λ,and then, we have(57)Redλdτ2λ=iω20∗-1=F51F53+F54F52F512+F522,where(58)F51=-Fω20∗2cos⁡ω20∗τ20∗+Gω20∗sin⁡ω20∗τ20∗,F52=Gω20∗cos⁡ω20∗τ20∗+Fω20∗2sin⁡ω20∗τ20∗,F53=-3ω20∗2+C+Dcos⁡ω20∗τ1-τ1Ecos⁡ω20∗τ1-τ1Dω20∗sin⁡ω20∗τ1+Fcos⁡ω20∗τ20∗,F54=2Bω20∗-Dsin⁡ω20∗τ1+τ1Esin⁡ω20∗τ1-τ1Dω20∗cos⁡ω20∗τ1-Fsin⁡ω20∗τ20∗.Supposing (H52): F51F53+F54F52≠0, then we have the following.Theorem 12.
For model (1), when τ1∈0,τ10 and τ2>0, if both H51 and H52 hold, then the positive equilibrium E∗ is locally asymptotically stable for τ2∈0,τ20∗ and Hopf bifurcation occurs at τ2=τ20∗.Case 5.
τ 1 > 0, τ2∈0,τ20.
Sinceτ1>0, τ2∈0,τ20, we consider τ1 as a parameter. The same as Case 4, letting λ=iω1∗ be the root of (23), we obtain:(59)R61cos⁡ω1∗τ1-R62sin⁡ω1∗τ1=R63,R61sin⁡ω1∗τ1+R62cos⁡ω1∗τ1=R64,where(60)R61=Dω1∗;R63=ω1∗3-Cω1∗-Fω1∗cos⁡ω1∗τ2+Gsin⁡ω1∗τ2;R62=E;R64=Bω1∗2-Fω1∗sin⁡ω1∗τ2-Gcos⁡ω1∗τ2.According to (59), the following holds:(61)G1ω1∗+G2ω1∗sin⁡ω1∗τ2+G3ω1∗cos⁡ω1∗τ2=0,where(62)G1ω1∗=ω1∗6+B2-2Cω1∗4+C2+F2-D2ω1∗2+G2-E2,G2ω1∗=2G-BFω1∗3-2CGω1∗,G3ω1∗=-2Fω1∗4+2CF-BGω2∗2.
Supposing (H61): (61) has finite positive root ω2k∗, k=1,2,…,l2, then we obtain (63)τ1k∗j=1ω1k∗arccos⁡R61·R63+R62·R64R612+R622+2πj,j=0,1,2…;k=1,2,…l2.
Assumingτ10∗=min⁡τ1k∗(0), (k=1,2,…l2) and differentiating left side of (23) with respect to τ1, therefore, the following is obtained:(64)dλdτ1-1=3λ2+2Bλ+C+F-Fλτ2-τ2Ge-λτ2+De-λτ1λFλ+Ge-λτ2-τ1λ,Hence, we have(65)Redλdτ1λ=iω10∗-1=F61F63+F64F62F612+F622,where(66)F61=-Dω10∗2cos⁡ω10∗τ10∗+Eω10∗sin⁡ω10∗τ10∗,F62=Eω10∗cos⁡ω10∗τ10∗+Dω10∗2sin⁡ω10∗τ10∗,F63=-3ω10∗2+C+Dcos⁡ω10∗τ10∗-τ2Gcos⁡ω10∗τ2-τ2Fω10∗sin⁡ω10∗τ2+Fcos⁡ω10∗τ2,F64=2Bω10∗-Dsin⁡ω10∗τ10∗+τ2Gsin⁡ω10∗τ2-τ2Fω10∗cos⁡ω10∗τ2-Fsin⁡ω10∗τ2,Supposing (H62): F61F63+F64F62≠0 holds, then we obtain the following theorem.Theorem 13.
For model (1), when τ1>0 and τ2∈0,τ20, if both H61 and H62 hold, then the positive equilibrium E∗ is locally asymptotically stable for τ1∈0,τ10∗ and Hopf bifurcation occurs at τ1=τ10∗.


## 4. Properties of Periodic Solution

In this section, we will discuss the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions under Case4 by using normal form method and center manifold theorem [43], and methods of other four cases are similar to Case 4. Assuming τ1∈0,τ10, τ20∗>τ1, and Hopf bifurcation occurs at N∗,P∗,A∗ in model (1) when τ=τ20∗.Letτ2=τ20∗+μ, t=sτ2, x(sτ2)=x^(s), y(sτ2)=y^(s), and z=(sτ2)=z^(s), and we also denote x^(s), y^(s), and z^(s) as x(s), y(s), and z(s). Then, model (1) could be rewritten as follows in C=C-1,0,R3:(67)u˙t=Lμut+fμ,ut,where u(t)=x(t),y(t),z(t)T∈R3. Lμ(ϕ):C→R3 and fμ,u(t) are given as follows: (68)Lμϕ=τ20∗+μAϕ0+Bϕ-τ1τ20∗+Cϕ-1,fμ,ϕ=τ20∗+μf1f2f3T,where(69)A=a110a130a22a23000,B=000000a3100,C=0000000a320,f1=-α1ϕ10ϕ30,f2=-α2ϕ20ϕ30,f3=β1α1ϕ1-τ1τ20∗ϕ30+β2α2ϕ2-1ϕ30.According to the Riesz representation theorem, we know that there exists a functionηθ,μ of bounded variation for θ∈-1,0 such that Lμϕ=∫-10dη(θ,μ)ϕ(θ), for all ϕ∈C-1,0,R3. Choosing (70)
η
θ
,
μ
=
τ
20
∗
+
μ
A
+
B
+
C
,
θ
=
0
τ
20
∗
+
μ
B
+
C
,
θ
∈
-
τ
1
τ
20
∗
,
0
τ
20
∗
+
μ
C
,
θ
∈
-
1
,
-
τ
1
τ
20
∗
0
,
θ
=
-
1Forϕ∈C1-1,0,R3, we define (71)Aμϕ=dϕθdθ,θ∈-1,0∫-10dημ,s·ϕs,θ=0Rμϕ=0,θ∈-1,0fμ,ϕ,θ=0Then, model (1) can be rewritten as(72)u˙t=Aμut+Rμut.For ψ∈C10,1,R3∗, the adjoint operator A∗ of A is defined as follows:(73)A∗ψs=-dψsds,s∈0,1,∫-10dηTt,0ψ-t,s=0.Associated with a bilinear inner product(74)ψs,ϕθ=ψ-0ϕ0-∫-10∫ξ=0θψ-ξ-θdηθϕξdξ,where ηθ=ηθ,0 and we know that ±iω20∗τ20∗ are the eigenvalues of A0 and A∗(0).Chooseqθ=1,q2,q3Teiω20∗τ20∗θ to be the eigenvector of A0 corresponding to the eigenvalue iω20∗τ20∗ and q∗s=D1,q2∗,q3∗eiω20∗τ20∗s to be the eigenvector of A∗(0) corresponding to the eigenvalue -iω20∗τ20∗. By computation, we obtain(75)q2=a23iω20∗-a11a13iω20∗-a22;q3=iω20∗-a11a13;q2∗=a32eiω20∗τ1a11+iω2∗a31eiω20∗τ1a22+iω2∗;q3∗=-a11+iω2∗a31eiω20∗τ1.Besides, from (74) we have (76)D-=11+q2q-2∗+q3q-3∗+q-3∗a31τ1e-iω20∗τ1+q2q-3∗a32τ20∗e-iω20∗τ20∗,such that q∗(s),q(θ)=1 and q∗(s),q¯(θ)=0.Then, according to [44], we obtain the following relevant parameter, which helps to determine the direction and stability of Hopf bifurcation:(77)g20=2D-τ20∗-α1q3-α2q-2∗q2q3+β1α1q3q-3∗e-iω20∗τ1+β2α2q2q3q-3∗e-iω20∗τ20∗,g11=D-τ20∗-α1q3+q-3-α2q-2∗q2q-3+q-2q3+β1α1q-3∗q-3e-iω20∗τ1+q3eiω20∗τ1+β2α2q-3∗q2q-3e-iω20∗τ20∗+q3q-2eiω20∗τ20∗,g02=2D-τ20∗-α1q-3-α2q-2∗q-2q-3+β1α1q-3q-3∗eiω20∗τ1+β2α2q-2q-3q-3∗eiω20∗τ20∗,g21=2D-τ20∗-α1W1130+12W2030+12q-3W2010+q3W1110-α2q-2∗q2W1130+q-212W2030+12q-3W2020+q3W1120+q-3∗β1α1q3W111-τ1τ20∗+q-312W201-τ1τ20∗+12W2030eiω20∗τ1+W1130e-iω20∗τ1+q-3∗β2α2q3W112-1+q-312W202-1+q-212W2030eiω20∗τ20∗+q2W1130e-iω20∗τ20∗,and (78)W20θ=ig20q0ω20∗τ20∗eiθω20∗τ20∗+g-02i3ω20∗τ20∗q-0e-iθω20∗τ20∗+E1e2iθω20∗τ20∗,W11θ=-ig11q0ω20∗τ20∗eiθω20∗τ20∗+g-11iω20∗τ20∗q-0e-iθω20∗τ20∗+E2,where E1and E2 can be determined by the following, respectively:(79)2iω20∗-a110-a1302iω20∗-a22-a23-a31e-2iω20∗τ1-a32e-2iω20∗τ20∗2iω20∗·E1=2K11K21K31,a110a130a22a23a31a320·E2=-K12K22K32,where(80)K11=-α1q3;K21=-α2q2q3;K31=β1α1q3e-iω20∗τ1+β2α2q2q3e-iω20∗τ20∗;K12=-α1q3+q-3;K22=-α2q2q-3+q-2q3;K32=β1α1q3eiω20∗τ1+q-3e-iω20∗τ1+β2α2q-2q3eiω20∗τ20∗+q2q-3e-iω20∗τ20∗.Then, we can compute the following values:(81)c10=i2ω20∗τ20∗g20g11-2g112-g0223+g212,μ=-Re⁡c10Re⁡λ′τ20∗,β=2Re⁡c10,T=-Imc10+μ2Imλ′τ20∗ω20∗τ20∗.This determine the properties of bifurcating periodic solutions and the Hopf bifurcation at τ=τ20∗. That is, (i)
μ determines the direction of the Hopf bifurcation. Specifically, when μ>0(<0), the Hopf bifurcation is supercritical (subcritical).(ii)
β determines the stability of the bifurcating periodic solutions; when β<0>0, the bifurcating periodic solution is stable (unstable).(iii)
T determines the period of the bifurcating periodic solutions; when T>0<0, the period of bifurcating periodic solution increases (decrease).


## 5. Numerical Simulations

Due to the complexity of model (1), we perform some numerical simulations in this section to investigate further how the delay influences dynamics in model (1). The following parameter values are used I1=I2=0.5, q1=q2=0.001, α1=α2=0.08, β1=0.3, and m=0.8. Other parameters are chosen as control parameters.According to the standard linear analysis, whenτ2 is equal to zero, the analysis reveals that the β2-τ1 parameter plane is divided into four parts (see Figure 1(a)). In Figure 1(a), before β2 reaching black dashed line, there exists τ10 in model (1) such that the unique positive equilibrium loses its stability when the condition, τ1>τ10, holds. When the locus of β2 is between black dashed line and green dashed line, the stability switches for positive equilibrium do not exist although (28) has two positive roots, which means that there exists τ10 in model (1) such that the unique positive equilibrium loses its stability when the condition, τ1>τ10, holds. However, the stability switches for positive equilibrium emerge when β2 is beyond green dashed line but it does not reach blue zone. When the locus of β2 is in blue zone, the unique positive equilibrium is always stable, which suggests that τ1 cannot influence the stability of the positive equilibrium. When τ1 equals zero, the similar results for β2-τ2 parameter plane are shown in Figure 1(b), but the sequence is reversed. Additionally, according to results in Section 4, we calculate the values of μ, β, and T at τ1=τ10 with β2∈(0,0.3), and the corresponding results are shown in Figure 1(c), where we can find that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable; especially, the period of the bifurcating periodic solutions increases as β2 increases. For other cases of τ1 and τ2, the same procedures with respect to calculations of μ, β, and T can be performed like Figure 1(c).Figure 1
(a) Bifurcation diagram withτ2=0 corresponding to τ1 V.S. β2, where solid, dashed, dash-dot, and dotted curves represent the critical values of τ1 in (30) for j=0,1,2,3, respectively, and the green dashed line denotes β2=0.417238. (b) Bifurcation diagram with τ1=0 corresponding to τ2 V.S. β2, where solid, dashed, dash-dot, and dotted curves represent the critical values of τ2 in (39) for j=0,1,2,3, respectively, and the green dashed line denotes β2=0.21. (c) Examples for μ, β, and T at τ1=τ10 with respect to β2, where τ2=0.
(a) (b) (c)As examples corresponding to stability of the positive equilibrium withβ2=0.2, taken τ1=1 and τ1=3 in Figure 1(a), respectively, the corresponding numerical solutions are shown in Figure 2. Obviously, the positive equilibrium is stable because τ1=1 is below τ10 (see Figure 2(a)). In contrast, due to 3=τ1 beyond τ10, a periodic solution exists (see Figure 2(b)). Furthermore, set β2=0.7, then we have τ1a0≈4.9050<τ1b0≈19.0615<τ1a1≈38.6683. Taken τ1=4, τ1=18 and τ1=21 in Figure 1(a), respectively, the corresponding numerical solutions are shown in Figure 3. Obviously, the positive equilibrium is stable when τ1=4 and τ1=21 (see Figures 3(a) and 3(c)), but the positive equilibrium is unstable when τ1=18 (see Figure 3(b)), which means that the positive equilibrium can gain its stability again for τ1>τ10. In Figure 3, the same initial values are applied, and other parameter values except for τ1 are also identical. Clearly, the delay is the principal factor giving rise to the difference among (a), (b), and (c).Figure 2
Numerical solutions of model (1) with τ2=0 and β2=0.2, (a) τ1=1; (b) τ1=3.
(a) (b)Figure 3
Numerical solutions of model (1) with τ2=0 and β2=0.7, (a) τ1=4; (b) τ1=18 (c) τ1=21.
(a) (b) (c)Numerical solutions in Figure3 suggest that the stability switches induced by delay may exist. Hence, the bifurcation diagram in τ1-τ2 parameter plane is given (see Figure 4(a)). For case τ2=0, there exists a τ1∗ such that the positive equilibrium with respect to τ1∈(0,τ1∗) is stable. For case τ1=0, there exists a τ2∗ such that the positive equilibrium with respect to τ2∈(0,τ2∗) is stable. Additionally, when τ2∈(0,τ2∗), Figure 4(a) shows that τ10 exists such that the positive equilibrium with respect to τ1∈(0,τ10) is stable. Likewise, when τ1∈(0,τ1∗), Figure 4(a) also display that τ20 exists such that the positive equilibrium with respect to τ2∈(0,τ20) is stable. Significantly, Figure 4(a) demonstrates that the stability switches for positive equilibrium with respect to τ1 emerge when τ2 below τ2∗ is fixed.Figure 4
(a) Bifurcation diagram withβ2=0.7 for τ1 V.S. τ2, where the symbol “S” denotes stable and the symbol “US” denotes unstable; (b) bifurcation diagram with β2=0.7 for the effect of τ2 on nutrient-phytoplankton dynamics, where dashed line denotes unstable; solid line denotes stable; the green solid square is Hopf bifurcation point, dot-dashed line corresponding to τ2∗ in (a), blue line represents equilibrium, and red line represents amplitude of periodic solutions. (c) Bifurcation diagram with β2=0.7 and τ2=1, where the yellow solid square is Hopf bifurcation point, the blue solid circle is bifurcation point for periodic-2 solution, the magenta zone indicates the existence of periodic-2 solutions, and the cyan solid diamond denotes the value of τ1 for a phase and a time-series in the inner of (c). (d) A periodic-3 solution for phytoplankton population, where β2=0.7, τ1=100, and τ2=2.
(a) (b) (c) (d)Figure4(b) depicts the dependence of stability of positive equilibrium on delay τ2 in the fully nonlinear regime for τ1 in sequence [0,2,10,30] and τ1=τ2, which is consistent with results in Figure 4(a). However, when τ2 is fixed, Figure 4(c) shows that the stability switches emerge with τ1 increases. Especially, periodic-2 solutions exist for some values of τ2 (see magenta zone in Figure 4(c)). As an example of periodic-2 solution existence, taking τ1=41, a phase and a time-series are given in the inner of Figure 4(c). Moreover, Figure 4(d) shows that there exist periodic-3 solutions in model (1). According to results in Section 2, the positive equilibrium is globally asymptotically stable in the model (1) without delay if it exists. Obviously, the results shown in Figure 4 are induced by delay.In Figure4(a), we can find that the number of intervals corresponding to stability of positive equilibrium for small values of τ2 equals 3. However, when τ2 is beyond dashed line (τ2∗∗), the number of intervals is 2. So the number of intervals for stability switches may be different for diverse τ2. Accordingly, we calculate the number of intervals with respect to parameter β2, as shown in Figure 5(a). Figure 5(b) shows that there exist 3 stable intervals when β2=0.7 and τ2=0, which is an example of Figure 5(a). Evidently, parameter β2 can remarkably influence the number of intervals for stability of positive equilibrium.Figure 5
(a) The number of intervals for stability switches of positive equilibrium with respect toβ2, where τ2=0; (b) based on (a), an example for comparison between results predicted by linear analysis and numerical results in the fully nonlinear regime with (β2,τ2) at (0,0.7), where top panel represents numerical results for model (1) in the fully nonlinear regime; bottom panel depicts the results predicted by linear analysis, and blue and red area denote stable and unstable, respectively.
(a) (b)


## 6. Conclusions

In this paper we proposed a nitrogen-phosphorus-phytoplankton model with multiple delays. The analysis focused on the effect of delay on nutrient-phytoplankton dynamics. In the absence of delay, theoretical analysis indicated that the unique positive equilibrium is globally asymptotically stable in model (1) if it exists. Deng et al. [48] also studied a nitrogen-phosphorus-phytoplankton model without delay, where Holling II function was employed to describe the nutrient uptake dynamics of phytoplankton. Although the function modelling the nutrient uptake dynamics of phytoplankton is different, they get the same results. These results mean that the nutrient-phytoplankton ecosystem will approach the stable equilibrium. However, it has been reported [18] that a constant population density may not exist because of the existence of some factors including noise, interval factors, and physical factors. And ecological studies [49, 50] also criticized this idea of “the balance of nature.” Actually, the existence of nutrient-plankton oscillations has been detected by laboratory experiments and field observation [19, 20]. Additionally, Benincà et al. [49] present the first experimental demonstration of chaos in a long-term experiment with a complex food web, where the food web was consisted of bacteria, several phytoplankton species, herbivorous and predatory zooplankton species, and detritivores. And they also find that the community moved back and forth between stabilizing and chaotic dynamics during the cyclic succession, and their findings provide a field demonstration of nonequilibrium coexistence of competing species through a cyclic succession at the edge of chaos [50]. These reports support that the nonequilibrium dynamics, such as oscillations and chaos, can exist in reality.In the present paper, we find that the unique positive equilibrium may lose its stability via Hopf bifurcation when delay appears, and then a periodic solution emerges, which means that nutrient-phytoplankton oscillation occurs. Obviously, the factor giving rise to nutrient- phytoplankton oscillation is delay in our studies. And the period and the stability of the bifurcating periodic solutions with respect to delay are discussed by using center manifold argument and normal form theory. In fact, instability induced by delay in nutrient-plankton model has been studied widely, and many studies indicate that the equilibrium is always unstable when delay is beyond a critical value [25, 29, 51, 52]. Yet, it should be emphasized in the present paper that the stability switches induced by delay can occur under some conditions.Moreover, numerical simulations showed how the delay influences nutrient-phytoplankton dynamics. Numerical results for model (1) in the fully nonlinear regime are consistent with the linear analysis. In numerical simulations, we found that delay indeed gives rise to the emergence of stability switches for the positive equilibrium. Yet, the numerical results show that the parameter intervals for stability switches may depend on other parameters as well, e.g., β2. Additionally, numerical results also indicated that periodic-2 solutions and periodic-3 solutions can emerge under some conditions for delay, which means that complex dynamics induced by delay exist in model (1). From biological viewpoint, the existence of periodic solutions implies that the fluctuations exist in density of phytoplankton population; that is, nutrient-phytoplankton oscillation emerges. Especially, by Li and York’s theory, periodic-3 solution implies chaos, which means that chaotic density fluctuations can display a variety of different periodicities and the long-term prediction of phytoplankton density can be fundamentally impossible. The chaotic density fluctuations do not contribute to the control of phytoplankton bloom. Consequently, the importance of the present paper is not the precision with which it predicts specific events for phytoplankton blooms but its contribution to the studies on how the delay influences nutrient-phytoplankton dynamics.



---

*Source: 3879626-2019-02-03.xml*