# Nonlinear Dynamics of a Nutrient-Phytoplankton Model with Time Delay

Authors: DeBing Mei; Min Zhao; Hengguo Yu; Chuanjun Dai; Yi Wang
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

We consider a nutrient-phytoplankton model with a Holling type II functional response and a time delay. By selecting the time delay used as a bifurcation parameter, we prove that the system is stable if the delay value is lower than the critical value but unstable when it is above this value. First, we investigate the existence and stability of the equilibria, as well as the existence of Hopf bifurcations. Second, we consider the direction, stability, and period of the periodic solutions from the steady state based on the normal form and the center manifold theory, thereby deriving explicit formulas. Finally, some numerical simulations are given to illustrate the main theoretical results.


---

## Body

## 1. Introduction

Ecological systems are characterized by the relationships between species and their natural environment. One of the key factors that affect population dynamics is predator. Due to its universal existence and importance in nature, the dynamic interactions between predators and their prey have remained one of the dominant themes in ecological population dynamics since the origin of this discipline [1–4].Phytoplankton has vital effects in aquatic ecosystems where it plays a significant role as the base of the food chain. Phytoplankton controls the global carbon cycle which has an important effect on climate regulation [5]. In some situations, marine waters, reservoir, and lakes may exhibit plankton or algal blooms [6, 7], where drinking water may be polluted and algal toxins can affect human health through the food chain, while the marine biodiversity may also be decreased.Nutrient-phytoplankton systems have a very important effect on aquatic ecosystems as part of predator-prey systems. Recently, many researchers have investigated a predator-prey model that involves nutrients, phytoplankton, and zooplankton [8–11], and various studies have simulated the dynamics of plankton systems.Huppert et al. [12] proposed a generic bottom-up nutrient-phytoplankton model:(1)N˙=a-bβtNP-eN,P˙=cβtNP-dP,where the forced NP compartment model comprises two variables: the nutrient concentration N and the phytoplankton biomass P. The initial conditions for this model are N(0)=N0>0 and P(0)=P0>0. a is an external source of nutrients that flows into the system from the environment, b and c are uptake rates, and e is the loss rate of nutrients. Seasonal environment conditions such as the water temperature, salinity, light, and thermocline depth are often important to different degrees depending on the system under investigation. The modulation of phytoplankton growth β can be represented as a periodic function β(t)=β(t+τ), where τ is the period of forcing, which is assumed to be annual in the present study. Analyses and numerical simulations can obtain insights into the mechanism of bloom recurrence, where modifications to the equations via the inclusion of appropriate functional forms can generate more realistic dynamics.Mäler [13] considered a better formulation of the dynamics of the eutrophication process:(2)dxdt=a-bx+x21+x2,where a is the runoff of nutrients into the lake and b is the constant rate of the natural removal of nutrients. After reducing the levels of green plants, the bottom sediments will be more vulnerable to winds, waves, and bottom-feeding fishes. Finally, nutrient sediments are released into the water and they contribute to the further growth of phytoplankton. The term x2/1+x2 represents the feedback from the bottom sediments. Depending on the numerical values for the parameters and specific simulations of the system, different negative externalities can affect the system in various ways.In this study, based on the idea of Huppert et al. [12] and Mäler [13], we consider a nutrient-phytoplankton model with delay:(3)dxdt=I-bx+x21+x2-axy,dydt=akxt-τy-my-fyyHy+y,where two variables x and y represent the concentrations of nutrients and phytoplankton, respectively. All of the parameters are positive constants, where τ is a positive delay, that is, the time required to convert nutrients into phytoplankton, I is a constant concentration, that is, the nutrient runoff from the environment, b is the loss rate of the nutrient concentration, a is the predation rate of nutrients by phytoplankton, k denotes the maximum nutrient intake rate of phytoplankton, m is the mortality of phytoplankton, fy represents the maximum predation rate of zooplankton on phytoplankton, and Hy is the half-saturation constant for phytoplankton, and fyy/Hy+yrepresents the predation of zooplankton on phytoplankton. The initial conditions are as follows: x(θ)=ϕ(θ)≥0, y(θ)=ψ(θ)≥0, θ∈-τ,0, ϕ(0)>0, and ψ(0)>0, where ϕ(θ), ψ(θ) are continuous bounded functions in the interval -τ,0.Huppert et al. [12] described the dynamics of the system above without considering the effect of delay. Time delay plays an important role in reflecting the real dynamical behaviors of biological systems and many studies have addressed this issue in recent years [14–23].In this study, in order to investigate the effects of a time delay on the system, we selected the delayτ as a bifurcation parameter. The remainder of this paper is organized as follows. In Section 2, we prove that positive equilibria exist under certain conditions and we then analyze the local stability of the boundary equilibrium. We also analyze the stability of the positive equilibria and the critical conditions when the Hopf bifurcation occurs. In Section 3, we derive an explicit algorithm and sufficient conditions for determining the direction of the Hopf bifurcation and the stability of the periodic solutions. In Section 4, we present some numerical simulations to illustrate our theoretical results. In Section 5, we give some brief conclusions.


## 2. Stability of the Positive Equilibria and Existence of the Hopf Bifurcation

### 2.1. The Existence of Positive Equilibria

