# Nonlinear Dynamics of a Toxin-Phytoplankton-Zooplankton System with Self- and Cross-Diffusion

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

---

## Abstract

A nonlinear system describing the interaction between toxin-producing phytoplankton and zooplankton was investigated analytically and numerically, where the system was represented by a couple of reaction-diffusion equations. We analyzed the effect of self- and cross-diffusion on the system. Some conditions for the local and global stability of the equilibrium were obtained based on the theoretical analysis. Furthermore, we found that the equilibrium lost its stability via Turing instability and patterns formation then occurred. In particular, the analysis indicated that cross-diffusion can play an important role in pattern formation. Subsequently, we performed a series of numerical simulations to further study the dynamics of the system, which demonstrated the rich dynamics induced by diffusion in the system. In addition, the numerical simulations indicated that the direction of cross-diffusion can influence the spatial distribution of the population and the population density. The numerical results agreed with the theoretical analysis. We hope that these results will prove useful in the study of toxic plankton systems.


---

## Body

## 1. Introduction

In marine ecosystems, most aquatic life relies on plankton, which comprises phytoplankton and zooplankton. Phytoplankton comprises most of the primary energy sources in aquatic food webs, and it accounts for a large proportion of the world’s fixed production. Phytoplankton is consumed by zooplankton, which provides food for fish and other aquatic animals. In fact, the phytoplankton can also render very useful service by producing a huge amount of oxygen for other animals after absorbing carbon dioxide from environments [1]. Thus, plankton forms the basis of all aquatic food chains and it has an essential role in the study of marine ecology [2, 3]. However, the biomass of phytoplankton population may be of rapid increase or almost equally rapid decrease, and this phenomenon of rapid change in phytoplankton population is called “bloom” at some fixed time. Because some phytoplankton can produce toxin and the accumulation of high biomass, some of these blooms are known as “harmful algal blooms” [4] which can have toxic effects on marine ecosystems or even human health, thereby causing great socioeconomic damage. Thus, research into bloom dynamics is widespread with a special emphasis on harmful algal blooms.In the past two decades, there have been major increases in harmful plankton blooms in aquatic ecosystems [5–7]. Studies have shown that there are at least eight different modes and mechanisms that allow harmful phytoplankton species to cause mortality, physiological impairment, or other negative in situ effects [8]. It is well known that the toxin-producing phytoplankton has important impacts on the growth of the zooplankton, and thus studies of marine plankton are ubiquitous and significant, while the dynamic behavior of interacting species in the marine is also a major topic [9]. However, researchers have paid less attention towards toxin-producing plankton blooms in recent years [10, 11], although some have recognized the role of toxin-producing phytoplankton in reducing the grazing pressure due to zooplankton [12–15]. In [14], a three-species model comprised a toxin-producing phytoplankton and zooplankton and fish population was studied by Upadhyay and Chattopadhyay, who employed Holling type I, II, and III functional responses to describe the liberation of the toxin, where they concluded that the toxin-producing phytoplankton made a stabilizing contribution to aquatic systems. In [15], the dynamical behaviors of toxin-producing phytoplankton and zooplankton were investigated, where the phytoplankton was divided into two groups, that is, susceptible and infected phytoplankton. Of course, given the global increase in harmful plankton blooms in aquatic ecosystems, the study of the effects of toxic phytoplankton has begun to become a growing concern in recent years [16–19].In marine ecosystems, the diffusion phenomenon of population exists widely, which can impact the distribution of population. Many factors can influence the diffusion, such as turbulence and foraging. Actually, it is indicated that the diffusion induced by mixing play an important role in the population dynamics of the phytoplankton [20]. And climate models predict that global warming will increase the stability of vertical stratification [21] and reduce vertical mixing [22], which means that climate change will affect the population dynamics. Hence, in order to understand how climate changes affect the population dynamics, the study of the diffusion mechanism of population becomes much more important. In this aspect, the reaction-diffusion equation may be a useful tool to study the spatiotemporal dynamics induced by diffusion [23, 24].The dynamics of interacting populations with self- and cross-diffusion have been studied widely [25–27]. In particular, Dubey et al. [25] analyzed a mathematical model of a predator-prey interaction with self- and cross-diffusion and obtained the criteria for local stability, instability, and global stability. Another study [26] considered a nutrient-plankton model of an aquatic environment in the context of phytoplankton blooms to analyze the diffusion-driven instability and stability, as well as cross-diffusion of zooplankton, under the influence of phytoplankton in the spatial model. The results [26] indicated the influence of cross-diffusion of zooplankton, and thus we can consider the effects of more complex cross-diffusion in marine ecosystems. In the previous work, some researchers have studied the dynamical properties of the toxic plankton system with diffusion [28–30]. Jang et al. studied a mathematic model which described the phytoplankton-zooplankton interactions with toxin-producing phytoplankton; the study revealed that passive diffusion of both populations can simplify the dynamics of the interactions and exhibit plankton patchiness [28]. In [29], authors studied a delay-diffusion model of marine plankton ecosystem, which exhibited cyclic nature of blooms. In [30], Roy considered a nontoxic phytoplankton-toxin-phytoplankton-zooplankton model which is described by a reaction-diffusion equation. The study demonstrated that spatial movements of planktonic systems in the presence of TPP generate and maintain inhomogeneous biomass distribution of competing phytoplankton, as well as grazer zooplankton. These works have studied the influence of diffusion on toxin plankton system and obtained some good results. However, the influence of cross-diffusion on the toxin plankton ecosystem was seldom considered. In order to study the influence of cross-diffusion in the toxic plankton ecosystem, we propose a toxin-phytoplankton-zooplankton system with self- and cross-diffusion.


## 2. Model and Linear Stability Analysis

### 2.1. Model Analysis

In this study, we consider a toxin-phytoplankton-zooplankton system with Holling type II response as follows:(1)
d
P
d
t
=
u
1
P
-
r
1
P
-
u
2
P
2
-
β
1
P
Z
K
+
P
,
d
Z
d
t
=
β
2
P
Z
K
+
P
-
v
1
Z
-
v
2
Z
2
-
β
3
P
Z
K
+
P
,where P
=
P
(
t
) is the density of the phytoplankton population and Z
=
Z
(
t
) is the density of the zooplankton population at time t. u
1 is the phytoplankton natural birth rate; r
1 is the phytoplankton natural mortality rate, and thus u
1
-
r
1 is the intrinsic growth rate of the phytoplankton. u
2 denotes phytoplankton intraspecies competition. P
Z
/
(
K
+
P
) is Holling type II response, K is the half-saturation constant for a Holling type II functional response, β
1 is the rate of predation for zooplankton, v
1 is the zooplankton’s natural mortality rate, v
2 is the zooplankton intraspecies competition, β
2 is the rate of biomass consumption by zooplankton during its growth, and β
3 denotes the rate of toxin liberation by the toxin-producing phytoplankton population. We make a basic assumption that β
2
>
β
3. The parameters u
1, r
1, u
2, K, β
1, β
2, v
1, v
2,  and β
3 are all positive constants.Considering the relationship between the climate and the diffusion of species and the fact on the existence of diffusion in population, system (1) is developed into a spatial system with diffusion. We expect to explore the effect of climate change phytoplankton population by studying the spatial dynamics of the diffusion system. The spatial system can be described as(2)
∂
P
∂
t
=
u
1
P
-
r
1
P
-
u
2
P
2
-
β
1
P
Z
K
+
P
+
D
11
Δ
P
+
D
12
Δ
Z
,
∂
Z
∂
t
=
β
2
P
Z
K
+
P
-
v
1
Z
-
v
2
Z
2
-
β
3
P
Z
K
+
P
+
D
21
Δ
P
+
D
22
Δ
Z
,where D
11 and D
22 are the constant diffusion coefficients for the phytoplankton and zooplankton, respectively, and D
12 and D
21 are the cross-diffusion coefficients for the phytoplankton and zooplankton, respectively. Biologically, cross-diffusion implies countertransport and it means that the prey exercised a self-defense mechanism to protect against attack by a predator [10], different from the self-diffusion, and the values of D
12 and D
21 may be positive or negative. We take the zero-flux boundary condition of the plankton populations P
(
s
,
t
) and Z
(
s
,
t
) in s
∈
Ω,  t
>
0 as follows:(3)
∂
P
∂
n
=
∂
Z
∂
n
=
0
,
s
∈
∂
Ω
,
t
>
0
,and the initial conditions are(4)
P
0
,
s
=
P
0
s
>
0
,
Z
0
,
s
=
Z
0
s
>
0
,
s
∈
Ω
,where Ω
⊆
R
2 is a bounded spatial domain with smooth boundary ∂
Ω.In addition, we take(5)
D
=
D
11
D
12
D
21
D
22as the diffusive matrix. The determinant of D is det
⁡
D
=
D
11
D
22
-
D
12
D
21. We assume that the self-diffusion is stronger than cross-diffusion; thus diffusive matrix D is positive definite.Letf
(
P
,
Z
) and g
(
P
,
Z
) be the right hand sides of the equations in (1). It is obvious that the equilibria E
0
=
(
0,0
) and E
1
=
(
(
u
1
-
r
1
)
/
u
2
,
0
) always exist.Byf
(
P
,
Z
)
=
0, we can obtain the vertical isoclines l
1
:
Z
=
-
u
2
/
β
1
P
2
+
K
-
u
1
-
r
1
/
u
2
P
-
u
1
-
r
1
K
/
u
2. This is a downward parabola, so the axis of symmetry is P
=
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, and the fellowship with the y-axis is Z
=
u
1
-
r
1
K
/
β
1. By g
(
P
,
Z
)
=
0, we can obtain the horizontal isoclines l
2
:
Z
=
(
β
2
-
β
3
)
P
/
v
2
(
K
+
P
)
-
v
1
/
v
2. It is easy to show that the line P
=
-
K is the asymptote of l
2. We can prove that l
2 increases monotonically when P
∈
-
K
,
+
∞, and the fellowship with the y-axis is Z
=
-
v
1
/
v
2. Based on the analysis above, we define(6)
F
P
=
-
u
2
β
1
P
2
+
K
-
u
1
-
r
1
u
2
P
-
u
1
-
r
1
u
2
K
,
G
P
=
β
2
-
β
3
P
v
2
K
+
P
-
v
1
v
2
,and the following conclusions can be established.Lemma 1.
(1) If(
u
1
-
r
1
)
/
u
2
>
K, u
2
(
u
1
-
r
1
/
u
2
+
K
)
2
/
4
β
1
<
(
β
2
-
β
3
)
(
(
u
1
-
r
1
)
/
u
2
-
K
)
/
v
2
(
(
u
1
-
r
1
)
/
u
2
+
K
)
-
v
1
/
v
2, then there exists one positive equilibrium point E
e
1
=
(
P
e
1
,
Z
e
1
) in system (1), where P
e
1
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K.
(2) If(
u
1
-
r
1
)
/
u
2
<
K
<
(
u
1
-
r
1
)
/
u
2
v
1
β
2
-
β
3
-
v
1, then there exists one positive equilibrium point E
e
2
=
(
P
e
2
,
Z
e
2
) in system (1); it is obvious that P
e
2
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K.
(3) If(
u
1
-
r
1
)
/
u
2
>
K and there at least exists P
0
∈
0
,
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K to guarantee F
(
P
0
)
<
G
(
P
0
) but F
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K
>
G
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then there exist three positive equilibria E
e
31
=
(
P
e
31
,
Z
e
31
), E
e
32
=
(
P
e
32
,
Z
e
32
), and E
e
33
=
(
P
e
33
,
Z
e
33
) in system (1), where P
e
31
<
P
e
32
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K,  P
e
33
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K.


