# Dynamical Analysis of a Nitrogen-Phosphorus-Phytoplankton Model

Authors: Yunli Deng; Min Zhao; Hengguo Yu; Yi Wang
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

This paper presents a nitrogen-phosphorus-phytoplankton model in a water ecosystem. The main aim of this research is to analyze the global system dynamics and to study the existence and stability of equilibria. It is shown that the phytoplankton-eradication equilibrium is globally asymptotically stable if the input nitrogen concentration is less than a certain threshold. However, the coexistence equilibrium is globally asymptotically stable as long as it exists. The system is uniformly persistent within threshold values of certain key parameters. Finally, to verify the results, numerical simulations are provided.


---

## Body

## 1. Introduction

In marine ecology, phytoplankton play a major role in nutrient cycling, primary production, and global carbon cycling. As is well known, nutrient validity is a necessary factor in phytoplankton population growth. However, the phytoplankton population will breed massively when nutrient input is excessive, eventually leading to eutrophication. Eutrophication can cause water-quality deterioration and fish killing and affect people’s health and recreational activities [1]. It is characterized by frequently recurring algal blooms and reducing species diversity in water bodies at all trophic levels [2]. Therefore, reducing nutrients and eutrophication of water ecosystems is a crucial environmental problem throughout the world.In 1932, Bertalanffy first proposed the use of mathematical models to study biological systems [3]. Then, some scientific researchers attempted to investigate the biological population dynamics use of mathematical models [4–6]. With the increasing prevalence of eutrophication and algae blooms, the history of mathematical modeling of plankton dynamics and biological eutrophication removal processes is already quite long and has been initiated by the biological sciences. Some approaches have been refined to provide more realistic descriptions of the development of biological natural populations. For instance, ecological models, including impulsive [7–9], diffusion [10, 11], and time delay [12–14], have been taken into account, which can explain certain phenomena in realistic world. In recent years, most efforts have focused on how to control eutrophication, how to predict algae outbreaks, and how to simulate algae spreading tendencies. In this context, many researchers have discussed the dynamic behavior of phytoplankton blooms [15–22]. Huppert et al. [23] presented a simple nutrient-phytoplankton model and explored the dynamics of phytoplankton blooms. Pei et al. [24] investigated a two-zooplankton and one-phytoplankton model with harvesting, which considered the impact of harvesting on the coexistence and competitive exclusion of competitive predators. Mukhopadhyay and Bhattacharyya [25] dealt with a nutrient-plankton model in an aquatic environment in the context of phytoplankton blooms. Zhang and Wang [26] considered a nutrient-phytoplankton-zooplankton model in an aquatic environment and analyzed its global dynamics. Fan et al. [27] proposed a new dynamic nutrient-plankton model and used it to study the relationship between nutritional enrichment and water-quality oscillations. However, researchers have paid little attention to modeling nitrogen-phosphorus-phytoplankton systems.The Sanyang wetland of Wenzhou is located in a subtropical area and draws on the Wenzhou Economic Development zone and the Longwan zone to the east, linking up with the Chasan Higher Education zone to the south and connecting with the Wenzhou city center area to the northwest; its total area is 13 square kilometers. A certain number of rivers are distributed in a crosswise pattern, forming more than 160 different sizes and shapes of small islands; the proportion of land to water is 1.1 : 1. 47% of the land area is used to grow Mandarin oranges and 15.2% for housing, with the remainder used as agricultural land and fallow land. For the development of Wenzhou, the Sanyang wetland has played an important role in providing water resources, climate regulation, water conservation, flood and drought control, degradation of pollutants, and protection of biodiversity. However, the Sanyang wetland is facing the threat of industrial pollution and sewage, and its formerly crystal-clear water has become a large contaminated area, with local water areas colored red or black and emitting foul odors. Judging from the results of the overall analysis, the quality of the water environment in the Sanyang wetland has been severely damaged, with indicators of nitrogen, phosphorus, and heavy metals seriously exceeding limits. This situation has resulted in frequent nuisance algal blooms, which cause clogging and blocking of filtration systems. An even more frightening threat is that, with the economic development of Wenzhou city, the land and water bodies in the Sanyang wetland are in danger of being heavily invaded by industrial and building land uses as well as agricultural reclamation projects. Therefore, researching on how to enhance the protection of Sanyang wetland natural ecosystems and how to achieve a significant ecological effect on this environment is particularly important and urgent.Nitrogen and phosphorus are necessary nutrient for plants to live. When small quantities of nutrients flow into a wetland, phytoplankton start to grow. If the process is allowed to proceed, blooms will break out. Laukkanen and Huhtala [1] stated that nitrogen and phosphorus are the primary factors limiting algae blooms, and therefore these two nutrients are considered in the present model. This paper presents a nitrogen-phosphorus-phytoplankton model and uses it to study the interaction between nutrient runoff and phytoplankton growth.This paper considers a nitrogen-phosphorus-phytoplankton model with Holling type II functional response. The basic model is described by the following ordinary differential equations:(1)
d
x
d
t
=
I
1
-
b
1
x
-
a
1
x
z
c
1
+
x
,
d
y
d
t
=
I
2
-
b
2
y
-
a
2
y
z
c
2
+
y
,
d
z
d
t
=
e
1
a
1
x
z
c
1
+
x
+
e
2
a
2
y
z
c
2
+
y
-
m
z
,
where x and y are, respectively, the density of total nitrogen (mg/L) and total phosphorus (mg/L) at time t, z is the biomass of the phytoplankton population (mg/L) at time t, I
1, I
2 are, respectively, the runoff of nitrogen and phosphorus into the wetland, b
1, b
2 are, respectively, the natural removal of nitrogen and phosphorus from the water, a
1, a
2 are the maximum uptake rates, e
1, e
2 are biomass conversion constants, m is the natural death rate of the phytoplankton population, c
1, c
2 are half-saturation constants, and the terms a
1
x
z
/
(
c
1
+
x
) and a
2
y
z
/
(
c
2
+
y
) represent the response function for nutrient uptake by phytoplankton. The system satisfies the following initial conditions: x
(
0
)
=
x
0
≥
0, y
(
0
)
=
y
0
≥
0, z
(
0
)
=
z
0
≥
0.This paper aims to obtain a theoretical result for which the values of the bifurcation parameter, the phytoplankton-eradication equilibrium point, and the coexistence equilibrium point are asymptotically stable. To verify the results, numerical simulation has been used to study the controlling relations on the bifurcation parameterI
1.This paper is organized as follows. In the next section, the theorem governing the positivity and boundedness of solutions is analyzed. In Section3, the conditions for existence and stability of equilibria are obtained. In Section 4, the uniform persistence of system (1) is examined. Finally, numerical simulations are described in Section 5, and discussion and conclusions are given in Section 6.


