# Stability and Hopf Bifurcation Analysis of a Nutrient-Phytoplankton Model with Delay Effect

Authors: Xinhong Pan; Min Zhao; Chuanjun Dai; Yapei Wang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A delay differential system is investigated based on a previously proposed nutrient-phytoplankton model. The time delay is regarded as a bifurcation parameter. Our aim is to determine how the time delay affects the system. First, we study the existence and local stability of two equilibria using the characteristic equation and identify the condition where a Hopf bifurcation can occur. Second, the formulae that determine the direction of the Hopf bifurcation and the stability of periodic solutions are obtained using the normal form and the center manifold theory. Furthermore, our main results are illustrated using numerical simulations.


---

## Body

## 1. Introduction

Phytoplankton plays a very important role as the first trophic level in aquatic ecosystems. To describe the complex dynamics of phytoplankton populations, the dynamic relationship between phytoplankton and nutrients has been investigated theoretically for a long time, as well as experimentally. Since the pioneering work of Riley et al. [1], various nutrient-phytoplankton models have been proposed and analyzed [2–9]. Several of these models have been shown to predict phytoplankton dynamics successfully in specific situations.The model proposed by Taylor et al. [10] in 1986 describes the nutrient-dependent growth of a single phytoplankton population by considering sinking and variable vertical mixing:
(1)
∂
N
∂
t
=
-
γ
(
α
(
I
)
ϕ
(
N
)
-
ε
m
)
P
+
∂
∂
z
(
K
(
z
,
t
)
∂
N
∂
z
)
,
∂
P
∂
t
=
(
α
(
I
)
ϕ
(
N
)
-
m
)
P
-
v
∂
P
∂
z
+
∂
∂
z
(
K
(
z
,
t
)
∂
P
∂
z
)
,
where N and P are the concentrations of the nutrient and phytoplankton, respectively; the specific growth rate of the phytoplankton is taken to be the product of two parts, a dependence (α) on incident light (I) and a function (ϕ) of the concentration of a single nutrient; m is the death rate of phytoplankton; ε
(
0
<
ε
≤
1
) is the regeneration efficiency; v is the sinking rate of phytoplankton; K
(
z
,
t
) is the turbulent diffusion coefficient. The diffusion coefficients for phytoplankton and nutrients are assumed to be the same for simplicity. It is well known that the abundance of the phytoplankton population is affected by many environmental factors, such as the water temperature, salinity, and sunlight intensity [11]. In system (1), I is the light intensity and numerical simulation indicates that the light intensity can affect the result.In this study, we consider an approximated model of system (1) with delay effect as a model of a layer of phytoplankton growing over a pool of nutrients:
(2)
d
N
d
t
=
-
γ
(
β
N
-
ε
m
)
P
+
k
h
(
N
0
-
N
)
,
d
P
d
t
=
(
β
N
(
t
-
τ
)
-
m
-
v
h
-
k
h
)
P
,
where β is the specific growth rate of phytoplankton; h is the thickness of the layer; k is the turbulent diffusion coefficient; τ is a positive delay, that is, the time required to convert nutrients into phytoplankton; N
0 and P
0 are the concentrations of nutrients and phytoplankton below the layer, respectively. We assume that all parameters are positive, except that P
0 is negligible and equal to zero.Taylor et al. [10] described the dynamics of the above system without considering the effect of delay, mainly using numerical methods. Pardo [12] conducted a mathematical study of the local and global stability of the equilibria in the same system. He proved the positivity and boundedness of the solutions, which made sure that the model is biologically sound. He found that the interior equilibrium point of the system is locally and globally asymptotically stable if it exists and that the boundary equilibrium point is also globally asymptotically stable if the system has only one equilibrium point.Time delay is known to play important roles in biological dynamical systems, which have been studied by many researchers in recent years [13–27]. Our aim is to investigate how time delay may affect the system (2). Thus, we use the delay τ as a bifurcation parameter.The remainder of this paper is organized as follows. In Section2, we consider the local stability of the equilibria and the condition where Hopf bifurcation can occur based on the characteristic equation. In Section 3, we derive an explicit algorithm to determine the direction of the Hopf bifurcation and the stability of the periodic solutions. The results of numerical simulations are presented to support the theoretical results in Section 4. The paper ends with a brief conclusion.


## 2. Stability Analysis of the Equilibria

In this section, we mainly consider the existence and stability of the nonnegative equilibria of system (2). The equations for the equilibria are as follows:
(3)
-
γ
(
β
N
-
ε
m
)
P
+
k
h
(
N
0
-
N
)
=
0
,
(
β
N
-
m
-
v
h
-
k
h
)
P
=
0
.There are two solutions:E
0
(
N
0
,
0
), E
*
(
N
*
,
P
*
) if P is not equal to zero and N
0
>
N
* exists, where
(4)
N
*
=
1
β
(
m
+
v
h
+
k
h
)
,
P
*
=
k
(
N
0
-
N
*
)
γ
h
(
β
N
*
-
ε
m
)
=
k
(
β
N
0
-
m
-
v
/
h
-
k
/
h
)
β
γ
h
(
m
+
v
/
h
+
k
/
h
-
ε
m
)
.