### 2.2. Linear Stability Analysis

ForE
0
=
(
0,0
) and E
1
=
(
u
1
-
r
1
)
/
u
2
,
0, we can find that E
0
=
(
0,0
) is always unstable and that E
1
=
(
u
1
-
r
1
)
/
u
2
,
0 is locally asymptotically stable when v
1
>
(
β
2
-
β
1
)
u
1
-
r
1
/
(
u
2
K
+
u
1
-
r
1
).Now, letE
¯
=
(
P
¯
,
Z
¯
) be an arbitrary positive equilibrium, and thus E
¯
=
(
P
¯
,
Z
¯
) satisfies the equations f
(
P
¯
,
Z
¯
)
=
0 and g
(
P
¯
,
Z
¯
)
=
0. Then, the Jacobian matrix for E
¯
=
(
P
¯
,
Z
¯
) is given by(7)
J
E
¯
=
J
11
J
12
J
21
J
22
=
P
¯
1
K
+
P
¯
u
1
-
r
1
-
u
2
P
¯
-
u
2
-
β
1
P
¯
K
+
P
¯
β
2
-
β
3
K
Z
¯
K
+
P
¯
2
-
v
2
Z
¯
.Note that(8)
tr
⁡
J
E
¯
=
P
¯
1
K
+
P
¯
u
1
-
r
1
-
u
2
P
¯
-
u
2
-
v
2
Z
¯
,
det
⁡
J
E
¯
=
P
¯
1
K
+
P
¯
u
1
-
r
1
-
u
2
P
¯
-
u
2
v
2
Z
¯
+
β
2
-
β
3
K
Z
¯
K
+
P
¯
2
β
1
P
¯
K
+
P
¯
.From J
E
¯, we can find that if P
¯
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then J
11
<
0, and thus tr
⁡
(
J
E
¯
)
<
0 and det
⁡
(
J
E
¯
)
>
0; in this case, the equilibrium E
¯
=
(
P
¯
,
Z
¯
) is always locally asymptotically stable. But if P
¯
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then J
11
>
0.Next, we analyze the stability and instability of system (2). Let 0
=
u
0
<
u
1
<
u
2
<
⋯
<
∞ (lim
i
→
∞
⁡
u
i
=
∞) be the eigenvalues of -
∇
2 on Ω with the zero-flux boundary condition. Set(9)
X
≔
U
=
P
,
Z
∈
C
1
Ω
2
:
∂
P
∂
n
=
∂
Z
∂
n
=
0
,and consider X
=
⨁
i
=
0
∞
X
i, where X
i is the eigenspace that corresponds to u
i.In the condition where(
u
1
-
r
1
)
/
u
2
<
K
<
(
u
1
-
r
1
)
/
u
2
v
1
β
2
-
β
3
-
v
1, we obtain the following conclusions.Lemma 2.
IfD
11
=
D
12
=
D
21
=
D
22
=
0, then the equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is always locally asymptotically stable.Theorem 3.
IfD
11
≠
0, D
22
≠
0, D
12
≠
0, and D
21
≠
0, then the equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is uniformly asymptotically stable when D
12
>
0, D
21
<
0.Proof.
The linearization of system (2) at the positive equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is(10)
∂
∂
t
P
Z
=
Γ
P
Z
+
f
1
P
-
P
e
2
,
Z
-
Z
e
2
f
2
P
-
P
e
2
,
Z
-
Z
e
2
,
Γ
=
P
e
2
1
K
+
P
e
2
u
1
-
r
1
-
u
2
P
e
2
-
u
2
+
D
11
Δ
-
β
1
P
e
2
K
+
P
e
2
+
D
12
Δ
β
2
-
β
3
K
Z
e
2
K
+
P
e
2
2
+
D
21
Δ
-
v
2
Z
e
2
+
D
22
Δ
=
Γ
11
+
D
11
Δ
Γ
12
+
D
12
Δ
Γ
21
+
D
21
Δ
Γ
22
+
D
22
Δ
,where f
n
ε
1
,
ε
2
=
O
ε
1
2
+
ε
2
2 (n
=
1,2).
For eachi (i
=
0,1
,
2
,
…), X
i is invariant under the operator Γ, and λ
i is an eigenvalue of Γ on X
i if and only if λ
i is an eigenvalue of the matrix(11)
H
i
=
Γ
11
-
D
11
u
i
Γ
12
-
D
12
u
i
Γ
21
-
D
21
u
i
Γ
22
-
D
22
u
i
.Note that(12)
tr
⁡
H
i
=
Γ
11
+
Γ
22
-
D
11
+
D
22
u
i
,
det
⁡
H
i
=
Γ
11
Γ
22
-
Γ
12
Γ
21
+
D
11
D
22
u
i
2
-
u
i
D
22
Γ
11
+
D
11
Γ
22
-
D
12
D
21
u
i
2
+
Γ
12
D
21
+
Γ
21
D
12
u
i
.When D
12
>
0
,
D
21
<
0, then tr
⁡
H
i
<
0, det
⁡
H
i
>
0, and the two eigenvalues λ
i
+ and λ
i
- have negative real parts. For any i
≥
0, we find the following.
Iftr
⁡
H
i
2
-
4
det
⁡
H
i
≤
0, then Re
⁡
λ
i
±
=
1
/
2
tr
⁡
H
i
≤
(
1
/
2
)
(
Γ
11
+
Γ
22
)
<
0, and if tr
⁡
H
i
2
-
4
det
⁡
H
i
>
0, since tr
⁡
(
H
i
)
<
0 and det
⁡
H
i
>
0, then(13)
Re
⁡
λ
i
-
=
tr
⁡
H
i
-
tr
⁡
H
i
2
-
4
det
⁡
H
i
2
≤
1
2
tr
⁡
H
i
≤
1
2
Γ
11
+
Γ
22
<
0
,
Re
⁡
λ
i
+
=
tr
⁡
H
i
+
tr
⁡
H
i
2
-
4
det
⁡
H
i
2
=
2
det
⁡
H
i
tr
⁡
H
i
-
tr
⁡
H
i
2
-
4
det
⁡
H
i
≤
det
⁡
H
i
tr
⁡
H
i
<
C
,where C is independent of i.
Thus, a negative constantC exists, which is independent of i, such that Re
⁡
λ
i
±
<
C for any i. By referring to [31], we can prove that the positive equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) of system (1) is uniformly asymptotically stable when D
12
>
0, D
21
<
0.Theorem 4.
(1) IfD
11
=
D
12
=
D
21
=
D
22
=
0, then the equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is globally asymptotically stable when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K.
(2) IfD
11
≠
0, D
22
≠
0, D
12
≠
0, D
21
≠
0, then the positive equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is globally asymptotically stable when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K and A
12
2
<
4
A
11
A
22, where(14)
A
11
=
D
11
P
¯
P
2
,
A
22
=
D
22
β
1
K
+
P
¯
β
2
-
β
3
K
Z
¯
Z
2
,
A
12
=
-
D
12
P
¯
P
2
-
D
21
β
1
K
+
P
¯
β
2
-
β
3
K
Z
¯
Z
2
.Proof.
(1) Define a Lyapunov function(15)
V
1
P
,
Z
=
∫
P
¯
P
X
-
P
¯
X
d
X
+
β
1
K
+
P
¯
β
2
-
β
3
K
∫
Z
¯
Z
Y
-
Z
¯
Y
d
Y
.We note that V
1
(
P
,
Z
) is nonnegative and V
1
(
P
,
Z
)
=
0 if and only if (
P
(
t
)
,
Z
(
t
)
)
=
(
P
¯
,
Z
¯
).
Furthermore,(16)
d
V
1
d
t
=
P
-
P
¯
P
d
P
d
t
+
β
1
K
+
P
¯
β
2
-
β
3
K
Z
-
Z
¯
Z
d
Z
d
t
.By substituting the expressions for d
P
/
d
t and d
Z
/
d
t from system (1), we obtain(17)
d
V
1
d
t
=
P
-
P
¯
u
1
-
r
1
-
u
2
P
-
β
1
Z
K
+
P
+
Z
+
β
1
K
+
P
β
2
-
β
3
K
Z
-
Z
¯
β
2
-
β
3
P
K
+
P
-
v
1
-
v
2
Z
.Using the condition that (18)
u
1
-
r
1
-
u
2
P
¯
-
β
1
Z
¯
K
+
P
¯
=
0
,
β
2
-
β
3
P
¯
K
+
P
¯
-
v
1
-
v
2
Z
¯
=
0
,(17) can be simplified as(19)
d
V
1
d
t
=
-
u
2
-
β
1
Z
¯
K
+
P
¯
K
+
P
P
-
P
¯
2
-
β
1
v
2
K
+
P
¯
β
2
-
β
3
K
Z
-
Z
¯
2
,
-
(
β
1
v
2
K
+
P
¯
/
β
2
-
β
3
K
)
(
Z
-
Z
¯
)
2
<
0, and -
(
u
2
-
β
1
Z
¯
/
(
K
+
P
¯
)
(
K
+
P
)
)
P
-
P
¯
2
<
0 always occurs when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K.
Therefore,d
V
1
/
d
t
<
0 if u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K holds.Proof.
(2) Then, by referring to [32, 33], we choose the following Lyapunov function:(20)
V
2
=
∬
Ω
V
1
P
,
Z
d
Ω
,and by differentiating V
2 with respect to time t along the solutions of system (2), we can obtain (21)
d
V
2
d
t
=
∬
Ω
d
V
1
d
t
d
Ω
+
∬
Ω
D
11
Δ
P
+
D
12
Δ
Z
∂
V
1
∂
P
+
D
21
Δ
P
+
D
22
Δ
Z
∂
V
1
∂
Z
d
Ω
.Using Green’s first identity in the plane, we obtain(22)
d
V
2
d
t
=
∬
Ω
d
V
1
d
t
d
Ω
-
D
11
∬
Ω
∂
2
V
1
∂
P
2
∇
P
2
d
Ω
-
D
22
∬
Ω
∂
2
V
1
∂
Z
2
∇
Z
2
d
Ω
-
D
21
∬
Ω
∂
2
V
1
∂
Z
2
∇
P
∇
Z
d
Ω
-
D
12
∬
Ω
∂
2
V
1
∂
P
2
∇
P
∇
Z
d
Ω
,where(23)
∂
2
V
1
∂
P
2
=
P
¯
P
2
>
0
,
∂
2
V
1
∂
Z
2
=
β
1
K
+
P
¯
β
2
-
β
3
K
Z
¯
Z
2
>
0
.Based on the analysis above,(24)
u
2
>
β
1
Z
¯
K
+
P
¯
K
.
In the case of increasingD
11, D
22 to sufficiently large values and if D
12, D
21 meet the conditions(25)
A
12
2
<
4
A
11
A
22then d
V
1
/
d
t
<
0, and thus d
V
2
/
d
t
<
0. Therefore E
e
2
=
(
P
e
2
,
Z
e
2
) of system (2) is globally asymptotically stable when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K and A
12
2
<
4
A
11
A
22.Under the condition of(26)
u
1
-
r
1
u
2
>
K
,
u
2
u
1
-
r
1
/
u
2
+
K
2
4
β
1
<
β
2
-
β
3
u
1
-
r
1
/
u
2
-
K
v
2
u
1
-
r
1
/
u
2
+
K
-
v
1
v
2
,we have the following conclusions.Lemma 5.
IfD
11
=
D
12
=
D
21
=
D
22
=
0, then the equilibrium E
e
1
=
(
P
e
1
,
Z
e
1
) is locally asymptotically stable when (27)
F
P
∗
-
G
P
∗
>
0
,
P
∗
=
u
1
-
r
1
-
u
2
K
-
β
2
+
β
3
+
v
1
+
u
1
-
r
1
-
u
2
K
-
β
2
+
β
3
+
v
1
2
+
8
u
2
v
1
K
4
u
2
,and the equilibrium E
e
1
=
(
P
e
1
,
Z
e
1
) is unstable when F
(
P
∗
)
-
G
(
P
∗
)
<
0.Proof.
Define(28)
H
P
=
F
P
-
G
P
,
0
<
P
<
u
1
-
r
1
/
u
2
-
K
2
,under the condition of Lemma 1, we can find that H
′
(
P
)
<
0
(
0
<
P
<
(
(
u
1
-
r
1
)
/
u
2
-
K
)
/
2
).
Thus(29)
β
2
-
β
3
K
v
2
K
+
P
2
<
-
u
2
β
1
2
P
+
K
-
u
1
-
r
1
u
2
,
0
<
P
<
u
1
-
r
1
/
u
2
-
K
2
,and combined with det
⁡
(
E
e
1
), we can find that det
⁡
(
E
e
1
)
>
0, and(30)
tr
⁡
E
e
1
=
-
2
u
2
P
e
1
2
+
u
1
-
r
1
-
u
2
K
-
β
2
+
β
3
+
v
1
P
e
1
+
v
1
K
K
+
P
e
1
.Hence, if F
(
P
∗
)
-
G
(
P
∗
)
>
0, then P
e
1
>
P
∗, and in this case tr
⁡
(
E
e
1
)
<
0; thus, E
e
1
=
(
P
e
1
,
Z
e
1
) is stable. If F
(
P
∗
)
-
G
(
P
∗
)
<
0, then P
e
1
<
P
∗ and in this case, tr
⁡
(
E
e
1
)
>
0; thus E
e
1
=
(
P
e
1
,
Z
e
1
) is unstable.
By analyzing Lemma5, we can find that F
(
P
∗
)
-
G
(
P
∗
)
=
0 is the value where the plankton system populations bifurcate into periodic oscillation. By F
(
P
∗
)
=
G
(
P
∗
), then we can obtain(31)
β
1
=
β
¯
=
-
u
2
v
2
K
+
P
∗
P
∗
2
+
K
-
u
1
-
r
1
/
u
2
P
∗
-
u
1
-
r
1
/
u
2
K
β
2
-
β
3
-
v
1
P
∗
-
v
1
K
,which is a bifurcation value of the parameter β
2. Now, the characteristic equation of the E
e
1 reduces to(32)
λ
2
+
det
⁡
E
e
1
=
0
,where det
⁡
(
E
e
1
)
>
0, and the equation has purely imaginary roots. We take λ
=
a
+
b
i. If the stability change occurs at β
1
=
β
¯, it is obvious that a
(
β
¯
)
=
0, b
(
β
¯
)
≠
0, and using a
(
β
¯
)
=
0, b
(
β
¯
)
≠
0, we get d
a
/
d
β
1
β
1
=
β
¯
=
(
(
1
/
2
)
d
tr
⁡
E
e
1
/
d
β
1
)
β
1
=
β
¯
>
0. Thus, the transversality condition for a Hopf bifurcation is satisfied.Theorem 6.
If(
u
1
-
r
1
)
/
u
2
>
K, u
2
u
1
-
r
1
/
u
2
+
K
2
/
4
β
1
<
(
β
2
-
β
3
)
(
(
u
1
-
r
1
)
/
u
2
-
K
)
/
v
2
(
(
u
1
-
r
1
)
/
u
2
+
K
)
-
v
1
/
v
2, and F
P
∗
-
G
P
∗
>
0, then the criterion for the Turing instability of system (2) satisfies the following condition a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
>
4
det
⁡
(
D
)
det
⁡
(
J
E
e
1
). Consider(33)
J
E
e
1
=
a
11
a
12
a
21
a
22
i
s
t
h
e
J
a
c
o
b
i
a
n
m
a
t
r
i
x
a
t
E
e
1
=
P
e
1
,
Z
e
1
.Proof.
The linearized form of system (2) corresponding to the equilibrium E
e
1
(
P
e
1
,
Z
e
1
) is given by(34)
∂
p
∂
t
=
a
11
p
+
a
12
z
+
D
11
∂
2
p
∂
s
2
+
D
12
∂
2
z
∂
s
2
,
∂
z
∂
t
=
a
21
p
+
a
22
z
+
D
21
∂
2
p
∂
s
2
+
D
22
∂
2
z
∂
s
2
,where P
=
P
e
1
+
p,  Z
=
Z
e
1
+
z, and p
,
z are small perturbations in P
,
Z about the equilibrium E
e
1
=
(
P
e
1
,
Z
e
1
). We expand the solution of system (34) into a Fourier series:(35)
p
z
=
∑
k
c
k
1
c
k
2
e
λ
k
t
+
i
k
s
,where k is the wave number of the solution. By combining (34) and (35), we obtain(36)
λ
k
c
k
1
c
k
2
=
a
11
-
k
2
D
11
a
12
-
k
2
D
12
a
21
-
k
2
D
21
a
22
-
k
2
D
22
c
k
1
c
k
2
.Hence, we can obtain the characteristic equation as follows:(37)
λ
2
-
t
r
k
λ
k
+
Δ
k
=
0
,where(38)
t
r
k
=
a
11
+
a
22
-
k
2
D
11
+
D
22
,
Δ
k
=
d
e
t
J
E
e
1
-
k
2
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
+
k
4
D
11
D
22
-
D
12
D
21
.
In order to allow the system to achieve Turing instability, at least one of the following conditions must be satisfied:t
r
k
>
0 and Δ
k
<
0. However, it is obvious that t
r
k
<
0. Thus only the condition that Δ
k
<
0 exists can give rise to instability. Therefore, a necessary condition for the system to be unstable is that(39)
det
⁡
J
E
e
1
-
k
2
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
+
k
4
D
11
D
22
-
D
12
D
21
<
0
.We define (40)
F
k
2
=
k
4
D
11
D
22
-
D
12
D
21
-
k
2
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
+
det
⁡
J
E
e
1
.If we let k
m
i
n
2 be the corresponding value of k
2 for the minimum value of F
(
k
2
), then (41)
k
m
i
n
2
=
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
D
11
D
22
-
D
12
D
21
>
0
.Thus, the corresponding minimum value of F
(
k
2
) is (42)
F
k
m
i
n
2
=
det
⁡
J
E
e
1
-
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
4
det
⁡
D
.In addition, if we let F
(
k
m
i
n
2
)
<
0, the sufficient condition for the system to be unstable reduces to(43)
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
>
4
det
⁡
D
det
⁡
J
E
e
1
.
Noting. If there exists three positive equilibria E
e
31
=
(
P
e
31
,
Z
e
31
), E
e
32
=
(
P
e
32
,
Z
e
32
), and E
e
33
=
(
P
e
33
,
Z
e
33
) in system (1), where P
e
31
<
P
e
32
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, P
e
33
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then a bistable system exists when(44)
P
e
31
>
u
1
-
r
1
-
u
2
K
-
β
1
+
v
1
+
u
1
-
r
1
-
u
2
K
-
β
1
+
v
1
2
+
8
u
2
v
1
K
4
u
2
.
We takeu
1
=
0.9, r
1
=
0.2, u
2
=
0.14, β
1
=
1.2, β
2
=
0.7, β
3
=
0.1, v
1
=
0.05, v
2
=
0.5. Table 1 show that the condition for Lemma 1(3) is achievable.Table 1
The condition for three positive equilibria exists in system.
K
Equilibrium
t
r
J
det
⁡
J
Stability
0.5
(0.7441, 0.6177)
<0
>0
Locally stable
(1, 0.7)
∖
<0
Unstable
(2.2559, 0.8823)
<0
>0
Locally stable