In this section, we consider the existence of the positive equilibria for system (3). We define f(x,y)=I-bx+x2/1+x2-axy, g(x,y)=akxy-my-fyy/Hy+y. The equations for the equilibria can be depicted as follows:(4)fx,y=0,gx,y=0.According tof(x,y)=0, a vertical isocline can be obtained, l1: y=1/ax(I-bx+x2/1+x2), which can be written as F(x)=1/ax(I-bx+x2/1+x2). For l1, when the value of x is close to positive infinite, the value of F(x) is close to -b/a. By contrast, when the value of x approaches zero from the right-hand side, the value of F(x) is approximately positive infinite. Due to the vertical isocline, l1 is continuous and derivable when x∈0,+∞ and x∈-∞,0; thus we can compute the derivative of F(x); that is, F′(x)=-(I+1)x4+(2I-1)x2+I/(1+x2)2ax2, and four roots exist when F′(x)=0, as follows:(5)x1=-2I+2-1+2I-1-8I,x2=--2I+2-1+2I-1-8I,x3=-2I+2-1+2I+1-8I,x4=--2I+2-1+2I+1-8I.From the analysis above, (i) if I≥1/8, we can obtain F′(x)≤0 for an arbitrary real number x; thus F(x) decreases monotonically when x∈0,+∞; (ii) if I<1/8, four real number roots exist when F′(x)=0. Obviously, x1>x3>0>x4>x2. According to the second derivative of F(x), we can find that x1, x4 are the points with the extreme maximum value and x2, x3 are the points with the extreme minimum value. We can find that F(x) decreases monotonically when x∈0,x3∪x1,+∞ and increases monotonically when x∈x3,x1.According tog(x,y)=0, the horizontal isocline l2: y=fy/akx-m-Hy and y=0, which can be written as G(x)=fy/akx-m-Hy and G(x)=0. When the horizontal isocline l2: G(x)≠0, that is, G(x)=fy/akx-m-Hy, it is easy to find that the line, x=m/ak, is the asymptote of G(x). Obviously, G(x) is continuous when x∈0,m/ak∪m/ak,+∞, and we can find that G(x) decreases monotonically when x∈0,m/ak∪m/ak,+∞. For G(x), G(x)<0 holds when x∈0,m/ak∪fy+mHy/akHy,+∞, whereas G(x)>0 holds when x∈m/ak,fy+mHy/akHy.First, when the system comprises only nutrients without phytoplankton in the system, according to the vertical isocline,l1: F(x)=1/ax(I-bx+x2/1+x2), and according to the horizontal isocline, l2: G(x)=0. From the analysis above, if I≥1/8, we find that there is always a boundary equilibrium. If I<1/8, one boundary equilibrium exists at least. The boundary equilibrium can be written as E0(x0,0), x0>0.In addition, according to the vertical isocline,l1: F(x)=1/ax(I-bx+x2/1+x2), and according to the horizontal isocline, l2: G(x)=fy/akx-m-Hy. From the discussion above, if I≥1/8, we find that a positive equilibrium exists in the system when the condition that fy+mHy/akHy<x0 holds. Thus, we can find that there is no positive equilibrium in the system when the condition that m/ak≥x0 holds. If I<1/8, one positive equilibrium exists at least when F(x1)≥G(x1)>0 holds. The positive equilibrium can be written as E∗(x∗,y∗).In the analysis above, we only prove that positive equilibria exist in the system when certain conditions are satisfied but we do not indicate that positive equilibria do not exist when the conditions are not satisfied.The Jacobian matrix of the system without time delay at the equilibriumE0(x0,0) is(6)T1=-b+2x01+x022-ax00aex0-m-fP1Hy.The index of equilibriumE0 is +1 when the condition 2x0<b(1+x02), aex0Hy<mHy+fy is satisfied, which is stable. According to the eigenvalues of T1, it is locally asymptotically stable when the above conditions are satisfied. In particular, the equilibrium E0(x0,0) is unstable, when 2x0>b(1+x02) or aex0Hy>mHy+fy.From the above discussion, we find a positive equilibrium in the system without time delay under some preconditions, defined byE∗(x∗,y∗). The Jacobian matrix of the system at the equilibrium E∗(x∗,y∗) is(7)T2=-b+2x∗1+x∗22-aPy∗-ax∗aey∗fyy∗Hy+y∗2.When2y∗(Hy+y∗)2+fyy∗(1+x∗2)2<(b+ay∗)(1+x∗2)2(Hy+y∗)2, the index of equilibrium N2 is +1 if a2ex∗y∗(Hy+y∗)2(1+x∗2)2>fyy∗(2x∗-b(1+x∗2)2-ax∗(1+x∗2)2). It is locally asymptotically stable using the Routh-Hurwitz criteria when the above conditions are satisfied. In particular, the equilibrium E∗(x∗,y∗) is unstable, when a2ex∗y∗(Hy+y∗)2(1+x∗2)2<fyy∗(2x∗-b(1+x∗2)2-ax∗(1+x∗2)2), which is a saddle.


### 2.2. Local Stability of EquilibriumE0

The linearization of (3) at the equilibrium E0 is(8)dxdt=-b+2x01+x02x-ax0y,dydt=akx0-m-fyHyy.Therefore, the characteristic matrix at equilibrium E0(x0,0) is(9)T0=λ--b+2x01+x02ax00λ-akx0-m-fyHy.Its characteristic equation is(10)λ2-b11+b22λ+b11b22=0,where b11=-b+2x0/1+(x0)2, b12=-ax0, b21=0, and b22=akx0-m-fy/Hy.Theorem 1.
If-b+2x0/1+(x0)2<0 and akx0-m-fy/Hy<0 hold, the boundary equilibrium E0 is locally asymptotically stable.Proof.
There are two eigenvalues where-b+2x0/1+(x0)2 and akx0-m-fy/Hy. When all the roots of the characteristic equation at E0 are negative, the boundary equilibrium E0 is a stable node. Thus, the boundary equilibrium E0 is locally asymptotically stable. The proof is complete.


