# Stability and Selective Harvesting of a Phytoplankton-Zooplankton System

Authors: Yong Wang; Hongbin Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Considering that some zooplankton can be harvested for food in some bodies of water, a phytoplankton-zooplankton model with continuous harvesting of zooplankton only is proposed and investigated. By using environmental carrying capacity as a parameter, possible dynamic behaviors, such as stability, global stability, Hopf bifurcation, and transcritical bifurcations, are analyzed. The optimal harvesting policy is disposed by imposing a tax per unit biomass of zooplankton. The problem of determining the optimal harvest policy is solved by using Pontryagin's maximum principle subject to the state equations and the control constraints, and the impact of tax is also discussed. Finally, some numerical simulations are performed to justify analytical findings.


---

## Body

## 1. Introduction

In marine ecology, the term plankton refers to the freely floating and weakly swimming organisms in some bodies of fresh water. There are two types of plankton: the plant species, commonly known as phytoplankton, which are unicellular and microscopic in size and use water and absorb carbon-dioxide from the air to grow, and the animal species, namely, zooplankton, which live on these phytoplankton. They are the basis of all aquatic food chain and eaten by many organisms, including mussels, fish, mammals, and people (see Figure1). In aquatic food webs, some zooplankton such as jellyfish, krill, andAcetes are harvested for food. Hence, the stocks of these tiny zooplankton play a significant role in marine reserves and fishery management.Figure 1
The ocean’s food chain.During the recent years, the problems of zooplankton-phytoplankton system have been discussed by many authors [1–9]. Chattopadhayay et al. [1, 2] investigated that toxin producing phytoplankton affected the growth of zooplankton population and had an impact on phytoplankton and zooplankton interaction. He and Ruan [3] considered plankton-nutrient interaction models and obtained some sufficient conditions for the global attractivity of the positive equilibrium. Das and Ray [4] considered phytoplankton and zooplankton interactions with delayed nutrient cycling from senescence and mortality of phytoplankton in Hooghly-Matla estuarine system. Roy [5] constructed a mathematical model for describing the interaction between a nontoxic and a toxic phytoplankton under a single nutrient. Gakkhar and Singh [6] proposed and analyzed an ecoepidemiological delay model for virally infected, toxin producing phytoplankton-zooplankton system, and so forth.But little attention has been paid to study the model about the effect of harvest on phytoplankton and zooplankton populations. The growing human need for more food and energy has led to increase in exploitation of several biological resources; during the past half century, the amount of the world’s fish has been greatly reduced. On the other hand, there is a global concern to protect the ecosystem at large. In the face of two opposing approaches, interest in renewable resources has increased greatly in recent years [10–12]. Determining socially acceptable harvesting policy is undoubtedly one of the most challenging and most controversial problems in the management of renewable resources. As is well known, the optimal harvesting problem results in a direct relationship to sustainable development. Taxation, lease of property rights, and seasonal harvesting are usually considered as possible governing instruments in fishery regulation. Economists are particularly attracted to taxation because an ecosystem can be better maintained under taxation rather than other regulatory methods. But little attention has been paid to study the dynamics of fishery resources using taxation as a control instrument; harvesting problems with tax have been studied by [13–16]. In this paper, in order to gain both the mathematical and biological generality, we use taxation as a control variable.The model we considered is based on the following plausible toxic phytoplankton and zooplankton system introduced by Chattopadhayay et al. [1] and Saha and Bandyopadhyay [7]:
(1)
d
P
d
t
=
r
P
(
1
-
P
k
)
-
β
P
Z
α
+
P
,
d
Z
d
t
=
β
1
P
Z
α
+
P
-
d
Z
-
ρ
P
Z
α
+
P
.
The following assumptions for model (1) are made: (1)
the variablesP and Z are the density of phytoplankton population and the density of zooplankton population at any instant of time t, respectively;
(2)
the parameterr is the intrinsic growth rate and k is the environmental carrying capacity of population. The constant β
(
>
0
) is the maximum uptake rate for zooplankton species; β
1
(
>
0
) denotes the ratio of biomass conversion (satisfying the obvious restriction 0
<
β
1
<
β) and d
(
>
0
) is the natural death rate of zooplankton;
(3)
the parameterρ
(
>
0
) denotes the rate of toxic substances produced per unit biomass of phytoplankton. It is assumed that β
1
<
ρ; that is, the ratio of biomass consumed by zooplankton is greater than the rate of toxic substance liberation by phytoplankton species;
(4)
the term(
β
P
Z
/
(
α
+
P
)
) represents the functional response for the grazing of phytoplankton by zooplankton and α is the half saturation constant for a Holling type II functional response. (
ρ
P
Z
/
(
α
+
P
)
) describes the distribution of toxic substance which ultimately contributes to the death of zooplankton populations.Now, we adapt model (1) and assume that zooplankton are subject to a harvesting effort governed by the differential equations
(2)
d
P
d
t
=
r
P
(
1
-
P
k
)
-
β
P
Z
α
+
P
,
d
Z
d
t
=
β
1
P
Z
α
+
P
-
d
Z
-
ρ
P
Z
α
+
P
-
c
1
E
Z
,
where the constant c
1 is the catchability coefficient of the zooplankton. E is harvesting effort.Here, we takeE as dynamic [9, 12] (i.e., time-dependent) variable governed by the equations
(3)
E
(
t
)
=
μ
1
Q
(
t
)
,
0
≤
μ
1
≤
1
,
d
Q
d
t
=
I
(
t
)
-
γ
Q
(
t
)
,
Q
(
0
)
=
Q
0
,
where I
(
t
) is the gross investment rate at time t. Q
(
t
) is the amount of capital invested in the fishery at time t; γ is constant rate of depreciation of capital. A regulatory agency controls exploitation of the fishery by imposing a tax τ
(
>
0
) per unit biomass of the landed fish. τ
<
0 denotes the subsidy given to the fisherman. The net economic revenue to the fisherman is {
c
1
(
p
-
τ
)
Z
-
C
}
E, where p is the constant price per unit biomass of the zooplankton species and C is the constant cost per unit of harvesting effort.We assume that the gross rate of investment of capital is proportional to the net economic revenue to the fisherman. Thus we have(4)
I
(
t
)
=
μ
2
{
c
1
(
p
-
τ
)
Z
-
C
}
E
(
t
)
,
0
≤
μ
2
≤
1
.
Equation (4) asserts that the maximum investment rate at any time equals the net economic revenue (for μ
2
=
1) at that time. By virtue of (4) and (3) yield the result
(5)
d
E
d
t
=
{
μ
1
μ
2
[
c
1
(
p
-
τ
)
Z
-
C
]
-
γ
}
E
.
Let m
=
μ
1
μ
2; therefore, we have the following system of equations:
(6)
d
P
d
t
=
r
P
(
1
-
P
k
)
-
β
P
Z
α
+
P
,
d
Z
d
t
=
β
1
P
Z
α
+
P
-
d
Z
-
ρ
P
Z
α
+
P
-
c
1
E
Z
,
d
E
d
t
=
{
m
[
c
1
(
p
-
τ
)
Z
-
C
]
-
γ
}
E
,
P
(
0
)
=
P
0
>
0
,
Z
(
0
)
=
Z
0
>
0
,
E
(
0
)
=
E
0
>
0
.In this paper, we choose tax as the management objective when discussing the impact of harvesting in the phytoplankton-zooplankton system and assume(
H
)
β
1
-
ρ
-
d
>
0, p
-
τ
>
0. The organization of this paper is as follows: to begin with, we construct and briefly describe our model. Then we study the local and global stabilities and bifurcations of the equilibria in the next section. In Section 3, we study the bifurcation phenomenon. In Section 4, by using Pontryagin’s maximum principle, we analyze the optimal tax policy. In order to illustrate our result, some numerical simulations are given in Section 5. We end the paper with a brief conclusion in Section 6.