## 2.1. Model Analysis

In this study, we consider a toxin-phytoplankton-zooplankton system with Holling type II response as follows:(1)
d
P
d
t
=
u
1
P
-
r
1
P
-
u
2
P
2
-
β
1
P
Z
K
+
P
,
d
Z
d
t
=
β
2
P
Z
K
+
P
-
v
1
Z
-
v
2
Z
2
-
β
3
P
Z
K
+
P
,where P
=
P
(
t
) is the density of the phytoplankton population and Z
=
Z
(
t
) is the density of the zooplankton population at time t. u
1 is the phytoplankton natural birth rate; r
1 is the phytoplankton natural mortality rate, and thus u
1
-
r
1 is the intrinsic growth rate of the phytoplankton. u
2 denotes phytoplankton intraspecies competition. P
Z
/
(
K
+
P
) is Holling type II response, K is the half-saturation constant for a Holling type II functional response, β
1 is the rate of predation for zooplankton, v
1 is the zooplankton’s natural mortality rate, v
2 is the zooplankton intraspecies competition, β
2 is the rate of biomass consumption by zooplankton during its growth, and β
3 denotes the rate of toxin liberation by the toxin-producing phytoplankton population. We make a basic assumption that β
2
>
β
3. The parameters u
1, r
1, u
2, K, β
1, β
2, v
1, v
2,  and β
3 are all positive constants.Considering the relationship between the climate and the diffusion of species and the fact on the existence of diffusion in population, system (1) is developed into a spatial system with diffusion. We expect to explore the effect of climate change phytoplankton population by studying the spatial dynamics of the diffusion system. The spatial system can be described as(2)
∂
P
∂
t
=
u
1
P
-
r
1
P
-
u
2
P
2
-
β
1
P
Z
K
+
P
+
D
11
Δ
P
+
D
12
Δ
Z
,
∂
Z
∂
t
=
β
2
P
Z
K
+
P
-
v
1
Z
-
v
2
Z
2
-
β
3
P
Z
K
+
P
+
D
21
Δ
P
+
D
22
Δ
Z
,where D
11 and D
22 are the constant diffusion coefficients for the phytoplankton and zooplankton, respectively, and D
12 and D
21 are the cross-diffusion coefficients for the phytoplankton and zooplankton, respectively. Biologically, cross-diffusion implies countertransport and it means that the prey exercised a self-defense mechanism to protect against attack by a predator [10], different from the self-diffusion, and the values of D
12 and D
21 may be positive or negative. We take the zero-flux boundary condition of the plankton populations P
(
s
,
t
) and Z
(
s
,
t
) in s
∈
Ω,  t
>
0 as follows:(3)
∂
P
∂
n
=
∂
Z
∂
n
=
0
,
s
∈
∂
Ω
,
t
>
0
,and the initial conditions are(4)
P
0
,
s
=
P
0
s
>
0
,
Z
0
,
s
=
Z
0
s
>
0
,
s
∈
Ω
,where Ω
⊆
R
2 is a bounded spatial domain with smooth boundary ∂
Ω.In addition, we take(5)
D
=
D
11
D
12
D
21
D
22as the diffusive matrix. The determinant of D is det
⁡
D
=
D
11
D
22
-
D
12
D
21. We assume that the self-diffusion is stronger than cross-diffusion; thus diffusive matrix D is positive definite.Letf
(
P
,
Z
) and g
(
P
,
Z
) be the right hand sides of the equations in (1). It is obvious that the equilibria E
0
=
(
0,0
) and E
1
=
(
(
u
1
-
r
1
)
/
u
2
,
0
) always exist.Byf
(
P
,
Z
)
=
0, we can obtain the vertical isoclines l
1
:
Z
=
-
u
2
/
β
1
P
2
+
K
-
u
1
-
r
1
/
u
2
P
-
u
1
-
r
1
K
/
u
2. This is a downward parabola, so the axis of symmetry is P
=
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, and the fellowship with the y-axis is Z
=
u
1
-
r
1
K
/
β
1. By g
(
P
,
Z
)
=
0, we can obtain the horizontal isoclines l
2
:
Z
=
(
β
2
-
β
3
)
P
/
v
2
(
K
+
P
)
-
v
1
/
v
2. It is easy to show that the line P
=
-
K is the asymptote of l
2. We can prove that l
2 increases monotonically when P
∈
-
K
,
+
∞, and the fellowship with the y-axis is Z
=
-
v
1
/
v
2. Based on the analysis above, we define(6)
F
P
=
-
u
2
β
1
P
2
+
K
-
u
1
-
r
1
u
2
P
-
u
1
-
r
1
u
2
K
,
G
P
=
β
2
-
β
3
P
v
2
K
+
P
-
v
1
v
2
,and the following conclusions can be established.Lemma 1.
(1) If(
u
1
-
r
1
)
/
u
2
>
K, u
2
(
u
1
-
r
1
/
u
2
+
K
)
2
/
4
β
1
<
(
β
2
-
β
3
)
(
(
u
1
-
r
1
)
/
u
2
-
K
)
/
v
2
(
(
u
1
-
r
1
)
/
u
2
+
K
)
-
v
1
/
v
2, then there exists one positive equilibrium point E
e
1
=
(
P
e
1
,
Z
e
1
) in system (1), where P
e
1
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K.
(2) If(
u
1
-
r
1
)
/
u
2
<
K
<
(
u
1
-
r
1
)
/
u
2
v
1
β
2
-
β
3
-
v
1, then there exists one positive equilibrium point E
e
2
=
(
P
e
2
,
Z
e
2
) in system (1); it is obvious that P
e
2
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K.
(3) If(
u
1
-
r
1
)
/
u
2
>
K and there at least exists P
0
∈
0
,
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K to guarantee F
(
P
0
)
<
G
(
P
0
) but F
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K
>
G
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then there exist three positive equilibria E
e
31
=
(
P
e
31
,
Z
e
31
), E
e
32
=
(
P
e
32
,
Z
e
32
), and E
e
33
=
(
P
e
33
,
Z
e
33
) in system (1), where P
e
31
<
P
e
32
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K,  P
e
33
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K.