### 2.3. Local Stability and the Hopf Bifurcation of EquilibriumE∗

In the following, we focus on the existence of a local Hopf bifurcation at positive equilibriumE∗(x∗,y∗). If we let P(x)=x2/1+x2, the linearization of (3) at the equilibrium E∗ is as follows:(11)dxdt=-b+P′x∗-ay∗x-ax∗y,dydt=aky∗xt-τ+fyy∗Hy+y∗2y.Thus, the characteristic matrix at equilibrium E∗(x∗,y∗) is(12)J∗=λ--b+P′x∗-ay∗ax∗-aky∗e-λτλ-fyy∗Hy+y∗2.Set a11=-b+2x∗/1+(x∗)2-ay∗,  a12=-ax∗,  a21=aky∗, and a22=fyy∗/(Hy+y∗)2.Its characteristic equation is(13)λ2-a11+a22λ+a11a22-a12a21e-λτ=0or(14)λ2+Aλ+Be-λτ+B0=0,where A=-(a11+a22), B=-a12a21, and B0=a11a22.Theorem 2.
For system (3), the positive equilibrium E∗(x∗,y∗) is locally asymptotically stable for all 0≤τ<τ0 and the system undergoes a Hopf bifurcation at the positive equilibrium E∗(x∗,y∗) for τ=τk, k=0,1,2,…,  when the conditions that -A<0 and B+B0>0 hold.Lemma 3.
Ifτ=0, -A<0, and B+B0>0 hold, the positive equilibrium E∗(x∗,y∗) is locally asymptotically stable.Proof.
Whenτ=0, the characteristic equation is(15)λ2+Aλ+B+B0=0.If -A<0 and B+B0>0, then all the roots of characteristic equation (15) have negative real parts. Thus, the positive equilibrium E∗(x∗,y∗) is locally asymptotically stable. The proof is complete.Lemma 4.
Let(16)w0=-A2-2B02+Δ2,-A<0,B+B0>0.If(17)τ=τk=1w0arccos⁡w02-B0B+2kπ,k=0,1,2,…,then (13) has a pair of purely imaginary roots ±iw0.Proof.
Whenτ≠0, we assume that λ=μ(τ)+iw(τ) is a root of (13). By substituting this into (13), we obtain(18)μ2τ-w2τ+2μτwτi+μτA+Awτi+Be-μττcos⁡τwτ-isin⁡τwτ+B0=0or(19)μ2τ-w2τ+μτA+B0=-Be-μττcos⁡τwτ,2μτwτ+Awτ=Be-μττsin⁡τwτ.If we let μ(τ)=0 and w(τ)=w0>0, (19) is equal to(20)w02-B0=Bcos⁡τw0,Aw0=Bsin⁡τw0and thus(21)w04+A2-2B0w02+B02-B2=0.If we let -A<0, B+B0>0, we can obtain 4((B0)2-B2)<0; therefore, Δ=(A2-2B0)2-4((B0)2-B2)>0, and then w0=-(A2-2B0)+Δ/2.
However, we can solve the first formula in (20) to obtain(22)τ=τk=1w0arccos⁡w02-B0B+2kπ,k=0,1,2,….Thus, when τ=τk, (13) has a pair of purely imaginary roots ±iw0. The proof is complete.Lemma 5.
For arbitraryτ=τk,  k=0,1,2,…, that satisfy the above conditions, one can find that the derivative of the real number part of the eigenvalue is greater than zero; that is,(23)dReλτdττ=τk>0,k=0,1,2,….We can find that the roots of (13) cross the imaginary axis in a transverse direction from left to right as τ increases.Proof.
Differentiating characteristic equation (13) with respect to τ, we obtain(24)2λτdλτdτ+Adλτdτ+Be-λττ-τdλτdτ-λτ=0.Thus, we have(25)dλτdτ=Bλτe-λττ2λτ+A-Bτe-λττ.From (20) and (21), we obtain(26)dλτdττ=τk=iBw0e-iw0τk2iw0+A-Bτke-iw0τk=iBw0cos⁡w0τk-isin⁡w0τk2iw0+A-Bτkcos⁡w0τk-isin⁡w0τk=Bw0sin⁡w0τk+iBw0cos⁡w0τkA-Bτkcos⁡w0τk+i2w0+Bτksin⁡w0τk=Aw02+iw0w02-B0A-τkw02-B0+i2w0+τkAw0,Redλτdττ=τk=Aw02A-τkw02-B0+w0w02-B02w0+τkAw0A-τkw02-B02+2w0+τkAw02=A2-2B0w02+2w04A-τkw02-B02+2w0+τkAw02=a112+a222w02+2w04A-τkw02-B02+2w0+τkAw02>0because of(27)dReλτdττ=τk=Redλτdττ=τkand thus(28)dReλτdττ=τk=a112+a222w02+2w04A-τkw02-B02+2w0+τkAw02>0.This shows that all of the roots cross the imaginary axis at iw from left to right as τ increases. The proof is complete.From Lemmas3–5 and the Hopf bifurcation theorem, the proof is complete.


## 2.1. The Existence of Positive Equilibria