## 2. Stability of the Equilibria

The system of system (6) has four feasible equilibria.(1)
The equilibrium pointsS
0
=
(
0,0
,
0
) and S
1
=
(
k
,
0,0
) exist for all permissible parameters.
(2)
Ifk
>
k
t
1
=
(
d
α
/
(
β
1
-
ρ
-
d
)
), the boundary equilibrium S
2
=
(
P
*
,
Z
*
,
0
) exists, where
(7)
P
*
=
d
α
β
1
-
ρ
-
d
,
Z
*
=
r
β
(
α
+
P
*
)
(
1
-
P
*
k
)
.
(3)
LetS
3
=
(
P
1
*
,
Z
1
*
,
E
1
*
) be the positive interior equilibrium, where
(8)
Z
1
*
=
m
C
+
γ
m
c
1
(
p
-
τ
)
,
E
1
*
=
(
β
1
-
ρ
)
P
1
*
c
1
(
α
+
P
1
*
)
-
d
c
1and P
1
* satisfies
(9)
(
P
1
*
)
2
+
(
α
-
k
)
P
1
*
-
k
α
+
k
β
(
m
C
+
γ
)
r
m
c
1
(
p
-
τ
)
=
0
.
Let x
1
,
x
2 be the roots of (9); we only consider that x
1
,
x
2 have only one positive root; then
(10)
x
1
x
2
=
-
k
α
+
k
β
(
m
C
+
γ
)
r
m
c
1
(
p
-
τ
)
<
0
kkkl
l
⟹
τ
<
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
α
m
r
c
1
,
l
l
l
l
Δ
=
(
α
-
k
)
2
+
4
(
k
α
-
k
β
(
m
C
+
γ
)
r
m
c
1
(
p
-
τ
)
)
>
0
,
and hence, P
1
* exists as a positive root:
(11)
P
1
*
=
1
2
[
k
-
α
+
Δ
]
.
We know that if k
>
k
t
2,  Z
1
*
>
0 and E
1
*
>
0, where k
t
2
=
(
(
d
2
α
2
/
(
β
1
-
ρ
-
d
)
2
)
+
(
d
α
2
/
(
β
1
-
ρ
-
d
)
)
)
/
(
α
+
(
d
α
/
(
β
1
-
ρ
-
d
)
)
-
(
β
(
m
C
+
γ
)
/
r
m
c
1
(
p
-
τ
)
)
).Therefore, whenk
>
k
t
2 together with τ
<
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1, the positive interior equilibrium S
3 exists.


### 2.1. Local Stability

In order to determine the stability of (6), we compute the Jacobian matrix of system (6):(12)
J
=
(
r
-
2
r
P
k
-
β
α
Z
(
α
+
P
)
2
-
β
P
α
+
P
0
(
β
1
-
ρ
)
α
Z
(
α
+
P
)
2
(
β
1
-
ρ
)
P
α
+
P
-
(
d
+
c
1
E
)
-
c
1
Z
0
m
c
1
(
p
-
τ
)
E
m
[
c
1
(
p
-
τ
)
Z
-
C
]
-
γ
)
.Theorem 1.
For system (6), we have the following: (1) the extinction equilibrium S
0 is unstable;
(2) when k
<
k
t
1, the equilibrium S
1 is locally asymptotically stable;
(3)
(i) if condition (10) is satisfied, then the equilibrium S
2
=
(
P
*
,
Z
*
,
0
) is locally asymptotically stable when k
t
1
<
k
<
k
t
2;
(ii) if condition (10) is not satisfied, then the equilibrium S
2
=
(
P
*
,
Z
*
,
0
) is locally asymptotically stable when 0
<
Z
*
<
(
m
C
+
γ
)
/
m
c
1
(
p
-
τ
), k
<
k
h
=
α
+
2
P
*;
(4) if system (6) has only one positive equilibrium S
3, then S
3 is locally asymptotically stable.Proof.
(
1
) For S
0, the eigenvalues of the matrix (12) are λ
1
=
r
>
0, λ
2
=
-
d
<
0, and λ
3
=
-
m
C
-
γ
<
0. It is seen that there are one unstable manifold and two stable manifolds. Therefore, the point S
0 is a saddle point.
(
2
) For S
1, the characteristic equation is
(13)
|
-
r
-
λ
-
β
k
α
+
k
0
0
(
β
1
-
ρ
)
k
α
+
k
-
d
-
λ
0
0
0
-
m
C
-
γ
-
λ
|
=
0
of which the roots are
(14)
λ
1
=
-
r
<
0
,
λ
2
=
(
β
1
-
ρ
)
k
α
+
k
-
d
,
λ
3
=
-
m
C
-
γ
<
0
.
Hence, the equilibrium S
1 is locally asymptotically stable, when k
<
k
t
1.
(
3
) For S
2, the characteristic equation of S
2 is given by
(15)
{
λ
-
[
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
]
}
×
{
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
λ
[
λ
-
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
]
K
K
K
K
+
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
}
=
0
.
Let λ
1
,
λ
2, and λ
3 be its eigenvalues; assume λ
1
,
λ
2, and λ
3 are all negative; hence,
(16)
λ
1
=
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
<
0
⟹
0
<
Z
*
<
m
C
+
γ
m
c
1
(
p
-
τ
)
.
(i)
Assuming that (10) is satisfied, then 0
<
Z
*
<
(
m
C
+
γ
)
/
m
c
1
(
p
-
τ
) is equivalent to 0
<
k
<
k
t
2. λ
2
,
λ
3 satisfy
(17)
λ
2
-
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
λ
+
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
=
0
.
Therefore(18)
λ
2
+
λ
3
=
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
<
0
kkkkkk
⟹
k
<
k
h
=
α
+
2
P
*
,
k
k
λ
2
λ
3
=
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
>
0
.
Obviously, we know thatk
t
1
<
k
t
2
<
k
h. Therefore,
(19)
λ
1
<
0
,
λ
2
<
0
,
λ
3
<
0
,
0
<
k
<
k
t
2
,
λ
1
=
0
,
λ
2
<
0
,
λ
3
<
0
,
k
=
k
t
2
,
λ
1
>
0
,
λ
2
<
0
,
λ
3
<
0
,
k
t
2
<
k
<
k
h
,
λ
1
>
0
,
λ
2
,
λ
3
=
±
i
ω
0
(
ω
0
>
0
)
,
k
=
k
h
,
λ
1
>
0
,
λ
2
>
0
,
λ
3
>
0
,
k
>
k
h
.
(ii)
Assuming that (10) is not satisfied, then λ
1
<
0, λ
2
<
0, and λ
3
<
0 are equivalent to 0
<
Z
*
<
(
m
C
+
γ
)
/
m
c
1
(
p
-
τ
), k
<
k
h. Therefore, from the paragraph above, combined with the condition of the existence for equilibrium S
2
=
(
P
*
,
Z
*
,
0
), one gets the result.
(
4
) For S
3, the characteristic equation of S
3 is the following:
(20)
λ
3
+
A
0
λ
2
+
A
1
λ
+
A
2
=
0
,
where
(21)
A
0
=
-
(
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
)
A
1
=
c
1
(
m
C
+
γ
)
E
1
*
+
(
β
1
-
ρ
)
α
β
Z
1
*
P
1
*
(
α
+
P
1
*
)
3
A
2
=
-
(
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
)
c
1
(
m
C
+
γ
)
E
1
*
.
According to Routh-Hurwitz criteria, the necessary and sufficient condition for local stability of equilibrium pointS
3 is
(22)
A
0
>
0
,
A
2
>
0
,
A
0
A
1
-
A
2
>
0
;
that is,
(23)
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
<
0
⟹
k
<
α
+
2
P
1
*
.
In fact, according to P
1
*
=
(
1
/
2
)
[
k
-
α
+
Δ
], k
<
α
+
2
P
1
* will always be satisfied. Hence, if τ
<
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1, then the positive equilibrium S
3 of system (6) is asymptotically stable provided that k
>
k
t
2.
The proof is completed.Remark 2.
(
1
) By Theorem 1
(
1
), the populations-extinction equilibrium S
0 is always unstable; that is to say, phytoplankton and zooplankton populations are not likely to be naturally extinct.
(
2
) Environmental carrying capacity can influence the stability of systems. With the increase of environmental carrying capacity, coexistence may occur.