## 2.2. Linear Stability Analysis

ForE
0
=
(
0,0
) and E
1
=
(
u
1
-
r
1
)
/
u
2
,
0, we can find that E
0
=
(
0,0
) is always unstable and that E
1
=
(
u
1
-
r
1
)
/
u
2
,
0 is locally asymptotically stable when v
1
>
(
β
2
-
β
1
)
u
1
-
r
1
/
(
u
2
K
+
u
1
-
r
1
).Now, letE
¯
=
(
P
¯
,
Z
¯
) be an arbitrary positive equilibrium, and thus E
¯
=
(
P
¯
,
Z
¯
) satisfies the equations f
(
P
¯
,
Z
¯
)
=
0 and g
(
P
¯
,
Z
¯
)
=
0. Then, the Jacobian matrix for E
¯
=
(
P
¯
,
Z
¯
) is given by(7)
J
E
¯
=
J
11
J
12
J
21
J
22
=
P
¯
1
K
+
P
¯
u
1
-
r
1
-
u
2
P
¯
-
u
2
-
β
1
P
¯
K
+
P
¯
β
2
-
β
3
K
Z
¯
K
+
P
¯
2
-
v
2
Z
¯
.Note that(8)
tr
⁡
J
E
¯
=
P
¯
1
K
+
P
¯
u
1
-
r
1
-
u
2
P
¯
-
u
2
-
v
2
Z
¯
,
det
⁡
J
E
¯
=
P
¯
1
K
+
P
¯
u
1
-
r
1
-
u
2
P
¯
-
u
2
v
2
Z
¯
+
β
2
-
β
3
K
Z
¯
K
+
P
¯
2
β
1
P
¯
K
+
P
¯
.From J
E
¯, we can find that if P
¯
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then J
11
<
0, and thus tr
⁡
(
J
E
¯
)
<
0 and det
⁡
(
J
E
¯
)
>
0; in this case, the equilibrium E
¯
=
(
P
¯
,
Z
¯
) is always locally asymptotically stable. But if P
¯
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then J
11
>
0.Next, we analyze the stability and instability of system (2). Let 0
=
u
0
<
u
1
<
u
2
<
⋯
<
∞ (lim
i
→
∞
⁡
u
i
=
∞) be the eigenvalues of -
∇
2 on Ω with the zero-flux boundary condition. Set(9)
X
≔
U
=
P
,
Z
∈
C
1
Ω
2
:
∂
P
∂
n
=
∂
Z
∂
n
=
0
,and consider X
=
⨁
i
=
0
∞
X
i, where X
i is the eigenspace that corresponds to u
i.In the condition where(
u
1
-
r
1
)
/
u
2
<
K
<
(
u
1
-
r
1
)
/
u
2
v
1
β
2
-
β
3
-
v
1, we obtain the following conclusions.Lemma 2.
IfD
11
=
D
12
=
D
21
=
D
22
=
0, then the equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is always locally asymptotically stable.Theorem 3.
IfD
11
≠
0, D
22
≠
0, D
12
≠
0, and D
21
≠
0, then the equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is uniformly asymptotically stable when D
12
>
0, D
21
<
0.Proof.
The linearization of system (2) at the positive equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is(10)
∂
∂
t
P
Z
=
Γ
P
Z
+
f
1
P
-
P
e
2
,
Z
-
Z
e
2
f
2
P
-
P
e
2
,
Z
-
Z
e
2
,
Γ
=
P
e
2
1
K
+
P
e
2
u
1
-
r
1
-
u
2
P
e
2
-
u
2
+
D
11
Δ
-
β
1
P
e
2
K
+
P
e
2
+
D
12
Δ
β
2
-
β
3
K
Z
e
2
K
+
P
e
2
2
+
D
21
Δ
-
v
2
Z
e
2
+
D
22
Δ
=
Γ
11
+
D
11
Δ
Γ
12
+
D
12
Δ
Γ
21
+
D
21
Δ
Γ
22
+
D
22
Δ
,where f
n
ε
1
,
ε
2
=
O
ε
1
2
+
ε
2
2 (n
=
1,2).
For eachi (i
=
0,1
,
2
,
…), X
i is invariant under the operator Γ, and λ
i is an eigenvalue of Γ on X
i if and only if λ
i is an eigenvalue of the matrix(11)
H
i
=
Γ
11
-
D
11
u
i
Γ
12
-
D
12
u
i
Γ
21
-
D
21
u
i
Γ
22
-
D
22
u
i
.Note that(12)
tr
⁡
H
i
=
Γ
11
+
Γ
22
-
D
11
+
D
22
u
i
,
det
⁡
H
i
=
Γ
11
Γ
22
-
Γ
12
Γ
21
+
D
11
D
22
u
i
2
-
u
i
D
22
Γ
11
+
D
11
Γ
22
-
D
12
D
21
u
i
2
+
Γ
12
D
21
+
Γ
21
D
12
u
i
.When D
12
>
0
,
D
21
<
0, then tr
⁡
H
i
<
0, det
⁡
H
i
>
0, and the two eigenvalues λ
i
+ and λ
i
- have negative real parts. For any i
≥
0, we find the following.
Iftr
⁡
H
i
2
-
4
det
⁡
H
i
≤
0, then Re
⁡
λ
i
±
=
1
/
2
tr
⁡
H
i
≤
(
1
/
2
)
(
Γ
11
+
Γ
22
)
<
0, and if tr
⁡
H
i
2
-
4
det
⁡
H
i
>
0, since tr
⁡
(
H
i
)
<
0 and det
⁡
H
i
>
0, then(13)
Re
⁡
λ
i
-
=
tr
⁡
H
i
-
tr
⁡
H
i
2
-
4
det
⁡
H
i
2
≤
1
2
tr
⁡
H
i
≤
1
2
Γ
11
+
Γ
22
<
0
,
Re
⁡
λ
i
+
=
tr
⁡
H
i
+
tr
⁡
H
i
2
-
4
det
⁡
H
i
2
=
2
det
⁡
H
i
tr
⁡
H
i
-
tr
⁡
H
i
2
-
4
det
⁡
H
i
≤
det
⁡
H
i
tr
⁡
H
i
<
C
,where C is independent of i.
Thus, a negative constantC exists, which is independent of i, such that Re
⁡
λ
i
±
<
C for any i. By referring to [31], we can prove that the positive equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) of system (1) is uniformly asymptotically stable when D
12
>
0, D
21
<
0.Theorem 4.
(1) IfD
11
=
D
12
=
D
21
=
D
22
=
0, then the equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is globally asymptotically stable when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K.
(2) IfD
11
≠
0, D
22
≠
0, D
12
≠
0, D
21
≠
0, then the positive equilibrium E
e
2
=
(
P
e
2
,
Z
e
2
) is globally asymptotically stable when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K and A
12
2
<
4
A
11
A
22, where(14)
A
11
=
D
11
P
¯
P
2
,
A
22
=
D
22
β
1
K
+
P
¯
β
2
-
β
3
K
Z
¯
Z
2
,
A
12
=
-
D
12
P
¯
P
2
-
D
21
β
1
K
+
P
¯
β
2
-
β
3
K
Z
¯
Z
2
.Proof.
(1) Define a Lyapunov function(15)
V
1
P
,
Z
=
∫
P
¯
P
X
-
P
¯
X
d
X
+
β
1
K
+
P
¯
β
2
-
β
3
K
∫
Z
¯
Z
Y
-
Z
¯
Y
d
Y
.We note that V
1
(
P
,
Z
) is nonnegative and V
1
(
P
,
Z
)
=
0 if and only if (
P
(
t
)
,
Z
(
t
)
)
=
(
P
¯
,
Z
¯
).
Furthermore,(16)
d
V
1
d
t
=
P
-
P
¯
P
d
P
d
t
+
β
1
K
+
P
¯
β
2
-
β
3
K
Z
-
Z
¯
Z
d
Z
d
t
.By substituting the expressions for d
P
/
d
t and d
Z
/
d
t from system (1), we obtain(17)
d
V
1
d
t
=
P
-
P
¯
u
1
-
r
1
-
u
2
P
-
β
1
Z
K
+
P
+
Z
+
β
1
K
+
P
β
2
-
β
3
K
Z
-
Z
¯
β
2
-
β
3
P
K
+
P
-
v
1
-
v
2
Z
.Using the condition that (18)
u
1
-
r
1
-
u
2
P
¯
-
β
1
Z
¯
K
+
P
¯
=
0
,
β
2
-
β
3
P
¯
K
+
P
¯
-
v
1
-
v
2
Z
¯
=
0
,(17) can be simplified as(19)
d
V
1
d
t
=
-
u
2
-
β
1
Z
¯
K
+
P
¯
K
+
P
P
-
P
¯
2
-
β
1
v
2
K
+
P
¯
β
2
-
β
3
K
Z
-
Z
¯
2
,
-
(
β
1
v
2
K
+
P
¯
/
β
2
-
β
3
K
)
(
Z
-
Z
¯
)
2
<
0, and -
(
u
2
-
β
1
Z
¯
/
(
K
+
P
¯
)
(
K
+
P
)
)
P
-
P
¯
2
<
0 always occurs when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K.
Therefore,d
V
1
/
d
t
<
0 if u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K holds.Proof.
(2) Then, by referring to [32, 33], we choose the following Lyapunov function:(20)
V
2
=
∬
Ω
V
1
P
,
Z
d
Ω
,and by differentiating V
2 with respect to time t along the solutions of system (2), we can obtain (21)
d
V
2
d
t
=
∬
Ω
d
V
1
d
t
d
Ω
+
∬
Ω
D
11
Δ
P
+
D
12
Δ
Z
∂
V
1
∂
P
+
D
21
Δ
P
+
D
22
Δ
Z
∂
V
1
∂
Z
d
Ω
.Using Green’s first identity in the plane, we obtain(22)
d
V
2
d
t
=
∬
Ω
d
V
1
d
t
d
Ω
-
D
11
∬
Ω
∂
2
V
1
∂
P
2
∇
P
2
d
Ω
-
D
22
∬
Ω
∂
2
V
1
∂
Z
2
∇
Z
2
d
Ω
-
D
21
∬
Ω
∂
2
V
1
∂
Z
2
∇
P
∇
Z
d
Ω
-
D
12
∬
Ω
∂
2
V
1
∂
P
2
∇
P
∇
Z
d
Ω
,where(23)
∂
2
V
1
∂
P
2
=
P
¯
P
2
>
0
,
∂
2
V
1
∂
Z
2
=
β
1
K
+
P
¯
β
2
-
β
3
K
Z
¯
Z
2
>
0
.Based on the analysis above,(24)
u
2
>
β
1
Z
¯
K
+
P
¯
K
.
In the case of increasingD
11, D
22 to sufficiently large values and if D
12, D
21 meet the conditions(25)
A
12
2
<
4
A
11
A
22then d
V
1
/
d
t
<
0, and thus d
V
2
/
d
t
<
0. Therefore E
e
2
=
(
P
e
2
,
Z
e
2
) of system (2) is globally asymptotically stable when u
2
>
β
1
Z
¯
/
(
K
+
P
¯
)
K and A
12
2
<
4
A
11
A
22.Under the condition of(26)
u
1
-
r
1
u
2
>
K
,
u
2
u
1
-
r
1
/
u
2
+
K
2
4
β
1
<
β
2
-
β
3
u
1
-
r
1
/
u
2
-
K
v
2
u
1
-
r
1
/
u
2
+
K
-
v
1
v
2
,we have the following conclusions.Lemma 5.
IfD
11
=
D
12
=
D
21
=
D
22
=
0, then the equilibrium E
e
1
=
(
P
e
1
,
Z
e
1
) is locally asymptotically stable when (27)
F
P
∗
-
G
P
∗
>
0
,
P
∗
=
u
1
-
r
1
-
u
2
K
-
β
2
+
β
3
+
v
1
+
u
1
-
r
1
-
u
2
K
-
β
2
+
β
3
+
v
1
2
+
8
u
2
v
1
K
4
u
2
,and the equilibrium E
e
1
=
(
P
e
1
,
Z
e
1
) is unstable when F
(
P
∗
)
-
G
(
P
∗
)
<
0.Proof.
Define(28)
H
P
=
F
P
-
G
P
,
0
<
P
<
u
1
-
r
1
/
u
2
-
K
2
,under the condition of Lemma 1, we can find that H
′
(
P
)
<
0
(
0
<
P
<
(
(
u
1
-
r
1
)
/
u
2
-
K
)
/
2
).
Thus(29)
β
2
-
β
3
K
v
2
K
+
P
2
<
-
u
2
β
1
2
P
+
K
-
u
1
-
r
1
u
2
,
0
<
P
<
u
1
-
r
1
/
u
2
-
K
2
,and combined with det
⁡
(
E
e
1
), we can find that det
⁡
(
E
e
1
)
>
0, and(30)
tr
⁡
E
e
1
=
-
2
u
2
P
e
1
2
+
u
1
-
r
1
-
u
2
K
-
β
2
+
β
3
+
v
1
P
e
1
+
v
1
K
K
+
P
e
1
.Hence, if F
(
P
∗
)
-
G
(
P
∗
)
>
0, then P
e
1
>
P
∗, and in this case tr
⁡
(
E
e
1
)
<
0; thus, E
e
1
=
(
P
e
1
,
Z
e
1
) is stable. If F
(
P
∗
)
-
G
(
P
∗
)
<
0, then P
e
1
<
P
∗ and in this case, tr
⁡
(
E
e
1
)
>
0; thus E
e
1
=
(
P
e
1
,
Z
e
1
) is unstable.
By analyzing Lemma5, we can find that F
(
P
∗
)
-
G
(
P
∗
)
=
0 is the value where the plankton system populations bifurcate into periodic oscillation. By F
(
P
∗
)
=
G
(
P
∗
), then we can obtain(31)
β
1
=
β
¯
=
-
u
2
v
2
K
+
P
∗
P
∗
2
+
K
-
u
1
-
r
1
/
u
2
P
∗
-
u
1
-
r
1
/
u
2
K
β
2
-
β
3
-
v
1
P
∗
-
v
1
K
,which is a bifurcation value of the parameter β
2. Now, the characteristic equation of the E
e
1 reduces to(32)
λ
2
+
det
⁡
E
e
1
=
0
,where det
⁡
(
E
e
1
)
>
0, and the equation has purely imaginary roots. We take λ
=
a
+
b
i. If the stability change occurs at β
1
=
β
¯, it is obvious that a
(
β
¯
)
=
0, b
(
β
¯
)
≠
0, and using a
(
β
¯
)
=
0, b
(
β
¯
)
≠
0, we get d
a
/
d
β
1
β
1
=
β
¯
=
(
(
1
/
2
)
d
tr
⁡
E
e
1
/
d
β
1
)
β
1
=
β
¯
>
0. Thus, the transversality condition for a Hopf bifurcation is satisfied.Theorem 6.
If(
u
1
-
r
1
)
/
u
2
>
K, u
2
u
1
-
r
1
/
u
2
+
K
2
/
4
β
1
<
(
β
2
-
β
3
)
(
(
u
1
-
r
1
)
/
u
2
-
K
)
/
v
2
(
(
u
1
-
r
1
)
/
u
2
+
K
)
-
v
1
/
v
2, and F
P
∗
-
G
P
∗
>
0, then the criterion for the Turing instability of system (2) satisfies the following condition a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
>
4
det
⁡
(
D
)
det
⁡
(
J
E
e
1
). Consider(33)
J
E
e
1
=
a
11
a
12
a
21
a
22
i
s
t
h
e
J
a
c
o
b
i
a
n
m
a
t
r
i
x
a
t
E
e
1
=
P
e
1
,
Z
e
1
.Proof.
The linearized form of system (2) corresponding to the equilibrium E
e
1
(
P
e
1
,
Z
e
1
) is given by(34)
∂
p
∂
t
=
a
11
p
+
a
12
z
+
D
11
∂
2
p
∂
s
2
+
D
12
∂
2
z
∂
s
2
,
∂
z
∂
t
=
a
21
p
+
a
22
z
+
D
21
∂
2
p
∂
s
2
+
D
22
∂
2
z
∂
s
2
,where P
=
P
e
1
+
p,  Z
=
Z
e
1
+
z, and p
,
z are small perturbations in P
,
Z about the equilibrium E
e
1
=
(
P
e
1
,
Z
e
1
). We expand the solution of system (34) into a Fourier series:(35)
p
z
=
∑
k
c
k
1
c
k
2
e
λ
k
t
+
i
k
s
,where k is the wave number of the solution. By combining (34) and (35), we obtain(36)
λ
k
c
k
1
c
k
2
=
a
11
-
k
2
D
11
a
12
-
k
2
D
12
a
21
-
k
2
D
21
a
22
-
k
2
D
22
c
k
1
c
k
2
.Hence, we can obtain the characteristic equation as follows:(37)
λ
2
-
t
r
k
λ
k
+
Δ
k
=
0
,where(38)
t
r
k
=
a
11
+
a
22
-
k
2
D
11
+
D
22
,
Δ
k
=
d
e
t
J
E
e
1
-
k
2
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
+
k
4
D
11
D
22
-
D
12
D
21
.
In order to allow the system to achieve Turing instability, at least one of the following conditions must be satisfied:t
r
k
>
0 and Δ
k
<
0. However, it is obvious that t
r
k
<
0. Thus only the condition that Δ
k
<
0 exists can give rise to instability. Therefore, a necessary condition for the system to be unstable is that(39)
det
⁡
J
E
e
1
-
k
2
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
+
k
4
D
11
D
22
-
D
12
D
21
<
0
.We define (40)
F
k
2
=
k
4
D
11
D
22
-
D
12
D
21
-
k
2
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
+
det
⁡
J
E
e
1
.If we let k
m
i
n
2 be the corresponding value of k
2 for the minimum value of F
(
k
2
), then (41)
k
m
i
n
2
=
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
D
11
D
22
-
D
12
D
21
>
0
.Thus, the corresponding minimum value of F
(
k
2
) is (42)
F
k
m
i
n
2
=
det
⁡
J
E
e
1
-
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
4
det
⁡
D
.In addition, if we let F
(
k
m
i
n
2
)
<
0, the sufficient condition for the system to be unstable reduces to(43)
a
11
D
22
+
a
22
D
11
-
D
12
a
21
-
D
21
a
12
2
>
4
det
⁡
D
det
⁡
J
E
e
1
.
Noting. If there exists three positive equilibria E
e
31
=
(
P
e
31
,
Z
e
31
), E
e
32
=
(
P
e
32
,
Z
e
32
), and E
e
33
=
(
P
e
33
,
Z
e
33
) in system (1), where P
e
31
<
P
e
32
<
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, P
e
33
>
(
1
/
2
)
(
u
1
-
r
1
)
/
u
2
-
K, then a bistable system exists when(44)
P
e
31
>
u
1
-
r
1
-
u
2
K
-
β
1
+
v
1
+
u
1
-
r
1
-
u
2
K
-
β
1
+
v
1
2
+
8
u
2
v
1
K
4
u
2
.
We takeu
1
=
0.9, r
1
=
0.2, u
2
=
0.14, β
1
=
1.2, β
2
=
0.7, β
3
=
0.1, v
1
=
0.05, v
2
=
0.5. Table 1 show that the condition for Lemma 1(3) is achievable.Table 1
The condition for three positive equilibria exists in system.
K
Equilibrium
t
r
J
det
⁡
J
Stability
0.5
(0.7441, 0.6177)
<0
>0
Locally stable
(1, 0.7)
∖
<0
Unstable
(2.2559, 0.8823)
<0
>0
Locally stable