In this section, we consider the existence of the positive equilibria for system (3). We define f(x,y)=I-bx+x2/1+x2-axy, g(x,y)=akxy-my-fyy/Hy+y. The equations for the equilibria can be depicted as follows:(4)fx,y=0,gx,y=0.According tof(x,y)=0, a vertical isocline can be obtained, l1: y=1/ax(I-bx+x2/1+x2), which can be written as F(x)=1/ax(I-bx+x2/1+x2). For l1, when the value of x is close to positive infinite, the value of F(x) is close to -b/a. By contrast, when the value of x approaches zero from the right-hand side, the value of F(x) is approximately positive infinite. Due to the vertical isocline, l1 is continuous and derivable when x∈0,+∞ and x∈-∞,0; thus we can compute the derivative of F(x); that is, F′(x)=-(I+1)x4+(2I-1)x2+I/(1+x2)2ax2, and four roots exist when F′(x)=0, as follows:(5)x1=-2I+2-1+2I-1-8I,x2=--2I+2-1+2I-1-8I,x3=-2I+2-1+2I+1-8I,x4=--2I+2-1+2I+1-8I.From the analysis above, (i) if I≥1/8, we can obtain F′(x)≤0 for an arbitrary real number x; thus F(x) decreases monotonically when x∈0,+∞; (ii) if I<1/8, four real number roots exist when F′(x)=0. Obviously, x1>x3>0>x4>x2. According to the second derivative of F(x), we can find that x1, x4 are the points with the extreme maximum value and x2, x3 are the points with the extreme minimum value. We can find that F(x) decreases monotonically when x∈0,x3∪x1,+∞ and increases monotonically when x∈x3,x1.According tog(x,y)=0, the horizontal isocline l2: y=fy/akx-m-Hy and y=0, which can be written as G(x)=fy/akx-m-Hy and G(x)=0. When the horizontal isocline l2: G(x)≠0, that is, G(x)=fy/akx-m-Hy, it is easy to find that the line, x=m/ak, is the asymptote of G(x). Obviously, G(x) is continuous when x∈0,m/ak∪m/ak,+∞, and we can find that G(x) decreases monotonically when x∈0,m/ak∪m/ak,+∞. For G(x), G(x)<0 holds when x∈0,m/ak∪fy+mHy/akHy,+∞, whereas G(x)>0 holds when x∈m/ak,fy+mHy/akHy.First, when the system comprises only nutrients without phytoplankton in the system, according to the vertical isocline,l1: F(x)=1/ax(I-bx+x2/1+x2), and according to the horizontal isocline, l2: G(x)=0. From the analysis above, if I≥1/8, we find that there is always a boundary equilibrium. If I<1/8, one boundary equilibrium exists at least. The boundary equilibrium can be written as E0(x0,0), x0>0.In addition, according to the vertical isocline,l1: F(x)=1/ax(I-bx+x2/1+x2), and according to the horizontal isocline, l2: G(x)=fy/akx-m-Hy. From the discussion above, if I≥1/8, we find that a positive equilibrium exists in the system when the condition that fy+mHy/akHy<x0 holds. Thus, we can find that there is no positive equilibrium in the system when the condition that m/ak≥x0 holds. If I<1/8, one positive equilibrium exists at least when F(x1)≥G(x1)>0 holds. The positive equilibrium can be written as E∗(x∗,y∗).In the analysis above, we only prove that positive equilibria exist in the system when certain conditions are satisfied but we do not indicate that positive equilibria do not exist when the conditions are not satisfied.The Jacobian matrix of the system without time delay at the equilibriumE0(x0,0) is(6)T1=-b+2x01+x022-ax00aex0-m-fP1Hy.The index of equilibriumE0 is +1 when the condition 2x0<b(1+x02), aex0Hy<mHy+fy is satisfied, which is stable. According to the eigenvalues of T1, it is locally asymptotically stable when the above conditions are satisfied. In particular, the equilibrium E0(x0,0) is unstable, when 2x0>b(1+x02) or aex0Hy>mHy+fy.From the above discussion, we find a positive equilibrium in the system without time delay under some preconditions, defined byE∗(x∗,y∗). The Jacobian matrix of the system at the equilibrium E∗(x∗,y∗) is(7)T2=-b+2x∗1+x∗22-aPy∗-ax∗aey∗fyy∗Hy+y∗2.When2y∗(Hy+y∗)2+fyy∗(1+x∗2)2<(b+ay∗)(1+x∗2)2(Hy+y∗)2, the index of equilibrium N2 is +1 if a2ex∗y∗(Hy+y∗)2(1+x∗2)2>fyy∗(2x∗-b(1+x∗2)2-ax∗(1+x∗2)2). It is locally asymptotically stable using the Routh-Hurwitz criteria when the above conditions are satisfied. In particular, the equilibrium E∗(x∗,y∗) is unstable, when a2ex∗y∗(Hy+y∗)2(1+x∗2)2<fyy∗(2x∗-b(1+x∗2)2-ax∗(1+x∗2)2), which is a saddle.


## 2.2. Local Stability of EquilibriumE0

The linearization of (3) at the equilibrium E0 is(8)dxdt=-b+2x01+x02x-ax0y,dydt=akx0-m-fyHyy.Therefore, the characteristic matrix at equilibrium E0(x0,0) is(9)T0=λ--b+2x01+x02ax00λ-akx0-m-fyHy.Its characteristic equation is(10)λ2-b11+b22λ+b11b22=0,where b11=-b+2x0/1+(x0)2, b12=-ax0, b21=0, and b22=akx0-m-fy/Hy.Theorem 1.
If-b+2x0/1+(x0)2<0 and akx0-m-fy/Hy<0 hold, the boundary equilibrium E0 is locally asymptotically stable.Proof.
There are two eigenvalues where-b+2x0/1+(x0)2 and akx0-m-fy/Hy. When all the roots of the characteristic equation at E0 are negative, the boundary equilibrium E0 is a stable node. Thus, the boundary equilibrium E0 is locally asymptotically stable. The proof is complete.