### 2.2. Global Stability

Now, we study the global behaviors of system (6). Firstly, we analyze the global stability of S
1.Theorem 3.
The equilibriumS
1 is globally asymptotically stable if condition k
≤
k
0
=
d
α
/
(
β
1
-
ρ
) holds.Proof.
Consider(24)
V
(
P
,
Z
,
E
)
=
P
-
k
-
k
ln
⁡
⁡
P
k
+
β
β
1
-
ρ
Z
+
β
m
(
p
-
τ
)
(
β
1
-
ρ
)
E
.
Its derivative along the solution of (1.7) is
(25)
V
˙
=
P
-
k
P
[
r
P
(
1
-
P
k
)
-
β
P
Z
α
+
P
]
+
β
β
1
-
ρ
[
β
1
P
Z
α
+
P
-
d
Z
-
ρ
P
Z
α
+
P
-
c
1
E
Z
]
+
β
m
(
p
-
τ
)
(
β
1
-
ρ
)
{
m
[
c
1
(
p
-
τ
)
Z
-
C
]
-
γ
}
E
=
-
r
(
P
-
k
)
2
k
-
β
P
Z
α
+
P
+
β
Z
α
+
P
k
+
β
P
Z
α
+
P
-
β
Z
β
1
-
ρ
d
-
β
c
1
β
1
-
ρ
E
Z
+
β
c
1
β
1
-
ρ
E
Z
-
β
(
m
C
+
γ
)
E
m
(
p
-
τ
)
(
β
1
-
ρ
)
=
-
r
(
P
-
k
)
2
k
+
β
P
Z
α
+
P
(
k
-
d
α
β
1
-
ρ
)
-
β
d
P
Z
(
β
1
-
ρ
)
(
α
+
P
)
-
β
(
m
C
+
γ
)
E
m
(
p
-
τ
)
(
β
1
-
ρ
)
;
d
V
/
d
t
≤
0 if k
≤
d
α
/
(
β
1
-
ρ
). At the same condition, the equilibrium S
1 is locally asymptotically stable.
By Theorem5.3
.
1 in [17], solution is limited to Φ, the largest invariant subset of {
V
′
(
t
)
=
0
}. Clearly, we see that V
′
(
t
)
=
0 if and only if P
=
k, Z
=
0, and E
=
0. Noting that Φ is invariant, for each element in Φ, we have P
(
t
)
=
k, Z
(
t
)
=
0, and E
(
t
)
=
0. Further, the Lasalle invariance principle implies that all solutions ultimately approach the equilibrium S
1. The proof is completed.Remark 4.
(
1
) By Theorem 3, we know that the zooplankton extinction equilibrium is globally asymptotically stable provided that environmental carrying capacity k
<
k
0. In biological terms, one knows that if k
<
k
0, in some bodies of fresh water, phytoplankton populations are overgrowth and can reach environmental carrying capacity and release toxins so that zooplankton can not survive and lead to extinction; populations adjust and settle down to a new equilibrium state.
(
2
) This phenomenon is called the Harmful Algal Blooms [18] which is a serious and increasing problem in marine waters; to avoid this harm, one can be achieved by controlling nutrient input [19] or biomanipulation [20].In the following, we consider the global stability ofS
3 by constructing a suitable Lyapunov function. The main result is the following.Theorem 5.
Ifτ
<
(
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1
), then the positive equilibrium S
3 of system (6) is globally asymptotically stable provided that k
>
k
t
2.Proof.
We define a Lyapunov function as(26)
V
(
P
,
Z
,
E
)
=
∫
P
1
*
P
ξ
-
P
1
*
ξ
d
ξ
+
d
1
∫
Z
1
*
Z
η
-
Z
1
*
η
d
η
+
d
2
∫
E
1
*
E
Γ
-
E
1
*
Γ
d
Γ
,
where d
1 and d
2 are positive constants to be chosen suitably in the subsequent steps.
It can be easily verified thatV
(
P
,
Z
,
E
) is zero at the equilibrium point and positive for all other positive values of (
P
,
Z
,
E
). Differentiating V with respect to the solutions of (6), a little algebraic manipulation yields(27)
V
˙
=
P
˙
P
(
P
-
P
1
*
)
+
d
1
Z
˙
Z
(
Z
-
Z
1
*
)
+
d
2
E
˙
E
(
E
-
E
1
*
)
=
[
r
(
1
-
P
k
)
-
β
Z
α
+
P
]
(
P
-
P
1
*
)
+
d
1
[
(
β
1
-
ρ
)
P
α
+
P
-
d
-
c
1
E
]
(
Z
-
Z
1
*
)
+
d
2
[
m
c
1
(
p
-
τ
)
Z
-
m
C
]
(
E
-
E
1
*
)
=
(
P
-
P
1
*
)
{
β
(
α
+
P
)
(
α
+
P
1
*
)
-
r
k
(
P
-
P
1
*
)
KKKKKKK
-
β
(
α
+
P
)
(
α
+
P
1
*
)
[
(
Z
-
Z
1
*
)
(
α
+
P
1
*
)
k
k
k
k
k
k
k
k
k
k
k
k
-
Z
*
(
P
-
P
1
*
)
]
-
β
(
α
+
P
)
(
α
+
P
1
*
)
}
+
d
1
[
(
β
1
-
ρ
)
P
α
+
P
-
d
-
c
1
E
]
(
Z
-
Z
1
*
)
+
d
2
[
m
c
1
(
p
-
τ
)
Z
-
m
C
]
(
E
-
E
1
*
)
=
[
-
r
k
+
β
Z
1
*
(
α
+
P
)
(
α
+
P
*
)
]
(
P
-
P
1
*
)
2
-
β
α
+
P
(
P
-
P
1
*
)
(
Z
-
Z
1
*
)
+
d
1
[
(
β
1
-
ρ
)
P
α
+
P
-
(
β
1
-
ρ
)
P
1
*
α
+
P
1
*
]
(
Z
-
Z
1
*
)
-
d
1
c
1
(
Z
-
Z
1
*
)
(
E
-
E
1
*
)
+
d
2
m
c
1
(
p
-
τ
)
(
Z
-
Z
1
*
)
(
E
-
E
1
*
)
.
Choose arbitrary constants d
1 and d
2 as
(28)
d
1
=
β
(
α
+
P
1
*
)
α
(
β
1
-
ρ
)
,
d
2
=
β
(
α
+
P
1
*
)
α
m
(
β
1
-
ρ
)
(
p
-
τ
)
then
(29)
d
V
(
P
,
Z
,
E
)
d
t
=
[
-
r
k
+
β
Z
1
*
(
α
+
P
)
(
α
+
P
1
*
)
]
(
P
-
P
1
*
)
2
.
We only need to consider the sign of-
(
r
/
k
)
+
(
β
Z
1
*
/
(
α
+
P
)
(
α
+
P
1
*
)
). Let M
=
-
(
r
/
k
)
+
(
β
Z
1
*
/
(
α
+
P
)
(
α
+
P
1
*
)
); then M
=
-
(
r
/
k
)
+
(
r
/
(
α
+
P
)
)
(
1
-
(
P
1
*
/
k
)
)
≤
(
r
/
α
k
)
(
k
-
α
-
P
1
*
). Since P
1
* satisfies (9), by the discussing of the last section, we get M
≤
0, and hence, V
′
(
t
)
≤
0.
By Theorem5.3
.
1 in [17], solution is limited to Φ, the largest invariant subset of {
V
′
(
t
)
=
0
}. Clearly, we see that V
′
(
t
)
=
0 if and only if P
=
P
1
*. Noting that Φ is invariant, for each element in Φ, we have P
(
t
)
=
P
1
*. It follows from the first equation of system (6) that 0
=
d
P
/
d
t
=
r
P
1
*
(
1
-
(
P
1
*
/
k
)
)
-
(
β
P
1
*
Z
/
(
α
+
P
1
*
)
), which yields Z
(
t
)
=
Z
1
*. From the third equation of system (6), E
(
t
)
=
E
1
*. Hence, V
′
(
t
)
=
0 if and only if P
=
P
1
*, Z
(
t
)
=
Z
1
*, and E
(
t
)
=
E
1
*. Accordingly, the global asymptotic stability of S
3 follows from the Lasalle invariance principle. The proof is completed.Remark 6.
(
1
) The results given in Theorems 5 can be applied in the context of biological control. If system (6) is exploited, we adjust the tax revenue to control parameters τ
<
(
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1
), k
>
k
t
2, and then interior equilibrium S
3 is globally asymptotically stable.
(
2
) In biological terms, in the case of harvesting zooplankton, by controlling parameter and taking the appropriate value, the phytoplankton and zooplankton populations can coexist and the system will asymptotically approach its equilibrium state.