## 3. Numerical Results

### 3.1. Numerical Analysis

In this section, system (2) is analyzed using the numerical technique to show its dynamic complexity. However, the existence of a positive equilibrium and its stability in system (1) are given first before this. Table 2 shows the existence of the equilibrium and its stability for different values of the parameter K, while the other parameters remain fixed: u
1
=
0.7
,
r
1
=
0.2,u
2
=
0.382,β
1
=
0.21,β
2
=
0.2428,β
3
=
0.1,v
1
=
0.09, and v
2
=
0.013
. Table 2 indicates that the positive equilibrium is unstable when K
=
0.1, but it is locally asymptotically stable when K
=
0.205 and global asymptotically stable when K
=
0.6.Table 2
The positive equilibrium and its stability.
K
Equilibrium
t
r
J
det
⁡
J
Stability
0.1
(0.2205, 0.6346)
>0
>0
Unstable
0.205
(0.5408, 1.0420)
<0
>0
Locally stable
0.6
(1.1915, 0.3826)
<0
>0
Global stableThe parameters are taken as follows:u
1
=
0.7, r
1
=
0.2, u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, and v
2
=
0.013.Figure2 presents a series of one-dimensional numerical solutions of system (2). Figure 2(a) shows that the solution of system (2) tends to a positive equilibrium, where D
21
=
D
12
=
0, which implies that positive equilibrium is stable. This means that the self-diffusion does not lead to the occurrence of instability. In fact, as shown in Figure 1, the self-diffusion does not cause instability in system (2) when D
21
=
D
12
=
0. However, instability may occur when the cross-diffusion coefficients satisfy some conditions.Figure 1
Bifurcation diagram withD
12 and D
21, where the parameters values are taken as K
=
0.205, u
1
=
0.7, r
1
=
0.2, u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, v
2
=
0.013, D
11
=
0.04, and D
22
=
3.5. The gray line represents D
11
a
22
+
D
22
a
11
-
D
12
a
21
-
D
21
a
12
=
0, the blue curve represents det
⁡
D
det
⁡
J
E
e
1
=
0, and the red curve represents D
11
a
22
+
D
22
a
11
-
D
12
a
21
-
D
21
a
12
-
2
det
⁡
D
det
⁡
J
E
e
1
=
0; thus I shows the zone where the instability induced by diffusion occurs. The remainder of this paper is organized as follows. In Section 2, we describe and analyze the toxin-phytoplankton-zooplankton system, and we derive the criteria for the stability and instability of a system with cross-diffusion. In Section 3, we present the results of numerical simulations and discuss their implications. Finally, we give our conclusions in Section 4.Figure 2
Numerical stability of system (2) with the following parameters: K
=
0.205, u
1
=
0.7, r
1
=
0.2, u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, v
2
=
0.013, D
11
=
0.04, and D
22
=
3.5. (a) D
21
=
0, D
12
=
0; (b): D
21
=
-
0.05, D
12
=
-
1.5; (c) D
21
=
0.5, D
12
=
0; (d) D
21
=
0.5, D
12
=
0.1; (e) D
21
=
1, D
12
=
-
0.5; (f) D
21
=
-
1.5, D
12
=
0.
(a)
(b)
(c)
(d)
(e)
(f)Therefore, different values ofD
12 and D
21 are employed to illustrate the instability induced by cross-diffusion. Figure 2(b) shows the phytoplankton density, where the values of D
12 and D
21 are negative, whereas the values of D
12 and D
21 are positive in Figure 2(d). It is obvious that the spatial distributions of phytoplankton in Figures 2(b) and 2(d) are markedly different. In ecology, positive cross-diffusion indicates that one species tends to move in the direction with a lower concentration of another species, whereas negative cross-diffusion denotes that the population tends to move in the direction with a higher concentration of another species [30]. In Figures 2(c) and 2(f), D
12
=
0, and D
21
=
0.5, -
1.5, respectively, which means that only the influence of zooplankton cross-diffusion is considered in the system. The spatial distributions of phytoplankton are not obviously different in Figures 2(c) and 2(f), but the amplitudes of the phytoplankton density differ, which means that the direction of zooplankton cross-diffusion can affect the phytoplankton density. In Figure 2(e), coefficient D
21 is positive and D
12 is negative, and it is not difficult to see that the spatial distribution of phytoplankton differs significantly from the other panels shown in Figure 2. These results demonstrate that cross-diffusion plays a significant role in the stability of the system, but they also show the effects of the various types of cross-diffusion that a plankton system may encounter.


