# Stability and Dynamical Analysis of a Biological System

Authors: Xinhong Pan; Min Zhao; Yapei Wang; Hengguo Yu; Zengling Ma; Qi Wang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This study considers the spatiotemporal dynamics of a reaction-diffusion phytoplankton-zooplankton system with a double Allee effect on prey under a homogeneous boundary condition. The qualitative properties are analyzed, including the local stability of all equilibria and the global asymptotic property of the unique positive equilibrium. We also discuss the Hopf bifurcation and the steady state bifurcation of the system. These results are expected to help understand the complexity of the Allee effect and the interaction between phytoplankton and zooplankton.


---

## Body

## 1. Introduction

The upper layer of the ocean contains large volumes of drifting plankton, which can be divided into phytoplankton and zooplankton. Phytoplankton is the autotrophic component of the plankton community, which is consumed by zooplankton, most of which are too small to observe individually with the naked eye. Zooplankton, which are heterotrophic organisms in oceans, are also mostly invisible to the naked eye. Therefore, it is difficult and expensive to quantify plankton directly. Plankton not only play an important role in the marine system because they are at the bottom level of the food chain that supports commercial fisheries, but also play important roles in the cycling of many chemical elements, such as carbon, which may affect climate change [1]. Furthermore, when plankton such as blue algae and dinoflagellates are present in large concentrations, the water appears to be discolored or murky, which is known as a red tide, and this can result in the death of marine and coastal species of fish, mammals, and other organisms [1]. Thus, analyzing the dynamics of plankton using mathematical models is beneficial for understanding the features of plankton populations, which have enormous economic and ecological value.However, the mechanism that leads to the occurrence of red tides is still an unsolved issue. Many models and theories have been proposed by mathematicians and ecologists to explain this phenomenon, but a general and correct explanation still remains a distant goal [2–5].The popular mathematical model called a Gause-type predator-prey model is used to consider the phytoplankton-zooplankton interaction in the following form:(1)
d
p
d
t
=
f
(
p
,
K
)
p
-
g
(
p
)
z
,
d
z
d
t
=
e
g
(
p
)
z
-
μ
z
,
where p
=
p
(
t
) and z
=
z
(
t
) are the population densities or biomass of phytoplankton and zooplankton at time t, respectively. f
(
p
,
K
) describes the intrinsic per capita growth rate of the phytoplankton, which may be a logistic growth function, exponential growth function, or other functions, where K is known as the environmental carrying capacity. g
(
p
) is the per unit-predator consumption rate of prey, which is commonly called the functional response. Some conventional forms of functional response include Holling types I, II, III, and IV and Ivlev type [6–8]. e
(
0
<
e
<
1
) is the conversion coefficient and μ represents the per capita predator mortality, which is assumed to have a linear form, although other forms are possible [9, 10]. The global dynamics of model (1) with a logistic growth rate have been studied during the last three decades based on theoretical analysis and numerical simulations, and many results have been reported [11–15].In recent years, the Allee effect has been the focus of increasing interest and it is recognized to be an important phenomenon in many fields of ecology and conservation biology by more and more people [16–21]. The Allee effect is named after W.C. Allee [22] and it describes a positive correlation between the density or number of population and individual fitness of population [16]. Standard population models assume that the fitness of population increases as the population density or size declines [11–15, 23, 24], whereas Allee effect states that when a population is below a critical density or size, the population cannot sustain itself and this leads to extinction. Thus, the Allee effect increases the likelihood of extinction [25]. Stephens et al. distinguished between a component Allee effect and a demographic Allee effect [17]. However, conservation biologists are usually more interested in the demographic Allee effect because it ultimately governs the probability of the extinction or recovery of populations with low abundances [16].Very recent ecological research has shown that two or more Allee effects can act on a single population simultaneously, which is known as the multiple (double) Allee effect [26, 27].There are many ways of describing the Allee effect [28], including the following differential equation:
(2)
d
p
d
t
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
,
where r describes the growth rate and n is an auxiliary parameter where n
>
0,  m
+
n
>
0. Indeed, it is considered that p represents the size of a fertile population and n is the nonfertile population, such as juvenile or oldest individuals [29]. In this case, m
(
-
K
<
m
<
K
) is called the Allee threshold because when the population density or size is below this threshold, the population is destined for extinction. When m
>
0, (2) describes a strong Allee effect [8, 30, 31]. In this case, the population growth rate decreases if the population size is below the threshold m and the population goes to extinction [29]. In addition, (2) describes a weak Allee effect [29, 31–33] for m
<
0.It is obvious that (2) is equal to
(3)
d
p
d
t
=
r
p
p
+
n
(
1
-
p
K
)
(
p
-
m
)
.
González-Olivares et al. [29] state that (3) describes a double Allee effect, that is, once in the factor m
(
p
)
=
p
-
m and the second time in the term r
(
p
)
=
r
p
/
(
p
+
n
) [34, 35].In a marine environment, the plankton populations tend to move in horizontal and vertical directions due to the strong water current. This movement is usually modeled by a reaction-diffusion equation. In this study, we consider the following reaction-diffusion model with constant diffusion coefficient as well as a strong Allee effect in different spatial locations within a fixed smooth bounded domainΩ
∈
R
n. We assume that the response function of the zooplankton follows the law of mass action [15]:
(4)
∂
p
(
x
,
t
)
∂
t
=
r
(
1
-
p
(
x
,
t
)
K
)
(
1
-
m
+
n
p
(
x
,
t
)
+
n
)
p
(
x
,
t
)
-
a
p
(
x
,
t
)
z
(
x
,
t
)
+
d
1
Δ
p
(
x
,
t
)
,
x
∈
Ω
,
t
>
0
,
∂
z
(
x
,
t
)
∂
t
=
e
a
p
(
x
,
t
)
z
(
x
,
t
)
-
μ
z
(
x
,
t
)
+
d
2
Δ
z
(
x
,
t
)
,
x
∈
Ω
,
t
>
0
,
∂
p
(
x
,
t
)
∂
n
=
∂
z
(
x
,
t
)
∂
n
=
0
,
x
∈
∂
Ω
,
t
>
0
,
p
(
x
,
0
)
=
p
0
(
x
)
≥
0
,
z
(
x
,
0
)
=
z
0
(
x
)
≥
0
,
x
∈
Ω
,
where m
(
0
<
m
<
K
) is the Allee threshold, d
1
,
d
2 are the diffusion coefficients of phytoplankton and zooplankton, respectively, Δ is the Laplacian operator, and x is the spatial habitat of two species, and we assume that the system is a close ecosystem and with a no-flux boundary condition.This paper is structured as follows. In Section2, we analyze the basic dynamics of (4) including estimates of the solution and the local and global stability of equilibria. In Section 3, we provide the analysis of the Hopf bifurcation and the steady state bifurcation. A brief discussion and summary are given in Section 4.


## 2. Main Results

### 2.1. Basic Dynamics