## 2. Positivity and Boundedness of the Solution

In this section, assuming the biologically meaningful initial conditionsx
(
0
)
=
x
0
≥
0
,
y
(
0
)
=
y
0
≥
0, z
(
0
)
=
z
0
≥
0, the theorem for the positivity and boundedness of system (1) is stated and proved.Theorem 1.
Under the given initial conditions, all solutions of system (1) are positive and uniformly bounded.Proof.
From the first equation of system (1),
(2)
d
x
d
t
≥
-
b
1
x
-
a
1
x
z
c
1
+
x
.
Hence, x
(
t
)
>
x
(
0
)
exp
⁡
[
-
∫
0
t
(
b
1
+
(
a
1
z
/
(
c
1
+
x
)
)
)
d
s
]
>
0.
Similarly, from the second equation of system (1),
(3)
y
(
t
)
>
y
0
exp
⁡
-
∫
0
t
b
2
+
a
2
z
c
2
+
y
d
s
>
0
.
From the third equation of system (1), it can be obtained that
(4)
z
(
t
)
=
z
0
exp
⁡
∫
0
t
e
1
a
1
x
c
1
+
x
+
e
2
a
2
y
c
2
+
y
-
m
d
s
>
0
.
Now let us define a Lyapunov function: M
(
t
)
=
e
1
x
+
e
2
y
+
z.
Letd
=
min
⁡
{
b
1
,
b
2
,
m
}.
Then,(5)
d
M
d
t
+
d
M
≤
e
1
I
1
+
e
2
I
2
=
L
>
0
.
The right-hand side of the inequality is bounded for
all
(
x
,
y
,
z
)
∈
R
+
3, and therefore
(6)
M
x
,
y
,
z
<
L
d
1
-
e
-
d
t
+
M
x
0
,
y
0
,
z
0
e
-
d
t
.
Moreover, as t
→
+
∞
,
0
<
M
<
L
/
d. Hence, by the definition of M
(
t
), there are three positive constants M
x
,
M
y, and M
z and T
1
>
0 such that x
(
t
)
≤
M
x
,
y
(
t
)
≤
M
y
,
z
(
t
)
≤
M
z, for t
≥
T
1. This completes the proof.