## 2.1. Local Stability

In order to determine the stability of (6), we compute the Jacobian matrix of system (6):(12)
J
=
(
r
-
2
r
P
k
-
β
α
Z
(
α
+
P
)
2
-
β
P
α
+
P
0
(
β
1
-
ρ
)
α
Z
(
α
+
P
)
2
(
β
1
-
ρ
)
P
α
+
P
-
(
d
+
c
1
E
)
-
c
1
Z
0
m
c
1
(
p
-
τ
)
E
m
[
c
1
(
p
-
τ
)
Z
-
C
]
-
γ
)
.Theorem 1.
For system (6), we have the following: (1) the extinction equilibrium S
0 is unstable;
(2) when k
<
k
t
1, the equilibrium S
1 is locally asymptotically stable;
(3)
(i) if condition (10) is satisfied, then the equilibrium S
2
=
(
P
*
,
Z
*
,
0
) is locally asymptotically stable when k
t
1
<
k
<
k
t
2;
(ii) if condition (10) is not satisfied, then the equilibrium S
2
=
(
P
*
,
Z
*
,
0
) is locally asymptotically stable when 0
<
Z
*
<
(
m
C
+
γ
)
/
m
c
1
(
p
-
τ
), k
<
k
h
=
α
+
2
P
*;
(4) if system (6) has only one positive equilibrium S
3, then S
3 is locally asymptotically stable.Proof.
(
1
) For S
0, the eigenvalues of the matrix (12) are λ
1
=
r
>
0, λ
2
=
-
d
<
0, and λ
3
=
-
m
C
-
γ
<
0. It is seen that there are one unstable manifold and two stable manifolds. Therefore, the point S
0 is a saddle point.
(
2
) For S
1, the characteristic equation is
(13)
|
-
r
-
λ
-
β
k
α
+
k
0
0
(
β
1
-
ρ
)
k
α
+
k
-
d
-
λ
0
0
0
-
m
C
-
γ
-
λ
|
=
0
of which the roots are
(14)
λ
1
=
-
r
<
0
,
λ
2
=
(
β
1
-
ρ
)
k
α
+
k
-
d
,
λ
3
=
-
m
C
-
γ
<
0
.
Hence, the equilibrium S
1 is locally asymptotically stable, when k
<
k
t
1.
(
3
) For S
2, the characteristic equation of S
2 is given by
(15)
{
λ
-
[
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
]
}
×
{
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
λ
[
λ
-
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
]
K
K
K
K
+
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
}
=
0
.
Let λ
1
,
λ
2, and λ
3 be its eigenvalues; assume λ
1
,
λ
2, and λ
3 are all negative; hence,
(16)
λ
1
=
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
<
0
⟹
0
<
Z
*
<
m
C
+
γ
m
c
1
(
p
-
τ
)
.
(i)
Assuming that (10) is satisfied, then 0
<
Z
*
<
(
m
C
+
γ
)
/
m
c
1
(
p
-
τ
) is equivalent to 0
<
k
<
k
t
2. λ
2
,
λ
3 satisfy
(17)
λ
2
-
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
λ
+
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
=
0
.
Therefore(18)
λ
2
+
λ
3
=
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
<
0
kkkkkk
⟹
k
<
k
h
=
α
+
2
P
*
,
k
k
λ
2
λ
3
=
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
>
0
.
Obviously, we know thatk
t
1
<
k
t
2
<
k
h. Therefore,
(19)
λ
1
<
0
,
λ
2
<
0
,
λ
3
<
0
,
0
<
k
<
k
t
2
,
λ
1
=
0
,
λ
2
<
0
,
λ
3
<
0
,
k
=
k
t
2
,
λ
1
>
0
,
λ
2
<
0
,
λ
3
<
0
,
k
t
2
<
k
<
k
h
,
λ
1
>
0
,
λ
2
,
λ
3
=
±
i
ω
0
(
ω
0
>
0
)
,
k
=
k
h
,
λ
1
>
0
,
λ
2
>
0
,
λ
3
>
0
,
k
>
k
h
.
(ii)
Assuming that (10) is not satisfied, then λ
1
<
0, λ
2
<
0, and λ
3
<
0 are equivalent to 0
<
Z
*
<
(
m
C
+
γ
)
/
m
c
1
(
p
-
τ
), k
<
k
h. Therefore, from the paragraph above, combined with the condition of the existence for equilibrium S
2
=
(
P
*
,
Z
*
,
0
), one gets the result.
(
4
) For S
3, the characteristic equation of S
3 is the following:
(20)
λ
3
+
A
0
λ
2
+
A
1
λ
+
A
2
=
0
,
where
(21)
A
0
=
-
(
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
)
A
1
=
c
1
(
m
C
+
γ
)
E
1
*
+
(
β
1
-
ρ
)
α
β
Z
1
*
P
1
*
(
α
+
P
1
*
)
3
A
2
=
-
(
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
)
c
1
(
m
C
+
γ
)
E
1
*
.
According to Routh-Hurwitz criteria, the necessary and sufficient condition for local stability of equilibrium pointS
3 is
(22)
A
0
>
0
,
A
2
>
0
,
A
0
A
1
-
A
2
>
0
;
that is,
(23)
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
<
0
⟹
k
<
α
+
2
P
1
*
.
In fact, according to P
1
*
=
(
1
/
2
)
[
k
-
α
+
Δ
], k
<
α
+
2
P
1
* will always be satisfied. Hence, if τ
<
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1, then the positive equilibrium S
3 of system (6) is asymptotically stable provided that k
>
k
t
2.
The proof is completed.Remark 2.
(
1
) By Theorem 1
(
1
), the populations-extinction equilibrium S
0 is always unstable; that is to say, phytoplankton and zooplankton populations are not likely to be naturally extinct.
(
2
) Environmental carrying capacity can influence the stability of systems. With the increase of environmental carrying capacity, coexistence may occur.


