# Dynamics Analysis of a Nutrient-Plankton Model with a Time Delay

Authors: Beibei Wang; Min Zhao; Chuanjun Dai; Hengguo Yu; Nan Wang; Pengfei Wang
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

We analyze a nutrient-plankton system with a time delay. We choose the time delay as a bifurcation parameter and investigate the stability of a positive equilibrium and the existence of Hopf bifurcations. By using the center manifold theorem and the normal form theory, the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions are researched. The theoretical results indicate that the time delay can induce a positive equilibrium to switch from a stable to an unstable to a stable state and so on. Numerical simulations show that the theoretical results are correct and feasible, and the system exhibits rich complex dynamics.


---

## Body

## 1. Introduction

Plankton is the basis of all aquatic food chains and its importance for marine ecosystems is widely recognized [1]. The dynamics of zooplankton-phytoplankton systems have been discussed by many authors [2–8]. A wealth of studies have shown that a time delay may have a complex effect on the dynamics of such systems, with effects that include stability switches for equilibria, the existence of Hopf bifurcation periodic solutions, and the direction and stability of Hopf bifurcation [9–18].Algal blooms are a common feature of marine ecosystems, and a high nutrient concentration has an important influence on algal blooms [19]. Plankton-nutrient interaction models can provide a better understanding of plankton dynamics. The model formulated by Huppert et al. [20] consists of two variables, nutrient level N and phytoplankton biomass P, in the following system:(1)
N
˙
=
a
-
b
β
t
N
P
-
e
N
,
P
˙
=
c
β
t
N
P
-
d
P
.Considering the biological significance, all the constants, a, b, c, d, and e, are assumed to be positive. The phytoplankton growth rate β can be represented as a periodic function β
(
t
)
=
β
(
t
+
T
), where T is the forcing period affected by seasonal conditions such as light, salinity, and water temperature.Harmful algal blooms have become a serious environmental problem worldwide and have been widely studied [8, 21]. On the basis of field studies and mathematical modeling, Chattopadhayay et al. [22] proposed the following phytoplankton-zooplankton model:(2)
P
˙
=
r
P
1
-
P
K
-
β
f
P
Z
,
Z
˙
=
β
1
f
P
Z
-
d
Z
-
ρ
g
P
Z
.The authors studied the existence and local stability of steady states and the existence of Hopf bifurcation for system (2) by taking various combinations of f
(
P
) and g
(
P
) [22]. The results showed that toxin-producing plankton is helpful in terminating planktonic blooms by decreasing the grazing pressure of zooplankton.To reflect maturation time, capture time, and other factors, a time delay is often included in mathematical models of population dynamics. Incorporation of a time delay can provide a better understanding of the dynamics of biological models. In recent years, many researchers have studied the fact that time delay plays important roles in biological dynamical systems [18, 23–28]. Following Huppert et al. [20] and Chattopadhayay et al. [22], we consider the following plausible nutrient-phytoplankton-zooplankton system with a time delay:(3)
N
˙
=
α
-
b
N
-
e
N
P
≜
F
1
N
,
P
,
Z
,
P
˙
=
β
N
t
-
τ
P
-
c
P
Z
h
+
P
-
m
P
-
r
P
2
≜
F
2
N
,
P
,
Z
,
Z
˙
=
d
P
Z
h
+
P
-
k
Z
-
ρ
P
Z
h
+
P
≜
F
3
N
,
P
,
Z
,where N is the nutrient concentration, P and Z are the density of the phytoplankton and zooplanktons population, respectively, α is the external source of nutrients flowing into the system, b is the small loss rate to reflect sinking of nutrients from the epilimnion down to the hypolimnion, e and c are the capture rates for phytoplankton on nutrient and zooplankton on phytoplankton, and β and d denote the rates of biomass conversion. h is the half-saturation constant, m and k are the death rates for phytoplankton and zooplankton, and r is the interspecies competition coefficient for phytoplankton, which reflects phytoplankton self-limitation. The parameter τ
≥
0 denotes the time delay, which arises because of the time required by phytoplankton to absorb nutrients.The remainder of the paper is organized as follows. In Section2, we determine the stability of the positive equilibrium and the existence of Hopf bifurcation. Section 3 illustrates the direction of Hopf bifurcation and the stability of bifurcating periodic solutions. To verify the theoretical analysis, we present numerical simulations in Section 4. Finally, Section 5 provides some conclusions.


## 2. Stability and Existence of Hopf Bifurcation