## 3. Existence and Stability of Equilibria

In this section, the existence of all possible nonnegative equilibria is first discussed.Obviously, the phytoplankton-eradication equilibriumE
1
I
1
/
b
1
,
I
2
/
b
2
,
0 exists in system (1). Next, in order to research the coexistence of three populations, let us consider the existence of the positive equilibrium.Theorem 2.
There exists a coexistence equilibriumE
2
(
x
*
,
y
*
,
z
*
) of system (1) if
(7)
0
<
m
<
a
1
e
1
I
1
b
1
c
1
+
I
1
+
a
2
e
2
I
2
b
2
c
2
+
I
2
.Proof.
If the coexistence equilibriumE
2
(
x
*
,
y
*
,
z
*
) exists, it must satisfy the following three equations:
(8)
I
1
-
b
1
x
*
-
a
1
x
*
z
*
c
1
+
x
*
=
0
,
I
2
-
b
2
y
*
-
a
2
y
*
z
*
c
2
+
y
*
=
0
,
e
1
a
1
x
*
c
1
+
x
*
+
e
2
a
2
y
*
c
2
+
y
*
-
m
=
0
.
From (8), we must have x
*
∈
(
0
,
I
1
/
b
1
) and y
*
∈
(
0
,
I
2
/
b
2
) if the coexistence equilibrium exists.
From the first and the second equations of system (8), the following second degree equation in respect of y
* is obtained:
(9)
a
1
b
2
x
*
y
*
2
+
a
2
I
1
-
b
1
x
*
c
1
+
x
*
-
a
1
x
*
I
2
-
b
2
c
2
y
*
-
a
1
I
2
c
2
x
*
=
0
.
It is easy to calculate that
(10)
Δ
=
a
2
I
1
-
b
1
x
*
c
1
+
x
*
-
a
1
x
*
I
2
-
b
2
c
2
2
+
4
a
1
2
b
2
c
2
I
2
x
*
2
>
0
.
Hence, there exist two roots in (9)  y
1 and y
2, as follows:
(11)
y
1
=
-
f
x
*
+
f
2
x
*
+
4
a
1
2
b
2
c
2
I
2
2
a
1
b
2
,
y
2
=
-
f
x
*
-
f
2
x
*
+
4
a
1
2
b
2
c
2
I
2
2
a
1
b
2
,
where f
(
x
*
)
=
(
a
2
(
I
1
-
b
1
x
*
)
(
c
1
+
x
*
)
-
a
1
(
I
2
-
b
2
c
2
)
x
*
)
/
x
*. Obviously, y
1
>
0 and y
2
<
0 in terms of expression (11). Thus, we just need to take into consideration y
1 because of y
*
>
0. Let y
1
=
φ
1
(
x
*
). Then,
(12)
φ
1
′
x
*
=
f
′
x
*
f
x
*
-
f
2
x
*
+
4
a
1
2
b
2
c
2
I
2
2
a
1
b
2
f
2
x
*
+
4
a
1
2
b
2
c
2
I
2
.
Due tof
′
(
x
*
)
=
-
a
2
(
I
1
c
1
/
x
*
2
)
+
b
1
<
0, obviously, φ
1
′
(
x
*
)
>
0. Therefore, φ
1
(
x
*
) is a strictly increasing and continuous function of x
*
∈
(
0
,
I
1
/
b
1
). Moreover, li
m
x
*
→
0
+
φ
1
(
x
*
)
=
0. From (9), we have y
*
=
0 when x
*
=
0; thus, y
1
=
φ
1
(
x
*
) is continuous at the point (
0,0
). In addition, φ
1
(
I
1
/
b
1
)
=
I
2
/
b
2; hence, φ
1
(
x
*
) satisfy (9) for x
*
∈
0
,
I
1
/
b
1.
From the third equation of system (8), it can be obtained that
(13)
y
*
=
c
2
m
-
e
1
a
1
x
*
+
m
c
1
c
2
e
1
a
1
+
e
2
a
2
-
m
x
*
+
c
1
e
2
a
2
-
m
=
Δ
φ
2
x
*
.
Becauseφ
2
′
(
x
*
)
<
0, φ
2
(
x
*
) is a strictly decreasing and continuous function of x
*
∈
(
-
∞
,
c
1
(
m
-
e
2
a
2
)
/
(
e
1
a
1
+
e
2
a
2
-
m
)
)
∪
(
c
1
(
m
-
e
2
a
2
)
/
(
e
1
a
1
+
e
2
a
2
-
m
)
,
+
∞
).
Furthermore,(14)
φ
2
0
=
m
c
2
e
2
a
2
-
m
,
φ
2
m
c
1
e
1
a
1
-
m
=
0
.
The vertical asymptote is
(15)
x
*
=
c
1
m
-
e
2
a
2
e
1
a
1
+
e
2
a
2
-
m
.
In addition,
(16)
lim
x
→
∞
⁡
φ
2
x
*
=
c
2
m
-
e
1
a
1
e
1
a
1
+
e
2
a
2
-
m
.
Based on the analysis of (14), (15), and (16), we know that y
*
=
φ
2
(
x
*
) and y
1
=
φ
1
(
x
*
) may intersect only if x
*
∈
(
c
1
(
m
-
e
2
a
2
)
/
(
e
1
a
1
+
e
2
a
2
-
m
)
,
+
∞
).
LetΨ
(
x
*
,
y
*
)
=
y
*
-
φ
2
(
x
*
), where x
*
∈
(
c
1
(
m
-
e
2
a
2
)
/
(
e
1
a
1
+
e
2
a
2
-
m
)
,
+
∞
). Because Ψ
(
0,0
)
=
-
m
<
0, when the condition Ψ
(
I
1
/
b
1
,
I
2
/
b
2
)
≥
0 holds, then (
x
*
,
y
*
)
∣
y
*
=
φ
2
(
x
*
)
,
x
*
∈
(
c
1
(
m
-
e
2
a
2
)
/
(
e
1
a
1
+
e
2
a
2
-
m
)
,
+
∞
)
∩
(
x
*
,
y
*
)
∣
y
1
=
φ
1
(
x
*
)
,
x
*
∈
[
0
,
I
1
/
b
1
]
≠
⌀. Whereas the condition Ψ
(
I
1
/
b
1
,
I
2
/
b
2
)
=
0 is inconsistent with the coexistence equilibrium E
2
(
x
*
,
y
*
,
z
*
), the condition Ψ
(
I
1
/
b
1
,
I
2
/
b
2
)
>
0 holds.
Obviously,Ψ
(
I
1
/
b
1
,
I
2
/
b
2
)
>
0 if (7) holds. Thus, there must be unique x
*
,
y
*, and z
* such that the expression (8) holds, because φ
1
(
x
*
) is a strictly increasing continuous function and φ
2
(
x
*
) is a strictly decreasing continuous function, where y
*
=
φ
2
(
x
*
), z
*
=
(
I
1
-
b
1
x
*
)
(
c
1
+
x
*
)
/
(
a
1
x
*
), and x
* is defined as the unique solution of φ
1
(
x
*
)
=
φ
2
(
x
*
).
Therefore, when0
<
m
<
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
) holds, there exists a unique coexistence equilibrium E
2
(
x
*
,
y
*
,
z
*
) in system (1). The proof is complete.Next, the stability of the phytoplankton-eradication equilibriumE
1
x
1
,
y
1
,
0 will be discussed, where x
1
=
I
1
/
b
1 and y
1
=
I
2
/
b
2.Theorem 3.
The phytoplankton-eradication equilibriumE
1
x
1
,
y
1
,
0 is locally asymptotically stable if m
>
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
).Proof.
ForE
1
x
1
,
y
1
,
0, the characteristic equation is
(17)
-
b
1
-
λ
0
-
a
1
I
1
b
1
c
1
+
I
1
0
-
b
2
-
λ
-
a
2
I
2
b
2
c
2
+
I
2
0
0
a
1
e
1
I
1
b
1
c
1
+
I
1
+
a
2
e
2
I
2
b
2
c
2
+
I
2
-
m
-
λ
=
0
.
The roots of the characteristics equation are
(18)
λ
1
=
-
b
1
<
0
,
λ
2
=
-
b
2
<
0
,
λ
3
=
a
1
e
1
I
1
b
1
c
1
+
I
1
+
a
2
e
2
I
2
b
2
c
2
+
I
2
-
m
.
It is easy to verify that λ
3
<
0 if m
>
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
).
Therefore, the phytoplankton-eradication equilibriumE
1
x
1
,
y
1
,
0 is locally asymptotically stable if m
>
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
). The proof is complete.Theorem 4.
The phytoplankton-eradication equilibriumE
1
x
1
,
y
1
,
0 is globally asymptotically stable if m
>
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
).Proof.
To investigate the global stability of system (1) at E
1
x
1
,
y
1
,
0, let us construct the following Lyapunov function:
(19)
V
x
,
y
,
z
=
w
1
∫
x
1
x
f
(
ξ
)
-
f
(
x
1
)
f
(
ξ
)
d
ξ
+
w
2
∫
y
1
y
g
(
ξ
)
-
g
(
y
1
)
g
(
ξ
)
d
ξ
+
w
3
z
,
where w
1
,
w
2
,
w
3 are constants to be determined in subsequent steps and f
(
x
)
=
a
1
x
/
(
c
1
+
x
), g
(
y
)
=
a
2
y
/
(
c
2
+
y
).
Along the trajectories of system (1),
(20)
d
V
d
t
=
w
1
f
(
x
)
-
f
x
1
f
x
d
x
d
t
+
w
2
g
(
y
)
-
g
y
1
g
y
d
y
d
t
+
w
3
d
z
d
t
=
-
w
1
b
1
f
x
(
f
(
x
)
-
f
(
x
1
)
)
x
-
x
1
-
w
1
f
x
-
f
x
1
z
-
w
2
b
2
g
y
g
y
-
g
y
1
y
-
y
1
-
w
2
g
y
-
g
y
1
z
+
w
3
e
1
f
x
z
+
e
2
g
y
z
-
m
z
.
Now, choosing w
1
=
e
1
w
3 and w
2
=
e
2
w
3,
(21)
d
V
d
t
=
-
e
1
w
3
b
1
f
x
f
x
-
f
x
1
x
-
x
1
-
e
2
w
3
b
2
g
y
g
y
-
g
y
1
y
-
y
1
+
w
3
e
1
f
x
1
+
e
2
g
y
1
-
m
z
.
Note thatf
(
x
) and g
(
y
) are increasing functions and that the solution of (1) is positive. Meanwhile, it can be determined that e
1
f
(
x
1
)
+
e
2
g
(
y
1
)
-
m
<
0 is equal to m
>
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
). Therefore, each item on the right-hand side of (21) is nonpositive; that is, d
V
/
d
t
≤
0 and d
V
/
d
t
=
0 if and only if x
=
x
1, y
=
y
1, andz
=
0. Based on the Lyapunov-LaSalle theorem, if m
>
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
), then the phytoplankton-eradication equilibrium E
1
x
1
,
y
1
,
0 is globally asymptotically stable. Hence, the proof is complete.Finally, the stability of the coexistence equilibriumE
2
=
(
x
*
,
y
*
,
z
*
) will be discussed.Theorem 5.
The coexistence equilibriumE
2
=
(
x
*
,
y
*
,
z
*
) is locally asymptotically stable if 0
<
m
<
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
).Proof.
ForE
2
=
(
x
*
,
y
*
,
z
*
), the characteristic equation is
(22)
λ
+
b
1
+
a
1
c
1
z
*
c
1
+
x
*
2
0
a
1
x
*
c
1
+
x
*
0
λ
+
b
2
+
a
2
c
2
z
*
c
2
+
y
*
2
a
2
y
*
c
2
+
y
*
-
e
1
a
1
c
1
z
*
c
1
+
x
*
2
-
e
2
a
2
c
2
z
*
c
2
+
y
*
2
λ
=
0
.
That is,
(23)
λ
3
+
A
1
λ
2
+
A
2
λ
+
A
3
=
0
,
where
(24)
A
1
=
b
1
+
b
2
+
a
1
c
1
z
*
c
1
+
x
*
2
+
a
2
c
2
z
*
c
2
+
y
*
2
>
0
,
A
2
=
b
1
+
a
1
c
1
z
*
c
1
+
x
*
2
b
2
+
a
2
c
2
z
*
c
2
+
y
*
2
+
a
1
2
e
1
c
1
x
*
z
*
c
1
+
x
*
3
+
a
2
2
e
2
c
2
y
*
z
*
c
2
+
y
*
3
>
0
,
A
3
=
a
1
2
e
1
c
1
x
*
z
*
c
1
+
x
*
3
b
2
+
a
2
c
2
z
*
c
2
+
y
*
2
+
a
2
2
e
2
c
2
y
*
z
*
c
2
+
y
*
3
b
1
+
a
1
c
1
z
*
c
1
+
x
*
2
>
0
.
Note that A
1 and A
3 are both positive. Based on the Routh-Hurwitz criterion, E
2 is locally asymptotically stable if and only if A
1
A
2
-
A
3
>
0. In fact, after some computations,
(25)
A
1
A
2
-
A
3
=
b
1
+
a
1
c
1
z
*
c
1
+
x
*
2
b
2
+
a
2
c
2
z
*
c
2
+
y
*
2
×
b
1
+
b
2
+
a
1
c
1
z
*
c
1
+
x
*
2
+
a
2
c
2
z
*
c
2
+
y
*
2
+
a
1
2
e
1
c
1
x
*
z
*
c
1
+
x
*
3
b
1
+
a
1
c
1
z
*
c
1
+
x
*
2
+
a
2
2
e
2
c
2
y
*
z
*
c
2
+
y
*
3
b
2
+
a
2
c
2
z
*
c
2
+
y
*
2
>
0
.
Hence, the coexistence equilibrium E
2
=
(
x
*
,
y
*
,
z
*
) is locally asymptotically stable if E
2 exists. The proof is complete.Theorem 6.
The coexistence equilibriumE
2
=
(
x
*
,
y
*
,
z
*
) is globally asymptotically stable if 0
<
m
<
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
).Proof.
To investigate the global stability of system (1), at E
2
=
(
x
*
,
y
*
,
z
*
), let us construct the following Lyapunov function:
(26)
V
x
,
y
,
z
=
w
1
∫
x
*
x
f
(
ξ
)
-
f
x
*
f
ξ
d
ξ
+
w
2
∫
y
*
y
g
(
ξ
)
-
g
y
*
g
ξ
d
ξ
+
w
3
∫
z
*
z
ξ
-
z
*
ξ
d
ξ
,
where w
1
,
w
2
,
w
3 are constants to be determined in subsequent steps and f
(
x
)
=
a
1
x
/
(
c
1
+
x
), g
(
y
)
=
a
2
y
/
(
c
2
+
y
).
Along the trajectories of system (1),
(27)
d
V
d
t
=
w
1
f
x
-
f
x
*
f
x
-
b
1
x
-
x
*
+
f
x
*
z
*
-
f
x
z
+
w
2
g
y
-
g
y
*
g
y
×
-
b
2
y
-
y
*
+
g
y
*
z
*
-
g
y
z
+
w
3
z
-
z
*
×
e
1
f
x
-
f
x
*
+
e
2
g
y
-
g
y
*
.
Because f
(
x
)
z
-
f
(
x
*
)
z
*
=
(
f
(
x
)
-
f
(
x
*
)
)
z
*
+
f
(
x
)
(
z
-
z
*
), g
(
y
)
z
-
g
(
y
*
)
z
*
=
(
g
(
y
)
-
g
(
y
*
)
)
z
*+g
(
y
)
(
z
-
z
*
) and choosing w
1
=
e
1
w
3 and w
2
=
e
2
w
3, then
(28)
d
V
d
t
=
-
e
1
w
3
b
1
f
x
f
x
-
f
x
*
x
-
x
*
-
e
1
w
3
f
(
x
)
f
x
-
f
x
*
2
z
*
-
e
2
w
3
b
2
g
y
g
y
-
g
y
*
y
-
y
*
-
e
2
w
3
g
y
g
y
-
g
y
*
2
z
*
.
The proof of this is similar to the proof of the global asymptotical stability ofE
1. It can be established that d
V
/
d
t
<
0. Based on the Lyapunov-LaSalle theorem, the coexistence equilibrium E
2
=
(
x
*
,
y
*
,
z
*
) is globally asymptotically stable if E
2 exists. Hence, the proof is complete.