### 3.2. Pattern Formation

In order to illustrate the spatial distribution of phytoplankton more clearly, we show the patterns formed by system (2) in a two-dimensional space in Figure 3.Figure 3
Patterns obtained with system (2) using the following parameters: K
=
0.205,u
1
=
0.7,r
1
=
0.2,u
2
=
0.382,β
1
=
0.21,β
2
=
0.2428,β
3
=
0.1,v
1
=
0.09,v
2
=
0.013,D
11
=
0.04, and D
22
=
3.5. (b) D
21
=
0.4, D
12
=
0.1; (c) D
21
=
0.91, D
12
=
-
0.5. Time steps: (a) 0, (b) 1000, and (c) 1000.
(a)
(b)
(c)Figure3 shows two different spatial distributions of phytoplankton. Figure 3(a) shows the initial spatial distributions of phytoplankton. Figures 3(b) and 3(c) illustrate different types of pattern formation for different values of D
12 and D
21 at t
=
1000. When D
21
=
0.4, D
12
=
0.1, Figure 3(b) shows that a spotted pattern and the spots appear to be cluttered with an irregular shape. However, when D
21
=
0.91, D
12
=
-
0.5, a circular spotted pattern appears, and the distribution of the spots is more regular, as shown in Figure 3(c). As discussed in Section 3.1, cross-diffusion will lead to instability in the system, where both positive cross-diffusion and negative cross-diffusion play important roles in the stability of the system. The spatial distributions of phytoplankton shown in Figures 2(b) and 2(c) agree well with the results given in the previous section.