### 2.1. Local Stability of EquilibriumE
0

We linearize (2) about E
0 to obtain the linear system
(5)
d
N
d
t
=
-
k
h
N
-
γ
β
(
N
0
-
ε
m
)
P
,
d
P
d
t
=
(
β
N
0
-
m
-
v
h
-
k
h
)
P
.There are two eigenvalues where-
k
/
h is always negative and β
N
0
-
m
-
v
/
h
-
k
/
h
>
0 only if N
0
>
N
*. Therefore, we have the following theorem.Theorem 1.
The equilibriumE
0
(
N
0
,
0
) is stable if β
N
0
-
m
-
v
/
h
-
k
/
h
<
0 and unstable if β
N
0
-
m
-
v
/
h
-
k
/
h
>
0.
We recall thatβ
N
0
-
m
-
v
/
h
-
k
/
h
>
0 is equivalent to N
0
>
N
*, which is the existing condition of the unique interior equilibrium. Thus, E
0 is stable only if E
* does not exist.


### 2.2. Local Stability and the Hopf Bifurcation of EquilibriumE
*

Similarly, we linearize (2) about E
* to obtain the linear system
(6)
d
N
d
t
=
(
-
γ
β
P
*
-
k
h
)
N
-
γ
(
β
N
*
-
ε
m
)
P
,
d
P
d
t
=
β
P
*
N
(
t
-
τ
)
.The corresponding characteristic equation is(7)
λ
2
+
(
γ
β
P
*
+
k
h
)
λ
+
γ
β
(
β
N
*
-
ε
m
)
P
*
e
-
λ
τ
=
0
or
(8)
λ
2
+
A
λ
+
B
e
-
λ
τ
=
0
,
where
(9)
A
=
(
γ
β
P
*
+
k
h
)
>
0
,
B
=
γ
β
(
β
N
*
-
ε
m
)
P
*
>
0
.Whenτ
=
0, the characteristic equation is
(10)
λ
2
+
A
λ
+
B
=
0
,
and the two eigenvalues satisfy
(11)
λ
1
+
λ
2
=
-
A
<
0
,
λ
1
·
λ
2
=
B
>
0
,
which indicates that E
* is locally asymptotically stable when τ
=
0.Whenτ
≠
0, we assume that λ
=
μ
(
τ
)
+
i
w
(
τ
) is a root of (8). Substituting this in (8) we have
(12)
μ
2
-
w
2
+
2
μ
w
i
+
μ
A
+
A
w
i
+
B
e
-
μ
τ
(
cos
⁡
τ
w
-
i
sin
τ
w
)
=
0
or
(13)
μ
2
-
w
2
+
μ
A
=
-
B
e
-
μ
τ
cos
⁡
τ
w
,
2
μ
w
+
A
w
=
B
e
-
μ
τ
sin
τ
w
.Letμ
(
τ
)
=
0 and w
(
τ
)
=
w
0
>
0, and hence (13) is equal to
(14)
w
0
2
=
B
cos
⁡
τ
w
0
,
A
w
0
=
B
sin
τ
w
0
.Hence,(15)
w
0
4
+
A
2
w
0
2
-
B
2
=
0
.LetΔ
=
A
4
+
4
B
2
>
0, and then w
0
=
(
-
A
2
+
Δ
)
/
2.However, we can solve the first formula in (13) to obtain
(16)
τ
=
τ
k
=
1
w
0
(
arccos
w
0
2
B
+
2
k
π
)
,
k
=
0,1
,
2
,
…
.Substitutingλ
(
τ
) in (8) and taking the derivatives with respect to τ, we have
(17)
2
λ
(
τ
)
d
λ
(
τ
)
d
τ
+
A
d
λ
(
τ
)
d
τ
+
B
e
-
λ
τ
(
-
τ
d
λ
(
τ
)
d
τ
-
λ
(
τ
)
)
=
0
.Forτ
=
τ
0,
(18)
sign
⁡
{
d
Re
(
λ
(
τ
)
)
d
τ
}
|
λ
=
i
w
0
=
sign
⁡
{
Re
(
d
λ
(
τ
)
d
τ
)
}
|
λ
=
i
w
0
=
sign
⁡
{
i
B
w
0
e
-
i
w
0
τ
0
2
i
w
0
+
A
-
B
τ
0
e
-
i
w
0
τ
0
}
=
sign
⁡
{
w
0
2
(
A
2
+
w
0
2
)
(
A
-
w
0
2
τ
0
)
2
+
w
0
2
(
A
τ
0
+
2
)
2
}
=
1
.This implies that all of the roots cross the imaginary axis ati
w from left to right as τ increases. Therefore, the transversality condition holds.Based on the above, we have the following theorem.Theorem 2.
The interior equilibriumE
*
(
N
*
,
P
*
) exists only if N
0
>
N
*, and it is stable for τ
∈
[
0
,
τ
0
) but unstable for τ
>
τ
0. Furthermore, system (2) undergoes a Hopf bifurcation at the positive equilibrium E
*  for τ
=
τ
k, k
=
0,1
,
2
,
….


## 2.1. Local Stability of EquilibriumE
0