## 2.2. Global Stability

Now, we study the global behaviors of system (6). Firstly, we analyze the global stability of S
1.Theorem 3.
The equilibriumS
1 is globally asymptotically stable if condition k
≤
k
0
=
d
α
/
(
β
1
-
ρ
) holds.Proof.
Consider(24)
V
(
P
,
Z
,
E
)
=
P
-
k
-
k
ln
⁡
⁡
P
k
+
β
β
1
-
ρ
Z
+
β
m
(
p
-
τ
)
(
β
1
-
ρ
)
E
.
Its derivative along the solution of (1.7) is
(25)
V
˙
=
P
-
k
P
[
r
P
(
1
-
P
k
)
-
β
P
Z
α
+
P
]
+
β
β
1
-
ρ
[
β
1
P
Z
α
+
P
-
d
Z
-
ρ
P
Z
α
+
P
-
c
1
E
Z
]
+
β
m
(
p
-
τ
)
(
β
1
-
ρ
)
{
m
[
c
1
(
p
-
τ
)
Z
-
C
]
-
γ
}
E
=
-
r
(
P
-
k
)
2
k
-
β
P
Z
α
+
P
+
β
Z
α
+
P
k
+
β
P
Z
α
+
P
-
β
Z
β
1
-
ρ
d
-
β
c
1
β
1
-
ρ
E
Z
+
β
c
1
β
1
-
ρ
E
Z
-
β
(
m
C
+
γ
)
E
m
(
p
-
τ
)
(
β
1
-
ρ
)
=
-
r
(
P
-
k
)
2
k
+
β
P
Z
α
+
P
(
k
-
d
α
β
1
-
ρ
)
-
β
d
P
Z
(
β
1
-
ρ
)
(
α
+
P
)
-
β
(
m
C
+
γ
)
E
m
(
p
-
τ
)
(
β
1
-
ρ
)
;
d
V
/
d
t
≤
0 if k
≤
d
α
/
(
β
1
-
ρ
). At the same condition, the equilibrium S
1 is locally asymptotically stable.
By Theorem5.3
.
1 in [17], solution is limited to Φ, the largest invariant subset of {
V
′
(
t
)
=
0
}. Clearly, we see that V
′
(
t
)
=
0 if and only if P
=
k, Z
=
0, and E
=
0. Noting that Φ is invariant, for each element in Φ, we have P
(
t
)
=
k, Z
(
t
)
=
0, and E
(
t
)
=
0. Further, the Lasalle invariance principle implies that all solutions ultimately approach the equilibrium S
1. The proof is completed.Remark 4.
(
1
) By Theorem 3, we know that the zooplankton extinction equilibrium is globally asymptotically stable provided that environmental carrying capacity k
<
k
0. In biological terms, one knows that if k
<
k
0, in some bodies of fresh water, phytoplankton populations are overgrowth and can reach environmental carrying capacity and release toxins so that zooplankton can not survive and lead to extinction; populations adjust and settle down to a new equilibrium state.
(
2
) This phenomenon is called the Harmful Algal Blooms [18] which is a serious and increasing problem in marine waters; to avoid this harm, one can be achieved by controlling nutrient input [19] or biomanipulation [20].In the following, we consider the global stability ofS
3 by constructing a suitable Lyapunov function. The main result is the following.Theorem 5.
Ifτ
<
(
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1
), then the positive equilibrium S
3 of system (6) is globally asymptotically stable provided that k
>
k
t
2.Proof.
We define a Lyapunov function as(26)
V
(
P
,
Z
,
E
)
=
∫
P
1
*
P
ξ
-
P
1
*
ξ
d
ξ
+
d
1
∫
Z
1
*
Z
η
-
Z
1
*
η
d
η
+
d
2
∫
E
1
*
E
Γ
-
E
1
*
Γ
d
Γ
,
where d
1 and d
2 are positive constants to be chosen suitably in the subsequent steps.
It can be easily verified thatV
(
P
,
Z
,
E
) is zero at the equilibrium point and positive for all other positive values of (
P
,
Z
,
E
). Differentiating V with respect to the solutions of (6), a little algebraic manipulation yields(27)
V
˙
=
P
˙
P
(
P
-
P
1
*
)
+
d
1
Z
˙
Z
(
Z
-
Z
1
*
)
+
d
2
E
˙
E
(
E
-
E
1
*
)
=
[
r
(
1
-
P
k
)
-
β
Z
α
+
P
]
(
P
-
P
1
*
)
+
d
1
[
(
β
1
-
ρ
)
P
α
+
P
-
d
-
c
1
E
]
(
Z
-
Z
1
*
)
+
d
2
[
m
c
1
(
p
-
τ
)
Z
-
m
C
]
(
E
-
E
1
*
)
=
(
P
-
P
1
*
)
{
β
(
α
+
P
)
(
α
+
P
1
*
)
-
r
k
(
P
-
P
1
*
)
KKKKKKK
-
β
(
α
+
P
)
(
α
+
P
1
*
)
[
(
Z
-
Z
1
*
)
(
α
+
P
1
*
)
k
k
k
k
k
k
k
k
k
k
k
k
-
Z
*
(
P
-
P
1
*
)
]
-
β
(
α
+
P
)
(
α
+
P
1
*
)
}
+
d
1
[
(
β
1
-
ρ
)
P
α
+
P
-
d
-
c
1
E
]
(
Z
-
Z
1
*
)
+
d
2
[
m
c
1
(
p
-
τ
)
Z
-
m
C
]
(
E
-
E
1
*
)
=
[
-
r
k
+
β
Z
1
*
(
α
+
P
)
(
α
+
P
*
)
]
(
P
-
P
1
*
)
2
-
β
α
+
P
(
P
-
P
1
*
)
(
Z
-
Z
1
*
)
+
d
1
[
(
β
1
-
ρ
)
P
α
+
P
-
(
β
1
-
ρ
)
P
1
*
α
+
P
1
*
]
(
Z
-
Z
1
*
)
-
d
1
c
1
(
Z
-
Z
1
*
)
(
E
-
E
1
*
)
+
d
2
m
c
1
(
p
-
τ
)
(
Z
-
Z
1
*
)
(
E
-
E
1
*
)
.
Choose arbitrary constants d
1 and d
2 as
(28)
d
1
=
β
(
α
+
P
1
*
)
α
(
β
1
-
ρ
)
,
d
2
=
β
(
α
+
P
1
*
)
α
m
(
β
1
-
ρ
)
(
p
-
τ
)
then
(29)
d
V
(
P
,
Z
,
E
)
d
t
=
[
-
r
k
+
β
Z
1
*
(
α
+
P
)
(
α
+
P
1
*
)
]
(
P
-
P
1
*
)
2
.
We only need to consider the sign of-
(
r
/
k
)
+
(
β
Z
1
*
/
(
α
+
P
)
(
α
+
P
1
*
)
). Let M
=
-
(
r
/
k
)
+
(
β
Z
1
*
/
(
α
+
P
)
(
α
+
P
1
*
)
); then M
=
-
(
r
/
k
)
+
(
r
/
(
α
+
P
)
)
(
1
-
(
P
1
*
/
k
)
)
≤
(
r
/
α
k
)
(
k
-
α
-
P
1
*
). Since P
1
* satisfies (9), by the discussing of the last section, we get M
≤
0, and hence, V
′
(
t
)
≤
0.
By Theorem5.3
.
1 in [17], solution is limited to Φ, the largest invariant subset of {
V
′
(
t
)
=
0
}. Clearly, we see that V
′
(
t
)
=
0 if and only if P
=
P
1
*. Noting that Φ is invariant, for each element in Φ, we have P
(
t
)
=
P
1
*. It follows from the first equation of system (6) that 0
=
d
P
/
d
t
=
r
P
1
*
(
1
-
(
P
1
*
/
k
)
)
-
(
β
P
1
*
Z
/
(
α
+
P
1
*
)
), which yields Z
(
t
)
=
Z
1
*. From the third equation of system (6), E
(
t
)
=
E
1
*. Hence, V
′
(
t
)
=
0 if and only if P
=
P
1
*, Z
(
t
)
=
Z
1
*, and E
(
t
)
=
E
1
*. Accordingly, the global asymptotic stability of S
3 follows from the Lasalle invariance principle. The proof is completed.Remark 6.
(
1
) The results given in Theorems 5 can be applied in the context of biological control. If system (6) is exploited, we adjust the tax revenue to control parameters τ
<
(
(
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
)
/
α
m
r
c
1
), k
>
k
t
2, and then interior equilibrium S
3 is globally asymptotically stable.
(
2
) In biological terms, in the case of harvesting zooplankton, by controlling parameter and taking the appropriate value, the phytoplankton and zooplankton populations can coexist and the system will asymptotically approach its equilibrium state.