## 3.1. Numerical Analysis

In this section, system (2) is analyzed using the numerical technique to show its dynamic complexity. However, the existence of a positive equilibrium and its stability in system (1) are given first before this. Table 2 shows the existence of the equilibrium and its stability for different values of the parameter K, while the other parameters remain fixed: u
1
=
0.7
,
r
1
=
0.2,u
2
=
0.382,β
1
=
0.21,β
2
=
0.2428,β
3
=
0.1,v
1
=
0.09, and v
2
=
0.013
. Table 2 indicates that the positive equilibrium is unstable when K
=
0.1, but it is locally asymptotically stable when K
=
0.205 and global asymptotically stable when K
=
0.6.Table 2
The positive equilibrium and its stability.
K
Equilibrium
t
r
J
det
⁡
J
Stability
0.1
(0.2205, 0.6346)
>0
>0
Unstable
0.205
(0.5408, 1.0420)
<0
>0
Locally stable
0.6
(1.1915, 0.3826)
<0
>0
Global stableThe parameters are taken as follows:u
1
=
0.7, r
1
=
0.2, u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, and v
2
=
0.013.Figure2 presents a series of one-dimensional numerical solutions of system (2). Figure 2(a) shows that the solution of system (2) tends to a positive equilibrium, where D
21
=
D
12
=
0, which implies that positive equilibrium is stable. This means that the self-diffusion does not lead to the occurrence of instability. In fact, as shown in Figure 1, the self-diffusion does not cause instability in system (2) when D
21
=
D
12
=
0. However, instability may occur when the cross-diffusion coefficients satisfy some conditions.Figure 1
Bifurcation diagram withD
12 and D
21, where the parameters values are taken as K
=
0.205, u
1
=
0.7, r
1
=
0.2, u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, v
2
=
0.013, D
11
=
0.04, and D
22
=
3.5. The gray line represents D
11
a
22
+
D
22
a
11
-
D
12
a
21
-
D
21
a
12
=
0, the blue curve represents det
⁡
D
det
⁡
J
E
e
1
=
0, and the red curve represents D
11
a
22
+
D
22
a
11
-
D
12
a
21
-
D
21
a
12
-
2
det
⁡
D
det
⁡
J
E
e
1
=
0; thus I shows the zone where the instability induced by diffusion occurs. The remainder of this paper is organized as follows. In Section 2, we describe and analyze the toxin-phytoplankton-zooplankton system, and we derive the criteria for the stability and instability of a system with cross-diffusion. In Section 3, we present the results of numerical simulations and discuss their implications. Finally, we give our conclusions in Section 4.Figure 2
Numerical stability of system (2) with the following parameters: K
=
0.205, u
1
=
0.7, r
1
=
0.2, u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, v
2
=
0.013, D
11
=
0.04, and D
22
=
3.5. (a) D
21
=
0, D
12
=
0; (b): D
21
=
-
0.05, D
12
=
-
1.5; (c) D
21
=
0.5, D
12
=
0; (d) D
21
=
0.5, D
12
=
0.1; (e) D
21
=
1, D
12
=
-
0.5; (f) D
21
=
-
1.5, D
12
=
0.
(a)
(b)
(c)
(d)
(e)
(f)Therefore, different values ofD
12 and D
21 are employed to illustrate the instability induced by cross-diffusion. Figure 2(b) shows the phytoplankton density, where the values of D
12 and D
21 are negative, whereas the values of D
12 and D
21 are positive in Figure 2(d). It is obvious that the spatial distributions of phytoplankton in Figures 2(b) and 2(d) are markedly different. In ecology, positive cross-diffusion indicates that one species tends to move in the direction with a lower concentration of another species, whereas negative cross-diffusion denotes that the population tends to move in the direction with a higher concentration of another species [30]. In Figures 2(c) and 2(f), D
12
=
0, and D
21
=
0.5, -
1.5, respectively, which means that only the influence of zooplankton cross-diffusion is considered in the system. The spatial distributions of phytoplankton are not obviously different in Figures 2(c) and 2(f), but the amplitudes of the phytoplankton density differ, which means that the direction of zooplankton cross-diffusion can affect the phytoplankton density. In Figure 2(e), coefficient D
21 is positive and D
12 is negative, and it is not difficult to see that the spatial distribution of phytoplankton differs significantly from the other panels shown in Figure 2. These results demonstrate that cross-diffusion plays a significant role in the stability of the system, but they also show the effects of the various types of cross-diffusion that a plankton system may encounter.