We linearize (2) about E
0 to obtain the linear system
(5)
d
N
d
t
=
-
k
h
N
-
γ
β
(
N
0
-
ε
m
)
P
,
d
P
d
t
=
(
β
N
0
-
m
-
v
h
-
k
h
)
P
.There are two eigenvalues where-
k
/
h is always negative and β
N
0
-
m
-
v
/
h
-
k
/
h
>
0 only if N
0
>
N
*. Therefore, we have the following theorem.Theorem 1.
The equilibriumE
0
(
N
0
,
0
) is stable if β
N
0
-
m
-
v
/
h
-
k
/
h
<
0 and unstable if β
N
0
-
m
-
v
/
h
-
k
/
h
>
0.
We recall thatβ
N
0
-
m
-
v
/
h
-
k
/
h
>
0 is equivalent to N
0
>
N
*, which is the existing condition of the unique interior equilibrium. Thus, E
0 is stable only if E
* does not exist.


## 2.2. Local Stability and the Hopf Bifurcation of EquilibriumE
*

Similarly, we linearize (2) about E
* to obtain the linear system
(6)
d
N
d
t
=
(
-
γ
β
P
*
-
k
h
)
N
-
γ
(
β
N
*
-
ε
m
)
P
,
d
P
d
t
=
β
P
*
N
(
t
-
τ
)
.The corresponding characteristic equation is(7)
λ
2
+
(
γ
β
P
*
+
k
h
)
λ
+
γ
β
(
β
N
*
-
ε
m
)
P
*
e
-
λ
τ
=
0
or
(8)
λ
2
+
A
λ
+
B
e
-
λ
τ
=
0
,
where
(9)
A
=
(
γ
β
P
*
+
k
h
)
>
0
,
B
=
γ
β
(
β
N
*
-
ε
m
)
P
*
>
0
.Whenτ
=
0, the characteristic equation is
(10)
λ
2
+
A
λ
+
B
=
0
,
and the two eigenvalues satisfy
(11)
λ
1
+
λ
2
=
-
A
<
0
,
λ
1
·
λ
2
=
B
>
0
,
which indicates that E
* is locally asymptotically stable when τ
=
0.Whenτ
≠
0, we assume that λ
=
μ
(
τ
)
+
i
w
(
τ
) is a root of (8). Substituting this in (8) we have
(12)
μ
2
-
w
2
+
2
μ
w
i
+
μ
A
+
A
w
i
+
B
e
-
μ
τ
(
cos
⁡
τ
w
-
i
sin
τ
w
)
=
0
or
(13)
μ
2
-
w
2
+
μ
A
=
-
B
e
-
μ
τ
cos
⁡
τ
w
,
2
μ
w
+
A
w
=
B
e
-
μ
τ
sin
τ
w
.Letμ
(
τ
)
=
0 and w
(
τ
)
=
w
0
>
0, and hence (13) is equal to
(14)
w
0
2
=
B
cos
⁡
τ
w
0
,
A
w
0
=
B
sin
τ
w
0
.Hence,(15)
w
0
4
+
A
2
w
0
2
-
B
2
=
0
.LetΔ
=
A
4
+
4
B
2
>
0, and then w
0
=
(
-
A
2
+
Δ
)
/
2.However, we can solve the first formula in (13) to obtain
(16)
τ
=
τ
k
=
1
w
0
(
arccos
w
0
2
B
+
2
k
π
)
,
k
=
0,1
,
2
,
…
.Substitutingλ
(
τ
) in (8) and taking the derivatives with respect to τ, we have
(17)
2
λ
(
τ
)
d
λ
(
τ
)
d
τ
+
A
d
λ
(
τ
)
d
τ
+
B
e
-
λ
τ
(
-
τ
d
λ
(
τ
)
d
τ
-
λ
(
τ
)
)
=
0
.Forτ
=
τ
0,
(18)
sign
⁡
{
d
Re
(
λ
(
τ
)
)
d
τ
}
|
λ
=
i
w
0
=
sign
⁡
{
Re
(
d
λ
(
τ
)
d
τ
)
}
|
λ
=
i
w
0
=
sign
⁡
{
i
B
w
0
e
-
i
w
0
τ
0
2
i
w
0
+
A
-
B
τ
0
e
-
i
w
0
τ
0
}
=
sign
⁡
{
w
0
2
(
A
2
+
w
0
2
)
(
A
-
w
0
2
τ
0
)
2
+
w
0
2
(
A
τ
0
+
2
)
2
}
=
1
.This implies that all of the roots cross the imaginary axis ati
w from left to right as τ increases. Therefore, the transversality condition holds.Based on the above, we have the following theorem.Theorem 2.
The interior equilibriumE
*
(
N
*
,
P
*
) exists only if N
0
>
N
*, and it is stable for τ
∈
[
0
,
τ
0
) but unstable for τ
>
τ
0. Furthermore, system (2) undergoes a Hopf bifurcation at the positive equilibrium E
*  for τ
=
τ
k, k
=
0,1
,
2
,
….


## 3. Direction and Stability of Hopf Bifurcations