## 4. Uniform Persistence

This section will discuss the uniform persistence of system (1). As discussed below, threshold expressions of some key parameters can be obtained under the condition of all-species persistence. The first step is to state the following lemma.Lemma 7.
Ifa
>
0, b
>
0, (
d
x
(
t
)
/
d
t
)
≤
(
≥
)
b
-
a
x
(
t
), and x
(
0
)
>
0, then x
(
t
)
≤
(
≥
)
(
b
/
a
)  +  (
x
(
0
)
-
(
b
/
a
)
)
e
-
a
t for t
>
0.Theorem 8.
IfU
>
0, then system (1) is uniformly persistent.Proof.
From Theorem1, there are positive constants M
x, M
y, and M
z and T
1
>
0 such that x
(
t
)
≤
M
x, y
(
t
)
≤
M
y, z
(
t
)
≤
M
z for t
≥
T
1.
On the other hand, from the first equation of system (1), it can be determined that
(29)
d
x
t
d
t
≥
I
1
-
b
1
+
a
1
M
z
c
1
x
.
Noting that Lemma 7 is satisfied, then
(30)
x
t
≥
I
1
c
1
b
1
c
1
+
a
1
M
z
+
x
0
-
I
1
c
1
b
1
c
1
+
a
1
M
z
e
-
(
b
1
+
(
a
1
M
z
/
c
1
)
)
t
for
t
>
0
.
Therefore, there are positive constants T
2, ε
1, and m
x such that
(31)
x
t
≥
I
1
c
1
b
1
c
1
+
a
1
M
z
-
ε
1
=
m
x
for
t
≥
T
2
.
From the second equation of system (1), it follows that (
d
y
(
t
)
/
d
t
)
≥
I
2
-
(
b
2
+
(
a
2
M
z
/
c
2
)
)
y. Therefore, there are positive constants T
3, ε
2, and m
y such that y
(
t
)
≥
(
I
2
c
2
/
(
b
2
c
2
+
a
2
M
z
)
)
-
ε
2
=
m
y for t
≥
T
3.
From the third equation of system (1), it follows that (
d
z
(
t
)
/
d
t
)
≥
(
(
e
1
a
1
m
x
/
(
c
1
+
M
x
)
)
+
(
e
2
a
2
m
y
/
(
c
2
+
M
y
)
)
-
m
)
z. Therefore, it is evident that if U
=
(
e
1
a
1
m
x
/
(
c
1
+
M
x
)
)
+
(
e
2
a
2
m
y
/
(
c
2
+
M
y
)
)
-
m
>
0, there are two positive constants T
4 and m
z such that z
(
t
)
≥
m
z for t
≥
T
4.
LetT
=
max
⁡
{
T
1
,
T
2
,
T
3
,
T
4
}; then, it can be obtained that m
x
≤
x
(
t
)
≤
M
x, m
y
≤
y
(
t
)
≤
M
y, and m
z
≤
z
(
t
)
≤
M
z for t
≥
T. Hence, the proof is complete.