## 2.3. Local Stability and the Hopf Bifurcation of EquilibriumE∗

In the following, we focus on the existence of a local Hopf bifurcation at positive equilibriumE∗(x∗,y∗). If we let P(x)=x2/1+x2, the linearization of (3) at the equilibrium E∗ is as follows:(11)dxdt=-b+P′x∗-ay∗x-ax∗y,dydt=aky∗xt-τ+fyy∗Hy+y∗2y.Thus, the characteristic matrix at equilibrium E∗(x∗,y∗) is(12)J∗=λ--b+P′x∗-ay∗ax∗-aky∗e-λτλ-fyy∗Hy+y∗2.Set a11=-b+2x∗/1+(x∗)2-ay∗,  a12=-ax∗,  a21=aky∗, and a22=fyy∗/(Hy+y∗)2.Its characteristic equation is(13)λ2-a11+a22λ+a11a22-a12a21e-λτ=0or(14)λ2+Aλ+Be-λτ+B0=0,where A=-(a11+a22), B=-a12a21, and B0=a11a22.Theorem 2.
For system (3), the positive equilibrium E∗(x∗,y∗) is locally asymptotically stable for all 0≤τ<τ0 and the system undergoes a Hopf bifurcation at the positive equilibrium E∗(x∗,y∗) for τ=τk, k=0,1,2,…,  when the conditions that -A<0 and B+B0>0 hold.Lemma 3.
Ifτ=0, -A<0, and B+B0>0 hold, the positive equilibrium E∗(x∗,y∗) is locally asymptotically stable.Proof.
Whenτ=0, the characteristic equation is(15)λ2+Aλ+B+B0=0.If -A<0 and B+B0>0, then all the roots of characteristic equation (15) have negative real parts. Thus, the positive equilibrium E∗(x∗,y∗) is locally asymptotically stable. The proof is complete.Lemma 4.
Let(16)w0=-A2-2B02+Δ2,-A<0,B+B0>0.If(17)τ=τk=1w0arccos⁡w02-B0B+2kπ,k=0,1,2,…,then (13) has a pair of purely imaginary roots ±iw0.Proof.
Whenτ≠0, we assume that λ=μ(τ)+iw(τ) is a root of (13). By substituting this into (13), we obtain(18)μ2τ-w2τ+2μτwτi+μτA+Awτi+Be-μττcos⁡τwτ-isin⁡τwτ+B0=0or(19)μ2τ-w2τ+μτA+B0=-Be-μττcos⁡τwτ,2μτwτ+Awτ=Be-μττsin⁡τwτ.If we let μ(τ)=0 and w(τ)=w0>0, (19) is equal to(20)w02-B0=Bcos⁡τw0,Aw0=Bsin⁡τw0and thus(21)w04+A2-2B0w02+B02-B2=0.If we let -A<0, B+B0>0, we can obtain 4((B0)2-B2)<0; therefore, Δ=(A2-2B0)2-4((B0)2-B2)>0, and then w0=-(A2-2B0)+Δ/2.
However, we can solve the first formula in (20) to obtain(22)τ=τk=1w0arccos⁡w02-B0B+2kπ,k=0,1,2,….Thus, when τ=τk, (13) has a pair of purely imaginary roots ±iw0. The proof is complete.Lemma 5.
For arbitraryτ=τk,  k=0,1,2,…, that satisfy the above conditions, one can find that the derivative of the real number part of the eigenvalue is greater than zero; that is,(23)dReλτdττ=τk>0,k=0,1,2,….We can find that the roots of (13) cross the imaginary axis in a transverse direction from left to right as τ increases.Proof.
Differentiating characteristic equation (13) with respect to τ, we obtain(24)2λτdλτdτ+Adλτdτ+Be-λττ-τdλτdτ-λτ=0.Thus, we have(25)dλτdτ=Bλτe-λττ2λτ+A-Bτe-λττ.From (20) and (21), we obtain(26)dλτdττ=τk=iBw0e-iw0τk2iw0+A-Bτke-iw0τk=iBw0cos⁡w0τk-isin⁡w0τk2iw0+A-Bτkcos⁡w0τk-isin⁡w0τk=Bw0sin⁡w0τk+iBw0cos⁡w0τkA-Bτkcos⁡w0τk+i2w0+Bτksin⁡w0τk=Aw02+iw0w02-B0A-τkw02-B0+i2w0+τkAw0,Redλτdττ=τk=Aw02A-τkw02-B0+w0w02-B02w0+τkAw0A-τkw02-B02+2w0+τkAw02=A2-2B0w02+2w04A-τkw02-B02+2w0+τkAw02=a112+a222w02+2w04A-τkw02-B02+2w0+τkAw02>0because of(27)dReλτdττ=τk=Redλτdττ=τkand thus(28)dReλτdττ=τk=a112+a222w02+2w04A-τkw02-B02+2w0+τkAw02>0.This shows that all of the roots cross the imaginary axis at iw from left to right as τ increases. The proof is complete.From Lemmas3–5 and the Hopf bifurcation theorem, the proof is complete.