In this section, we study the local stability of the positive equilibrium of system (3). The positive equilibrium E
∗
(
N
∗
,
P
∗
,
Z
∗
) corresponds to the coexistence of the three species, where(4)
P
∗
=
h
k
d
-
k
-
ρ
,
N
∗
=
α
b
+
e
P
∗
,
Z
∗
=
β
N
∗
-
m
-
r
P
∗
h
+
P
∗
c
.It is obvious that E
∗ exists if and only if(5)
α
>
m
d
-
k
-
ρ
+
h
k
r
b
d
-
k
-
ρ
+
e
h
k
β
d
-
k
-
ρ
2
,
d
>
k
+
ρ
.From (5), we know that the positive equilibrium E
∗ of system (3) exists if external nutrients flowing into the system exceed a certain critical value, and the ratio of biomass consumed by zooplankton is greater than the sum of the zooplankton mortality rate and the rate of toxic substance production by phytoplankton. Hereafter we assume that the conditions stated in (5) always hold.Using the translations(6)
u
1
t
=
N
t
-
N
∗
,
u
2
t
=
P
t
-
P
∗
,
u
3
t
=
Z
t
-
Z
∗
,we translate the positive equilibrium E
∗ to the origin. Then the linearized system for system (3) near (
N
∗
,
P
∗
,
Z
∗
) is(7)
u
˙
1
=
a
11
u
1
+
a
12
u
2
,
u
˙
2
=
a
21
u
1
t
-
τ
+
a
22
u
2
+
a
23
u
3
,
u
˙
3
=
a
32
u
2
,where(8)
a
11
=
-
b
-
e
P
∗
,
a
12
=
-
e
N
∗
,
a
32
=
d
-
ρ
h
Z
∗
h
+
P
∗
2
,
a
21
=
β
P
∗
,
a
22
=
c
P
∗
Z
∗
h
+
P
∗
2
-
r
P
∗
,
a
23
=
-
c
P
∗
h
+
P
∗
.The associated characteristic equation of system (7) is(9)
λ
3
+
A
λ
2
+
B
λ
+
C
+
D
λ
e
-
λ
τ
=
0
,where(10)
A
=
-
a
11
+
a
22
,
B
=
a
11
a
22
-
a
23
a
32
,
C
=
a
11
a
23
a
32
>
0
,
D
=
-
a
12
a
21
<
0
.It is well known that the positive equilibrium E
∗ of system (3) is locally asymptotically stable if all roots of (9) have negative real parts and is unstable if (9) has a root with positive real parts. Now, we shall discuss the distribution of the roots of (9). First, when τ
=
0, the characteristic equation becomes(11)
λ
3
+
A
λ
2
+
B
+
D
λ
+
C
=
0
.According to the Routh-Hurwitz criterion, E
∗ is locally asymptotically stable if and only if(12)
A
B
+
D
-
C
>
0
,
A
>
0
.By Corollary 2.4 in the paper of Ruan and Wei [29], we know that if instability occurs for a particular value of time delay, a characteristic root of (9) must intersect the imaginary axis. Hence, suppose that i
ω
(
ω
>
0
) is a purely imaginary root of (9). Inserting this into (9) and separating the real and imaginary parts, we obtain(13)
-
ω
3
+
B
ω
+
D
ω
cos
⁡
ω
τ
=
0
,
-
A
ω
2
+
C
+
D
ω
sin
⁡
ω
τ
=
0
.Squaring and adding these two equations, we have(14)
ω
6
+
M
1
ω
4
+
M
2
ω
2
+
C
2
=
0
,where M
1
=
A
2
-
2
B and M
2
=
B
2
-
2
A
C
-
D
2. Let χ
=
ω
2; then (14) becomes(15)
χ
3
+
M
1
χ
2
+
M
2
χ
+
C
2
=
0
.Form (15), we define a function f
(
χ
)
=
χ
3
+
M
1
χ
2
+
M
2
χ
+
C
2. It is obvious that if Δ
1
=
M
1
2
-
3
M
2
≤
0, then f
′
(
χ
)
=
3
χ
2
+
2
M
1
χ
+
M
2
≥
0
. When Δ
1
>
0, the equation f
′
(
χ
)
=
0 has two real roots:(16)
χ
1
∗
=
-
M
1
+
Δ
1
3
,
χ
2
∗
=
-
M
1
-
Δ
1
3
.Noticing that lim
χ
→
+
∞
⁡
f
(
χ
)
=
+
∞ and f
(
0
)
=
C
2
>
0, then we introduce the following results, which have been proved by many authors [30, 31].Lemma 1.
For (15), we obtain the following results:(i)
IfΔ
1
=
M
1
2
-
3
M
2
≤
0, then (15) has no positive roots.
(ii)
IfΔ
1
=
M
1
2
-
3
M
2
>
0, then (15) has two positive real roots if and only if χ
1
∗
>
0 and f
(
χ
1
∗
)
<
0.
Hence, if the condition,(17)
Δ
1
=
M
1
2
-
3
M
2
>
0
,
χ
0
≔
-
M
1
+
Δ
1
3
>
0
,
f
χ
0
<
0
,is satiated, then (15) has two positive roots, denoted by χ
1, χ
2. Then (14) has two positive roots, namely, ω
1
=
χ
1 and ω
2
=
χ
2
. We define(18)
τ
j
k
=
1
ω
k
φ
k
+
2
j
π
,
k
=
1,2
;
j
=
0,1
,
2
,
…
,where φ
k
∈
(
0,2
π
] satisfies(19)
φ
k
=
2
π
-
arccos
⁡
ω
k
2
-
B
D
,
i
f
A
ω
k
2
-
C
D
ω
k
<
0
,
arccos
⁡
ω
k
2
-
B
D
,
i
f
A
ω
k
2
-
C
D
ω
k
≥
0
.For convenience, let τ
j
1
<
τ
j
2
. Let λ
(
τ
)
=
μ
(
τ
)
+
i
ω
(
τ
) denote the root of (9) such that(20)
μ
τ
j
k
=
0
,
ω
τ
j
k
=
ω
k
k
=
1,2
,
j
=
0,1
,
2
,
…
.Then we have the following lemma.Lemma 2.
If (17) holds, then (21)
d
Re
⁡
λ
τ
d
τ
τ
=
τ
j
1
>
0
,
d
Re
⁡
λ
τ
d
τ
τ
=
τ
j
2
<
0
,
j
=
0,1
,
2
,
…
.Proof.
Differentiating both sides of (9) with respect to τ, we have(22)
d
λ
d
τ
=
D
λ
2
e
-
λ
τ
3
λ
2
+
2
A
λ
+
B
+
D
-
D
λ
τ
e
-
λ
τ
.Then(23)
d
R
e
λ
τ
j
k
d
τ
=
ω
k
2
Δ
2
3
ω
k
4
+
2
A
2
-
2
B
ω
k
2
+
B
2
-
2
A
C
-
D
2
,where(24)
Δ
2
=
B
-
3
ω
k
2
+
D
cos
⁡
ω
k
τ
j
k
-
D
ω
k
τ
j
k
sin
⁡
ω
k
τ
j
k
2
+
2
A
ω
k
-
D
sin
⁡
ω
k
τ
j
k
-
D
ω
k
τ
j
k
cos
⁡
ω
k
τ
j
k
2
.As χ
k
=
ω
k
2, then d
Re
⁡
(
λ
(
τ
j
k
)
)
/
d
τ
=
(
ω
k
2
/
Δ
2
)
f
′
(
χ
k
)
.
If (17) holds, (15) has two positive roots, denoted by χ
1, χ
2, and f
(
χ
0
) is the local minimum value. Assuming τ
j
1
<
τ
j
2, we obtain f
′
(
χ
1
)
>
0 and f
′
(
χ
2
)
<
0
. Hence, (
d
Re
⁡
λ
(
τ
)
/
d
τ
)
τ
=
τ
j
1
>
0, (
d
Re
⁡
λ
(
τ
)
/
d
τ
)
τ
=
τ
j
2
<
0, j
=
0,1
,
2
,
…. This completes the proof.On the basis of the above analysis, we have the following theorem.Theorem 3.
Suppose that (5) and (12) hold; then we obtain the following:(i)
IfΔ
1
=
M
1
2
-
3
M
2
≤
0, then the positive equilibrium E
∗ of system (3) is locally asymptotically stable for all values of τ
≥
0.
(ii)
If (17) holds, there exists a nonnegative integer n, such that the positive equilibrium E
∗ is locally asymptotically stable whenever τ
∈
[
0
,
τ
0
1
)
∪
(
τ
0
2
,
τ
1
1
)
∪
⋯
∪
(
τ
n
-
1
2
,
τ
n
1
) and is unstable whenever τ
∈
(
τ
0
1
,
τ
0
2
)
∪
(
τ
1
1
,
τ
1
2
)
∪
⋯
∪
(
τ
n
-
1
1
,
τ
n
-
1
2
)
∪
(
τ
n
1
,
+
∞
)
. Then system (3) undergoes Hopf bifurcation around E
∗ for every τ
=
τ
j
1,2, j
=
0,1
,
2
…
.