In this section, we consider the direction, stability, and period of the periodic solutions from the steady state using the method introduced by Hassard et al. [28]. Let x
(
t
)
=
N
(
τ
t
)
-
N
*, y
(
t
)
=
P
(
τ
t
)
-
P
*, τ
=
μ
+
τ
k, u
=
(
x
,
y
)
T, and u
t
(
θ
)
=
u
(
t
+
θ
), θ
∈
[
-
1,0
], where μ is a real number.Then (2) can be written as
(19)
u
′
(
t
)
=
L
(
μ
)
u
t
+
f
(
μ
)
u
t
;
that is,
(20)
(
x
˙
(
t
)
y
˙
(
t
)
)
=
(
μ
+
τ
k
)
(
L
1
(
x
(
t
)
y
(
t
)
)
+
L
2
(
x
(
t
-
1
)
y
(
t
-
1
)
)
+
f
)
,
where
(21)
L
1
=
(
-
γ
β
P
*
-
k
h
-
γ
(
β
N
*
-
ε
m
)
0
0
)
,
L
2
=
(
0
0
β
P
*
0
)
,
f
=
f
(
μ
)
u
t
=
(
-
γ
β
x
(
t
)
y
(
t
)
β
x
(
t
-
1
)
y
(
t
)
)
.According to the Riesz representation theorem, there exists a matrixη
(
θ
,
μ
), which is bounded for θ
∈
[
-
1,0
] such that
(22)
L
(
μ
)
ϕ
=
∫
-
1
0
d
η
(
θ
,
μ
)
ϕ
(
θ
)
for
ϕ
∈
C
(
[
-
1,0
]
,
R
2
)
.In particular, we can select(23)
η
(
θ
,
μ
)
=
(
μ
+
τ
k
)
(
-
γ
β
P
*
-
k
h
-
γ
(
β
N
*
-
ε
m
)
0
0
)
δ
(
θ
)
η
(
θ
,
μ
)
=
-
(
μ
+
τ
k
)
(
0
0
β
P
*
0
)
δ
(
θ
+
1
)
,
where
(24)
δ
(
θ
)
=
{
1
,
θ
=
0
0
,
θ
≠
0
.Forϕ
∈
C
1
(
[
0,1
]
,
R
2
), we define
(25)
A
(
μ
)
ϕ
=
{
d
ϕ
d
θ
,
θ
∈
[
-
1,0
)
L
(
μ
)
ϕ
,
θ
=
0
,
F
(
μ
)
ϕ
=
{
0
,
θ
∈
[
-
1,0
)
f
(
μ
)
ϕ
,
θ
=
0
.Then, (19) can be rewritten as
(26)
u
˙
(
t
)
=
A
(
μ
)
u
t
+
F
(
μ
)
u
t
.Forψ
∈
C
1
(
[
0,1
]
,
(
R
2
)
*
), we define the adjoint operator A
* as
(27)
A
*
ψ
(
s
)
=
{
-
d
ψ
d
s
,
s
∈
(
0,1
]
∫
-
1
0
ψ
(
-
θ
)
d
η
(
θ
,
0
)
,
s
=
0
.Forϕ
∈
C
1
(
[
0,1
]
,
R
2
) and ψ
∈
C
1
(
[
0,1
]
,
(
R
2
)
*
), we define
(28)
〈
ψ
,
ϕ
〉
=
ψ
-
(
0
)
ϕ
(
0
)
-
∫
-
1
0
∫
ξ
=
0
θ
ψ
(
ξ
-
θ
)
d
η
(
θ
,
0
)
ϕ
(
ξ
)
d
ξ
.Suppose thatq
(
θ
) and q
*
(
s
) are the eigenvectors that correspond to i
w
0
τ
k and -
i
w
0
τ
k, respectively. Let
(29)
q
(
θ
)
=
(
1
C
)
e
i
w
0
τ
k
θ
.With the conditionA
q
(
0
)
=
i
w
0
τ
k
q
(
0
), that is,
(30)
τ
k
(
i
w
0
+
γ
β
P
*
+
k
h
γ
(
β
N
*
-
ε
m
)
-
β
P
*
e
-
i
w
0
τ
k
i
w
0
)
(
1
C
)
=
(
0
0
)
,
we can obtain
(31)
C
=
-
i
w
0
+
γ
β
P
*
+
k
/
h
γ
(
β
N
*
-
ε
m
)
.Similarly, let(32)
q
*
(
s
)
=
E
(
1
D
)
e
i
w
0
τ
k
s
and with
(33)
A
*
q
*
(
0
)
=
-
i
w
0
τ
k
q
*
(
0
)
,
that is,
(34)
τ
k
(
-
i
w
0
+
γ
β
P
*
+
k
h
-
β
P
*
e
-
i
w
0
τ
k
γ
(
β
N
*
-
ε
m
)
-
i
w
0
)
(
1
D
)
=
(
0
0
)
,
we can obtain
(35)
D
=
γ
(
β
N
*
-
ε
m
)
i
w
0
.Given the condition〈
q
*
,
q
〉
=
1, 〈
q
*
,
q
-
〉
=
0 and since
(36)
〈
q
*
,
q
〉
=
q
*
¯
(
0
)
q
(
0
)
〈
q
*
,
q
〉
=
-
∫
-
1
0
∫
ξ
=
0
θ
q
*
¯
(
ξ
-
θ
)
d
η
(
θ
)
q
(
ξ
)
d
ξ
〈
q
*
,
q
〉
=
E
-
(
1
D
-
)
(
1
C
)
〈
q
*
,
q
〉
=
-
∫
-
1
0
∫
ξ
=
0
θ
E
-
(
1
D
-
)
e
-
i
w
0
τ
k
(
ξ
-
θ
)
d
η
(
θ
)
(
1
C
)
e
i
w
0
τ
k
ξ
d
ξ
〈
q
*
,
q
〉
=
E
-
{
1
+
D
-
C
-
∫
-
1
0
θ
(
1
D
-
)
d
η
(
θ
)
(
1
C
)
e
i
w
0
τ
k
θ
}
〈
q
*
,
q
〉
=
E
-
{
1
+
D
-
C
+
τ
k
D
-
β
P
*
e
-
i
w
0
τ
k
}
,
we obtain
(37)
E
=
1
1
+
D
C
-
+
τ
k
D
β
P
*
e
i
w
0
τ
k
.In the following, we compute the coordinates to describe the center manifoldC
0 at μ
=
0 according to [28]. We define
(38)
z
(
t
)
=
〈
q
*
,
u
t
〉
,
W
(
t
,
θ
)
=
u
t
(
θ
)
-
2
Re
{
z
(
t
)
q
(
θ
)
}
.On the center manifoldC
0, we have
(39)
W
(
t
,
θ
)
=
W
(
z
(
t
)
,
z
-
(
t
)
,
θ
)
,
where
(40)
W
(
z
,
z
-
,
θ
)
=
W
20
(
θ
)
z
2
2
+
W
11
(
θ
)
z
z
-
+
W
02
(
θ
)
z
-
2
2
W
(
z
,
z
-
,
θ
)
=
+
W
30
(
θ
)
z
3
6
+
⋯
.Then,(41)
z
˙
(
t
)
=
〈
q
*
,
u
˙
(
t
)
〉
=
i
w
0
τ
k
z
(
t
)
z
˙
(
t
)
=
+
〈
q
*
(
θ
)
,
F
(
0
)
(
W
(
z
,
z
-
,
θ
)
+
2
Re
{
z
(
t
)
q
(
θ
)
}
)
〉
z
˙
(
t
)
=
i
w
0
τ
k
z
(
t
)
+
q
*
¯
(
0
)
f
(
0
)
z
˙
(
t
)
=
×
(
W
(
z
,
z
-
,
0
)
+
2
Re
{
z
(
t
)
q
(
0
)
}
)
z
˙
(
t
)
=
Δ
i
w
0
τ
k
z
(
t
)
+
g
(
z
,
z
-
)
,
where
(42)
g
(
z
,
z
-
)
=
q
*
¯
(
0
)
f
(
0
)
(
W
(
z
,
z
-
,
0
)
+
2
Re
{
z
(
t
)
q
(
0
)
}
)
g
(
z
,
z
-
)
=
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
g
21
z
2
z
-
2
+
⋯
.
Based on (26) and (38), we have
(43)
W
˙
=
u
˙
t
-
2
Re
{
z
˙
q
}
=
u
˙
t
-
z
˙
q
-
z
-
q
-
˙
W
˙
=
A
u
t
+
F
u
t
-
z
˙
q
-
z
-
˙
q
-
W
˙
=
A
(
W
+
2
Re
{
z
q
}
)
+
F
(
W
+
2
Re
{
z
q
}
)
W
˙
=
-
2
Re
{
g
(
z
,
z
-
)
q
(
θ
)
}
W
˙
=
{
A
W
-
2
Re
{
q
*
-
(
0
)
F
(
0
)
q
(
θ
)
}
,
-
1
≤
θ
<
0
A
W
-
2
Re
{
q
*
-
(
0
)
F
(
0
)
q
(
θ
)
}
+
F
(
0
)
,
θ
=
0
W
˙
=
Δ
A
W
+
H
(
z
,
z
-
,
θ
)
,
where
(44)
H
(
z
,
z
-
,
θ
)
=
H
20
(
θ
)
z
2
2
+
H
11
(
θ
)
z
z
-
+
H
02
(
θ
)
z
-
2
2
+
⋯
.Comparing the coefficients, we have(45)
(
2
i
w
0
τ
k
I
-
A
(
0
)
)
W
20
(
θ
)
=
H
20
(
θ
)
,
(46)
-
A
(
0
)
W
11
(
θ
)
=
H
11
(
θ
)
,
(47)
(
A
(
0
)
+
2
i
w
0
τ
k
I
)
W
02
(
θ
)
=
-
H
02
(
θ
)
.Note that(48)
q
(
0
)
=
(
1
C
)
,
q
*
(
0
)
=
E
(
1
D
)
,
x
(
t
)
=
z
+
z
-
+
W
20
(
1
)
(
0
)
z
2
2
+
W
11
(
1
)
(
0
)
z
z
-
x
(
t
)
=
+
W
02
(
1
)
(
0
)
z
-
2
2
+
⋯
,
y
(
t
)
=
C
z
+
C
-
z
-
+
W
20
(
2
)
(
0
)
z
2
2
+
W
11
(
2
)
(
0
)
z
z
-
y
(
t
)
=
+
W
02
(
2
)
(
0
)
z
-
2
2
+
⋯
,
x
(
t
-
1
)
=
e
-
i
w
0
τ
k
z
+
e
i
w
0
τ
k
z
-
+
W
20
(
1
)
(
-
1
)
z
2
2
x
(
t
-
1
)
=
+
W
11
(
1
)
(
-
1
)
z
z
-
+
W
02
(
1
)
(
-
1
)
z
-
2
2
+
⋯
,
F
(
0
)
=
(
-
τ
γ
β
x
(
t
)
y
(
t
)
τ
β
x
(
t
-
1
)
y
(
t
)
)
.Thus, according to the definition ofg
(
z
,
z
-
), we can easily obtain
(49)
g
20
=
2
τ
k
E
-
[
-
γ
β
C
+
D
-
β
C
e
-
i
w
0
τ
k
]
,
g
11
=
τ
k
E
-
[
-
γ
β
(
C
+
C
-
)
+
D
-
β
(
C
-
e
-
i
w
0
τ
k
+
C
e
i
w
0
τ
k
)
]
,
g
02
=
2
τ
k
E
-
[
-
γ
β
C
-
+
D
-
β
C
-
e
i
w
0
τ
k
]
,
g
21
=
2
τ
k
E
-
[
-
γ
β
(
W
11
(
2
)
(
0
)
+
W
11
(
1
)
(
0
)
2
+
C
W
11
(
1
)
(
0
)
1111111111111111
+
C
-
W
20
(
1
)
(
0
)
2
)
g
21
=
2
τ
k
E
-
1
+
D
-
β
(
e
-
i
w
0
τ
k
W
11
(
2
)
(
0
)
+
e
i
w
0
τ
k
W
20
(
2
)
(
0
)
2
g
21
=
2
τ
k
E
-
1111111
+
C
W
11
(
1
)
(
-
1
)
+
C
-
W
20
(
1
)
(
-
1
)
2
)
]
.In the final formula,W
20
(
0
), W
20
(
-
1
), W
11
(
0
), and W
11
(
-
1
) are still unknown. Thus, we need to compute W
20
(
θ
) and W
11
(
θ
) accurately, as follows.Consider(50)
H
(
z
,
z
-
,
θ
)
=
-
2
Re
{
q
-
*
(
0
)
F
(
0
)
q
(
θ
)
}
=
-
g
q
(
θ
)
-
g
-
q
-
(
θ
)
H
(
z
,
z
-
,
θ
)
=
-
(
g
20
z
2
2
+
g
11
z
z
-
+
g
02
z
-
2
2
+
⋯
)
q
(
θ
)
H
(
z
,
z
-
,
θ
)
=
-
(
g
-
20
z
-
2
2
+
g
-
11
z
z
-
+
g
-
02
z
2
2
+
⋯
)
q
-
(
θ
)
.Comparing the coefficients, we have(51)
H
20
(
θ
)
=
-
g
20
q
(
θ
)
-
g
-
02
q
-
(
θ
)
,
H
11
(
θ
)
=
-
g
11
q
(
θ
)
-
g
-
11
q
-
(
θ
)
.From (45), we can obtain
(52)
W
˙
20
(
θ
)
=
2
i
w
0
τ
k
W
20
(
θ
)
+
g
20
q
(
0
)
+
g
-
02
q
-
(
0
)
W
˙
20
(
θ
)
=
2
i
w
0
τ
k
W
20
(
θ
)
+
g
20
q
(
0
)
e
i
w
0
τ
k
θ
W
˙
20
(
θ
)
=
+
g
-
02
q
-
(
0
)
e
-
i
w
0
τ
k
θ
.Solving the formula, we have(53)
W
20
(
θ
)
=
i
g
20
q
(
0
)
w
0
τ
k
e
i
w
0
τ
k
θ
+
i
g
02
q
-
(
0
)
3
w
0
τ
k
e
-
i
w
0
τ
k
θ
+
E
1
e
2
i
w
0
τ
k
θ
.Using the same method, we can obtain(54)
W
11
(
θ
)
=
-
i
g
11
q
(
0
)
w
0
τ
k
e
i
w
0
τ
k
θ
+
i
g
-
11
q
-
(
0
)
w
0
τ
k
e
-
i
w
0
τ
k
θ
+
E
2
,
where E
1 and E
2 are both two-dimensional vectors.Forθ
=
0, H
(
z
,
z
-
,
0
)
=
-
2
Re
{
q
-
*
(
0
)
F
(
0
)
q
(
0
)
}
+
F
(
0
); thus
(55)
H
20
(
0
)
=
-
g
20
q
(
0
)
-
g
-
02
q
-
(
0
)
+
τ
β
(
-
γ
C
C
e
-
i
w
0
τ
k
)
,
H
11
(
0
)
=
-
g
11
q
(
0
)
-
g
-
11
q
-
(
0
)
+
τ
β
(
-
γ
C
+
C
-
C
e
i
w
0
τ
k
+
C
-
e
-
i
w
0
τ
k
)
.According to (45) and (46), we have
(56)
τ
k
L
1
W
20
(
0
)
-
τ
k
L
2
W
20
(
-
1
)
=
2
i
w
0
τ
k
W
20
(
0
)
-
H
20
(
0
)
,
(57)
τ
k
L
1
W
11
(
0
)
-
τ
k
L
2
W
11
(
-
1
)
=
-
H
11
(
0
)
.Substituting (56) and (57) into (53) and (54), respectively, we can obtain
(58)
E
1
=
1
τ
(
2
i
w
0
I
-
L
1
-
L
2
e
-
2
i
w
0
τ
k
)
-
1
E
1
=
×
i
g
-
02
e
i
w
0
τ
k
L
2
q
-
(
0
)
3
w
0
[
H
20
(
0
)
+
2
g
20
q
(
0
)
+
g
-
02
q
-
(
0
)
3
+
i
g
20
L
1
q
(
0
)
w
0
a
a
a
a
a
a
a
a
+
i
g
-
02
L
1
q
-
(
0
)
3
w
0
+
i
g
20
e
-
i
w
0
τ
k
L
2
q
(
0
)
w
0
a
a
a
a
a
a
a
a
+
i
g
-
02
e
i
w
0
τ
k
L
2
q
-
(
0
)
3
w
0
]
,
E
2
=
1
τ
(
-
L
1
-
L
2
)
-
1
E
2
=
×
[
H
11
(
0
)
-
i
g
11
L
1
q
(
0
)
w
0
+
i
g
-
11
L
1
q
-
(
0
)
w
0
111111111
-
i
g
11
e
-
i
w
0
τ
k
L
2
q
(
0
)
w
0
+
i
g
-
11
e
i
w
0
τ
k
L
2
q
-
(
0
)
w
0
]
.From the above, we already knowW
20
(
θ
) and W
11
(
θ
), and g
21 is also expressed using these parameters. Thus, we can compute the following values:
(59)
c
1
(
0
)
=
i
2
w
0
τ
k
(
g
20
g
11
-
2
|
g
11
|
2
-
|
g
02
|
2
3
)
+
g
21
2
,
μ
2
=
-
Re
(
c
1
(
0
)
)
Re
λ
′
(
τ
k
)
,
β
2
=
2
Re
(
c
1
(
0
)
)
,
T
2
=
-
Im
⁡
{
c
1
(
0
)
}
+
μ
2
Im
⁡
{
λ
′
(
τ
k
)
}
w
0
τ
k
.Theorem 3.
μ
2, β
2, and T
2 are defined above, so 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)
ifT
2
>
0
(
T
2
<
0
), then the period of the bifurcating periodic solution of system (2) increases (decreases).


## 4. Numerical Simulation

According to Pardo [12], if there exists an interior equilibrium in system (2) without a delay effect, then the equilibrium is always globally asymptotically stable. According to our study, however, the equilibrium is not always stable due to the delay effect. That is, an oscillation will occur in system (2) under some conditions. These conditions are obtained using an analytical technique. To verify the validity of the results, a series of numerical results is given in this section. In the following, the parameters are fixed, excluding N
0, h, and τ. The three parameters are taken as control parameters because we want to know how the concentration of nutrients below the layer and the thickness of the layer, but mainly the delay, affect the system.Consider(60)
γ
=
0.503
,
β
=
0.18
,
m
=
0.0087
,
v
=
0.036
,
k
=
1.152
,
ε
=
0.2
.In this study, we are interested in the interior equilibrium, though the interior does not always exist. Thus, the field where there exists an interior equilibrium is given whenN
0, h change, as shown in Figure 1(a). In Figure 1(a), both the interior equilibrium and boundary equilibrium exist in zone I (yellow), whereas there is only a boundary equilibrium in zone II (green), which is separated by the blue line. Based on the theoretical results, an oscillation will occur when the delay τ is larger than τ
0 (the critical value). The critical value τ
0 is determined by the parameters of system (2) except τ, and the relationship of τ
0 with respect to N
0, h is shown in Figure 1(b). In Figure 1(b), when h is fixed, τ
0 will decrease with the increase of N
0. If N
0 is fixed, τ
0 decreases first with the increase in h (0
<
h
<
30), before it decreases. It should be noted that τ
0 does not always exist when 0
<
h
<
30 and N
0 is much smaller, because of the reason illustrated in Figure 1(a).(a) The field where two equilibria exist; (b) the relationship ofτ
0 with N
0, h.
(a)
(b)Based on Figure1, the parameter N
0 is taken as 10 and h is set as 1, so τ
0 is obtained according to τ
0
=
(
1
/
w
0
) arccos (
w
0
2
/
B
). To investigate the effect of a delay on system (2), we set τ as 0 and 2 initially. The numerical solutions for phytoplankton are shown in Figure 2(a). In Figure 2(a), the solution eventually converges to the interior equilibrium, but the solution with the delay oscillates initially. The interior equilibrium is globally asymptotically stable in system (2) without delay, so the solution must converge to the interior equilibrium. If the solution with a delay converges to the interior equilibrium because τ
<
τ
0
=
3.45, an oscillation does not appear. In the following, we set τ as 0 and 4. The numerical solutions for phytoplankton are shown in Figure 2(b). Hopf bifurcation occurs at τ
=
τ
0, so an oscillation appears because τ
>
τ
0
=
3.45. From Figure 2(b), it is obvious that an oscillation occurs.Solution for phytoplankton forN
0
=
10, h
=
1, and (a) τ
=
2; (b) τ
=
4.
(a)
(b)Furthermore, to investigate the relationship between nutrients and phytoplankton, the numerical solutions for nutrients and phytoplankton are shown in Figure3, where N
0
=
10, h
=
1, and τ
=
3.6. Figure 3(a) shows that the solutions oscillate. Thus, the relationship between nutrients and phytoplankton is mutually constrained. The phase of system (2) is shown in Figure 3(b).(a) The solutions for nutrients and phytoplankton forN
0
=
10, h
=
1, and τ
=
3.6; (b) phase diagram.
(a)
(b)To determine howN
0, h, and τ affect system (2), the bifurcation of the stable state is considered and Figure 4 shows the results. In Figure 4(a), the parameters N
0, h are fixed, and the delay τ varies. The solid line denotes the stable state, the symbol “*” is the Hopf bifurcation point, and the dashed line represents the unstable state. There is only one stable state when τ is smaller than τ
0. However, the stable state becomes unstable when τ is larger than τ
0. Next, a periodic solution arises from the stable state, which is stable when τ
0
<
τ
<
5. It is obvious that a delay destroys the stability of the equilibrium. In addition, the delay τ is fixed, where τ
=
6. We consider the effects of N
0, h on system (2), as shown in Figure 4(b). The solid line denotes the stable state, the symbol “*” is the Hopf bifurcation point, and the dashed line represents the unstable state. The red line and the blue line both denote the equilibrium (phytoplankton), while the magenta line represents the minimum amount of phytoplankton. There are two Hopf bifurcation points when N
0
=
7 and 1
<
h
<
7, but there is only one when N
0
=
10 and N
0
=
12. When N
0
=
7, the equilibrium is stable at first, but the stable state becomes unstable when h reaches h
* and a Hopf bifurcation occurs, before a stable periodic solution emerges. When h reaches h
*
*, another Hopf bifurcation occurs and the stable periodic solution disappears, before the equilibrium becomes stable again. In particular, the minimum periodic solution decreases initially with the increase in h, before it increases. When N
0
=
10 and N
0
=
12, the minimum periodic solution increases with the increase in h and the unstable equilibrium becomes stable via a Hopf bifurcation; that is, the oscillation disappears.The stable state bifurcation: (a)N
0
=
10 and h
=
1 are fixed; (b) τ
=
6 is fixed.
(a)
(b)


## 5. Conclusion

In this study, we investigate a biological system and consider the effect of time delay. Our results show that the time delay plays a vital role in system (2). Model (1), which was proposed by Taylor et al. [10], includes nutrient recycling without any nutrient loss and vertical variation. Using a numerical simulation, Taylor discussed how the light intensity, diurnal variation in light, and diurnal variation in turbulence affected the system. This is a very realistic model and it is useful for studying the dynamic relationship between nutrients and phytoplankton. However, because it considered too many realistic factors, the model dynamics were very complicated and difficult to discern in the full model. Therefore, an approximated model was proposed, which considered a layer of phytoplankton growing over a pool of nutrients. The model was the same as (2) but it lacked the time delay. Similar to the analysis of system (1), the condition where an oscillation occurred was considered based on a theoretical analysis, including how factors such as the phytoplankton loss rate affected the system based on a numerical simulation. Most of this analysis was based on the results obtained using a numerical method.Pardo [12] focused on the approximated model in 2000. In a mathematical study, Pardo proved that if the interior equilibrium point exists, it is always globally asymptotically stable. The boundary equilibrium point is also globally asymptotically stable if the system has only one equilibrium point.In our study, we consider the model with delay as (2). Our theoretical analysis shows that the boundary equilibrium is still stable if N
0
<
N
*, but the unique positive equilibrium will lose its stability if the time delay exceeds a critical value and an oscillation occurs. Furthermore, we also prove that system (2) undergoes a Hopf bifurcation with a specific delay. These results are verified using numerical simulation based on published parameters [10] to ensure that they are realistic. Using a computer simulation, we find that the concentration of phytoplankton will increase to a large value within a short period of time because of the effect of delay, which is very similar to the phytoplankton blooms found in nature. Thus, the time delay may be a reason for phytoplankton blooms, but more research is needed to confirm this hypothesis.The limitations of our study are that we only consider a discrete delay in the phytoplankton increase and we assumed that the recycling of nutrients is instantaneous. A model with delayed recycling would be more complicated but more similar to reality. Moreover, the model we considered is only an approximate model based on (1) and we consider that β is a constant. Thus, we did not consider the influence of the light intensity, although the light intensity is depth independent. Therefore, further research is required to consider these aspects fully.



---

*Source: 471507-2014-03-12.xml*