## 3. Direction and Stability of the Hopf Bifurcations

In the previous section, we obtained some conditions for the occurrence of Hopf bifurcations. In this section, we consider the direction, stability, and period of the periodic solutions from the steady state, as well as deriving the explicit formulae that determine these factors at the critical valueτ=τk, k=0,1,2,…, using the method introduced by Hassard et al. [24].If we letu1(t)=x(τt)-x∗,  u2(t)=y(τt)-y∗,  τ=μ+τk, and Q(x)=fyy/Hy+y, system (3) can be written as(29)u′t=Lμut+fμut,where u=(u1,u2)T, ut(θ)=u(t+θ), θ∈[-1,0], for ϕ=(ϕ1,ϕ2)∈C([-1,0],R2), and(30)Lμϕ=μ+τk-b+P′x∗-ay∗-ax∗0Q′y∗ϕ10ϕ20+00aky∗0ϕ1-1ϕ2-1,(31)fμϕ=μ+τk·-aϕ10ϕ20+P′′x∗2!ϕ120+P′′′x∗3!ϕ130akϕ1-1ϕ20-Q′′x∗2!ϕ220-Q′′′x∗3!ϕ130+H.O.T.,where H.O.T. denotes “higher order terms.” By the Riesz representation theorem, a function η(θ,μ) of bounded variation exists for θ∈[-1,0] such that(32)Lμϕ=∫-10dηθ,μϕθforϕ∈C-1,0,R2.In fact, we can choose(33)ηθ,μ=μ+τk-b+P′x∗-ay∗-ax∗0Q′y∗δθ-μ+τk00aky∗0δθ+1,where(34)δθ=1,θ=00,θ≠0.In the following discussion, we omit theH.O.T. term in (31) because determining the direction and stability of the Hopf bifurcation only requires up to the third term. For ϕ∈C1([-1,0],R2), we define(35)Aμϕ=dϕdθ,θ∈-1,0Lμϕ,θ=0,(36)Fμϕ=0,θ∈-1,0fμϕ,θ=0.Then, (29) can be rewritten as(37)u˙t=Aμut+Fμut.For ψ∈C1([0,1],(R2)∗), we define the adjoint operator A∗ as follows:(38)A∗ψs=-dψds,s∈0,1∫-10ψ-θdηθ,0,s=0.For ψ∈C1([0,1],(R2)∗) and ϕ∈C1([0,1],R2), we define the bilinear form:(39)ψ,ϕ=ψ-0ϕ0-∫-10∫ξ=0θψξ-θdηθϕξdξ,where η(θ)=η(θ,0). In Section 2, we found that iw0τk and -iw0τk are eigenvalues of A(0) and A∗. Thus, we need to calculate the eigenvectors of A(0) and A∗ that correspond to iw0τk and -iw0τk, respectively. Suppose that q(θ) and q∗(s) are the eigenvectors that correspond to iw0τk and -iw0τk, respectively. Let(40)qθ=1Ceiw0τkθ.Using the condition that A(0)q(θ)=iw0τkq(θ), that is,(41)τkiw0+b-P′x∗+ay∗ax∗-aky∗e-iw0τkiw0-Q′y∗1C=00,thus we can obtain(42)C=-iw0+b-P′x∗+ay∗ax∗.Similarly, suppose that q∗(s)=E(1,D)eiw0τks is the eigenvector of A∗ that corresponds to -iw0τk; then A∗q∗(s)=-iw0τkq∗(s), and from the definition of A∗, we have(43)τk-iw0+b-P′x∗+ay∗-aky∗eiw0τkax∗-iw0-Q′y∗1D=00.Therefore, we obtain(44)D=ax∗iw0+Q′y∗.Let q∗,q=1, q∗,q-=0 and due to(45)q∗,q=q∗¯0q0-∫-10∫ξ=0θq∗¯ξ-θdηθ·qξdξ=E-1D-1C-∫-10∫ξ=0θE-1D-e-iw0τkξ-θdηθ1Ceiw0τkξdξ=E-1+D-C-∫-10θ1D-dηθ1Ceiw0τkθ=E-1+D-C-1D-τk-00aky∗0·1Ce-iw0τkθ=E-1+D-C+τkD-aky∗e-iw0τkwe can obtain(46)E=11+DC-+τkDaky∗eiw0τk.In the following, we first construct the coordinates to describe the center manifoldC0 at μ=0 (i.e., τ=τk). We define(47)zt=q∗,ut,Wt,θ=utθ-2Reztqθ.On the center manifold C0, we have(48)Wt,θ=Wzt,z-t,θ,where(49)Wz,z-,θ=W20θz22+W11θzz-+W02θz-22+W30θz36+⋯z and z- are local coordinates for the center manifold C0 at μ=0. From (36), (37), (38), (39), and (47), we can obtain(50)z˙t=q∗,u˙t=iw0τkzt+q∗θ,F0Wz,z-,θ+2Reztqθ=iw0τkzt+q∗¯0f0Wz,z-,0+2Reztq0≜iw0τkzt+gz,z-,where(51)gz,z-=q∗¯0f0Wz,z-,0+2Reztq0=g20z22+g11zz-+g02z-22+g21z2z-2+⋯.Based on (37) and (47), we have(52)W˙=u˙t-2Rez˙q=u˙t-z˙q-z-˙q-=Aut+Fut-z˙q-z-˙q-=AW+2Rezq+FW+2Rezq-2Regz,z-qθ=AW-2Req∗¯0F0qθ,-1≤θ<0AW-2Req∗¯0F0qθ+F0,θ=0≜AW+Hz,z-,θ,where(53)Hz,z-,θ=H20θz22+H11θzz-+H02θz-22+⋯.Near the origin on the center manifold C0, by (49)–(51), we have(54)W˙z,z-=Wzz˙+Wz-z-˙=zW20θz˙+W11θz-z˙+z-W02θz-˙+W11θzz-˙+Oz,z-3=zW20θiw0τkz+gz,z-+W11θz-iw0τkz+gz,z-+z-W02θ-iw0τkz+g-z,z-+W11θz-iw0τkz+g-z,z-+Oz,z-3=iw0τkW20θz2-iw0τkW02θz-2+Oz,z-3.By (49), (52), (53), and (54) and by comparing the coefficients, we have(55)2iw0τkI-A0W20θ=H20θ,(56)-A0W11θ=H11θ,(57)A0+2iw0τkIW02θ=-H02θ.We consider(58)q0=1C,q∗0=E1Dand because of(59)utθ=u1tθ,u2tθT=Wt,θ+2Reztqθ=Wt,θ+zqθ+z-q-θwe obtain(60)utθ=Wz,z-,θ=z1Ceiw0τkθ+z-1C-e-iw0τkθ+W20θz22+W11θzz-+W02θz-22+W30θz36+⋯,u1t0=z+z-+W2010z22+W1110zz-+W0210z-22+⋯,u2t0=Cz+C-z-+W2020z22+W1120zz-+W0220z-22+⋯,u1t-1=e-iw0τkz+eiw0τkz-+W201-1z22+W111-1zz-+W021-1z-22+⋯,F0=τk-au1t0u2t0+P′′x∗2!u1t20+P′′′x∗3!u1t30P′′′x∗3!2aku1t-1u2t0-Q′′y∗2!u2t20-Q′′′y∗3!u2t30,gz,z-=q∗¯0f0Wz,z-,0+2Reztq0=E-τk1,D-·-au1t0u2t0+P′′x∗2!u1t20+P′′′x∗3!u1t30P′′′x∗3!2aku1t-1u2t0-Q′′y∗2!u2t20-Q′′′y∗3!u2t30=E-τk-az+z-+W2010z22+W1110zz-+W0210z-22+⋯·Cz+C-z-+W2020z22+W1120zz-+W0220z-22+⋯+P′′x∗2!z+z-+W2010z22+W1110zz-+W0210z-22+⋯2+P′′′x∗3!z+z-+W2010z22+W1110zz-+W0210z-22+⋯3+D-ake-iw0τkz+eiw0τkz-+W201-1z22+W111-1zz-+W021-1z-22+⋯·Cz+C-z-+W2020z22+W1120zz-+W0220z-22+⋯-D-·Q′′y∗2!Cz+C-z-+W2020z22+W1120zz-+W0220z-22+⋯2-D-·Q′′′y∗3!Cz+C-z-+W2020z22+W1120zz-+W0220z-22+⋯3.Thus, by comparing the coefficients with (51), we obtain(61)g20=2τkE--aC+P′′x∗2!+D-akCe-iw0τk-Q′′x∗2!C2,g11=τkE--aC+C-+P′′x∗+D-akCeiw0τk+C-e-iw0τk-Q′′y∗CC-,g02=2τkE--aC-+P′′x∗2!+D-akC-eiw0τk-P′′x∗2!C-2,g21=2τkE--aW1120+W20202+W20102C-+W1110C+P′′x∗2!W2010+2W1110+P′′′x∗3!×3+D-ake-iw0τkW1120+eiw0τkW20202+CW111-1+C-W201-12-D-Q′′y∗2!W2020C-+2W1120C-D-·Q′′′y∗3!×3C2C-.W20(θ) and W11(θ) are still unknown; thus, we need to compute them. For θ∈-1,0, we have(62)Hz,z-,θ=-2Req-∗0F0qθ=-gqθ-g-q-θ=-g20z22+g11zz-+g02z-22+⋯qθ-g-20z-22+g-11zz-+g-02z22+⋯q-θ.By comparing the coefficients, we have(63)H20θ=-g20qθ-g-02q-θ,H11θ=-g11qθ-g-11q-θ.It follows from (55) that(64)W˙20θ=2iw0τkW20θ+g20q0+g-02q-0=2iw0τkW20θ+g20q0eiw0τkθ+g-02q-0e-iw0τkθ.Solving for W20(θ), we obtain(65)W20θ=ig20q0w0τkeiw0τkθ+ig02q-03w0τke-iw0τkθ+E1e2iw0τkθ.Similarly, we obtain(66)W11θ=-ig11q0w0τkeiw0τkθ+ig-11q-0w0τke-iw0τkθ+E2,where E1 and E2 are both two-dimensional vectors, which can be determined by setting θ=0 in H. For θ=0, we have H(z,z-,0)=-2Re{q-∗(0)F(0)q(0)}+F(0); thus(67)H200=-g20q0-g-02q-0+τk-aC+P′′x∗2!akCe-iw0τk-Q′′x∗2!C2,H110=-g11q0-g-11q-0+τk-aC+C-+P′′x∗2!×2akCeiw0τk+C-e-iw0τk-P′′x∗2!-2CC-.According to (55) and (56), we have(68)τkL1W200+τkL2W20-1=2iw0τkW200-H200,(69)τkL1W110+τkL2W11-1=-H110.By substituting (68) and (69) into (65) and (66), respectively, we obtain(70)E1=1τk2iw0I-L1-L2e-2iw0τk-1H200+2g20q0+2g-02q-03+ig20L1q0w0+ig-02L1q-03w0+ig20e-iw0τkL2q0w0+ig-02eiw0τkL2q-03w0,E2=1τk-L1-L2-1H110-ig11L1q0w0+ig-11L1q-0w0-ig11e-iw0τkL2q0w0+ig-11eiw0τkL2q-0w0,where I is the 2×2 identity matrix. From the above, we already know W20(θ) and W11(θ), and g21 can be expressed in terms of the parameters and delay τ; thus, we can compute the values as follows:(71)c10=i2w0τkg20g11-2g112-g0223+g212,μ2=-Rec10Reλ′τk,β2=2Rec10,T2=-Im⁡c10+μ2Im⁡λ′τkw0τk.Theorem 6.
μ
2, β2, and T2 are defined above; thus the following are true:(i)
ifμ2>0  (μ2<0), then the Hopf bifurcation is supercritical (subcritical);(ii)
ifβ2<0  (β2>0), then the periodic solutions are stable (unstable);(iii)
ifT2>0  (T2<0), then the period of the bifurcating periodic solution of system (2) increases (decreases).