## 3. Direction and Stability of Hopf Bifurcation

From Theorem3 (ii), we have obtained the conditions for the occurrence of Hopf bifurcation when τ
=
τ
j
1,2, j
=
0,1
,
2
,
…
. For convenience, we define τ
0
=
τ
j
1,2, j
=
0,1
,
2
,
…, and ω
0
=
ω
k and k
=
1,2
. In this section, we consider the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions by using the center manifold and normal form theories presented by Hassard et al. [32]. Let x
(
t
)
=
N
(
τ
t
)
-
N
∗, y
(
t
)
=
P
(
τ
t
)
-
P
∗, z
(
t
)
=
Z
(
τ
t
)
-
Z
∗, τ
=
δ
+
τ
0, ϕ
=
(
x
,
y
,
z
)
T, and ϕ
t
(
θ
)
=
ϕ
(
t
+
θ
), where θ
∈
[
-
1,0
]. Using Taylor series expansion about E
∗, system (3) can be written as(25)
x
˙
t
y
˙
t
z
˙
t
=
τ
K
1
x
t
y
t
z
t
+
K
2
x
t
-
1
y
t
-
1
z
t
-
1
+
f
,where (26)
K
1
=
a
11
a
12
0
0
a
22
a
23
0
a
32
0
,
K
2
=
0
0
0
a
21
0
0
0
0
0
,
f
=
f
δ
ϕ
t
=
f
1
f
2
f
3
,
f
1
=
-
e
x
t
y
t
,
f
2
=
β
x
t
-
1
y
t
+
c
h
Z
∗
h
+
P
∗
-
3
-
r
y
2
t
-
c
h
h
+
P
∗
-
2
y
t
z
t
-
c
h
Z
∗
h
+
P
∗
-
4
y
3
t
+
⋯
,
f
3
=
ρ
-
d
h
Z
∗
h
+
P
∗
-
3
y
2
t
+
d
-
ρ
h
h
+
P
∗
-
2
y
t
z
t
+
d
-
ρ
h
Z
∗
h
+
P
∗
-
4
y
3
t
+
⋯
.Then δ
=
0 is Hopf bifurcation value of system (3).LetL
δ
(
ϕ
)
=
(
δ
+
τ
0
)
K
1
ϕ
(
0
)
+
K
2
ϕ
(
-
1
)
. According to the Riesz representation theorem, there exists a 3
×
3 matrix η
(
θ
,
δ
) of the bounded variation for θ
∈
[
-
1,0
] such that(27)
L
δ
ϕ
=
∫
-
1
0
d
η
θ
,
δ
ϕ
θ
,
for
ϕ
∈
C
-
1,0
,
R
3
.In fact, we can select(28)
η
θ
,
δ
=
τ
0
+
δ
K
1
υ
θ
-
τ
0
+
δ
K
2
υ
θ
+
1
,where υ denote the Dirac delta function:(29)
υ
θ
=
0
,
θ
≠
0
,
1
,
θ
=
0
.For ϕ
∈
C
(
[
-
1,0
]
,
R
3
), we define(30)
Γ
δ
ϕ
θ
=
d
ϕ
θ
d
θ
,
θ
∈
-
1,0
,
∫
-
1
0
d
η
s
,
δ
ϕ
s
,
θ
=
0
,
R
δ
ϕ
θ
=
0
,
θ
∈
-
1,0
,
f
δ
ϕ
,
θ
=
0
.Then system (25) can be rewritten as(31)
u
˙
t
=
Γ
δ
u
t
+
R
δ
u
t
,where u
=
(
x
,
y
,
z
)
T and u
t
=
u
(
t
+
θ
) for θ
∈
[
-
1,0
]
.Forψ
∈
C
1
(
[
0,1
]
,
(
R
3
)
∗
), we define the adjoint operator Γ
∗ as(32)
Γ
∗
ψ
s
=
-
d
ϕ
s
d
s
,
s
∈
0,1
,
∫
-
1
0
ψ
-
ξ
d
η
ξ
,
0
,
s
=
0
,and a bilinear inner product(33)
ψ
s
,
ϕ
θ
=
ψ
-
0
ϕ
0
-
∫
-
1
0
∫
ξ
=
0
θ
ψ
-
ξ
-
θ
d
η
θ
ϕ
ξ
d
ξ
,where η
(
θ
)
=
η
(
θ
,
0
). From the results in the previous section, we know that ±
i
ω
0
τ
0 are the eigenvalues of Γ
(
0
). Thus, ±
i
ω
0
τ
0 are also the eigenvalues of Γ
∗.Letq
(
θ
)
=
(
1
,
γ
1
,
γ
2
)
T
e
i
ω
0
τ
0
θ be the eigenvector of Γ
(
0
) corresponding to the eigenvalue i
ω
0
τ
0, and let q
∗
(
s
)
=
D
(
1
,
γ
3
,
γ
4
)
e
i
ω
0
τ
0
s be the eigenvector of Γ
∗ corresponding to the eigenvalue -
i
ω
0
τ
0. With the condition Γ
(
0
)
q
(
θ
)
=
i
ω
0
τ
0
q
(
θ
), we obtain(34)
i
w
0
I
-
K
1
+
K
2
e
-
i
ω
0
τ
0
q
0
=
0
,where I is an identity matrix of order 3; that is,(35)
i
ω
0
-
a
11
-
a
12
0
-
a
21
e
-
i
ω
0
τ
0
i
ω
0
-
a
22
-
a
23
0
-
a
32
i
ω
0
1
γ
1
γ
2
=
0
0
0
.Then we can obtain(36)
γ
1
=
i
ω
0
-
a
11
a
12
,
γ
2
=
i
ω
0
-
a
11
a
32
i
ω
0
a
12
.Similarly, with Γ
∗
q
∗
(
0
)
=
-
i
ω
0
τ
0
q
∗
(
0
) or [
-
i
w
0
I
-
(
K
1
T
+
K
2
T
e
i
ω
0
τ
0
)
]
(
q
∗
(
0
)
)
T
=
0, we have(37)
i
ω
0
+
a
11
a
21
e
-
i
ω
0
τ
0
0
a
12
i
ω
0
+
a
22
a
32
0
a
23
i
ω
0
1
γ
3
γ
4
=
0
0
0
.Solving the above equations, we obtain(38)
γ
3
=
-
i
ω
0
+
a
11
a
21
e
-
i
ω
0
τ
0
,
γ
4
=
a
23
i
ω
0
+
a
11
i
ω
0
a
21
e
-
i
ω
0
τ
0
.Hence, (39)
q
∗
s
,
q
θ
=
q
-
∗
0
q
0
-
∫
-
1
0
∫
ξ
=
0
θ
q
-
∗
ξ
-
θ
d
η
θ
q
ξ
d
ξ
=
D
-
D
1
,
γ
-
3
,
γ
-
4
1
,
γ
1
,
γ
2
T
-
∫
-
1
0
∫
ξ
=
0
θ
D
-
1
,
γ
-
3
,
γ
-
4
e
-
i
ω
0
τ
0
ξ
-
θ
d
η
θ
1
,
γ
1
,
γ
2
T
e
i
ω
0
τ
0
ξ
d
ξ
=
D
-
1
+
γ
1
γ
-
3
+
γ
2
γ
-
4
-
∫
-
1
0
θ
D
-
1
,
γ
-
3
,
γ
-
4
e
i
ω
0
τ
0
θ
d
η
θ
1
,
γ
1
,
γ
2
T
=
D
-
1
+
γ
1
γ
-
3
+
γ
2
γ
-
4
+
τ
0
D
-
e
-
i
ω
0
τ
0
1
,
γ
-
3
,
γ
-
4
K
2
1
,
γ
1
,
γ
2
T
=
D
-
1
+
γ
1
γ
-
3
+
γ
2
γ
-
4
+
τ
0
γ
-
3
β
P
∗
e
-
i
ω
0
τ
0
.Then we can choose(40)
D
=
1
+
γ
-
1
γ
3
+
γ
-
2
γ
4
+
τ
0
γ
3
β
P
∗
e
i
ω
0
τ
0
-
1such that q
∗
(
s
)
,
q
(
θ
)
=
1 and q
∗
(
s
)
,
q
-
(
θ
)
=
0
.In the following, we first compute the coordinates to describe the center manifoldC
0 at δ
=
0 according to the approach of Hassard et al. [32]. Let u
t be the solution of (31) when δ
=
0. We define(41)
z
t
=
q
∗
,
u
t
,
W
t
,
θ
=
u
t
θ
-
2
R
e
z
t
q
θ
.On the center manifoldC
0, we have(42)
W
t
,
θ
=
W
z
,
z
-
,
θ
=
W
20
θ
z
2
2
+
W
11
θ
z
z
-
+
W
02
θ
z
-
2
2
+
W
30
θ
z
3
6
+
⋯
,where z and z
- are local coordinates for the center manifold C
0 in the direction of q
∗ and q
-
∗, respectively. In this paper, we only consider real solutions as W is real if u
t is real. Since δ
=
0, for the solution u
t
∈
C
0 of (31), we obtain(43)
z
˙
t
=
q
∗
,
u
˙
t
=
i
ω
0
τ
0
z
t
+
q
∗
θ
,
R
0
W
z
,
z
-
,
θ
+
2
R
e
z
t
q
θ
=
i
ω
0
τ
0
z
t
+
q
-
∗
0
f
0
W
z
,
z
-
,
0
+
2
R
e
z
t
q
0
≜
i
ω
0
τ
0
z
t
+
g
z
,
z
-
,where(44)
g
z
,
z
-
=
q
-
∗
0
f
0
W
z
,
z
-
,
0
+
2
R
e
z
t
q
0
=
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
g
21
z
2
z
-
2
+
⋯
.From (31) and (41), we have(45)
W
˙
=
u
˙
t
-
2
R
e
z
˙
q
=
u
t
-
z
˙
q
-
z
-
˙
q
-
=
Γ
u
t
+
R
u
t
-
z
˙
q
-
z
-
˙
q
-
=
Γ
W
+
2
R
e
z
q
+
R
W
+
2
R
e
z
q
-
2
R
e
g
z
,
z
-
q
θ
=
Γ
W
-
2
R
e
q
-
∗
θ
R
0
q
θ
,
θ
∈
-
1,0
Γ
W
-
2
R
e
q
-
∗
0
R
0
q
0
+
R
0
,
θ
=
0
≜
Γ
W
+
H
z
,
z
-
,
θ
,where(46)
H
z
,
z
-
,
θ
=
H
20
θ
z
2
2
+
H
11
θ
z
z
-
+
H
02
θ
z
-
2
2
+
⋯
.On the center manifold C
0, from (42)–(44), we obtain(47)
W
˙
z
,
z
-
=
W
z
z
˙
+
W
z
-
z
-
˙
=
W
20
θ
z
z
˙
+
W
11
θ
z
-
z
˙
+
W
02
θ
z
-
z
-
˙
+
W
11
θ
z
z
-
˙
+
O
z
,
z
-
3
=
W
20
θ
z
i
w
0
τ
0
z
+
g
z
,
z
-
+
W
11
θ
z
-
i
w
0
τ
0
z
+
g
z
,
z
-
+
W
02
θ
z
-
-
i
w
0
τ
0
z
-
+
g
-
z
,
z
-
+
W
11
θ
z
-
i
w
0
τ
0
z
-
+
g
-
z
,
z
-
+
O
z
,
z
-
3
=
i
w
0
τ
0
W
20
θ
z
2
-
i
w
0
τ
0
W
02
θ
z
-
2
+
O
z
,
z
-
3
.Combining (45)–(47) and (42) and comparing the coefficients, we have(48)
2
i
ω
0
τ
0
I
-
Γ
0
W
20
θ
=
H
20
θ
,
Γ
0
W
11
θ
=
-
H
11
θ
,
2
i
ω
0
τ
0
I
+
Γ
0
W
02
θ
=
-
H
02
θ
.Since(49)
u
t
θ
=
x
t
θ
,
y
t
θ
,
z
t
θ
T
=
2
Re
⁡
z
t
q
θ
+
W
t
,
θ
=
z
q
θ
+
z
-
q
-
θ
+
W
t
,
θ
,noting that q
(
0
)
=
(
1
,
γ
1
,
γ
2
)
T, we have(50)
x
t
=
z
+
z
-
+
W
20
1
0
z
2
2
+
W
11
1
0
z
z
-
+
W
02
1
0
z
-
2
2
+
⋯
,
y
t
=
γ
1
z
+
γ
-
1
z
-
+
W
20
2
0
z
2
2
+
W
11
2
0
z
z
-
+
W
02
2
0
z
-
2
2
+
⋯
,
z
t
=
γ
2
z
+
γ
-
2
z
-
+
W
20
3
0
z
2
2
+
W
11
3
0
z
z
-
+
W
02
3
0
z
-
2
2
+
⋯
,
x
t
-
1
=
e
-
i
w
0
τ
0
z
+
e
i
w
0
τ
0
z
-
+
W
20
1
-
1
z
2
2
+
W
11
1
-
1
z
z
-
+
W
02
1
-
1
z
-
2
2
+
⋯
.According to the definition of g
(
z
,
z
-
), we obtain(51)
g
z
,
z
-
=
q
-
∗
0
f
0
W
z
,
z
-
,
0
+
2
Re
⁡
z
t
q
0
=
D
-
1
,
γ
-
1
,
γ
-
4
τ
0
f
1
,
f
2
,
f
3
T
.Comparison of the coefficients with (44) yields(52)
g
20
=
2
τ
0
D
-
γ
1
γ
1
σ
1
+
γ
2
σ
2
+
β
γ
-
3
e
-
i
ω
0
τ
0
-
e
,
g
11
=
τ
0
D
-
γ
1
γ
-
1
σ
1
+
γ
-
2
σ
2
+
β
γ
-
3
e
i
ω
0
τ
0
-
e
+
γ
-
1
γ
1
σ
1
+
γ
2
σ
2
+
β
γ
-
3
e
-
i
ω
0
τ
0
-
e
,
g
02
=
2
τ
0
D
-
γ
-
1
γ
-
1
σ
1
+
γ
-
2
σ
2
+
β
γ
-
3
e
i
ω
0
τ
0
-
e
,
g
21
=
2
τ
0
D
-
γ
1
σ
1
W
11
2
0
+
σ
2
W
11
3
0
+
β
γ
-
3
W
11
1
-
1
-
e
W
11
1
0
+
γ
-
1
2
σ
1
W
20
2
0
+
σ
2
W
20
3
0
+
β
γ
-
3
W
20
1
-
1
-
e
W
20
1
0
+
W
20
2
0
2
σ
1
γ
-
1
+
σ
2
γ
-
2
+
β
γ
-
3
e
i
ω
0
τ
0
-
e
+
W
11
2
0
σ
1
γ
1
+
σ
2
γ
2
+
β
γ
-
3
e
-
i
ω
0
τ
0
-
e
,where σ
1
=
c
γ
-
3
+
γ
-
4
(
ρ
-
d
)
h
Z
∗
/
(
h
+
P
∗
)
3
-
γ
-
3
r, σ
2
=
(
γ
-
4
(
d
-
ρ
)
h
Z
∗
-
c
h
γ
-
3
)
/
(
h
+
P
∗
)
2, and W
20
(
0
), W
20
(
-
1
)
,
W
11
(
0
), and W
11
(
-
1
) are still unknown. Next, we compute W
20
(
θ
) and W
11
(
θ
) accurately.Forθ
∈
[
-
1,0
), we have(53)
H
z
,
z
-
,
θ
=
-
2
Re
⁡
q
-
∗
θ
R
0
q
θ
=
-
g
z
,
z
-
q
θ
-
g
-
z
,
z
-
q
-
θ
=
-
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
⋯
q
θ
-
g
-
20
z
-
2
2
+
g
-
11
z
z
-
+
g
-
02
z
2
2
+
⋯
q
-
θ
.Comparing the coefficients with (46), we obtain(54)
H
20
θ
=
-
g
20
q
θ
-
g
-
02
q
-
θ
,
H
11
θ
=
-
g
11
q
θ
-
g
-
11
q
-
θ
.From the definition of Γ
(
0
) and (48) and (54), we have(55)
W
˙
20
θ
=
2
i
w
0
τ
0
W
20
θ
+
g
20
q
θ
+
g
-
02
q
-
θ
.For q
(
θ
)
=
(
1
,
γ
1
,
γ
2
)
T
e
i
ω
0
τ
0
θ, solving (55) yields(56)
W
20
θ
=
i
g
20
q
0
w
0
τ
0
e
i
w
0
τ
0
θ
+
i
g
-
02
q
-
0
3
w
0
τ
0
e
-
i
w
0
τ
0
θ
+
E
1
e
2
i
w
0
τ
0
θ
.Similarly, using the same method, we obtain(57)
W
11
θ
=
-
i
g
11
q
0
w
0
τ
0
e
i
w
0
τ
0
θ
+
i
g
-
11
q
-
0
w
0
τ
0
e
-
i
w
0
τ
0
θ
+
E
2
.
E
1 and E
2 are three-dimensional vectors that can be determined by setting θ
=
0 in H
(
z
,
z
-
,
θ
)
.For (45), when θ
=
0, H
(
z
,
z
-
,
0
)
=
-
2
Re
⁡
{
q
-
∗
(
0
)
R
(
0
)
q
(
0
)
}
+
R
(
0
). Thus,(58)
H
20
0
=
-
g
20
q
0
-
g
-
02
q
-
0
+
τ
0
-
e
γ
1
,
γ
1
Θ
1
,
γ
1
Θ
3
T
,
H
11
0
=
-
g
11
q
0
-
g
-
11
q
-
0
+
τ
0
-
e
γ
1
+
γ
-
1
,
γ
1
Θ
2
+
γ
-
1
Θ
1
,
γ
-
1
Θ
3
+
γ
1
Θ
4
T
,where(59)
Θ
1
=
β
e
-
i
ω
0
τ
0
-
r
γ
1
+
γ
1
c
h
Z
∗
h
+
P
∗
3
-
γ
2
c
h
h
+
P
∗
2
,
Θ
2
=
β
e
i
ω
0
τ
0
-
r
γ
-
1
+
γ
-
1
c
h
Z
∗
h
+
P
∗
3
-
γ
-
2
c
h
h
+
P
∗
2
,
Θ
3
=
h
d
-
ρ
γ
2
h
+
γ
2
P
∗
-
γ
1
Z
∗
h
+
P
∗
3
,
Θ
4
=
h
d
-
ρ
γ
-
2
h
+
γ
-
2
P
∗
-
γ
-
1
Z
∗
h
+
P
∗
3
.According to (48), we obtain(60)
2
i
w
0
τ
0
W
20
0
-
τ
0
K
1
W
20
0
-
τ
0
K
2
W
20
-
1
=
H
20
0
,
τ
0
K
1
W
11
0
+
τ
0
K
2
W
11
-
1
=
-
H
11
0
.From (56), (57), and (60), we have(61)
E
1
=
1
τ
0
2
i
w
0
I
-
K
1
-
K
2
e
-
2
i
w
0
τ
0
-
1
H
20
0
+
2
g
20
q
0
+
2
g
-
02
q
-
0
3
+
i
g
20
K
1
q
0
w
0
+
i
g
-
02
K
1
q
-
0
3
w
0
+
i
g
20
e
-
i
w
0
τ
0
K
2
q
0
w
0
+
i
g
-
02
e
i
w
0
τ
0
K
2
q
-
0
3
w
0
,
E
2
=
1
τ
0
-
L
1
-
L
2
-
1
H
11
0
-
i
g
11
K
1
q
0
w
0
+
i
g
-
11
K
1
q
-
0
w
0
-
i
g
11
e
-
i
w
0
τ
0
L
2
q
0
w
0
+
i
g
-
11
e
i
w
0
τ
0
K
2
q
-
0
w
0
.Now, all g
i
j can be expressed in terms of parameters. Therefore, we can evaluate the following values:(62)
c
1
0
=
i
2
w
0
τ
0
g
20
g
11
-
2
g
11
2
-
g
02
2
3
+
g
21
2
,
μ
2
=
-
Re
⁡
c
1
0
Re
⁡
λ
′
τ
0
,
β
2
=
Re
⁡
c
1
0
,
T
2
=
-
I
m
c
1
0
+
μ
2
I
m
λ
′
τ
0
w
0
τ
0
.According to the above analysis, we can obtain the following theorem about the properties of Hopf bifurcation.Theorem 4.
Forμ
2
,
β
2, and T
2 defined as above, the properties of Hopf bifurcation at the critical value τ
=
τ
0 are as follows:(i)
Ifμ
2
>
0
(
<
0
), Hopf bifurcation is supercritical (subcritical).
(ii)
Ifβ
2
<
0
(
>
0
), the Hopf periodic solutions are stable (unstable).
(iii)
IfT
2
>
0
(
<
0
), the period of the bifurcation periodic solution of system (3) increases (decreases).