## 3. Bifurcation Phenomenon

Theorem 7.
Ifk
=
k
t
1, then the system (6) exhibits a transcritical bifurcation about S
2 which is branched out from S
1 at k
=
k
t
1.Proof.
By analysing Section2, one knows that the system (6) only has two feasible equilibria S
0
,
S
1, when k
<
k
t
1. In fact, the system (6) has another equilibrium S
2
=
(
P
*
,
Z
*
,
0
); in this case, we know that Z
*
<
0. Based on the biological significance, zooplankton populations must be positive, so we abandon it. Similar to the proof of Theorem 1, we easily know that when k
<
k
t
1, equilibrium S
2 is a saddle point and equilibrium S
1 is stable. Equilibria S
1 and S
2 coalesce into equilibrium S
1 which becomes a nonhyperbolic equilibrium at k
=
k
t
1; when k
>
k
t
1, the system (6) appears as a boundary equilibrium S
2 which is stable by Theorem 1, and equilibrium S
1 becomes a saddle point by (13). Thus, an exchange of stability has occurred at k
=
k
t
1. This type of bifurcation is called a transcritical bifurcation.Theorem 8.
Ifk
=
k
t
2, then the system (6) exhibits a transcritical bifurcation about S
3 which is branched out from S
2 at k
=
k
t
2.Proof.
Obviously,
k
t
1
<
k
t
2
<
k
h. When k
t
1
<
k
<
k
t
2, the system (6) has equilibria S
0
,
S
1
,
S
2, and  S
3; in this case, E
1
*
<
0. Based on the biological significance, we abandon equilibrium S
3. For equilibrium S
3, it is a saddle point as A
2
<
0 in (20), and equilibrium S
2 is stable by (19). Equilibria S
2 and S
3 coalesce into equilibrium S
2 which becomes a nonhyperbolic equilibrium at k
=
k
t
2. When k
>
k
t
2, the system (6) appears as a positive interior equilibrium S
3 which is stable by Theorem 1, and equilibrium S
2 becomes a saddle point by (19). Thus, an exchange of stability has occurred at k
=
k
t
2. This type of bifurcation is called a transcritical bifurcation.When condition (10) is satisfied, accordingly by (19), the characteristic equation of S
2 has a pair of pure imaginary roots at k
=
k
h, and the system (6) exhibits Hopf bifurcation and can be stated as follows.Theorem 9.
When condition (10) is satisfied, the model system (6) exhibits a Hopf bifurcation around S
2
=
(
P
*
,
Z
*
,
0
) and it is unstable at k
=
k
h.Proof.
The eigenvalues ofJ
(
S
2
) can be expressed as the solutions of the characteristic equation
(30)
λ
3
+
a
0
λ
2
+
a
1
λ
+
a
2
=
0
,
where
(31)
a
0
=
-
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
-
[
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
]
,
a
1
=
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
[
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
]
+
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
,
a
2
=
[
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
]
(
β
1
-
ρ
)
α
β
(
α
+
P
*
)
3
Z
*
P
*
.
The theorem will be proved if we show that the conditions for Hopf bifurcation are satisfied. By the discussion of (19), the characteristic equation of S
2 has a pair of pure imaginary roots at k
=
k
h; obviously, we only need to verify the transversality conditions
(32)
[
d
d
k
(
Re
⁡
λ
)
]
k
=
k
h
≠
0
,
where λ
(
k
)
=
μ
(
k
)
±
i
ω
(
k
), with μ
(
k
h
)
=
0, ω
(
k
h
)
=
ω
0.
In this case, the eigenvalues areλ
1
=
[
m
c
1
(
p
-
τ
)
Z
*
-
m
C
-
γ
] and λ
2,3
=
λ
(
k
)
=
μ
(
k
)
±
i
ω
(
k
).
Hence,(33)
λ
3
+
a
0
λ
2
+
a
1
λ
+
a
2
=
(
λ
-
λ
1
)
(
λ
2
-
2
μ
(
k
)
λ
+
μ
(
k
)
2
+
ω
(
k
2
)
)
.
Comparative coefficient
(34)
μ
(
k
)
=
1
2
(
r
-
2
r
P
*
k
-
β
α
Z
*
(
α
+
P
*
)
2
)
.
Therefore,
(35)
[
d
d
k
(
Re
⁡
λ
)
]
k
=
k
h
=
[
d
d
k
(
μ
(
k
)
)
]
k
h
=
r
α
P
*
+
2
r
(
P
*
)
2
2
(
k
*
)
2
(
α
+
P
*
)
≠
0
.
Hence, all the conditions for Hopf bifurcation are satisfied. By (19), we know λ
1
>
0 at k
=
k
h; therefore, Hopf bifurcation of S
2 is unstable.Remark 10.
By Theorems1 and 9, when environmental carrying capacity k
=
k
h, the system is unstable and exhibits a Hopf bifurcation around S
2; when condition (10) is satisfied, the Hopf bifurcation is unstable; when condition (10) is not satisfied, the Hopf bifurcation may be stable. Because we care more about dynamical behavior of positive interior equilibrium S
3, we do not specifically analyze the stability of Hopf bifurcation and only give some numerical simulation in Section 5.In order to understand the problem more comfortably, we draw a graph in Figure2.Figure 2
(
1
) Equilibria S
0 and S
1 always exist, where S
0 is unstable. (
2
) When 0
<
k
<
k
t
1, S
1 is stable; when 0
<
k
<
k
0, S
1 is globally asymptotically stable; when k
>
k
t
1, S
1 is unstable. (
3
) When k
=
k
t
1, equilibrium point S
2 is branched out from S
1; k
=
k
t
1 is a transcritical bifurcation value. (
4
) When k
1
<
k
<
k
t
2, equilibrium point S
2 is asymptotically stable. (
5
) When k
=
k
t
2, equilibrium point S
3 is branched out from S
2; k
=
k
t
2 is a transcritical bifurcation value. (
6
) When k
>
k
t
2, S
3 is globally asymptotically stable. (
7
)
k
=
k
h is a Hopf bifurcation value of S
2.