Suppose thatn
,
a
,
μ
,
d
1
,
d
2
>
0, 0
<
e
<
1, 0
<
m
<
K, and Ω is a bounded domain; then we obtain the following results.Theorem 1.
The system (4) has a unique solution and the solution is bounded. Furthermore, the solution (
p
,
z
) of (4) satisfies
(5)
lim
⁡
t
→
∞
sup
⁡
p
(
x
,
t
)
≤
K
,
lim
⁡
t
→
∞
sup
⁡
∫
Ω
z
(
x
,
t
)
d
x
≤
K
[
e
+
e
r
K
n
μ
(
K
-
m
2
)
2
]
|
Ω
|
if d
1
=
d
2, lim
⁡
t
→
∞
sup
⁡
z
(
x
,
t
)
≤
K
[
e
+
(
e
r
/
K
n
μ
)
(
(
K
-
m
)
/
2
)
2
].Proof.
Let(6)
M
(
p
,
z
)
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
,
N
(
p
,
z
)
=
e
a
p
z
-
μ
z
;
then ∂
M
/
∂
z
=
-
a
p
≤
0, ∂
N
/
∂
p
=
e
a
z
≥
0 for {
p
≥
0
,
z
≥
0
}, which implies that (4) is a mixed quasi-monotone system [36]. We define (
p
_
(
x
,
t
)
,
z
_
(
x
,
t
)
)
=
(
0,0
), (
p
¯
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
=
(
p
*
(
t
)
,
z
*
(
t
)
), where (
p
*
(
t
)
,
z
*
(
t
)
) is the solution of the system
(7)
d
p
d
t
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
,
d
z
d
t
=
e
a
p
z
-
μ
z
,
p
(
0
)
=
p
*
,
z
(
0
)
=
z
*
,
where p
*
=
su
p
Ω
¯
p
0
(
x
),  z
*
=
su
p
Ω
¯
z
0
(
x
). We find that
(8)
∂
p
¯
(
x
,
t
)
∂
t
-
Δ
p
¯
(
x
,
t
)
-
M
(
p
¯
(
x
,
t
)
,
z
_
(
x
,
t
)
)
=
0
,
∂
p
_
(
x
,
t
)
∂
t
-
Δ
p
_
(
x
,
t
)
-
M
(
p
_
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
=
0
.
Therefore,
(9)
∂
p
¯
(
x
,
t
)
∂
t
-
Δ
p
¯
(
x
,
t
)
-
M
(
p
¯
(
x
,
t
)
,
z
_
(
x
,
t
)
)
≥
∂
p
_
(
x
,
t
)
∂
t
-
Δ
p
_
(
x
,
t
)
-
M
(
p
_
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
.
Similarly,
(10)
∂
z
¯
(
x
,
t
)
∂
t
-
Δ
z
¯
(
x
,
t
)
-
N
(
p
¯
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
≥
∂
z
_
(
x
,
t
)
∂
t
-
Δ
z
_
(
x
,
t
)
-
N
(
p
_
(
x
,
t
)
,
z
_
(
x
,
t
)
)
.
This implies that (
p
_
(
x
,
t
)
,
z
_
(
x
,
t
)
)
=
(
0,0
) and (
p
¯
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
=
(
p
*
(
t
)
,
z
*
(
t
)
) are the lower solution and upper solution to (4), respectively (Figure 2). Therefore, Theorem 5.3
.
3 in [36] shows that system (4) has a unique solution (
(
p
(
x
,
t
)
,
z
(
x
,
t
)
) that satisfies
(11)
0
≤
p
(
x
,
t
)
≤
p
*
(
t
)
,
0
≤
z
(
x
,
t
)
≤
z
*
(
t
)
.
According to the strong maximum principle and the boundary condition, p
(
x
,
t
)
>
0, z
(
x
,
t
)
>
0 for t
>
0 and x
∈
Ω
¯.
From the first equation of (7), we can see that p
*
(
t
)
→
0 for p
*
<
m and p
*
(
t
)
→
K for p
*
>
m. Thus, lim
⁡
t
→
∞
⁡
sup
⁡
p
(
x
,
t
)
≤
K. To estimate z
(
x
,
t
), let ∫
Ω
p
(
x
,
t
)
d
x
=
P
(
t
), ∫
Ω
z
(
x
,
t
)
d
x
=
Z
(
t
); then
(12)
d
P
d
t
=
∫
Ω
p
t
d
x
=
∫
Ω
[
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
]
d
x
+
∫
Ω
d
1
Δ
p
d
x
,
(13)
d
Z
d
t
=
∫
Ω
e
a
p
z
d
x
-
μ
Z
+
∫
Ω
d
2
Δ
z
d
x
,
(12)  *
e
+ (13) and with boundary conditions, we obtain
(14)
(
e
P
+
Z
)
t
=
e
∫
Ω
r
(
1
-
p
K
)
p
-
m
p
+
n
p
d
x
-
μ
Z
=
e
∫
Ω
r
(
1
-
p
K
)
p
-
m
p
+
n
p
d
x
-
μ
(
e
P
+
Z
)
-
μ
e
P
≤
-
μ
(
e
P
+
Z
)
+
[
μ
e
+
e
r
K
n
(
K
-
m
2
)
2
]
P
.
From the proof above, we know that lim
⁡
t
→
∞
sup
⁡
P
(
t
)
≤
K
·
|
Ω
|. Thus, for any ε
>
0, ∃
T
1
>
0, when t
>
T
1
(15)
(
e
P
+
Z
)
t
≤
-
μ
(
e
P
+
Z
)
+
[
μ
e
+
e
r
K
n
(
K
-
m
2
)
2
]
(
K
+
ε
)
|
Ω
|
.
By integrating (15), ∃
T
2
>
T
1, we have
(16)
∫
Ω
z
(
x
,
t
)
d
x
=
Z
(
t
)
<
e
P
(
t
)
+
Z
(
t
)
≤
K
+
ε
μ
[
μ
e
+
e
r
K
n
(
K
-
m
2
)
2
]
|
Ω
|
+
ε
for
t
>
T
2
.
This implies that lim
⁡
t
→
∞
sup
⁡
∫
Ω
z
(
x
,
t
)
d
x
≤
K
[
e
+
(
e
r
/
K
n
μ
)
(
(
K
-
m
)
/
2
)
2
]
|
Ω
|.
Ifd
1
=
d
2, we can add the two equations in (4) and we have
(17)
∂
w
∂
t
=
e
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
μ
z
+
d
1
Δ
w
,
x
∈
Ω
,
t
>
T
,
∂
w
∂
n
=
0
,
x
∈
∂
Ω
,
t
>
T
,
w
(
x
,
T
)
=
e
p
(
x
,
T
)
+
z
(
x
,
T
)
,
x
∈
Ω
.
Since
(18)
e
r
(
1
-
p
K
)
p
-
m
p
+
n
p
-
μ
z
=
e
r
(
1
-
p
K
)
p
-
m
p
+
n
p
-
μ
w
+
μ
e
p
≤
[
e
r
K
n
(
K
-
m
2
)
2
+
μ
e
]
p
-
μ
w
≤
[
e
r
K
n
(
K
-
m
2
)
2
+
μ
e
]
(
K
+
ε
)
-
μ
w
,
we know that the solution of the equation
(19)
∂
Φ
∂
t
=
[
e
r
K
n
(
K
-
m
2
)
2
+
μ
e
]
(
K
+
ε
)
-
μ
v
+
d
1
Δ
v
,
x
∈
Ω
,
t
>
T
,
∂
Φ
∂
n
=
0
,
x
∈
∂
Ω
,
t
>
T
is Φ
(
x
,
t
)
→
(
1
/
μ
)
[
(
e
r
/
K
n
)
(
(
K
-
m
)
/
2
)
2
+
μ
e
]
(
K
+
ε
) for t
→
∞. The comparison argument implies that lim
⁡
t
→
∞
sup
⁡
z
(
x
,
t
)
≤
lim
⁡
t
→
∞
sup
⁡
w
(
x
,
t
)
≤
(
1
/
μ
)
[
(
e
r
/
K
n
)
(
(
K
-
m
)
/
2
)
2
+
μ
e
]
(
K
+
ε
).Time evolution of (4) around (
0,0
). Phytoplankton is in (a) and zooplankton is in (b). The parameter values are d
1
=
0.001, d
2
=
0.004, e
=
0.5, a
=
1, n
=
0.5, m
=
0.5, μ
=
0.4, r
=
2, and K
=
10.
(a)
(b)Solution to(
p
(
x
,
t
)
,
z
(
x
,
t
)
) for (4). Phytoplankton is in (a) and zooplankton is in (b). The parameter values are the same as those given in Figure 1.
(a)
(b)Theorem 2.
Ifp
0
(
x
)
≤
m, then (
p
(
x
,
t
)
,
z
(
x
,
t
)
)
→
(
0,0
) or (
p
(
x
,
t
)
,
z
(
x
,
t
)
)
→
(
m
,
0
) as t
→
∞.Proof.
From the proof of Theorem1, if p
0
(
x
)
≤
p
*
<
m, then p
*
(
t
)
→
0 and consequently z
*
(
t
)
→
0 as t
→
∞. This completes the proof.From a biological viewpoint, this implies that if the initial population density is below the thresholdm, the phytoplankton become extinct so the zooplankton would become extinct.


### 2.2. Local and Global Stability of Equilibria

System (4) has four nonnegative steady state solutions: (
0,0
), (
m
,
0
), (
K
,
0
), and (
u
,
v
u
), where u
=
μ
/
e
a,  v
u
=
(
r
/
a
)
(
1
-
(
u
/
K
)
)
(
(
u
-
m
)
/
(
u
+
n
)
).The local stability of the steady state solutions can be analyzed as follows.Theorem 3.
(1)(
0,0
) is locally asymptotically stable.
(2)(
K
,
0
) is unstable.
(3)(
m
,
0
) is locally asymptotically stable when e
a
K
-
μ
<
0 and is unstable for e
a
K
-
μ
>
0.
(4)(
u
,
v
u
) is locally asymptotically stable for u
¯
<
u
<
K and is unstable for u
¯
<
u
<
K, where
(20)
u
¯
=
2
n
-
4
n
2
+
4
(
m
n
+
m
K
+
n
K
)
-
2
.Proof.
The linearization of (4) at solution (
p
,
z
) can be expressed as
(21)
∂
U
∂
t
=
(
D
Δ
+
J
(
p
,
z
)
)
U
,
where
(22)
U
=
(
p
(
x
,
t
)
,
z
(
x
,
t
)
)
T
,
D
=
diag
⁡
(
d
1
,
d
2
)
,
J
(
p
,
z
)
=
(
A
(
p
,
z
)
B
(
p
,
z
)
C
(
p
,
z
)
D
(
p
,
z
)
)
,
A
(
p
,
z
)
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
-
r
K
p
(
1
-
m
+
n
p
+
n
)
+
r
(
1
-
p
K
)
p
m
+
n
(
p
+
n
)
2
-
a
z
,
B
(
p
,
z
)
=
-
a
p
,
C
(
p
,
z
)
=
e
a
z
,
D
(
p
,
z
)
=
e
a
p
-
μ
.
According to Theorems 5.1
.
1 and 5.1
.
3 from [37], we know that if all the eigenvalues of the operator L have negative real parts, then the solution (
p
,
z
) is asymptotically stable; if there is at least one eigenvalue with a positive real part, then the solution (
p
,
z
) is unstable; if some eigenvalues have zero real parts, then the stability cannot be determined using this method.
Letλ
i
(
i
=
0,1
,
2
,
…
) be the eigenvalues of -
Δ on Ω under a homogeneous Neumann boundary condition and 0
=
λ
0
<
λ
1
≤
λ
2
≤
⋯, lim
⁡
i
→
∞
λ
i
=
∞. Thus, it is known that ξ is an eigenvalue of L if and only if ξ is the eigenvalue of the matrix J
i
=
-
λ
i
D
+
J
(
p
,
z
) for some i
≥
0.
Consider the characteristic equation(23)
det
⁡
(
ξ
I
-
J
i
)
=
ξ
2
-
trace
J
i
ξ
+
det
⁡
J
i
,
where
(24)
trace
J
i
=
-
λ
i
(
d
1
+
d
2
)
+
A
(
p
,
z
)
+
D
(
p
,
z
)
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
-
(
A
(
p
,
z
)
d
2
+
D
(
p
,
z
)
d
1
)
λ
i
+
det
⁡
J
(
p
,
z
)
.
(1) If(
p
,
z
)
=
(
0,0
), then J
(
0,0
)
=
(
r
·
(
-
(
m
/
n
)
)
0
0
-
μ
)
(25)
trace
J
i
=
-
λ
i
(
d
1
+
d
2
)
-
(
r
m
n
+
μ
)
<
0
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
+
(
r
m
n
d
2
+
μ
d
1
)
λ
i
+
r
m
μ
n
>
0
.
This implies that (
0,0
) is locally asymptotically stable.
(2) If(
p
,
z
)
=
(
m
,
0
), then J
(
m
,
0
)
=
(
r
·
(
1
-
(
m
/
K
)
)
(
m
/
(
m
+
n
)
)
-
a
m
0
e
a
m
-
μ
).
Fori
=
0, one of the eigenvalues is r
(
1
-
(
m
/
K
)
)
(
m
/
(
m
+
n
)
)
>
0, which implies that (
m
,
0
) is an unstable point.
(3) If(
p
,
z
)
=
(
K
,
0
), then J
(
K
,
0
)
=
(
-
r
(
(
K
-
m
)
/
(
K
+
n
)
)
-
a
K
0
e
a
K
-
μ
).
Whene
a
K
-
μ
<
0,
(26)
trace
J
i
=
-
λ
i
(
d
1
+
d
2
)
-
r
K
-
m
K
+
n
+
e
a
K
-
μ
<
0
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
+
(
-
r
K
-
m
K
+
n
d
2
+
(
e
a
K
-
μ
)
d
1
)
λ
i
-
r
K
-
m
K
+
n
(
e
a
K
-
μ
)
>
0
.
This implies that (
K
,
0
) is locally asymptotically stable.
Whene
a
K
-
μ
>
0, for i
=
0, det
⁡
J
i
=
-
r
(
(
K
-
m
)
/
(
K
+
n
)
)
(
e
a
K
-
μ
)
<
0, which implies that J
i has at least one eigenvalue with positive real part. This implies that (
K
,
0
) is unstable.
(4) If(
p
,
z
)
=
(
u
,
v
u
), then J
(
u
,
v
u
)
=
(
A
(
u
)
B
(
u
)
C
(
u
)
0
), where
(27)
A
(
u
)
=
r
(
1
-
u
K
)
u
m
+
n
(
u
+
n
)
2
-
r
K
u
u
-
m
u
+
n
,
B
(
u
)
=
-
μ
e
,
C
(
u
)
=
e
r
(
1
-
u
K
)
u
-
m
u
+
n
and v
>
0 implies that m
<
u
<
K. Let u
¯ be the largest root of A
(
u
)
=
0.
Whenu
¯
<
u
<
K, then A
(
u
)
<
0 and
(28)
trace
(
J
i
)
=
-
λ
i
(
d
1
+
d
2
)
+
A
(
u
)
<
0
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
-
A
(
u
)
d
2
λ
i
-
B
(
u
)
C
(
u
)
>
0
.
This implies that (
u
,
v
u
) is locally asymptotically stable.
Whenm
<
u
<
u
¯, then A
(
u
)
>
0 and trace
(
J
i
)
=
A
(
u
)
>
0 for i
=
0, which implies that J
i has at least one eigenvalue with a positive real part. This implies that (
u
,
v
u
) is unstable.Theorem 4.
IfK
m
+
m
n
+
K
n
-
n
(
μ
/
e
a
)
<
0, then the positive constant steady state (
u
,
v
u
) is globally asymptotically stable.Proof.
Let us consider a Lyapunov functionV as
(29)
V
=
∫
Ω
W
(
p
(
x
,
t
)
,
z
(
x
,
t
)
)
d
x
,
where
(30)
W
(
p
(
x
,
t
)
,
z
(
x
,
t
)
)
=
c
1
∫
u
p
p
-
u
p
d
p
+
c
2
∫
v
u
z
z
-
v
u
z
d
z
and c
i
>
0
(
i
=
1,2
) will be determined next.
By taking the time derivative ofV, we have
(31)
d
V
d
t
=
∫
Ω
(
W
p
·
p
t
+
W
z
·
z
t
)
d
x
=
∫
Ω
(
c
1
p
-
u
p
∂
p
∂
t
+
c
2
z
-
v
u
z
∂
z
∂
t
)
d
x
=
∫
Ω
(
c
1
(
p
-
u
)
[
r
(
1
-
x
K
)
p
-
m
p
+
n
-
a
z
]
+
c
2
(
z
-
v
u
)
e
a
(
p
-
u
)
[
r
(
1
-
x
K
)
p
-
m
p
+
n
-
a
z
]
)
d
x
+
∫
Ω
(
c
1
p
-
u
p
d
1
Δ
p
+
c
2
z
-
v
u
z
d
2
Δ
z
)
d
x
=
I
1
(
t
)
+
I
2
(
t
)
.
Due to the Neumann boundary condition, it can easily be derived that
(32)
I
2
(
t
)
=
-
∫
Ω
(
c
1
d
1
u
p
2
|
∇
p
|
2
+
c
2
d
2
v
u
z
2
|
∇
z
|
2
)
d
x
≤
0
.
Further,
(33)
I
1
(
t
)
=
∫
Ω
(
c
1
(
p
-
u
)
[
r
(
1
-
p
K
)
p
-
m
p
+
n
-
r
(
1
-
u
K
)
u
-
m
u
+
n
+
a
v
u
-
a
z
(
1
-
p
K
)
p
-
m
p
+
n
]
+
c
2
e
a
(
p
-
u
)
(
z
-
v
u
)
(
1
-
p
K
)
p
-
m
p
+
n
)
d
x
=
∫
Ω
(
c
1
(
p
-
u
)
2
r
K
[
-
p
u
+
K
m
-
n
(
p
+
u
)
+
n
(
m
+
K
)
(
p
+
n
)
(
u
+
n
)
]
+
(
e
a
c
2
-
a
c
1
)
(
p
-
u
)
(
z
-
v
u
)
-
p
u
+
K
m
-
n
(
p
+
u
)
+
n
(
m
+
K
)
(
p
+
n
)
(
u
+
n
)
)
d
x
.
If we choose c
2
>
0 arbitrarily and c
1
=
e
c
2, then we can obtain
(34)
d
V
d
t
=
I
1
(
t
)
+
I
2
(
t
)
≤
I
1
(
t
)
=
∫
Ω
c
1
r
K
[
-
p
u
+
K
m
+
n
(
m
+
K
)
-
n
(
p
+
u
)
(
p
+
n
)
(
u
+
n
)
]
×
(
p
-
u
)
2
d
x
≤
∫
Ω
c
1
r
K
K
m
+
m
n
+
K
n
-
n
u
n
(
u
+
n
)
(
p
-
u
)
2
d
x
.
Therefore, if K
m
+
m
n
+
K
n
-
n
u
<
0, then d
V
/
d
t
≤
0 and d
V
/
d
t
=
0 if p
=
u
,
z
=
v
u. This completes the proof.


## 2.1. Basic Dynamics

Suppose thatn
,
a
,
μ
,
d
1
,
d
2
>
0, 0
<
e
<
1, 0
<
m
<
K, and Ω is a bounded domain; then we obtain the following results.Theorem 1.
The system (4) has a unique solution and the solution is bounded. Furthermore, the solution (
p
,
z
) of (4) satisfies
(5)
lim
⁡
t
→
∞
sup
⁡
p
(
x
,
t
)
≤
K
,
lim
⁡
t
→
∞
sup
⁡
∫
Ω
z
(
x
,
t
)
d
x
≤
K
[
e
+
e
r
K
n
μ
(
K
-
m
2
)
2
]
|
Ω
|
if d
1
=
d
2, lim
⁡
t
→
∞
sup
⁡
z
(
x
,
t
)
≤
K
[
e
+
(
e
r
/
K
n
μ
)
(
(
K
-
m
)
/
2
)
2
].Proof.
Let(6)
M
(
p
,
z
)
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
,
N
(
p
,
z
)
=
e
a
p
z
-
μ
z
;
then ∂
M
/
∂
z
=
-
a
p
≤
0, ∂
N
/
∂
p
=
e
a
z
≥
0 for {
p
≥
0
,
z
≥
0
}, which implies that (4) is a mixed quasi-monotone system [36]. We define (
p
_
(
x
,
t
)
,
z
_
(
x
,
t
)
)
=
(
0,0
), (
p
¯
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
=
(
p
*
(
t
)
,
z
*
(
t
)
), where (
p
*
(
t
)
,
z
*
(
t
)
) is the solution of the system
(7)
d
p
d
t
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
,
d
z
d
t
=
e
a
p
z
-
μ
z
,
p
(
0
)
=
p
*
,
z
(
0
)
=
z
*
,
where p
*
=
su
p
Ω
¯
p
0
(
x
),  z
*
=
su
p
Ω
¯
z
0
(
x
). We find that
(8)
∂
p
¯
(
x
,
t
)
∂
t
-
Δ
p
¯
(
x
,
t
)
-
M
(
p
¯
(
x
,
t
)
,
z
_
(
x
,
t
)
)
=
0
,
∂
p
_
(
x
,
t
)
∂
t
-
Δ
p
_
(
x
,
t
)
-
M
(
p
_
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
=
0
.
Therefore,
(9)
∂
p
¯
(
x
,
t
)
∂
t
-
Δ
p
¯
(
x
,
t
)
-
M
(
p
¯
(
x
,
t
)
,
z
_
(
x
,
t
)
)
≥
∂
p
_
(
x
,
t
)
∂
t
-
Δ
p
_
(
x
,
t
)
-
M
(
p
_
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
.
Similarly,
(10)
∂
z
¯
(
x
,
t
)
∂
t
-
Δ
z
¯
(
x
,
t
)
-
N
(
p
¯
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
≥
∂
z
_
(
x
,
t
)
∂
t
-
Δ
z
_
(
x
,
t
)
-
N
(
p
_
(
x
,
t
)
,
z
_
(
x
,
t
)
)
.
This implies that (
p
_
(
x
,
t
)
,
z
_
(
x
,
t
)
)
=
(
0,0
) and (
p
¯
(
x
,
t
)
,
z
¯
(
x
,
t
)
)
=
(
p
*
(
t
)
,
z
*
(
t
)
) are the lower solution and upper solution to (4), respectively (Figure 2). Therefore, Theorem 5.3
.
3 in [36] shows that system (4) has a unique solution (
(
p
(
x
,
t
)
,
z
(
x
,
t
)
) that satisfies
(11)
0
≤
p
(
x
,
t
)
≤
p
*
(
t
)
,
0
≤
z
(
x
,
t
)
≤
z
*
(
t
)
.
According to the strong maximum principle and the boundary condition, p
(
x
,
t
)
>
0, z
(
x
,
t
)
>
0 for t
>
0 and x
∈
Ω
¯.
From the first equation of (7), we can see that p
*
(
t
)
→
0 for p
*
<
m and p
*
(
t
)
→
K for p
*
>
m. Thus, lim
⁡
t
→
∞
⁡
sup
⁡
p
(
x
,
t
)
≤
K. To estimate z
(
x
,
t
), let ∫
Ω
p
(
x
,
t
)
d
x
=
P
(
t
), ∫
Ω
z
(
x
,
t
)
d
x
=
Z
(
t
); then
(12)
d
P
d
t
=
∫
Ω
p
t
d
x
=
∫
Ω
[
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
]
d
x
+
∫
Ω
d
1
Δ
p
d
x
,
(13)
d
Z
d
t
=
∫
Ω
e
a
p
z
d
x
-
μ
Z
+
∫
Ω
d
2
Δ
z
d
x
,
(12)  *
e
+ (13) and with boundary conditions, we obtain
(14)
(
e
P
+
Z
)
t
=
e
∫
Ω
r
(
1
-
p
K
)
p
-
m
p
+
n
p
d
x
-
μ
Z
=
e
∫
Ω
r
(
1
-
p
K
)
p
-
m
p
+
n
p
d
x
-
μ
(
e
P
+
Z
)
-
μ
e
P
≤
-
μ
(
e
P
+
Z
)
+
[
μ
e
+
e
r
K
n
(
K
-
m
2
)
2
]
P
.
From the proof above, we know that lim
⁡
t
→
∞
sup
⁡
P
(
t
)
≤
K
·
|
Ω
|. Thus, for any ε
>
0, ∃
T
1
>
0, when t
>
T
1
(15)
(
e
P
+
Z
)
t
≤
-
μ
(
e
P
+
Z
)
+
[
μ
e
+
e
r
K
n
(
K
-
m
2
)
2
]
(
K
+
ε
)
|
Ω
|
.
By integrating (15), ∃
T
2
>
T
1, we have
(16)
∫
Ω
z
(
x
,
t
)
d
x
=
Z
(
t
)
<
e
P
(
t
)
+
Z
(
t
)
≤
K
+
ε
μ
[
μ
e
+
e
r
K
n
(
K
-
m
2
)
2
]
|
Ω
|
+
ε
for
t
>
T
2
.
This implies that lim
⁡
t
→
∞
sup
⁡
∫
Ω
z
(
x
,
t
)
d
x
≤
K
[
e
+
(
e
r
/
K
n
μ
)
(
(
K
-
m
)
/
2
)
2
]
|
Ω
|.
Ifd
1
=
d
2, we can add the two equations in (4) and we have
(17)
∂
w
∂
t
=
e
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
μ
z
+
d
1
Δ
w
,
x
∈
Ω
,
t
>
T
,
∂
w
∂
n
=
0
,
x
∈
∂
Ω
,
t
>
T
,
w
(
x
,
T
)
=
e
p
(
x
,
T
)
+
z
(
x
,
T
)
,
x
∈
Ω
.
Since
(18)
e
r
(
1
-
p
K
)
p
-
m
p
+
n
p
-
μ
z
=
e
r
(
1
-
p
K
)
p
-
m
p
+
n
p
-
μ
w
+
μ
e
p
≤
[
e
r
K
n
(
K
-
m
2
)
2
+
μ
e
]
p
-
μ
w
≤
[
e
r
K
n
(
K
-
m
2
)
2
+
μ
e
]
(
K
+
ε
)
-
μ
w
,
we know that the solution of the equation
(19)
∂
Φ
∂
t
=
[
e
r
K
n
(
K
-
m
2
)
2
+
μ
e
]
(
K
+
ε
)
-
μ
v
+
d
1
Δ
v
,
x
∈
Ω
,
t
>
T
,
∂
Φ
∂
n
=
0
,
x
∈
∂
Ω
,
t
>
T
is Φ
(
x
,
t
)
→
(
1
/
μ
)
[
(
e
r
/
K
n
)
(
(
K
-
m
)
/
2
)
2
+
μ
e
]
(
K
+
ε
) for t
→
∞. The comparison argument implies that lim
⁡
t
→
∞
sup
⁡
z
(
x
,
t
)
≤
lim
⁡
t
→
∞
sup
⁡
w
(
x
,
t
)
≤
(
1
/
μ
)
[
(
e
r
/
K
n
)
(
(
K
-
m
)
/
2
)
2
+
μ
e
]
(
K
+
ε
).Time evolution of (4) around (
0,0
). Phytoplankton is in (a) and zooplankton is in (b). The parameter values are d
1
=
0.001, d
2
=
0.004, e
=
0.5, a
=
1, n
=
0.5, m
=
0.5, μ
=
0.4, r
=
2, and K
=
10.
(a)
(b)Solution to(
p
(
x
,
t
)
,
z
(
x
,
t
)
) for (4). Phytoplankton is in (a) and zooplankton is in (b). The parameter values are the same as those given in Figure 1.
(a)
(b)Theorem 2.
Ifp
0
(
x
)
≤
m, then (
p
(
x
,
t
)
,
z
(
x
,
t
)
)
→
(
0,0
) or (
p
(
x
,
t
)
,
z
(
x
,
t
)
)
→
(
m
,
0
) as t
→
∞.Proof.
From the proof of Theorem1, if p
0
(
x
)
≤
p
*
<
m, then p
*
(
t
)
→
0 and consequently z
*
(
t
)
→
0 as t
→
∞. This completes the proof.From a biological viewpoint, this implies that if the initial population density is below the thresholdm, the phytoplankton become extinct so the zooplankton would become extinct.


## 2.2. Local and Global Stability of Equilibria

System (4) has four nonnegative steady state solutions: (
0,0
), (
m
,
0
), (
K
,
0
), and (
u
,
v
u
), where u
=
μ
/
e
a,  v
u
=
(
r
/
a
)
(
1
-
(
u
/
K
)
)
(
(
u
-
m
)
/
(
u
+
n
)
).The local stability of the steady state solutions can be analyzed as follows.Theorem 3.
(1)(
0,0
) is locally asymptotically stable.
(2)(
K
,
0
) is unstable.
(3)(
m
,
0
) is locally asymptotically stable when e
a
K
-
μ
<
0 and is unstable for e
a
K
-
μ
>
0.
(4)(
u
,
v
u
) is locally asymptotically stable for u
¯
<
u
<
K and is unstable for u
¯
<
u
<
K, where
(20)
u
¯
=
2
n
-
4
n
2
+
4
(
m
n
+
m
K
+
n
K
)
-
2
.Proof.
The linearization of (4) at solution (
p
,
z
) can be expressed as
(21)
∂
U
∂
t
=
(
D
Δ
+
J
(
p
,
z
)
)
U
,
where
(22)
U
=
(
p
(
x
,
t
)
,
z
(
x
,
t
)
)
T
,
D
=
diag
⁡
(
d
1
,
d
2
)
,
J
(
p
,
z
)
=
(
A
(
p
,
z
)
B
(
p
,
z
)
C
(
p
,
z
)
D
(
p
,
z
)
)
,
A
(
p
,
z
)
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
-
r
K
p
(
1
-
m
+
n
p
+
n
)
+
r
(
1
-
p
K
)
p
m
+
n
(
p
+
n
)
2
-
a
z
,
B
(
p
,
z
)
=
-
a
p
,
C
(
p
,
z
)
=
e
a
z
,
D
(
p
,
z
)
=
e
a
p
-
μ
.
According to Theorems 5.1
.
1 and 5.1
.
3 from [37], we know that if all the eigenvalues of the operator L have negative real parts, then the solution (
p
,
z
) is asymptotically stable; if there is at least one eigenvalue with a positive real part, then the solution (
p
,
z
) is unstable; if some eigenvalues have zero real parts, then the stability cannot be determined using this method.
Letλ
i
(
i
=
0,1
,
2
,
…
) be the eigenvalues of -
Δ on Ω under a homogeneous Neumann boundary condition and 0
=
λ
0
<
λ
1
≤
λ
2
≤
⋯, lim
⁡
i
→
∞
λ
i
=
∞. Thus, it is known that ξ is an eigenvalue of L if and only if ξ is the eigenvalue of the matrix J
i
=
-
λ
i
D
+
J
(
p
,
z
) for some i
≥
0.
Consider the characteristic equation(23)
det
⁡
(
ξ
I
-
J
i
)
=
ξ
2
-
trace
J
i
ξ
+
det
⁡
J
i
,
where
(24)
trace
J
i
=
-
λ
i
(
d
1
+
d
2
)
+
A
(
p
,
z
)
+
D
(
p
,
z
)
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
-
(
A
(
p
,
z
)
d
2
+
D
(
p
,
z
)
d
1
)
λ
i
+
det
⁡
J
(
p
,
z
)
.
(1) If(
p
,
z
)
=
(
0,0
), then J
(
0,0
)
=
(
r
·
(
-
(
m
/
n
)
)
0
0
-
μ
)
(25)
trace
J
i
=
-
λ
i
(
d
1
+
d
2
)
-
(
r
m
n
+
μ
)
<
0
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
+
(
r
m
n
d
2
+
μ
d
1
)
λ
i
+
r
m
μ
n
>
0
.
This implies that (
0,0
) is locally asymptotically stable.
(2) If(
p
,
z
)
=
(
m
,
0
), then J
(
m
,
0
)
=
(
r
·
(
1
-
(
m
/
K
)
)
(
m
/
(
m
+
n
)
)
-
a
m
0
e
a
m
-
μ
).
Fori
=
0, one of the eigenvalues is r
(
1
-
(
m
/
K
)
)
(
m
/
(
m
+
n
)
)
>
0, which implies that (
m
,
0
) is an unstable point.
(3) If(
p
,
z
)
=
(
K
,
0
), then J
(
K
,
0
)
=
(
-
r
(
(
K
-
m
)
/
(
K
+
n
)
)
-
a
K
0
e
a
K
-
μ
).
Whene
a
K
-
μ
<
0,
(26)
trace
J
i
=
-
λ
i
(
d
1
+
d
2
)
-
r
K
-
m
K
+
n
+
e
a
K
-
μ
<
0
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
+
(
-
r
K
-
m
K
+
n
d
2
+
(
e
a
K
-
μ
)
d
1
)
λ
i
-
r
K
-
m
K
+
n
(
e
a
K
-
μ
)
>
0
.
This implies that (
K
,
0
) is locally asymptotically stable.
Whene
a
K
-
μ
>
0, for i
=
0, det
⁡
J
i
=
-
r
(
(
K
-
m
)
/
(
K
+
n
)
)
(
e
a
K
-
μ
)
<
0, which implies that J
i has at least one eigenvalue with positive real part. This implies that (
K
,
0
) is unstable.
(4) If(
p
,
z
)
=
(
u
,
v
u
), then J
(
u
,
v
u
)
=
(
A
(
u
)
B
(
u
)
C
(
u
)
0
), where
(27)
A
(
u
)
=
r
(
1
-
u
K
)
u
m
+
n
(
u
+
n
)
2
-
r
K
u
u
-
m
u
+
n
,
B
(
u
)
=
-
μ
e
,
C
(
u
)
=
e
r
(
1
-
u
K
)
u
-
m
u
+
n
and v
>
0 implies that m
<
u
<
K. Let u
¯ be the largest root of A
(
u
)
=
0.
Whenu
¯
<
u
<
K, then A
(
u
)
<
0 and
(28)
trace
(
J
i
)
=
-
λ
i
(
d
1
+
d
2
)
+
A
(
u
)
<
0
,
det
⁡
J
i
=
d
1
d
2
λ
i
2
-
A
(
u
)
d
2
λ
i
-
B
(
u
)
C
(
u
)
>
0
.
This implies that (
u
,
v
u
) is locally asymptotically stable.
Whenm
<
u
<
u
¯, then A
(
u
)
>
0 and trace
(
J
i
)
=
A
(
u
)
>
0 for i
=
0, which implies that J
i has at least one eigenvalue with a positive real part. This implies that (
u
,
v
u
) is unstable.Theorem 4.
IfK
m
+
m
n
+
K
n
-
n
(
μ
/
e
a
)
<
0, then the positive constant steady state (
u
,
v
u
) is globally asymptotically stable.Proof.
Let us consider a Lyapunov functionV as
(29)
V
=
∫
Ω
W
(
p
(
x
,
t
)
,
z
(
x
,
t
)
)
d
x
,
where
(30)
W
(
p
(
x
,
t
)
,
z
(
x
,
t
)
)
=
c
1
∫
u
p
p
-
u
p
d
p
+
c
2
∫
v
u
z
z
-
v
u
z
d
z
and c
i
>
0
(
i
=
1,2
) will be determined next.
By taking the time derivative ofV, we have
(31)
d
V
d
t
=
∫
Ω
(
W
p
·
p
t
+
W
z
·
z
t
)
d
x
=
∫
Ω
(
c
1
p
-
u
p
∂
p
∂
t
+
c
2
z
-
v
u
z
∂
z
∂
t
)
d
x
=
∫
Ω
(
c
1
(
p
-
u
)
[
r
(
1
-
x
K
)
p
-
m
p
+
n
-
a
z
]
+
c
2
(
z
-
v
u
)
e
a
(
p
-
u
)
[
r
(
1
-
x
K
)
p
-
m
p
+
n
-
a
z
]
)
d
x
+
∫
Ω
(
c
1
p
-
u
p
d
1
Δ
p
+
c
2
z
-
v
u
z
d
2
Δ
z
)
d
x
=
I
1
(
t
)
+
I
2
(
t
)
.
Due to the Neumann boundary condition, it can easily be derived that
(32)
I
2
(
t
)
=
-
∫
Ω
(
c
1
d
1
u
p
2
|
∇
p
|
2
+
c
2
d
2
v
u
z
2
|
∇
z
|
2
)
d
x
≤
0
.
Further,
(33)
I
1
(
t
)
=
∫
Ω
(
c
1
(
p
-
u
)
[
r
(
1
-
p
K
)
p
-
m
p
+
n
-
r
(
1
-
u
K
)
u
-
m
u
+
n
+
a
v
u
-
a
z
(
1
-
p
K
)
p
-
m
p
+
n
]
+
c
2
e
a
(
p
-
u
)
(
z
-
v
u
)
(
1
-
p
K
)
p
-
m
p
+
n
)
d
x
=
∫
Ω
(
c
1
(
p
-
u
)
2
r
K
[
-
p
u
+
K
m
-
n
(
p
+
u
)
+
n
(
m
+
K
)
(
p
+
n
)
(
u
+
n
)
]
+
(
e
a
c
2
-
a
c
1
)
(
p
-
u
)
(
z
-
v
u
)
-
p
u
+
K
m
-
n
(
p
+
u
)
+
n
(
m
+
K
)
(
p
+
n
)
(
u
+
n
)
)
d
x
.
If we choose c
2
>
0 arbitrarily and c
1
=
e
c
2, then we can obtain
(34)
d
V
d
t
=
I
1
(
t
)
+
I
2
(
t
)
≤
I
1
(
t
)
=
∫
Ω
c
1
r
K
[
-
p
u
+
K
m
+
n
(
m
+
K
)
-
n
(
p
+
u
)
(
p
+
n
)
(
u
+
n
)
]
×
(
p
-
u
)
2
d
x
≤
∫
Ω
c
1
r
K
K
m
+
m
n
+
K
n
-
n
u
n
(
u
+
n
)
(
p
-
u
)
2
d
x
.
Therefore, if K
m
+
m
n
+
K
n
-
n
u
<
0, then d
V
/
d
t
≤
0 and d
V
/
d
t
=
0 if p
=
u
,
z
=
v
u. This completes the proof.


## 3. Bifurcation Analysis

In this section, we mainly analyze the stability of the steady state(
u
,
v
u
) and take u as the bifurcation parameter (or equivalently take a as a parameter). In particular, we assume that all of the eigenvalues of -
Δ are simple.We know from the proof of Theorem2 that the stability of (
u
,
v
u
) is determined by the trace and determinant of J
i. Let
(35)
T
(
u
,
λ
)
=
A
(
u
)
-
λ
(
d
1
+
d
2
)
,
D
(
u
,
λ
)
=
d
1
d
2
λ
2
-
d
2
A
(
u
)
λ
-
B
(
u
)
C
(
u
)
.
We refer to {
(
u
,
λ
)
∈
R
+
2
:
T
(
u
,
λ
)
=
0
} as the Hopf bifurcation curve and {
(
u
,
λ
)
∈
R
+
2
:
D
(
u
,
λ
)
=
0
} as the steady state bifurcation curve [38] (Figures 4 and 5).First,T
(
u
,
λ
)
=
0 is equal to λ
=
(
A
(
u
)
/
(
d
1
+
d
2
)
). We can summarize the properties of A
(
u
) as follows, which are easy to prove so the proof is omitted.Lemma 5.
ConsiderA
(
u
)
=
(
r
u
/
K
(
u
+
n
)
2
)
[
-
u
2
-
2
n
u
+
m
n
+
m
K
+
n
K
]; then 0
<
u
*
<
u
¯
<
K exists such that the following hold.(1)
Ifm
n
+
m
K
+
n
K
≥
0, then A
(
u
)
>
0 in (
0
,
u
¯
) and A
(
0
)
=
A
(
u
¯
)
=
0.
(2)
Ifm
n
+
m
K
+
n
K
<
0, then ∃
u
c
∈
(
0
,
u
*
) such that A
(
u
)
>
0 in (
u
c
,
u
¯
), A
(
u
)
<
0 in (
0
,
u
c
), and A
(
0
)
=
A
(
u
c
)
=
A
(
u
¯
)
=
0.
(3)
A
′
(
u
)
>
0 in (
max
⁡
{
0
,
u
c
}
,
u
*
), A
′
(
u
)
>
0 in (
u
*
,
u
¯
), and A
′
(
u
*
)
=
0.
Second,D
(
u
,
λ
)
=
0 is equal to λ
±
(
u
)
=
(
d
2
A
(
u
)
±
d
2
2
A
2
(
u
)
+
4
d
1
d
2
B
(
u
)
C
(
u
)
)
/
2
d
1
d
2. Let
(36)
S
(
u
)
=
d
2
2
A
2
(
u
)
+
4
d
1
d
1
B
(
u
)
C
(
u
)
.
Since
(37)
S
(
m
)
=
d
2
2
A
2
(
u
)
>
0
,
S
(
u
¯
)
=
4
d
1
d
2
B
(
u
¯
)
C
(
u
¯
)
<
0
according to the continuity of S
(
u
), there exists a u
s
∈
(
m
,
u
¯
), such that S
(
u
s
)
=
0. Thus, we can summarize the properties of λ as follows.Lemma 6.
Ifu
*
≤
m, then λ
+
(
u
) is decreasing and λ
-
(
u
) is increasing in u
∈
(
m
,
λ
s
).Proof.
DifferentiatingD
(
u
,
λ
)
=
0 with respect to u, we have
(38)
2
d
1
d
2
λ
(
u
)
λ
′
(
u
)
-
d
2
A
′
(
u
)
λ
(
u
)
-
d
2
A
(
u
)
λ
′
(
u
)
-
B
(
u
)
C
′
(
u
)
=
0
.
Therefore, λ
′
(
u
)
=
(
B
(
u
)
C
′
(
u
)
+
d
2
A
′
(
u
)
λ
(
u
)
)
/
d
2
(
2
d
1
λ
(
u
)
-
A
(
u
)
). From the definition of λ
±
(
u
), when u
∈
(
m
,
u
s
), d
2
(
2
d
1
λ
+
(
u
)
-
A
(
u
)
)
>
0 and d
2
(
2
d
1
λ
-
(
u
)
-
A
(
u
)
)
<
0.
However,C
(
u
)
=
e
r
(
1
-
(
u
/
K
)
)
(
1
-
(
(
m
+
n
)
/
(
u
+
n
)
)
), C
′
(
u
)
=
e
r
(
-
u
2
-
2
n
u
+
m
n
+
m
K
+
n
K
)
/
K
(
u
+
n
)
2, and C
′
(
u
)
=
0 imply that u
=
u
~
=
(
2
n
±
4
n
2
+
4
(
m
n
+
m
K
+
n
K
)
)
/
-
2, and
(39)
A
(
u
~
)
=
r
(
1
-
u
~
K
)
u
~
m
+
n
(
u
~
+
n
)
2
-
r
K
u
~
u
~
-
m
u
~
+
n
=
r
u
~
K
(
u
~
+
n
)
2
[
-
u
~
2
-
2
n
u
~
+
m
n
+
n
K
+
m
K
]
=
0
.
Since A
(
m
)
>
0 and m
<
u
~
<
K, then u
~
=
u
¯. Therefore, B
(
u
)
<
0, C
′
(
u
)
>
0, and λ
±
(
u
)
>
0 when u
∈
(
m
,
λ
s
). If u
*
≤
m, then A
′
(
u
)
≤
0 for u
∈
(
m
,
λ
s
). Thus, B
(
u
)
C
′
(
u
)
+
d
2
A
′
(
u
)
λ
±
(
u
)
<
0. This implies that λ
+
′
(
u
)
<
0.


### 3.1. Hopf Bifurcation Analysis

In this section, we mainly analyze the properties of the Hopf bifurcation for system (4). According to [38], we know that a Hopf bifurcation point u must satisfy the following conditions.(A)
There existsi
(
i
=
0,1
,
2
,
…
) such that
(40)
T
i
(
u
)
=
0
,
D
i
(
u
)
>
0
,
T
j
(
u
)
≠
0
,
D
j
(
u
)
≠
0
for
j
≠
i
and the unique pair of complex eigenvaluesα
(
u
)
±
i
ω
(
u
), α
(
u
)
=
0, ω
(
u
)
>
0 exist and are continuously differentiable in u, with α
′
(
u
)
≠
0.Theorem 7.
Ifu
*
≤
m holds, then k
∈
N exists such that the system (4) undergoes a Hopf bifurcation at u
=
u
i
(
0
≤
i
≤
k
) and a smooth curve of positive periodic orbits of (4) bifurcates from (
u
,
p
,
z
)
=
(
u
i
,
u
i
,
v
u
i
). The bifurcating periodic orbits from u
=
u
¯ are spatially homogeneous and the Hopf bifurcation at u
=
u
¯ is supercritical and backward if n
>
u
¯.Proof.
It can be verified thatT
i
(
u
)
<
0 and D
i
(
u
)
>
0 for u
∈
(
u
¯
,
K
), which implies that the potential Hopf bifurcation point must be in (
m
,
u
¯
]. However, T
0
(
u
¯
)
=
0, D
0
(
u
¯
)
>
0 and T
j
(
u
¯
)
=
-
(
d
1
+
d
2
)
λ
j
<
0,  D
j
(
u
¯
)
>
0 for j
≥
1. Therefore u
0
=
u
¯ is a Hopf bifurcation point. If u
*
≤
m holds, we know that A
(
u
) decreases strictly in (
m
,
u
¯
). For every i
>
0, let u
i be the solution of A
(
u
)
=
(
d
1
+
d
2
)
λ
i, so we have
(41)
m
<
u
h
<
u
h
-
1
<
⋯
<
u
1
<
u
¯
,
where λ
h is the largest eigenvalue for λ
i
<
A
(
u
)
/
(
d
1
+
d
2
). Therefore, T
i
(
u
i
)
=
0 and T
j
(
u
i
)
≠
0 for j
≠
i, 1
≤
i
≤
h. Geometrically, from T
i
(
u
)
=
0 we can determine that λ is a parabola of u with λ
(
u
¯
)
=
λ
(
u
c
)
=
λ
(
0
)
=
0 in (
u
,
λ
) coordinate system. However, T
(
m
,
λ
)
=
0 implies that λ
(
m
)
=
A
(
m
)
/
(
d
1
+
d
2
)
=
r
(
1
-
(
m
/
K
)
)
(
m
/
(
m
+
n
)
)
>
0, D
(
m
,
λ
)
=
0 implies that λ
=
0 or λ
(
m
)
=
A
(
m
)
/
d
1
>
A
(
m
)
/
(
d
1
+
d
2
), and Lemma 6 implies that λ
+
(
u
) is decreasing and λ
-
(
u
) is increasing in u
∈
(
m
,
λ
s
) if u
*
≤
m. Thus, the curves T
(
m
,
λ
)
=
0 and D
(
m
,
λ
)
=
0 have only one intersection point, which is noted as (
u
l
,
λ
l
) in u
∈
(
m
,
u
¯
). Then, D
i
(
u
i
)
>
0 if u
i
>
u
b and D
i
(
u
i
)
<
0 if u
i
<
u
b, 1
≤
i
≤
h. According to Theorem 2.1 in [13], a smooth curve of positive periodic orbits of (4) bifurcates from (
u
,
p
,
z
)
=
(
u
i
,
u
i
,
v
u
i
). According to Theorem 3.1 in [8],
(42)
a
(
u
¯
)
=
u
¯
f
′′
(
u
¯
)
16
(
f
′′′
(
u
¯
)
f
′′
(
u
¯
)
+
2
u
¯
)
=
-
u
¯
(
n
-
u
¯
)
(
2
n
2
+
m
n
+
K
n
+
K
m
)
8
(
u
¯
+
n
)
5
<
0
if n
>
u
¯ holds.


### 3.2. Steady State Bifurcation Analysis

In this section, we consider the steady state bifurcation of system (4). The nonnegative steady state solutions of (4) satisfy the following system:
(43)
-
d
1
Δ
p
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
,
x
∈
Ω
,
-
d
2
Δ
z
=
e
a
p
z
-
μ
z
,
x
∈
Ω
,
∂
p
∂
n
=
∂
z
∂
n
=
0
,
x
∈
∂
Ω
.
Apparently (43) has spatially homogeneous solutions (
0,0
), (
K
,
0
), (
m
,
0
), (
u
,
v
u
). First, we discuss the nonnegative steady state solutions of (4). Recall the maximum principle [13].Lemma 8.
LetΩ be a bounded Lipschitz domain in R
n, and let g
∈
C
(
Ω
¯
×
R
). If y
∈
W
1,2
(
Ω
) is a weak solution of the inequalities
(44)
Δ
y
+
g
(
x
,
y
(
x
)
)
≥
0
i
n
Ω
,
∂
y
(
x
)
∂
n
≤
0
o
n
∂
Ω
,
and if there is a constant M such that g
(
x
,
y
)
<
0 for y
>
M, then y
≤
M a.e. in Ω.From Lemma8, it can easily be derived that all nontrivial solutions of equation
(45)
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
+
d
1
Δ
p
=
0
,
x
∈
Ω
,
∂
p
∂
n
=
0
,
x
∈
∂
Ω
satisfy 0
≤
p
(
x
)
≤
K.Theorem 9.
The solutions of system (43) are in the form of either (
p
(
x
)
,
0
) or (
p
(
x
)
,
z
(
x
)
) satisfying
(46)
0
<
p
(
x
)
<
K
,
0
<
z
(
x
)
<
e
r
n
μ
(
K
-
m
2
)
2
-
d
1
K
2
≜
L
.Proof.
Ifx
0
∈
Ω
¯ exists such that z
(
x
0
)
=
0, from the strong maximum principle and the boundary condition ∂
z
/
∂
n
=
0, we can derive z
(
x
)
≡
0 in Ω
¯, so p
(
x
) satisfies (45). Similarly, if x
0
∈
Ω
¯ exists such that p
(
x
0
)
=
0, we can derive p
(
x
)
≡
0 and consequently z
(
x
)
≡
0. If, for all x
∈
Ω
¯, p
(
x
)
>
0 and z
(
x
)
>
0 hold, according to the above discussion and the strong maximum principle we have 0
<
p
(
x
)
<
K for all x
∈
Ω
¯. After adding the two equations in (43), we have
(47)
-
(
e
d
1
Δ
p
+
d
2
Δ
z
)
=
e
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
μ
z
=
e
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
μ
d
2
(
d
1
p
+
d
2
z
)
+
μ
d
1
d
2
p
≤
[
e
r
K
n
(
K
-
m
2
)
2
-
μ
d
1
d
2
]
K
-
μ
d
2
(
d
1
p
+
d
2
z
)
.
The maximum principle and Green formula imply that
(48)
d
2
z
<
d
1
p
+
d
2
z
<
d
2
μ
[
e
r
K
n
(
K
-
m
2
)
2
-
μ
d
1
2
]
K
.
Therefore, we can show the nonexistence of positive steady state solutions when the diffusion coefficients are large.Theorem 10.
LetD
1
>
(
K
2
+
m
K
-
m
n
)
K
/
n
2
λ
1 be a fixed constant. Then, another constant D
2 exists such that if d
1
≥
D
1 and d
2
≥
D
2 hold, (43) has no nonconstant positive solution.Proof.
Assume that(
p
(
x
)
,
z
(
x
)
) is a positive solution of (43). For convenience, we denote
(49)
f
1
(
p
,
z
)
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
,
f
2
(
p
,
z
)
=
e
a
p
z
-
μ
z
.
Let
(50)
p
¯
=
1
|
Ω
|
∫
Ω
p
(
x
)
d
x
,
z
¯
=
1
|
Ω
|
∫
Ω
z
(
x
)
d
x
;
then ∫
Ω
(
p
-
p
¯
)
d
x
=
0, ∫
Ω
(
z
-
z
¯
)
d
x
=
0. By multiplying (
p
-
p
¯
) by the first equation in (43) and then integrating on Ω, we have
(51)
∫
Ω
d
1
|
∇
p
|
2
d
x
=
∫
Ω
f
1
(
p
,
z
)
(
p
-
p
¯
)
d
x
=
∫
Ω
(
f
1
(
p
,
z
)
-
f
1
(
p
¯
,
z
¯
)
)
(
p
-
p
¯
)
d
x
=
∫
Ω
(
(
-
p
p
¯
(
p
+
p
¯
)
(
p
-
p
¯
)
2
+
(
K
+
m
)
p
p
¯
(
p
-
p
¯
)
2
-
K
m
n
(
p
-
p
¯
)
)
×
(
(
p
+
n
)
(
p
¯
+
n
)
)
-
1
-
a
p
(
p
-
p
¯
)
(
z
-
z
¯
)
-
a
z
¯
(
z
-
z
¯
)
2
(
p
-
p
¯
)
2
)
d
x
≤
∫
Ω
(
(
K
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
+
a
K
|
p
-
p
¯
|
|
z
-
z
¯
|
(
K
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
)
d
x
.
Similarly,
(52)
∫
Ω
d
2
|
∇
z
|
2
d
x
=
∫
Ω
f
2
(
p
,
z
)
(
z
-
z
¯
)
d
x
=
∫
Ω
(
f
2
(
p
,
z
)
-
f
2
(
p
¯
,
z
¯
)
)
(
z
-
z
¯
)
d
x
=
∫
Ω
(
(
e
a
p
-
μ
)
(
z
-
z
¯
)
2
-
e
a
z
¯
(
p
-
p
¯
)
(
z
-
z
¯
)
)
d
x
≤
∫
Ω
(
(
e
a
K
-
μ
)
(
z
-
z
¯
)
2
+
e
a
L
|
p
-
p
¯
|
|
z
-
z
¯
|
)
d
x
,
where L is defined in Theorem 9.
Therefore,(53)
∫
Ω
(
d
1
|
∇
p
|
2
+
d
2
|
∇
z
|
2
)
d
x
≤
∫
Ω
(
(
K
2
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
+
(
a
K
+
e
a
L
)
|
p
-
p
¯
|
|
z
-
z
¯
|
+
(
e
a
K
-
μ
)
(
z
-
z
¯
)
2
(
K
2
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
)
d
x
≤
∫
Ω
(
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
(
p
-
p
¯
)
2
+
(
e
a
L
+
a
K
2
ε
+
e
a
K
-
μ
)
(
z
-
z
¯
)
2
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
)
d
x
.
Using the Poincáre inequality, we can obtain
(54)
∫
Ω
(
d
1
λ
1
|
p
-
p
¯
|
2
+
d
2
λ
1
|
z
-
z
¯
|
2
)
d
x
≤
∫
Ω
(
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
(
p
-
p
¯
)
2
+
(
e
a
L
+
a
K
2
ε
+
e
a
K
-
μ
)
(
z
-
z
¯
)
2
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
)
d
x
.
Under the assumption d
1
λ
1
>
(
K
2
+
m
K
-
m
n
)
K
/
n
2, a sufficiently small ε
>
0 exists such that d
1
λ
1
≥
(
K
2
+
m
K
-
m
n
)
K
/
n
2
+
(
(
a
K
+
e
a
L
)
/
2
)
ε. Let D
2
=
(
(
e
a
L
+
a
K
)
/
2
ε
+
e
a
K
-
μ
)
(
1
/
λ
1
); then we can have p
=
p
¯, z
=
z
¯. This completes the proof.We assume that all of the eigenvalues of-
Δ are simple and the corresponding eigenfunctions are denoted by ϕ
i
(
x
). Reference [13] gives an example of n
=
1 with λ
i
=
i
2
/
l
2 and ϕ
i
(
x
)
=
cos
⁡
(
i
x
/
l
). Let
(55)
T
i
(
u
)
=
-
λ
i
(
d
1
+
d
2
)
+
A
(
u
)
,
D
i
(
u
)
=
d
1
d
2
λ
i
2
-
A
(
u
)
d
2
λ
i
-
B
(
u
)
C
(
u
)
.
According to [13], we know that a steady state bifurcation point u must satisfy the following conditions.(B)
i
(
i
=
0,1
,
2
,
…
) exists such that
(56)
D
i
(
u
)
=
0
,
T
i
(
u
)
≠
0
,
D
j
(
u
)
≠
0
,
T
j
(
u
)
≠
0
,
for
j
≠
i
,
d
D
i
(
u
)
d
u
≠
0
.Theorem 11.
k
′
∈
N exists such that system (43) undergoes a steady state bifurcation at u
=
u
i
(
1
≤
i
≤
k
′
) if u
*
≤
m holds.Proof.
Apparently (B) is not established fori
=
0. It is known that D
i
(
u
) is a degree 3 polynomial of u and there are at most three u
i for D
i
(
u
i
)
=
D
(
u
i
,
λ
i
)
=
0. In particular, if the parameters are selected such that  u
*
≤
m holds, then, for each i
>
0, there exists a unique u
i
∈
(
m
,
u
s
) such that D
i
(
u
i
)
=
D
(
u
i
,
λ
i
)
=
0. Thus, there is at most one bifurcation point. However, D
i
′
(
u
)
=
-
A
(
u
)
d
2
λ
i
-
B
(
u
)
C
′
(
u
). From the proof of Lemma 6 we know that B
(
u
)
C
′
(
u
)
+
d
2
A
′
(
u
)
λ
±
(
u
)
<
0 for u
∈
(
m
,
λ
s
). Therefore, d
D
i
(
u
)
/
d
u
≠
0 holds if u
*
≤
m.


## 3.1. Hopf Bifurcation Analysis

In this section, we mainly analyze the properties of the Hopf bifurcation for system (4). According to [38], we know that a Hopf bifurcation point u must satisfy the following conditions.(A)
There existsi
(
i
=
0,1
,
2
,
…
) such that
(40)
T
i
(
u
)
=
0
,
D
i
(
u
)
>
0
,
T
j
(
u
)
≠
0
,
D
j
(
u
)
≠
0
for
j
≠
i
and the unique pair of complex eigenvaluesα
(
u
)
±
i
ω
(
u
), α
(
u
)
=
0, ω
(
u
)
>
0 exist and are continuously differentiable in u, with α
′
(
u
)
≠
0.Theorem 7.
Ifu
*
≤
m holds, then k
∈
N exists such that the system (4) undergoes a Hopf bifurcation at u
=
u
i
(
0
≤
i
≤
k
) and a smooth curve of positive periodic orbits of (4) bifurcates from (
u
,
p
,
z
)
=
(
u
i
,
u
i
,
v
u
i
). The bifurcating periodic orbits from u
=
u
¯ are spatially homogeneous and the Hopf bifurcation at u
=
u
¯ is supercritical and backward if n
>
u
¯.Proof.
It can be verified thatT
i
(
u
)
<
0 and D
i
(
u
)
>
0 for u
∈
(
u
¯
,
K
), which implies that the potential Hopf bifurcation point must be in (
m
,
u
¯
]. However, T
0
(
u
¯
)
=
0, D
0
(
u
¯
)
>
0 and T
j
(
u
¯
)
=
-
(
d
1
+
d
2
)
λ
j
<
0,  D
j
(
u
¯
)
>
0 for j
≥
1. Therefore u
0
=
u
¯ is a Hopf bifurcation point. If u
*
≤
m holds, we know that A
(
u
) decreases strictly in (
m
,
u
¯
). For every i
>
0, let u
i be the solution of A
(
u
)
=
(
d
1
+
d
2
)
λ
i, so we have
(41)
m
<
u
h
<
u
h
-
1
<
⋯
<
u
1
<
u
¯
,
where λ
h is the largest eigenvalue for λ
i
<
A
(
u
)
/
(
d
1
+
d
2
). Therefore, T
i
(
u
i
)
=
0 and T
j
(
u
i
)
≠
0 for j
≠
i, 1
≤
i
≤
h. Geometrically, from T
i
(
u
)
=
0 we can determine that λ is a parabola of u with λ
(
u
¯
)
=
λ
(
u
c
)
=
λ
(
0
)
=
0 in (
u
,
λ
) coordinate system. However, T
(
m
,
λ
)
=
0 implies that λ
(
m
)
=
A
(
m
)
/
(
d
1
+
d
2
)
=
r
(
1
-
(
m
/
K
)
)
(
m
/
(
m
+
n
)
)
>
0, D
(
m
,
λ
)
=
0 implies that λ
=
0 or λ
(
m
)
=
A
(
m
)
/
d
1
>
A
(
m
)
/
(
d
1
+
d
2
), and Lemma 6 implies that λ
+
(
u
) is decreasing and λ
-
(
u
) is increasing in u
∈
(
m
,
λ
s
) if u
*
≤
m. Thus, the curves T
(
m
,
λ
)
=
0 and D
(
m
,
λ
)
=
0 have only one intersection point, which is noted as (
u
l
,
λ
l
) in u
∈
(
m
,
u
¯
). Then, D
i
(
u
i
)
>
0 if u
i
>
u
b and D
i
(
u
i
)
<
0 if u
i
<
u
b, 1
≤
i
≤
h. According to Theorem 2.1 in [13], a smooth curve of positive periodic orbits of (4) bifurcates from (
u
,
p
,
z
)
=
(
u
i
,
u
i
,
v
u
i
). According to Theorem 3.1 in [8],
(42)
a
(
u
¯
)
=
u
¯
f
′′
(
u
¯
)
16
(
f
′′′
(
u
¯
)
f
′′
(
u
¯
)
+
2
u
¯
)
=
-
u
¯
(
n
-
u
¯
)
(
2
n
2
+
m
n
+
K
n
+
K
m
)
8
(
u
¯
+
n
)
5
<
0
if n
>
u
¯ holds.


## 3.2. Steady State Bifurcation Analysis

In this section, we consider the steady state bifurcation of system (4). The nonnegative steady state solutions of (4) satisfy the following system:
(43)
-
d
1
Δ
p
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
,
x
∈
Ω
,
-
d
2
Δ
z
=
e
a
p
z
-
μ
z
,
x
∈
Ω
,
∂
p
∂
n
=
∂
z
∂
n
=
0
,
x
∈
∂
Ω
.
Apparently (43) has spatially homogeneous solutions (
0,0
), (
K
,
0
), (
m
,
0
), (
u
,
v
u
). First, we discuss the nonnegative steady state solutions of (4). Recall the maximum principle [13].Lemma 8.
LetΩ be a bounded Lipschitz domain in R
n, and let g
∈
C
(
Ω
¯
×
R
). If y
∈
W
1,2
(
Ω
) is a weak solution of the inequalities
(44)
Δ
y
+
g
(
x
,
y
(
x
)
)
≥
0
i
n
Ω
,
∂
y
(
x
)
∂
n
≤
0
o
n
∂
Ω
,
and if there is a constant M such that g
(
x
,
y
)
<
0 for y
>
M, then y
≤
M a.e. in Ω.From Lemma8, it can easily be derived that all nontrivial solutions of equation
(45)
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
+
d
1
Δ
p
=
0
,
x
∈
Ω
,
∂
p
∂
n
=
0
,
x
∈
∂
Ω
satisfy 0
≤
p
(
x
)
≤
K.Theorem 9.
The solutions of system (43) are in the form of either (
p
(
x
)
,
0
) or (
p
(
x
)
,
z
(
x
)
) satisfying
(46)
0
<
p
(
x
)
<
K
,
0
<
z
(
x
)
<
e
r
n
μ
(
K
-
m
2
)
2
-
d
1
K
2
≜
L
.Proof.
Ifx
0
∈
Ω
¯ exists such that z
(
x
0
)
=
0, from the strong maximum principle and the boundary condition ∂
z
/
∂
n
=
0, we can derive z
(
x
)
≡
0 in Ω
¯, so p
(
x
) satisfies (45). Similarly, if x
0
∈
Ω
¯ exists such that p
(
x
0
)
=
0, we can derive p
(
x
)
≡
0 and consequently z
(
x
)
≡
0. If, for all x
∈
Ω
¯, p
(
x
)
>
0 and z
(
x
)
>
0 hold, according to the above discussion and the strong maximum principle we have 0
<
p
(
x
)
<
K for all x
∈
Ω
¯. After adding the two equations in (43), we have
(47)
-
(
e
d
1
Δ
p
+
d
2
Δ
z
)
=
e
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
μ
z
=
e
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
μ
d
2
(
d
1
p
+
d
2
z
)
+
μ
d
1
d
2
p
≤
[
e
r
K
n
(
K
-
m
2
)
2
-
μ
d
1
d
2
]
K
-
μ
d
2
(
d
1
p
+
d
2
z
)
.
The maximum principle and Green formula imply that
(48)
d
2
z
<
d
1
p
+
d
2
z
<
d
2
μ
[
e
r
K
n
(
K
-
m
2
)
2
-
μ
d
1
2
]
K
.
Therefore, we can show the nonexistence of positive steady state solutions when the diffusion coefficients are large.Theorem 10.
LetD
1
>
(
K
2
+
m
K
-
m
n
)
K
/
n
2
λ
1 be a fixed constant. Then, another constant D
2 exists such that if d
1
≥
D
1 and d
2
≥
D
2 hold, (43) has no nonconstant positive solution.Proof.
Assume that(
p
(
x
)
,
z
(
x
)
) is a positive solution of (43). For convenience, we denote
(49)
f
1
(
p
,
z
)
=
r
(
1
-
p
K
)
(
1
-
m
+
n
p
+
n
)
p
-
a
p
z
,
f
2
(
p
,
z
)
=
e
a
p
z
-
μ
z
.
Let
(50)
p
¯
=
1
|
Ω
|
∫
Ω
p
(
x
)
d
x
,
z
¯
=
1
|
Ω
|
∫
Ω
z
(
x
)
d
x
;
then ∫
Ω
(
p
-
p
¯
)
d
x
=
0, ∫
Ω
(
z
-
z
¯
)
d
x
=
0. By multiplying (
p
-
p
¯
) by the first equation in (43) and then integrating on Ω, we have
(51)
∫
Ω
d
1
|
∇
p
|
2
d
x
=
∫
Ω
f
1
(
p
,
z
)
(
p
-
p
¯
)
d
x
=
∫
Ω
(
f
1
(
p
,
z
)
-
f
1
(
p
¯
,
z
¯
)
)
(
p
-
p
¯
)
d
x
=
∫
Ω
(
(
-
p
p
¯
(
p
+
p
¯
)
(
p
-
p
¯
)
2
+
(
K
+
m
)
p
p
¯
(
p
-
p
¯
)
2
-
K
m
n
(
p
-
p
¯
)
)
×
(
(
p
+
n
)
(
p
¯
+
n
)
)
-
1
-
a
p
(
p
-
p
¯
)
(
z
-
z
¯
)
-
a
z
¯
(
z
-
z
¯
)
2
(
p
-
p
¯
)
2
)
d
x
≤
∫
Ω
(
(
K
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
+
a
K
|
p
-
p
¯
|
|
z
-
z
¯
|
(
K
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
)
d
x
.
Similarly,
(52)
∫
Ω
d
2
|
∇
z
|
2
d
x
=
∫
Ω
f
2
(
p
,
z
)
(
z
-
z
¯
)
d
x
=
∫
Ω
(
f
2
(
p
,
z
)
-
f
2
(
p
¯
,
z
¯
)
)
(
z
-
z
¯
)
d
x
=
∫
Ω
(
(
e
a
p
-
μ
)
(
z
-
z
¯
)
2
-
e
a
z
¯
(
p
-
p
¯
)
(
z
-
z
¯
)
)
d
x
≤
∫
Ω
(
(
e
a
K
-
μ
)
(
z
-
z
¯
)
2
+
e
a
L
|
p
-
p
¯
|
|
z
-
z
¯
|
)
d
x
,
where L is defined in Theorem 9.
Therefore,(53)
∫
Ω
(
d
1
|
∇
p
|
2
+
d
2
|
∇
z
|
2
)
d
x
≤
∫
Ω
(
(
K
2
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
+
(
a
K
+
e
a
L
)
|
p
-
p
¯
|
|
z
-
z
¯
|
+
(
e
a
K
-
μ
)
(
z
-
z
¯
)
2
(
K
2
+
m
K
-
m
n
)
K
n
2
(
p
-
p
¯
)
2
)
d
x
≤
∫
Ω
(
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
(
p
-
p
¯
)
2
+
(
e
a
L
+
a
K
2
ε
+
e
a
K
-
μ
)
(
z
-
z
¯
)
2
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
)
d
x
.
Using the Poincáre inequality, we can obtain
(54)
∫
Ω
(
d
1
λ
1
|
p
-
p
¯
|
2
+
d
2
λ
1
|
z
-
z
¯
|
2
)
d
x
≤
∫
Ω
(
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
(
p
-
p
¯
)
2
+
(
e
a
L
+
a
K
2
ε
+
e
a
K
-
μ
)
(
z
-
z
¯
)
2
(
(
K
2
+
m
K
-
m
n
)
K
n
2
+
a
K
+
e
a
L
2
ε
)
)
d
x
.
Under the assumption d
1
λ
1
>
(
K
2
+
m
K
-
m
n
)
K
/
n
2, a sufficiently small ε
>
0 exists such that d
1
λ
1
≥
(
K
2
+
m
K
-
m
n
)
K
/
n
2
+
(
(
a
K
+
e
a
L
)
/
2
)
ε. Let D
2
=
(
(
e
a
L
+
a
K
)
/
2
ε
+
e
a
K
-
μ
)
(
1
/
λ
1
); then we can have p
=
p
¯, z
=
z
¯. This completes the proof.We assume that all of the eigenvalues of-
Δ are simple and the corresponding eigenfunctions are denoted by ϕ
i
(
x
). Reference [13] gives an example of n
=
1 with λ
i
=
i
2
/
l
2 and ϕ
i
(
x
)
=
cos
⁡
(
i
x
/
l
). Let
(55)
T
i
(
u
)
=
-
λ
i
(
d
1
+
d
2
)
+
A
(
u
)
,
D
i
(
u
)
=
d
1
d
2
λ
i
2
-
A
(
u
)
d
2
λ
i
-
B
(
u
)
C
(
u
)
.
According to [13], we know that a steady state bifurcation point u must satisfy the following conditions.(B)
i
(
i
=
0,1
,
2
,
…
) exists such that
(56)
D
i
(
u
)
=
0
,
T
i
(
u
)
≠
0
,
D
j
(
u
)
≠
0
,
T
j
(
u
)
≠
0
,
for
j
≠
i
,
d
D
i
(
u
)
d
u
≠
0
.Theorem 11.
k
′
∈
N exists such that system (43) undergoes a steady state bifurcation at u
=
u
i
(
1
≤
i
≤
k
′
) if u
*
≤
m holds.Proof.
Apparently (B) is not established fori
=
0. It is known that D
i
(
u
) is a degree 3 polynomial of u and there are at most three u
i for D
i
(
u
i
)
=
D
(
u
i
,
λ
i
)
=
0. In particular, if the parameters are selected such that  u
*
≤
m holds, then, for each i
>
0, there exists a unique u
i
∈
(
m
,
u
s
) such that D
i
(
u
i
)
=
D
(
u
i
,
λ
i
)
=
0. Thus, there is at most one bifurcation point. However, D
i
′
(
u
)
=
-
A
(
u
)
d
2
λ
i
-
B
(
u
)
C
′
(
u
). From the proof of Lemma 6 we know that B
(
u
)
C
′
(
u
)
+
d
2
A
′
(
u
)
λ
±
(
u
)
<
0 for u
∈
(
m
,
λ
s
). Therefore, d
D
i
(
u
)
/
d
u
≠
0 holds if u
*
≤
m.


## 4. Discussion

Reaction-diffusion phytoplankton-zooplankton models with Allee effects have been studied extensively in recent years. In this study, we rigorously considered a Gause-type predator-prey model with a double Allee effect on prey, which was formulated as (4). It is known that the predator-prey model with the most usual form of Allee effect has a unique limit cycle, but the existence of two limit cycles was proved by González-Olivares et al. [29] with a double Allee effect. Thus, the double Allee effect produces different results with different mathematical expressions.The paper [15] found that the system without Allee effect was always stable and without fluctuations, but in this paper the results of the stability of the equilibrium and the bifurcation analysis based on a rigorous theoretical analysis show that this system has complex spatiotemporal dynamics: for p
0
(
x
)
≤
m, the phytoplankton is destined to become extinct and leads to the extinction of zooplankton; after considering the strong Allee effect in phytoplankton, extinction for both species is always a locally stable equilibrium. But for u
>
m, which is the condition in which the interior equilibrium exists, the interior equilibrium is globally stable in some case and there always exist some other spatiotemporal patterns in other cases (Figure 3).Time evolution of (4) around the interior equilibrium (
8,2.30769
). Phytoplankton is in (a) and zooplankton is in (b). The parameter values are d
1
=
0.001, d
2
=
0.004, e
=
0.5, a
=
0.1, n
=
5, m
=
0.5, μ
=
0.4, r
=
2, and K
=
10. The figure shows that both types of plankton exhibit stable behavior.
(a)
(b)Figure 4
The graph ofT
(
u
,
λ
)
=
0 and D
(
u
,
λ
)
=
0, where d
1
=
0.001, d
2
=
0.004, e
=
0.5, a
=
1, n
=
5, m
=
5, μ
=
0.4, r
=
2, and K
=
10. In this case, u
*
≤
m.Figure 5
The graph ofT
(
u
,
λ
)
=
0 and D
(
u
,
λ
)
=
0, where d
1
=
0.001, d
2
=
0.004, e
=
0.5, a
=
1, n
=
5, m
=
1, μ
=
0.4, r
=
2, and K
=
10. In this case, u
*
>
m.Overall, our results indicate that the impact of the Allee effect increases the spatiotemporal complexity of the system. The mathematical form which expresses the double Allee effect has a strong impact of the dynamics of system. Thus, we think it is important for ecologists to be aware of the difference of the selection on the forms of Allee effect.The limitations of our study are that we only consider a simplest phytoplankton-zooplankton interaction and the special formalization to describe the Allee effect. What is more is that, compared with the ODE dynamics, the results shown here are still coarse. Therefore, further research is still needed to elaborate a general theory on the influence of this ecological phenomenon.



---

*Source: 952840-2014-07-22.xml*