## 4. Numerical Simulation

According to the analysis above, the positive equilibrium does not always exist. Thus, based on the numerical technology, we can obtain the positive equilibria that exist under certain conditions. When some parameters are set, we selectb=0.3, a=0.5, k=0.75, m=0.2, fy=0.6, and Hy=4, and we can obtain the different results shown in Figure 1. Figure 1(a) shows that different nutrient concentrations correspond to various vertical isoclines and we can determine the different intersection points (i.e., equilibria). Figure 1(b) shows the stability with the change in the nutrient concentration when the other parameters are fixed, where the blue solid line shows that the positive equilibrium is unstable, whereas the red solid line shows that the positive equilibrium is stable.Figure 1
(a) The intersections between the vertical isocline and horizontal isocline with different nutrient concentrationsI, where the red solid line represents the horizontal isocline, the grey dashed line represents the asymptote of the horizontal isocline, and the other lines with different colors represent the vertical isocline that corresponds to different nutrient concentrations I. (b) The stability of the positive equilibrium with the change in the nutrient concentration I.
(a)
(b)Based on the two previous sections, we know that the stability of the boundary equilibriaE0 is relatively simple. However, the stability of the positive equilibria E∗ is complex. We studied the effect of a time delay τ on the positive equilibria of system (3) and our results showed that the time delay has a vital role in system. The critical value τ0 can be obtained by τ0=1/w0arccos⁡w02-B0/B. We set the parameter I=1.3 and the others are the same as the above. By combining the analysis above, we can obtain the critical value: τ0=3.91. The numerical solutions for nutrients are shown in Figures 2(a) and 2(b). In order to study the effect of a time delay, we set τ as 0 and 1 initially, and Figure 2(a) shows that the solution converged to a positive equilibrium in the end. The positive equilibrium is locally asymptotically stable in system (3) without delay. The oscillation does not occur when the system had a delay and the parameter τ<τ0=3.91. The Hopf bifurcation occurred at the critical value of τ0. Thus, the oscillation will occur if τ>τ0=3.91. The numerical solutions for nutrients are shown in Figure 2(b), and we set τ as 0 and 4. Figure 2(b) shows that the oscillation occurs when τ=4.Figure 2
The solution for nutrients whereI=1.3: (a) τ=0 (red dashed line), τ=1 (blue solid line); (b) τ=0 (red dashed line), τ=4 (blue solid line).
(a)
(b)Furthermore, in order to study the relationship between nutrients and phytoplankton, we set the parameterτ as 5 and the numerical solutions for nutrients and phytoplankton are shown in Figure 3. Figure 3(a) shows the relationship between the spatial distributions of nutrients and phytoplankton, which indicates that the density of phytoplankton reaches its maximal value when the concentration of nutrients reaches its minimal value. The density of phytoplankton affects the concentration of nutrients. As the density of phytoplankton increases, the density of nutrients decreases. Therefore, the density of nutrients and phytoplankton are mutually constrained. The phase of system (3) is shown in Figure 3(b).Figure 3
(a) Solutions for nutrients and phytoplankton whereI=1.3, τ=5. (b) Phase diagram.
(a)
(b)


## 5. Conclusion

In this study, we considered a biological system that comprises nutrients and phytoplankton, where we focused on the effects of a delay on the system dynamics.In Section2, we showed that the positive equilibrium E∗ is locally asymptotically stable if there is no time delay. However, as the time delay increases above the critical value, the equilibrium state becomes unstable and a Hopf bifurcation occurs. Using the method introduced by Hassard et al. [24], we derived the direction and stability of the Hopf bifurcation.Based on the computer simulation, we showed that the relationship between nutrients and phytoplankton is mutually constrained. Thus, an abundance of nutrients leads to a major increase of phytoplankton, which depleted the nutrients via consumption. Due to the effect of time delay, the system shows oscillation. It reveals that time delay has a vital effect on the system.The main limitations of this study are that the main theoretical results were confirmed by numerical simulations with hypothetical parameter values. Thus, we aim to obtain some real data to confirm the validity of our system. Furthermore, we did not consider the diffusion and advection of nutrients and phytoplankton; therefore, these aspects should be considered comprehensively in our future research.



---

*Source: 939187-2015-10-11.xml*