## 4. Optimal Taxation Policy

Our focus so far has been on the dynamic behaviors of the system (6). Biologically, in the presence of harvesting, in order to maintain the survival of both species, particularly we care more about the positive interior equilibrium. The objective is to maximize the monetary social benefit as well as conservation of the ecosystem.The fisherman and regulatory agency are actually two different components of the society at large. Hence, the revenues earned by them are the revenues accrued to the society through the fishery. The net economic revenue to the society is(36)
(
c
1
p
Z
-
C
)
E
=
[
c
1
(
p
-
τ
)
Z
-
C
]
E
+
τ
c
1
Z
E
which equals the net economic revenue to the fisherman (perceived rent) plus the economic revenue to the regulatory agency. Note that
(37)
π
(
P
,
Z
,
E
)
=
(
c
1
p
Z
-
C
)
E
.For optimal harvest policy, this objective amounts to maximizing the present valueJ of a continuous time steam of revenues given by
(38)
J
=
∫
0
∞
e
-
δ
t
(
c
1
p
Z
-
C
)
E
d
t
,
where δ denotes the instantaneous annual rate of discount [21].Our objective is to determine a tax policyτ
=
τ
(
t
) to maximize J subject to the state equation of (6) by invoking Pontryagin’s maximum principle [22]; the control variable τ
(
t
) is subjected to the constraints τ
min
⁡
≤
τ
≤
τ
max
⁡. The case in which τ
min
⁡
<
0 allows us to consider subsidies, which in this case would have the effect of increasing the rate of expansion of the fishery.The Hamiltonian of this control problem is(39)
H
=
e
-
δ
t
(
c
1
p
Z
-
C
)
E
+
λ
1
d
P
d
t
+
λ
2
d
Z
d
t
+
λ
3
d
E
d
t
,
where λ
1
,
λ
2, and λ
3 are the adjoint variables.Hamiltonian (39) must be maximized for τ
∈
[
τ
min
⁡
,
τ
max
⁡
]. Assuming that the constraints are not binding (i.e., the optimal solution does not occur at τ
=
τ
min
⁡ or τ
max
⁡), we have singular control given by
(40)
∂
H
∂
τ
=
-
λ
3
m
c
1
Z
E
=
0
⟹
λ
3
=
0
.
According to Pontryagin’s maximum principle, the adjoint equations are
(41)
d
λ
1
d
t
=
-
∂
H
∂
P
,
d
λ
2
d
t
=
-
∂
H
∂
Z
,
d
λ
3
d
t
=
-
∂
H
∂
E
=
0
.
Substitution and simplification yield
(42)
d
λ
1
d
t
=
-
∂
H
∂
P
=
-
λ
1
[
r
-
2
r
P
k
-
β
α
Z
(
α
+
P
)
2
]
kkkkkkl
-
λ
2
(
β
1
-
ρ
)
α
Z
(
α
+
P
)
2
,
d
λ
2
d
t
=
-
∂
H
∂
Z
=
-
c
1
p
e
-
δ
t
E
+
λ
1
β
P
α
+
P
kkkkkkl
-
λ
2
[
(
β
1
-
ρ
)
P
α
+
P
-
(
d
+
c
1
E
)
]
,
(43)
d
λ
3
d
t
=
-
∂
H
∂
E
=
-
(
c
1
p
Z
-
C
)
e
-
δ
t
+
λ
2
c
1
Z
=
0
.To obtain an optimal equilibrium solution, the solution of (43) is discussed in the following by considering the interior equilibrium as
(44)
λ
2
=
e
-
δ
t
(
p
-
C
c
1
Z
1
*
)
.
Let
(45)
A
1
=
r
-
2
r
P
1
*
k
-
β
α
Z
1
*
(
α
+
P
1
*
)
2
A
2
=
(
β
1
-
ρ
)
α
Z
1
*
(
α
+
P
1
*
)
2
(
p
-
C
c
1
Z
1
*
)
A
3
=
-
c
1
p
E
1
*
+
A
2
A
1
+
δ
β
P
1
*
α
+
P
1
*
.
We will rewrite (42) by considering the interior equilibrium as
(46)
d
λ
1
d
t
=
-
A
1
λ
1
-
A
2
e
-
δ
t
,
d
λ
2
d
t
=
A
3
e
-
δ
t
.
Solutions of the above linear differential equation are
(47)
λ
1
(
t
)
=
A
2
A
1
+
δ
e
-
δ
t
,
λ
2
(
t
)
=
A
3
δ
e
-
δ
t
.
Substituting the value of λ
2
(
t
) from (44) into (47), we get
(48)
(
p
-
C
c
1
Z
1
*
)
=
A
3
δ
.Now using the value ofP
1
*
,
Z
1
*, and E
1
* from Section 2 into (48), we get an equation for τ; let τ
* be a solution of this equation. Using this value of τ
=
τ
*, we get the optimal equilibrium solutions P
=
P
(
τ
*
), Z
=
Z
(
τ
*
), and E
=
E
(
τ
*
). Thus, we have established the existence of an optimal equilibrium solution that satisfies the necessary conditions of the maximum principle.From the above analysis carried out in this section, we observe the following.(1)
From (43), we get
(49)
λ
2
c
1
Z
=
(
c
1
p
Z
-
C
)
e
-
δ
t
=
e
-
δ
t
∂
π
∂
E
.
Putting the value ofλ
2
(
t
) into (48), we get
(50)
c
1
p
Z
-
A
3
δ
c
1
Z
=
C
.
Whenδ
→
∞, (50) leads to the result c
1
p
Z
=
C which implies that the economic rent is completely dissipated.
(2)
By (48), we get the optimal equilibrium populations P
=
P
(
τ
*
), Z
=
Z
(
τ
*
), and E
=
E
(
τ
*
); hence, we have
(51)
π
=
(
c
1
p
Z
-
C
)
E
=
A
3
δ
c
1
Z
E
.
Thus,π is a decreasing function of δ; we, therefore, conclude that π leads to maximization when δ leads to 0.