## 3.2. Pattern Formation

In order to illustrate the spatial distribution of phytoplankton more clearly, we show the patterns formed by system (2) in a two-dimensional space in Figure 3.Figure 3
Patterns obtained with system (2) using the following parameters: K
=
0.205,u
1
=
0.7,r
1
=
0.2,u
2
=
0.382,β
1
=
0.21,β
2
=
0.2428,β
3
=
0.1,v
1
=
0.09,v
2
=
0.013,D
11
=
0.04, and D
22
=
3.5. (b) D
21
=
0.4, D
12
=
0.1; (c) D
21
=
0.91, D
12
=
-
0.5. Time steps: (a) 0, (b) 1000, and (c) 1000.
(a)
(b)
(c)Figure3 shows two different spatial distributions of phytoplankton. Figure 3(a) shows the initial spatial distributions of phytoplankton. Figures 3(b) and 3(c) illustrate different types of pattern formation for different values of D
12 and D
21 at t
=
1000. When D
21
=
0.4, D
12
=
0.1, Figure 3(b) shows that a spotted pattern and the spots appear to be cluttered with an irregular shape. However, when D
21
=
0.91, D
12
=
-
0.5, a circular spotted pattern appears, and the distribution of the spots is more regular, as shown in Figure 3(c). As discussed in Section 3.1, cross-diffusion will lead to instability in the system, where both positive cross-diffusion and negative cross-diffusion play important roles in the stability of the system. The spatial distributions of phytoplankton shown in Figures 2(b) and 2(c) agree well with the results given in the previous section.


## 4. Conclusion and Discussion

The system is simple because it is only an abstract description of a marine plankton system. However, the system exhibited some rich features of the marine plankton system. The results show that the spatial distribution of phytoplankton is relatively uniform and stable under some conditions if we only consider the self-diffusion in system (2). When the influence of cross-diffusion is considered, however, the distribution of phytoplankton will change significantly, as shown in Figure 2. Our results demonstrate that the cross-diffusion can greatly affect the dynamic behavior of plankton system.In Section2, we showed that complex equilibria exist in the system and that a bistable system exists when the equilibrium satisfies certain conditions. The results shown in Table 1 verified these conclusions. We also found the condition for Hopf bifurcation. In Table 2, we take parameters values of u
1
=
0.7, r
1
=
0.2,  u
2
=
0.382, β
1
=
0.21, β
2
=
0.2428, β
3
=
0.1, v
1
=
0.09, and v
2
=
0.013, and we illustrated the changes in stability by varying the half-saturation constant for Holling type II functional response in system (1). Although the system is simple, the behavior of the system was found to be very rich.In contrast to previous studies of marine plankton systems, diffusion system (2) considered the mutual effects of self-diffusion and cross-diffusion. Furthermore, various effects of the saturation of cross-diffusion were considered in our study and we analyzed the influence of differences in real saturation. In the paper, theoretical proof of the stability of the system, we considered positive or negative or zero cross-diffusion and we found that the conditions for system instability were caused by cross-diffusion. In our numerical analysis, we focused on the influence of cross-diffusion on the system. Figure 2 shows the influence of cross-diffusion on the system, where we considered the complexity of the actual situation, and the results are illustrated for different values of D
12 and D
21. In Figures 2(a)–2(f), the results of our theoretical derivation are verified, and they also reflect the complex form of the system under the influence of cross-diffusion. Turing first proposed the reaction-diffusion equation and described the concept of Turing instability in 1952 [34], and pattern formation has now become an important component when studying the reaction-diffusion equation. Thus, based on Figure 2, we analyzed the patterns formed to further illustrate the instability of the system as in Figure 3.In addition, it is known that climate change can affect the diffusion of population in marine. Our studies show that the diffusion of phytoplankton and zooplankton can affect the spatial distribution of population under some conditions. Hence, climate change may have a significant effect on the population in marine. We hope that this finding can be further studied and proved to be useful in the study of toxic plankton systems.



---

*Source: 4893451-2016-05-04.xml*