## 4. Numerical Simulations

In this section, we verify the theoretical results proved in previous sections using numerical simulations for the following parameter values:(63)
b
=
0.4
,
e
=
0.1
,
β
=
1.2
,
c
=
0.1
,
h
=
3.5
,
m
=
0.2
,
r
=
0.2
,
d
=
1.2
,
k
=
0.2
,
ρ
=
0.1
.According to the analysis above, the existence and stability of the interior equilibrium for system (3) both change with the value of the external source of nutrients α, as illustrated in Figure 1. In Figure 1(a), n (the conversion frequency of stability switches for E
∗) changes with α. From (5), the condition for the existence of E
∗ is α
>
0
.
1416. Region I in Figure 1(a) indicates that system (3) has no positive equilibrium if the external source of nutrients flowing into the system is less than the critical value, 0.1416. In region IV, E
∗ is unstable when condition (12) is not satisfied. This shows that when α is greater than a critical value, 0.755, E
∗ is always unstable and a time delay will no longer affect system (3). When the value of α is in region II, condition (17) is not satisfied and (14) has no positive roots. Hence, E
∗ is stable only when τ
=
0. In region III, Theorem 3 (ii) holds and there exists a nonnegative integer n such that E
∗ is locally asymptotically stable when (64)
τ
∈
0
,
τ
0
1
∪
τ
0
2
,
τ
1
1
∪
⋯
∪
τ
n
-
1
2
,
τ
n
1
,and it is unstable when (65)
τ
∈
τ
0
1
,
τ
0
2
∪
τ
1
1
,
τ
1
2
∪
⋯
∪
τ
n
-
1
1
,
τ
n
-
1
2
∪
τ
n
1
,
+
∞
.To simulate clear stability switches, Figure 1(b) shows interval values for stability switches and Hopf bifurcation about E
∗ with the values of α selected as 0.355, 0.2991, 0.2948, 0.2819, 0.2776, and 0.2733 and τ
j
k
(
k
=
1,2
;
j
=
0,1
,
…
,
n
). For example, when α
=
0.2819, Theorem 3 (ii) is satisfied. Then (14) has two positive roots, ω
1
=
0.2625 andω
2
=
0.2157. From (18), we obtain(66)
τ
0
1
=
8.3830
<
τ
0
2
=
14.9808
<
τ
1
1
=
32.3148
<
τ
1
2
=
44.1079
<
τ
2
1
=
56.2466
<
τ
2
2
=
73.2349
<
τ
3
1
=
80.1783
<
τ
3
2
=
102.362
<
τ
4
1
=
104.1101
<
τ
5
1
=
128.0419
<
τ
4
2
=
131.4890
.In order to more clearly show the interval values for stability switches, we use the rectangles to represent the stable and unstable interval of positive equilibrium E
∗. Therefore, in Figure 1(b), it is easy to see that the stability switches of positive equilibrium E
∗ occur with the change of time delay τ.Figure 1
(a) Relationship between the conversion frequencyn of stability switches for positive equilibrium E
∗ and the external source of nutrients α, where the blue spots represent the value of n, the black line represents the boundary of the existence of E
∗, the green line represents the boundary of the appearance of n, and the red line represents the disappearance of n. (b) Interval values for stability switches and Hopf bifurcation about E
∗, where the green diamonds and the blue dots represent τ
j
1
(
0
≤
j
≤
n
) and τ
j
2
(
0
≤
j
≤
n
-
1
), respectively, and the honeydew and deep pink rectangles, respectively, represent the stable and unstable interval of E
∗.
(a)
(b)On the basis of the above analysis, we present examples in Figures2–5. The equilibrium point E
∗ is locally asymptotically stable in Figure 2. When τ passes through critical values, E
∗ loses its stability, and bifurcating periodic solutions occur and are stable, as shown in Figures 3-4. In Figure 5, the bifurcating solutions disappear and chaos occurs when τ passes through critical values.Figure 2
Numerical solutions of system (3) with α
=
0.2819, τ
=
5, and initial point (0.7, 1.2, 18). (a) shows the time series of N (black line), P (red line), and Z (blue line) and (b) shows the associated phase diagram of system (3) (magenta line). The equilibrium E
∗ is locally asymptotically stable.
(a)
(b)Figure 3
Numerical solutions of system (3) with α
=
0.2819, τ
=
10, and initial point (0.6, 0.7, 15). (a) shows the time series of N (black line), P (red line), and Z (blue line) and (b) shows the associated phase diagram of system (3) (magenta line). There exists a positive period-1 solution of system (3).
(a)
(b)Figure 4
Numerical solutions of system (3) with α
=
0.4, τ
=
50, and initial point (0.8, 1, 18). (a) shows the time series of N (black line), P (red line), and Z (blue line) and (b) shows the associated phase diagram of system (3) (magenta line). There exists a positive period-3 solution of system (3).
(a)
(b)Figure 5
Numerical solutions of system (3) with α
=
0.4, τ
=
70, and initial point (0.8, 1, 18). (a) shows the time series of N (black line), P (red line), and Z (blue line) and (b) shows the associated phase diagram of system (3) (magenta line). There exists a chaos attractor of system (3).
(a)
(b)The numerical results in Figures2–5 show that oscillating solutions are strongly affected by τ. To generalize the dynamical behavior of system (3) influenced by parameter τ, we make numerical simulation using a wide range of values of parameter τ to show the bifurcation diagram (see Figure 6(a)). Results for the nutrient density and phytoplankton density as a function of τ are similar and are not shown here. From Figure 6(a), the order-2 periodic solution appears in the system (3) by increasing parameter τ. Keeping increasing τ, chaos occurs. Subsequently, an order-3 periodic solution bifurcates from chaos. Finally, chaos occurs again with increase of τ. To investigate the influence of parameter τ on amplitudes of the oscillations observed when the positive equilibrium undergoes Hopf bifurcations, we simulate the stability switches of positive equilibrium in the component Z (see Figure 6(b)). In Figure 6(b), the green points are τ
0
1
=
6.7822
,
τ
0
2
=
15.5408
,
τ
1
1
=
28.6033
,
τ
1
2
=
44.0317, and τ
2
1
=
50.4243, respectively. The black solid lines denote that the positive equilibrium component Z
∗
=
16.92598 is stable when τ
∈
[
0
,
τ
0
1
)
∪
(
τ
0
2
,
τ
1
1
)
∪
(
τ
1
2
,
τ
2
1
). The blue and red curves denote the maximum and minimum Z values, respectively. The black dashed lines indicate that when τ
∈
(
τ
0
1
,
τ
0
2
)
∪
(
τ
1
1
,
τ
1
2
)
∪
(
τ
2
1
,
+
∞
), the equilibrium loses its stability and oscillations occur. Therefore, the numerical simulations agree with the theoretical predictions.Figure 6
(a) Bifurcation diagram of system (3) with α
=
0.5 and initial point (0.6, 0.7, 16); (b) the stability switches of positive equilibrium in the component Z with respect to time delay τ for α
=
0.2
991 and Z
∗
=
16.92598, where the green point denotes the Hopf bifurcation point, the black solid line denotes the stable positive equilibrium, the black dashed line denotes the unstable positive equilibrium, and the blue and red curves denote the maximum and minimum values of zooplankton density Z, respectively.
(a)
(b)


## 5. Conclusion

In this paper, we have studied the dynamics of a nutrient-plankton system with a time delay. We have investigated the stability of the positive equilibrium and the existence of Hopf bifurcation. By using center manifold theory and the normal form method, we determined the direction of Hopf bifurcation and the stability of bifurcating periodic solutions.In detail, we have found that if some conditions are satisfied, the phenomenon of stability switches arises from system (3). When the time delay passes through some critical values, there exists a nonnegative integer n such that the positive equilibrium switches n times from stability to instability to stability and so on, and Hopf bifurcations occur at the positive equilibrium. The numerical simulations in Figure 1 indicate that the theoretical results are correct, and the number of stability switches changes with the value of the external source of nutrients α. This result indicates that the external source of nutrients flowing into the system has important influence on the complex dynamics of system (3).Numerical simulations also showed that as the time delay further increases, the periodic solutions disappear and chaos appears. All these results not only will help in further investigating the dynamics of pelagic ecosystem in theory but also are very useful to understand the complex phenomena really happening in marine ecosystem.



---

*Source: 9797624-2016-01-27.xml*