## 5. Numerical Simulation

We have considered dynamic behaviors and optimal taxation policy of the system (6). To facilitate the interpretation of our mathematical results in model (6), we perform some numerical simulations.(1)
Letr
=
5, k
=
6, α
=
9, β
=
0.6, β
1
=
0.43, ρ
=
0.2, c
1
=
0.64, m
=
0.53, p
=
0.7, γ
=
0.1, d
=
0.2, τ
=
0.01, and C
=
0.3. It is easy to verify that 6
=
k
<
d
α
/
(
β
1
-
ρ
)
=
7.826
<
d
α
/
(
β
1
-
ρ
-
d
); then the condition of Theorem 5 is satisfied. Hence, the equilibrium S
1
(
6,0
,
0
) is globally asymptotically stable, which is shown in Figure 3.
(2)
We choose a set of parameters as follows:r
=
1, k
=
5.5, α
=
4, β
=
1
,
β
1
=
0.7, ρ
=
0.1, c
1
=
0.1, m
=
0.5, p
=
0.7, γ
=
0.1, d
=
0.1, τ
=
0.2, and C
=
0.1. It easy to compute that
(52)
P
*
=
0.8
,
Z
*
=
4.10
,
Z
*
<
m
C
+
γ
m
c
1
(
p
-
τ
)
=
6
,
P
*
<
k
<
α
+
2
P
*
=
5.6
.So condition (
3
)
(
ii
) of Theorem 1 holds. The equilibrium S
2
(
0.8,4.10,0
) is locally asymptotically stable, which is shown in Figure 4(a).(3)
Choosek
=
7 and other parameters are the same as case (
2
); then
(53)
P
*
=
0.8
,
Z
*
=
4.25
,
Z
*
<
m
C
+
γ
m
c
1
(
p
-
τ
)
=
6
,
k
>
α
+
2
P
*
=
5.6
.Therefore, the equilibrium S
2
(
0.8,4.25,0
) is unstable and exhibits a Hopf bifurcation. From the numerical simulation, we note that the system (6) has a cycle in Figure 4(b).(4)
Consider the following choice of parametric values:r
=
16, k
=
1.6, α
=
1, β
=
1, β
1
=
0.9, ρ
=
0.1, c
1
=
0.1, m
=
0.5, p
=
0.7, γ
=
0.01, d
=
0.1, τ
=
0.2, andC
=
0.1. It is easy to verify that
(54)
k
=
1.6
>
k
t
2
=
0.1615
,
τ
=
0.2
<
min
⁡
{
p
,
α
m
r
c
1
p
-
β
(
m
C
+
γ
)
α
m
r
c
1
}
=
0.625
.So, the conditions of Theorem 5 are satisfied; then interior equilibrium S
3
=
(
1.5,2.4,3.8
) is globally asymptotically stable, which is shown in Figure 5.Figure 3
The stability of the equilibriumS
1.The stability of the equilibriumS
2. (a) The equilibrium S
2
(
0.8,4.10,0
) is locally asymptotically stable when k
=
5.5
<
5.6. (b) The equilibrium S
2
(
0.8,4.25,0
) exhibits a Hopf bifurcation when k
=
7
>
5.6.
(a)
(b)The globally asymptotical stability of positive equilibrium point. (a) Time series portrait. (b) The phase portrait with different initial values.
(a)
(b)From the numerical examples discussed above, we may note the following points: the equilibriumS
2 exists in absence of harvesting effort but at a lower population level (
P
*
,
Z
*
,
0
)
=
(
0.8,4.10,0
) compared to equilibrium S
3
=
(
P
1
*
,
Z
1
*
,
E
1
*
)
=
(
1.5,2.4,3.8
) (in presence of harvesting effort) for the phytoplankton and at a higher population level for the zooplankton. This agrees with the facts.(5)
Letr
=
4, k
=
1.6, α
=
1, β
=
1, β
1
=
0.9, ρ
=
0.1, c
1
=
0.1, m
=
0.5, p
=
0.7, γ
=
0.01, d
=
0.1, andC
=
0.1, in appropriate units. According to the dissuasion of Section 3, then, for the above values of the parameter, optimal tax becomes τ
*
=
0.39, and corresponding stable optimal equilibrium is (0.676, 3.871, 2.227). The three-dimensional phase space trajectories corresponding to the optimal tax τ
*
=
0.39, beginning with different initial levels, are depicted in Figure 6. For this optimal tax τ
*
=
0.39, trajectories clearly indicate that the optimal equilibrium S
*
(
0.676,3.871,2.227
) is found to be stable because condition of Theorem 1 holds in this case.Figure 6
The optimal equilibriumS
*
(
0.676,3.871,2.227
) is stable with different initial levels.In Figures7 and 8, variations of phytoplankton, zooplankton, and harvesting effort against time are plotted for different tax levels. From these plots, we observe that as the rate of tax increases zooplankton populations increase, while phytoplankton populations and harvesting effort decrease as expected.(a) Variation of phytoplankton population with time for different tax levels. (b) Variation of zooplankton population with time for different tax levels.
(a)
(b)Figure 8
Variation of harvesting effort with time for different tax levels.


## 6. Conclusion

In this work, we have elaborated a phytoplankton-zooplankton system model, in which zooplankton are assumed to undergo commercial exploitation. The model is realistic because we force the fishing effort to remain under control by imposing a tax to keep the ecological balance. The most important feature of the present model is that it assumes a fully dynamic interaction between fishing effort and the net economic revenue to the fisherman in the case of a phytoplankton-zooplankton fishery.First, stability criteria of the model are analyzed both from local and global points of view. Deserving to be mentioned, we investigate the global stability of the equilibria and give the corresponding parameter regions. The consequence of global stability shows that exploitation will not irreversibly change the system, as long as the tax keeps a threshold value; that is, the zooplankton are not excessive exploitation, and the system is able to recover. Furthermore, existence of Hopf bifurcation around the equilibriumS
2 has been established with carrying capacity as the bifurcation parameter. That is, when carrying capacity k crosses a threshold value k
h, the system undergoes a Hopf bifurcation. We also get that the system exhibits transcritical bifurcations.In addition, the optimal harvesting policy for harvesting zooplankton is studied by imposing a tax. The monetary social benefit is maximized by using Pontryagin’s maximum principle. We discuss the case of the optimal equilibrium solution and have established the optimal equilibrium solution by using taxτ
=
τ
*. It is established that the zero discounting leads to the maximization of economic revenue and that an infinite discount rate leads to complete dissipation of economic rent.



---

*Source: 684790-2014-05-29.xml*