## 5. Simulation Analysis and Results

The stability of the phytoplankton-eradication equilibrium point and the coexistence equilibrium point of system (1) has been demonstrated in the previous section. In this section, numerical simulations will be used to analyze the dynamic behavior of system (1). I
1 is chosen as the bifurcation parameter. Table 1 provides the values of the other fixed parameters and their units, which have been obtained from previous studies [27–29].Table 1
Parameter values of system (1).
Symbol
Value
Units
a
1
0.4
day−1
a
2
0.5
day−1
b
1
0.4
day−1
b
2
0.5
day−1
c
1
2
mg/L
c
2
3
mg/L
e
1
0.6
Dimensionless
e
2
0.5
Dimensionless
m
0.1
day−1
I
2
0.3
mg/LBased on the previous analysis, we know that a unique coexistence equilibrium exists in system (1) if 0
<
m
<
(
a
1
e
1
I
1
/
(
b
1
c
1
+
I
1
)
)
+
(
a
2
e
2
I
2
/
(
b
2
c
2
+
I
2
)
). According to the parameter values in Table 1, it is easy to compute that
(32)
a
1
e
1
I
1
b
1
c
1
+
I
1
+
a
2
e
2
I
2
b
2
c
2
+
I
2
≈
0.12
.From Figure1, it is easy to see that the curves y
=
φ
2
(
x
) and y
=
φ
1
(
x
) intersect when m
=
0.05 and 0.1, while, if m
=
0.15 and 0.2, the curves y
=
φ
2
(
x
) and y
=
φ
1
(
x
) do not intersect which is shown in Figure 1. Therefore, Figure 1 verifies the correctness of the theoretical results.Figure 1
The illustration on the existence of the coexistence equilibrium. The gray closed area denotes the regionΩ
=
[
0
,
I
1
/
b
1
]
×
[
0
,
I
2
/
b
2
], the blue curve denotes the curve y
=
φ
1
(
x
), the green curves represent the curves y
=
φ
2
(
x
) with the change of m, m
=
0.05, m
=
0.1, m
=
0.15, and m
=
0.2, and the red diamond implies the crossover point of y
=
φ
2
(
x
) and y
=
φ
1
(
x
), where I
1
=
0.4 and the other parameters are the same as in Table 1.From Theorem3, the phytoplankton-eradication equilibrium E
1 is locally asymptotically stable if
(33)
m
>
a
1
e
1
I
1
b
1
c
1
+
I
1
+
a
2
e
2
I
2
b
2
c
2
+
I
2
.
Therefore, a condition has been obtained which is equivalent to inequality (33):
(34)
I
1
<
a
1
b
1
c
1
e
1
b
2
c
2
+
I
2
a
1
e
1
-
m
b
2
c
2
+
I
2
+
a
2
e
2
I
2
-
b
1
c
1
.
Based on the parameter values in Table 1, it is easy to obtain that
(35)
a
1
b
1
c
1
e
1
b
2
c
2
+
I
2
a
1
e
1
-
m
b
2
c
2
+
I
2
+
a
2
e
2
I
2
-
b
1
c
1
≈
0.2569
.It can be determined thatE
1
(
0.25,0.6,0
) is locally asymptotically stable if the input nitrogen concentration is 0.1 mg/L, as shown in Figure 2(a). Therefore, the conditions of Theorem 3 are satisfied. When the input nitrogen concentration reaches 0.3 mg/L, E
2
(
0.68,0.56,0.28
) exists and is locally asymptotically stable, as shown in Figure 2(b). Therefore, the conditions of Theorem 5 are satisfied.Phase diagram of system (1) under initial conditions (
x
(
0
)
,
y
(
0
)
,
z
(
0
)
)
=
(
0.1,0.1,0.1
). (a) shows that E
1
(
0.25,0.6,0
) is locally asymptotically stable if I
1
=
0.1 mg/L; (b) shows that E
2
(
0.68,0.56,0.28
) exists and that E
2 is locally asymptotically stable if I
1
=
0.3 mg/L. The other parameters are the same as in Table 1.
(a)
(b)In addition, Figure3 shows the bifurcation diagram in the I
1
-
z plane. This figure shows the dynamic behavior of system (1) for increasing input nitrogen concentration. A threshold for input nitrogen concentration can be seen to exist. When I
1 is less than this threshold, that is, I
1
<
0.2569 mg/L, only E
1 exists, and E
1 is globally asymptotically stable. When I
1 is greater than this threshold, that is, I
1
>
0.2569 mg/L, E
1 and E
2 exist, but E
1 is unstable and E
2 is globally asymptotically stable. These simulation results are in line with Theorems 4 and 6. Therefore, it is possible to speculate that E
2 is globally asymptotically stable if and only if E
1 is unstable.Figure 3
Bifurcation diagram of phytoplankton population with increasing input nitrogen concentration. The other parameters are the same as in Table1. The red solid line indicates that E
1 exists and is globally asymptotically stable. The blue dashed line implies that E
1 exists, but is unstable. The green solid line indicates that E
2 exists and is globally asymptotically stable.


## 6. Discussion

This paper has discussed the nitrogen-phosphorus-phytoplankton model. The model is simple because it assumes that phytoplankton is not affected by other environmental factors. The model is only an abstraction of real ecological phenomena; however, it generates many characteristics of these phenomena. Certain conditions for boundedness of solutions and existence and stability of equilibria have been obtained. Global stability of the system has been proved by constructing the Lyapunov function. Furthermore, the uniform persistence of the system has been analyzed.According to the theoretical analysis and numerical simulation results, the concentrations of input nitrogen and phosphorus are an important factor in the system dynamics. Using the input nitrogen concentrationI
1 as a bifurcation parameter, the relationship between phytoplankton and nitrogen input has been examined. Figures 2 and 3 show that if I
1 is less than a certain threshold, E
1 is globally asymptotically stable; when I
1 is greater than this threshold, E
1 is unstable and E
2 is globally asymptotically stable. Hence, it can be conjectured that E
2 is globally asymptotically stable if and only if E
1 is unstable. All these results are expected to be of significance in the study of dynamic ecosystem complexity.



---

*Source: 823026-2015-02-22.xml*