Stability and Selective Harvesting of a Phytoplankton-Zooplankton System

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.

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 Figure 1). In aquatic food webs, some zooplankton such as jellyfish, krill, and Acetes 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−𝑘−𝛽𝑃𝑍,𝛼+𝑃𝑑𝑍=𝛽𝑑𝑡¹𝑃𝑍𝛼+𝑃−𝑑𝑍−𝜌𝑃𝑍.𝛼+𝑃 (1) The following assumptions for model (1) are made: (1)the variables 𝑃 and 𝑍 are the density of phytoplankton population and the density of zooplankton population at any instant of time 𝑡, respectively;(2)the parameter 𝑟 is the intrinsic growth rate and 𝑘 is the environmental carrying capacity of population. The constant 𝛽(>0) is the maximum uptake rate for zooplankton species; 𝛽₁(>0) denotes the ratio of biomass conversion (satisfying the obvious restriction 0<𝛽₁<𝛽) and 𝑑(>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 𝛽₁<𝜌; that is, the ratio of biomass consumed by zooplankton is greater than the rate of toxic substance liberation by phytoplankton species;(4)the term (𝛽𝑃𝑍/(𝛼+𝑃)) represents the functional response for the grazing of phytoplankton by zooplankton and 𝛼 is the half saturation constant for a Holling type II functional response. (𝜌𝑃𝑍/(𝛼+𝑃)) 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 𝑑𝑃𝑃𝑑𝑡=𝑟𝑃1−𝑘−𝛽𝑃𝑍,𝛼+𝑃𝑑𝑍=𝛽𝑑𝑡¹𝑃𝑍𝛼+𝑃−𝑑𝑍−𝜌𝑃𝑍𝛼+𝑃−𝑐₁𝐸𝑍, (2) where the constant 𝑐₁ is the catchability coefficient of the zooplankton. 𝐸 is harvesting effort.

Here, we take 𝐸 as dynamic [9, 12] (i.e., time-dependent) variable governed by the equations 𝐸(𝑡)=𝜇₁𝑄(𝑡),0≤𝜇₁≤1,𝑑𝑄(𝑡)(𝑡)(0)𝑑𝑡=𝐼−𝛾𝑄,𝑄=𝑄₀, (3) where 𝐼(𝑡) is the gross investment rate at time 𝑡. 𝑄(𝑡) is the amount of capital invested in the fishery at time 𝑡; 𝛾 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 {𝑐₁(𝑝−𝜏)𝑍−𝐶}𝐸, where 𝑝 is the constant price per unit biomass of the zooplankton species and 𝐶 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 𝐼(𝑡)=𝜇₂𝑐₁𝑝−𝜏𝑍−𝐶𝐸(𝑡),0≤𝜇₂≤1. (4) Equation (4) asserts that the maximum investment rate at any time equals the net economic revenue (for 𝜇₂=1) at that time. By virtue of (4) and (3) yield the result 𝑑𝐸=𝜇𝑑𝑡¹𝜇₂𝑐₁𝑝−𝜏𝑍−𝐶−𝛾𝐸. (5) Let 𝑚=𝜇₁𝜇₂; therefore, we have the following system of equations: 𝑑𝑃𝑃𝑑𝑡=𝑟𝑃1−𝑘−𝛽𝑃𝑍,𝛼+𝑃𝑑𝑍=𝛽𝑑𝑡¹𝑃𝑍𝛼+𝑃−𝑑𝑍−𝜌𝑃𝑍𝛼+𝑃−𝑐₁𝐸𝑍,𝑑𝐸=𝑚𝑐𝑑𝑡¹𝑃(0)𝑝−𝜏𝑍−𝐶−𝛾𝐸,=𝑃₀>0,𝑍(0)=𝑍₀𝐸(0)>0,=𝐸₀>0. (6)

In this paper, we choose tax as the management objective when discussing the impact of harvesting in the phytoplankton-zooplankton system and assume (𝐇)𝛽₁−𝜌−𝑑>0, 𝑝−𝜏>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 points 𝑆₀=(0,0,0) and 𝑆₁=(𝑘,0,0) exist for all permissible parameters.(2)If 𝑘>𝑘_𝑡1=(𝑑𝛼/(𝛽₁−𝜌−𝑑)), the boundary equilibrium 𝑆₂=(𝑃^∗,𝑍^∗,0) exists, where 𝑃^∗=𝑑𝛼𝛽₁−𝜌−𝑑,𝑍^∗=𝑟𝛽𝛼+𝑃^∗𝑃1−^∗𝑘. (7)(3)Let 𝑆₃=(𝑃^∗₁,𝑍^∗₁,𝐸^∗₁) be the positive interior equilibrium, where 𝑍^∗₁=𝑚𝐶+𝛾𝑚𝑐₁𝑝−𝜏,𝐸^∗₁=𝛽₁−𝜌𝑃^∗₁𝑐₁𝛼+𝑃^∗₁−𝑑𝑐₁ (8)and 𝑃^∗₁ satisfies 𝑃^∗₁²+(𝛼−𝑘)𝑃^∗₁−𝑘𝛼+𝑘𝛽𝑚𝐶+𝛾𝑟𝑚𝑐₁𝑝−𝜏=0. (9) Let 𝑥₁,𝑥₂ be the roots of (9); we only consider that 𝑥₁,𝑥₂ have only one positive root; then 𝑥₁𝑥₂=−𝑘𝛼+𝑘𝛽𝑚𝐶+𝛾𝑟𝑚𝑐₁𝑝−𝜏<0⟹𝜏<𝛼𝑚𝑟𝑐₁𝑝−𝛽𝑚𝐶+𝛾𝛼𝑚𝑟𝑐₁,Δ=(𝛼−𝑘)²+4𝑘𝛼−𝑘𝛽𝑚𝐶+𝛾𝑟𝑚𝑐₁𝑝−𝜏>0, (10) and hence, 𝑃^∗₁ exists as a positive root: 𝑃^∗₁=12𝑘−𝛼+√Δ. (11) We know that if 𝑘>𝑘_𝑡2,  𝑍^∗₁>0 and 𝐸^∗₁>0, where 𝑘_𝑡2=((𝑑²𝛼²/(𝛽₁−𝜌−𝑑)²)+(𝑑𝛼²/(𝛽₁−𝜌−𝑑)))/(𝛼+(𝑑𝛼/(𝛽₁−𝜌−𝑑))−(𝛽(𝑚𝐶+𝛾)/𝑟𝑚𝑐₁(𝑝−𝜏))).

Therefore, when 𝑘>𝑘_𝑡2 together with 𝜏<(𝛼𝑚𝑟𝑐₁𝑝−𝛽(𝑚𝐶+𝛾))/𝛼𝑚𝑟𝑐₁, the positive interior equilibrium 𝑆₃ exists.

2.1. Local Stability

In order to determine the stability of (6), we compute the Jacobian matrix of system (6):𝐽=𝑟−2𝑟𝑃𝑘−𝛽𝛼𝑍(𝛼+𝑃)²−𝛽𝑃0𝛽𝛼+𝑃₁−𝜌𝛼𝑍(𝛼+𝑃)²𝛽₁−𝜌𝑃−𝛼+𝑃𝑑+𝑐₁𝐸−𝑐₁𝑍0𝑚𝑐₁𝑝−𝜏𝐸𝑚𝑐₁.𝑝−𝜏𝑍−𝐶−𝛾 (12)

Theorem 1. For system (6), we have the following:  (1) the extinction equilibrium 𝑆₀ is unstable; (2) when 𝑘<𝑘_𝑡1, the equilibrium 𝑆₁ is locally asymptotically stable; (3)  (i) if condition (10) is satisfied, then the equilibrium 𝑆₂=(𝑃^∗,𝑍^∗,0) is locally asymptotically stable when 𝑘_𝑡1<𝑘<𝑘_𝑡2; (ii) if condition (10) is not satisfied, then the equilibrium 𝑆₂=(𝑃^∗,𝑍^∗,0) is locally asymptotically stable when 0<𝑍^∗<(𝑚𝐶+𝛾)/𝑚𝑐₁(𝑝−𝜏), 𝑘<𝑘_ℎ=𝛼+2𝑃^∗; (4) if system (6) has only one positive equilibrium 𝑆₃, then 𝑆₃ is locally asymptotically stable.

Proof. (1) For 𝑆₀, the eigenvalues of the matrix (12) are 𝜆₁=𝑟>0, 𝜆₂=−𝑑<0, and 𝜆₃=−𝑚𝐶−𝛾<0. It is seen that there are one unstable manifold and two stable manifolds. Therefore, the point 𝑆₀ is a saddle point.
(2) For 𝑆₁, the characteristic equation is −𝑟−𝜆−𝛽𝑘00𝛽𝛼+𝑘₁−𝜌𝑘𝛼+𝑘−𝑑−𝜆000−𝑚𝐶−𝛾−𝜆=0 (13) of which the roots are 𝜆₁=−𝑟<0,𝜆₂=𝛽₁−𝜌𝑘𝜆𝛼+𝑘−𝑑,₃=−𝑚𝐶−𝛾<0. (14) Hence, the equilibrium 𝑆₁ is locally asymptotically stable, when 𝑘<𝑘_𝑡1.
(3) For 𝑆₂, the characteristic equation of 𝑆₂ is given by 𝜆−𝑚𝑐₁𝑝−𝜏𝑍^∗×𝜆−𝑚𝐶−𝛾𝜆−𝑟−2𝑟𝑃^∗𝑘−𝛽𝛼𝑍^∗(𝛼+𝑃^∗)²+𝛽₁−𝜌𝛼𝛽(𝛼+𝑃^∗)³𝑍^∗𝑃^∗=0. (15) Let 𝜆₁,𝜆₂, and 𝜆₃ be its eigenvalues; assume 𝜆₁,𝜆₂, and 𝜆₃ are all negative; hence, 𝜆₁=𝑚𝑐₁𝑝−𝜏𝑍^∗−𝑚𝐶−𝛾<0⟹0<𝑍^∗<𝑚𝐶+𝛾𝑚𝑐₁.𝑝−𝜏 (16)(i)Assuming that (10) is satisfied, then 0<𝑍^∗<(𝑚𝐶+𝛾)/𝑚𝑐₁(𝑝−𝜏) is equivalent to 0<𝑘<𝑘_𝑡2. 𝜆₂,𝜆₃ satisfy 𝜆²−𝑟−2𝑟𝑃^∗𝑘−𝛽𝛼𝑍^∗(𝛼+𝑃^∗)²𝛽𝜆+¹−𝜌𝛼𝛽(𝛼+𝑃^∗)³𝑍^∗𝑃^∗=0. (17) Therefore 𝜆₂+𝜆₃=𝑟−2𝑟𝑃^∗𝑘−𝛽𝛼𝑍^∗(𝛼+𝑃^∗)²<0⟹𝑘<𝑘_ℎ=𝛼+2𝑃^∗,𝜆₂𝜆₃=𝛽₁−𝜌𝛼𝛽(𝛼+𝑃^∗)³𝑍^∗𝑃^∗>0. (18) Obviously, we know that 𝑘_𝑡1<𝑘_𝑡2<𝑘_ℎ. Therefore, 𝜆₁<0,𝜆₂<0,𝜆₃<0,0<𝑘<𝑘_𝑡2,𝜆₁=0,𝜆₂<0,𝜆₃<0,𝑘=𝑘_𝑡2,𝜆₁>0,𝜆₂<0,𝜆₃<0,𝑘_𝑡2<𝑘<𝑘_ℎ,𝜆₁>0,𝜆₂,𝜆₃=±𝑖𝜔₀𝜔₀>0,𝑘=𝑘_ℎ,𝜆₁>0,𝜆₂>0,𝜆₃>0,𝑘>𝑘_ℎ. (19)(ii)Assuming that (10) is not satisfied, then 𝜆₁<0, 𝜆₂<0, and 𝜆₃<0 are equivalent to 0<𝑍^∗<(𝑚𝐶+𝛾)/𝑚𝑐₁(𝑝−𝜏), 𝑘<𝑘_ℎ. Therefore, from the paragraph above, combined with the condition of the existence for equilibrium 𝑆₂=(𝑃^∗,𝑍^∗,0), one gets the result.
(4) For 𝑆₃, the characteristic equation of 𝑆₃ is the following: 𝜆³+𝐴₀𝜆²+𝐴₁𝜆+𝐴₂=0, (20) where 𝐴₀=−𝑟−2𝑟𝑃^∗₁𝑘−𝛽𝛼𝑍^∗₁𝛼+𝑃^∗₁²𝐴₁=𝑐₁𝑚𝐶+𝛾𝐸^∗₁+𝛽₁−𝜌𝛼𝛽𝑍^∗₁𝑃^∗₁𝛼+𝑃^∗₁³𝐴₂=−𝑟−2𝑟𝑃^∗₁𝑘−𝛽𝛼𝑍^∗₁𝛼+𝑃^∗₁²𝑐₁𝑚𝐶+𝛾𝐸^∗₁. (21)
According to Routh-Hurwitz criteria, the necessary and sufficient condition for local stability of equilibrium point 𝑆₃ is 𝐴₀>0,𝐴₂>0,𝐴₀𝐴₁−𝐴₂>0; (22) that is, 𝑟−2𝑟𝑃^∗₁𝑘−𝛽𝛼𝑍^∗₁𝛼+𝑃^∗₁²<0⟹𝑘<𝛼+2𝑃^∗₁. (23) In fact, according to 𝑃^∗₁=(1/2)[𝑘−𝛼+√Δ], 𝑘<𝛼+2𝑃^∗₁ will always be satisfied. Hence, if 𝜏<(𝛼𝑚𝑟𝑐₁𝑝−𝛽(𝑚𝐶+𝛾))/𝛼𝑚𝑟𝑐₁, then the positive equilibrium 𝑆₃ of system (6) is asymptotically stable provided that 𝑘>𝑘_𝑡2.
The proof is completed.

Remark 2. (1) By Theorem 1(1), the populations-extinction equilibrium 𝑆₀ 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 𝑆₁.

Theorem 3. The equilibrium 𝑆₁ is globally asymptotically stable if condition 𝑘≤𝑘₀=𝑑𝛼/(𝛽₁−𝜌) holds.

Proof. Consider 𝑉()𝑃𝑃,𝑍,𝐸=𝑃−𝑘−𝑘ln𝑘+𝛽𝛽₁𝛽−𝜌𝑍+𝑚𝑝−𝜏𝛽₁−𝜌𝐸. (24) Its derivative along the solution of (1.7) is 𝑉̇=𝑃−𝑘𝑃𝑃𝑟𝑃1−𝑘−𝛽𝑃𝑍+𝛽𝛼+𝑃𝛽₁𝛽−𝜌¹𝑃𝑍𝛼+𝑃−𝑑𝑍−𝜌𝑃𝑍𝛼+𝑃−𝑐₁+𝛽𝐸𝑍𝑚𝑝−𝜏𝛽₁𝑚𝑐−𝜌¹𝐸()𝑝−𝜏𝑍−𝐶−𝛾=−𝑟𝑃−𝑘²𝑘−𝛽𝑃𝑍+𝛼+𝑃𝛽𝑍𝑘+𝛼+𝑃𝛽𝑃𝑍−𝛼+𝑃𝛽𝑍𝛽₁−𝜌𝑑−𝛽𝑐₁𝛽₁−𝜌𝐸𝑍+𝛽𝑐₁𝛽₁−𝛽𝐸−𝜌𝐸𝑍𝑚𝐶+𝛾𝑚𝑝−𝜏𝛽₁()−𝜌=−𝑟𝑃−𝑘²𝑘+𝛽𝑃𝑍𝛼+𝑃𝑘−𝑑𝛼𝛽₁−−𝜌𝛽𝑑𝑃𝑍𝛽₁()−𝛽𝐸−𝜌𝛼+𝑃𝑚𝐶+𝛾𝑚𝑝−𝜏𝛽₁;−𝜌 (25)𝑑𝑉/𝑑𝑡≤0 if 𝑘≤𝑑𝛼/(𝛽₁−𝜌). At the same condition, the equilibrium 𝑆₁ is locally asymptotically stable.
By Theorem 5.3.1 in [17], solution is limited to Φ, the largest invariant subset of {𝑉^′(𝑡)=0}. Clearly, we see that 𝑉^′(𝑡)=0 if and only if 𝑃=𝑘, 𝑍=0, and 𝐸=0. Noting that Φ is invariant, for each element in Φ, we have 𝑃(𝑡)=𝑘, 𝑍(𝑡)=0, and 𝐸(𝑡)=0. Further, the Lasalle invariance principle implies that all solutions ultimately approach the equilibrium 𝑆₁. 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 𝑘<𝑘₀. In biological terms, one knows that if 𝑘<𝑘₀, 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 of 𝑆₃ by constructing a suitable Lyapunov function. The main result is the following.

Theorem 5. If 𝜏<((𝛼𝑚𝑟𝑐₁𝑝−𝛽(𝑚𝐶+𝛾))/𝛼𝑚𝑟𝑐₁), then the positive equilibrium 𝑆₃ of system (6) is globally asymptotically stable provided that 𝑘>𝑘_𝑡2.

Proof. We define a Lyapunov function as 𝑉(𝑃,𝑍,𝐸)=^𝑃_𝑃∗1𝜉−𝑃^∗₁𝜉𝑑𝜉+𝑑₁^𝑍_𝑍∗1𝜂−𝑍^∗₁𝜂𝑑𝜂+𝑑₂^𝐸_𝐸∗1Γ−𝐸^∗₁Γ𝑑Γ, (26) where 𝑑₁ and 𝑑₂ are positive constants to be chosen suitably in the subsequent steps.
It can be easily verified that 𝑉(𝑃,𝑍,𝐸) is zero at the equilibrium point and positive for all other positive values of (𝑃,𝑍,𝐸). Differentiating 𝑉 with respect to the solutions of (6), a little algebraic manipulation yieldṡ𝑃̇𝑉=𝑃𝑃−𝑃^∗₁+𝑑₁𝑍̇𝑍𝑍−𝑍^∗₁+𝑑₂𝐸̇𝐸𝐸−𝐸^∗₁=𝑟𝑃1−𝑘−𝛽𝑍𝛼+𝑃𝑃−𝑃^∗₁+𝑑₁𝛽₁−𝜌𝑃𝛼+𝑃−𝑑−𝑐₁𝐸𝑍−𝑍^∗₁+𝑑₂𝑚𝑐₁𝑝−𝜏𝑍−𝑚𝐶𝐸−𝐸^∗₁=𝑃−𝑃^∗₁−𝑟𝑘𝑃−𝑃^∗₁−𝛽(𝛼+𝑃)𝛼+𝑃^∗₁𝑍−𝑍^∗₁𝛼+𝑃^∗₁−𝑍^∗𝑃−𝑃^∗₁+𝑑₁𝛽₁−𝜌𝑃𝛼+𝑃−𝑑−𝑐₁𝐸𝑍−𝑍^∗₁+𝑑₂𝑚𝑐₁𝑝−𝜏𝑍−𝑚𝐶𝐸−𝐸^∗₁=−𝑟𝑘+𝛽𝑍^∗₁(𝛼+𝑃)(𝛼+𝑃^∗)𝑃−𝑃^∗₁²−𝛽𝛼+𝑃𝑃−𝑃^∗₁𝑍−𝑍^∗₁+𝑑₁𝛽₁−𝜌𝑃−𝛽𝛼+𝑃¹−𝜌𝑃^∗₁𝛼+𝑃^∗₁𝑍−𝑍^∗₁−𝑑₁𝑐₁𝑍−𝑍^∗₁𝐸−𝐸^∗₁+𝑑₂𝑚𝑐₁𝑝−𝜏𝑍−𝑍^∗₁𝐸−𝐸^∗₁. (27) Choose arbitrary constants 𝑑₁ and 𝑑₂ as 𝑑₁=𝛽𝛼+𝑃^∗₁𝛼𝛽₁−𝜌,𝑑₂=𝛽𝛼+𝑃^∗₁𝛼𝑚𝛽₁−𝜌𝑝−𝜏 (28) then 𝑑𝑉(𝑃,𝑍,𝐸)=−𝑟𝑑𝑡𝑘+𝛽𝑍^∗₁(𝛼+𝑃)𝛼+𝑃^∗₁𝑃−𝑃^∗₁². (29)
We only need to consider the sign of −(𝑟/𝑘)+(𝛽𝑍^∗₁/(𝛼+𝑃)(𝛼+𝑃^∗₁)). Let 𝑀=−(𝑟/𝑘)+(𝛽𝑍^∗₁/(𝛼+𝑃)(𝛼+𝑃^∗₁)); then 𝑀=−(𝑟/𝑘)+(𝑟/(𝛼+𝑃))(1−(𝑃^∗₁/𝑘))≤(𝑟/𝛼𝑘)(𝑘−𝛼−𝑃^∗₁). Since 𝑃^∗₁ satisfies (9), by the discussing of the last section, we get 𝑀≤0, and hence, 𝑉^′(𝑡)≤0.
By Theorem 5.3.1 in [17], solution is limited to Φ, the largest invariant subset of {𝑉^′(𝑡)=0}. Clearly, we see that 𝑉^′(𝑡)=0 if and only if 𝑃=𝑃^∗₁. Noting that Φ is invariant, for each element in Φ, we have 𝑃(𝑡)=𝑃^∗₁. It follows from the first equation of system (6) that 0=𝑑𝑃/𝑑𝑡=𝑟𝑃^∗₁(1−(𝑃^∗₁/𝑘))−(𝛽𝑃^∗₁𝑍/(𝛼+𝑃^∗₁)), which yields 𝑍(𝑡)=𝑍^∗₁. From the third equation of system (6), 𝐸(𝑡)=𝐸^∗₁. Hence, 𝑉^′(𝑡)=0 if and only if 𝑃=𝑃^∗₁, 𝑍(𝑡)=𝑍^∗₁, and 𝐸(𝑡)=𝐸^∗₁. Accordingly, the global asymptotic stability of 𝑆₃ 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 𝜏<((𝛼𝑚𝑟𝑐₁𝑝−𝛽(𝑚𝐶+𝛾))/𝛼𝑚𝑟𝑐₁), 𝑘>𝑘_𝑡2, and then interior equilibrium 𝑆₃ 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. If 𝑘=𝑘_𝑡1, then the system (6) exhibits a transcritical bifurcation about 𝑆₂ which is branched out from 𝑆₁ at 𝑘=𝑘_𝑡1.

Proof. By analysing Section 2, one knows that the system (6) only has two feasible equilibria 𝑆₀,𝑆₁, when 𝑘<𝑘_𝑡1. In fact, the system (6) has another equilibrium 𝑆₂=(𝑃^∗,𝑍^∗,0); in this case, we know that 𝑍^∗<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 𝑘<𝑘_𝑡1, equilibrium 𝑆₂ is a saddle point and equilibrium 𝑆₁ is stable. Equilibria 𝑆₁ and 𝑆₂ coalesce into equilibrium 𝑆₁ which becomes a nonhyperbolic equilibrium at 𝑘=𝑘_𝑡1; when 𝑘>𝑘_𝑡1, the system (6) appears as a boundary equilibrium 𝑆₂ which is stable by Theorem 1, and equilibrium 𝑆₁ becomes a saddle point by (13). Thus, an exchange of stability has occurred at 𝑘=𝑘_𝑡1. This type of bifurcation is called a transcritical bifurcation.

Theorem 8. If 𝑘=𝑘_𝑡2, then the system (6) exhibits a transcritical bifurcation about 𝑆₃ which is branched out from 𝑆₂ at 𝑘=𝑘_𝑡2.

Proof. Obviously,𝑘_𝑡1<𝑘_𝑡2<𝑘_ℎ. When 𝑘_𝑡1<𝑘<𝑘_𝑡2, the system (6) has equilibria 𝑆₀,𝑆₁,𝑆₂, and  𝑆₃; in this case, 𝐸^∗₁<0. Based on the biological significance, we abandon equilibrium 𝑆₃. For equilibrium 𝑆₃, it is a saddle point as 𝐴₂<0 in (20), and equilibrium 𝑆₂ is stable by (19). Equilibria 𝑆₂ and 𝑆₃ coalesce into equilibrium 𝑆₂ which becomes a nonhyperbolic equilibrium at 𝑘=𝑘_𝑡2. When 𝑘>𝑘_𝑡2, the system (6) appears as a positive interior equilibrium 𝑆₃ which is stable by Theorem 1, and equilibrium 𝑆₂ becomes a saddle point by (19). Thus, an exchange of stability has occurred at 𝑘=𝑘_𝑡2. This type of bifurcation is called a transcritical bifurcation.

When condition (10) is satisfied, accordingly by (19), the characteristic equation of 𝑆₂ has a pair of pure imaginary roots at 𝑘=𝑘_ℎ, 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 𝑆₂=(𝑃^∗,𝑍^∗,0) and it is unstable at 𝑘=𝑘_ℎ.

Proof. The eigenvalues of 𝐽(𝑆₂) can be expressed as the solutions of the characteristic equation 𝜆³+𝑎₀𝜆²+𝑎₁𝜆+𝑎₂=0, (30) where 𝑎₀=−𝑟−2𝑟𝑃^∗𝑘−𝛽𝛼𝑍^∗(𝛼+𝑃^∗)²−𝑚𝑐₁𝑝−𝜏𝑍^∗,𝑎−𝑚𝐶−𝛾₁=𝑟−2𝑟𝑃^∗𝑘−𝛽𝛼𝑍^∗(𝛼+𝑃^∗)²𝑚𝑐₁𝑝−𝜏𝑍^∗+𝛽−𝑚𝐶−𝛾₁−𝜌𝛼𝛽(𝛼+𝑃^∗)³𝑍^∗𝑃^∗,𝑎₂=𝑚𝑐₁𝑝−𝜏𝑍^∗𝛽−𝑚𝐶−𝛾¹−𝜌𝛼𝛽(𝛼+𝑃^∗)³𝑍^∗𝑃^∗. (31)
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 𝑆₂ has a pair of pure imaginary roots at 𝑘=𝑘_ℎ; obviously, we only need to verify the transversality conditions 𝑑𝑑𝑘(Re𝜆)_𝑘=𝑘ℎ≠0, (32) where 𝜆(𝑘)=𝜇(𝑘)±𝑖𝜔(𝑘), with 𝜇(𝑘_ℎ)=0, 𝜔(𝑘_ℎ)=𝜔₀.
In this case, the eigenvalues are 𝜆₁=[𝑚𝑐₁(𝑝−𝜏)𝑍^∗−𝑚𝐶−𝛾] and 𝜆_2,3=𝜆(𝑘)=𝜇(𝑘)±𝑖𝜔(𝑘).
Hence, 𝜆³+𝑎₀𝜆²+𝑎₁𝜆+𝑎₂=𝜆−𝜆₁𝜆²−2𝜇(𝑘)𝜆+𝜇(𝑘)²+𝜔𝑘². (33) Comparative coefficient 𝜇(𝑘)=12𝑟−2𝑟𝑃^∗𝑘−𝛽𝛼𝑍^∗(𝛼+𝑃^∗)². (34) Therefore, 𝑑𝑑𝑘(Re𝜆)_𝑘=𝑘ℎ=𝑑𝑑𝑘(𝜇(𝑘))_𝑘ℎ=𝑟𝛼𝑃^∗+2𝑟𝑃^∗22𝑘^∗2(𝛼+𝑃^∗)≠0. (35) Hence, all the conditions for Hopf bifurcation are satisfied. By (19), we know 𝜆₁>0 at 𝑘=𝑘_ℎ; therefore, Hopf bifurcation of 𝑆₂ is unstable.

Remark 10. By Theorems 1 and 9, when environmental carrying capacity 𝑘=𝑘_ℎ, the system is unstable and exhibits a Hopf bifurcation around 𝑆₂; 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 𝑆₃, 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 Figure 2.

Figure 2 

(1) Equilibria 𝑆0 and 𝑆1 always exist, where 𝑆0 is unstable. (2) When 0<𝑘<𝑘𝑡1, 𝑆1 is stable; when 0<𝑘<𝑘0, 𝑆1 is globally asymptotically stable; when 𝑘>𝑘𝑡1, 𝑆1 is unstable. (3) When 𝑘=𝑘𝑡1, equilibrium point 𝑆2 is branched out from 𝑆1; 𝑘=𝑘𝑡1 is a transcritical bifurcation value. (4) When 𝑘1<𝑘<𝑘𝑡2, equilibrium point 𝑆2 is asymptotically stable. (5) When 𝑘=𝑘𝑡2, equilibrium point 𝑆3 is branched out from 𝑆2; 𝑘=𝑘𝑡2 is a transcritical bifurcation value. (6) When 𝑘>𝑘𝑡2, 𝑆3 is globally asymptotically stable. (7)𝑘=𝑘ℎ is a Hopf bifurcation value of 𝑆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) which equals the net economic revenue to the fisherman (perceived rent) plus the economic revenue to the regulatory agency. Note that 𝜋(𝑃,𝑍,𝐸)=𝑐₁𝑝𝑍−𝐶𝐸. (37)

For optimal harvest policy, this objective amounts to maximizing the present value 𝐽 of a continuous time steam of revenues given by 𝐽=^∞₀𝑒^−𝛿𝑡𝑐₁𝑝𝑍−𝐶𝐸𝑑𝑡, (38) where 𝛿 denotes the instantaneous annual rate of discount [21].

Our objective is to determine a tax policy 𝜏=𝜏(𝑡) to maximize 𝐽 subject to the state equation of (6) by invoking Pontryagin’s maximum principle [22]; the control variable 𝜏(𝑡) 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) where 𝜆₁,𝜆₂, and 𝜆₃ 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 𝜕𝐻𝜕𝜏=−𝜆₃𝑚𝑐₁𝑍𝐸=0⟹𝜆₃=0. (40) According to Pontryagin’s maximum principle, the adjoint equations are 𝑑𝜆₁𝑑𝑡=−𝜕𝐻,𝜕𝑃𝑑𝜆₂𝑑𝑡=−𝜕𝐻,𝜕𝑍𝑑𝜆₃𝑑𝑡=−𝜕𝐻𝜕𝐸=0. (41) Substitution and simplification yield 𝑑𝜆₁𝑑𝑡=−𝜕𝐻𝜕𝑃=−𝜆₁𝑟−2𝑟𝑃𝑘−𝛽𝛼𝑍(𝛼+𝑃)²−𝜆₂𝛽₁−𝜌𝛼𝑍(𝛼+𝑃)²,𝑑𝜆₂𝑑𝑡=−𝜕𝐻𝜕𝑍=−𝑐₁𝑝𝑒^−𝛿𝑡𝐸+𝜆₁𝛽𝑃𝛼+𝑃−𝜆₂𝛽₁−𝜌𝑃−𝛼+𝑃𝑑+𝑐₁𝐸,𝑑𝜆₃𝑑𝑡=−𝜕𝐻𝑐𝜕𝐸=−₁𝑝𝑍−𝐶𝑒^−𝛿𝑡+𝜆₂𝑐₁𝑍=0. (42)(43)

To obtain an optimal equilibrium solution, the solution of (43) is discussed in the following by considering the interior equilibrium as 𝜆₂=𝑒^−𝛿𝑡𝐶𝑝−𝑐₁𝑍^∗₁. (44) Let 𝐴₁=𝑟−2𝑟𝑃^∗₁𝑘−𝛽𝛼𝑍^∗₁𝛼+𝑃^∗₁²𝐴₂=𝛽₁−𝜌𝛼𝑍^∗₁𝛼+𝑃^∗₁²𝐶𝑝−𝑐₁𝑍^∗₁𝐴₃=−𝑐₁𝑝𝐸^∗₁+𝐴₂𝐴₁+𝛿𝛽𝑃^∗₁𝛼+𝑃^∗₁. (45) We will rewrite (42) by considering the interior equilibrium as 𝑑𝜆₁𝑑𝑡=−𝐴₁𝜆₁−𝐴₂𝑒^−𝛿𝑡,𝑑𝜆₂𝑑𝑡=𝐴₃𝑒^−𝛿𝑡. (46) Solutions of the above linear differential equation are 𝜆₁(𝑡)=𝐴₂𝐴₁𝑒+𝛿^−𝛿𝑡,𝜆₂(𝑡)=𝐴₃𝛿𝑒^−𝛿𝑡. (47) Substituting the value of 𝜆₂(𝑡) from (44) into (47), we get 𝐶𝑝−𝑐₁𝑍^∗₁=𝐴₃𝛿. (48)

Now using the value of 𝑃^∗₁,𝑍^∗₁, and 𝐸^∗₁ 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 𝑃=𝑃(𝜏^∗), 𝑍=𝑍(𝜏^∗), and 𝐸=𝐸(𝜏^∗). 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) Putting the value of 𝜆₂(𝑡) into (48), we get 𝑐₁𝐴𝑝𝑍−³𝛿𝑐₁𝑍=𝐶. (50) When 𝛿→∞, (50) leads to the result 𝑐₁𝑝𝑍=𝐶 which implies that the economic rent is completely dissipated.(2)By (48), we get the optimal equilibrium populations 𝑃=𝑃(𝜏^∗), 𝑍=𝑍(𝜏^∗), and 𝐸=𝐸(𝜏^∗); hence, we have 𝜋=𝑐₁𝐴𝑝𝑍−𝐶𝐸=³𝛿𝑐₁𝑍𝐸. (51) 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)Let 𝑟=5, 𝑘=6, 𝛼=9, 𝛽=0.6, 𝛽₁=0.43, 𝜌=0.2, 𝑐₁=0.64, 𝑚=0.53, 𝑝=0.7, 𝛾=0.1, 𝑑=0.2, 𝜏=0.01, and 𝐶=0.3. It is easy to verify that 6=𝑘<𝑑𝛼/(𝛽₁−𝜌)=7.826<𝑑𝛼/(𝛽₁−𝜌−𝑑); then the condition of Theorem 5 is satisfied. Hence, the equilibrium 𝑆₁(6,0,0) is globally asymptotically stable, which is shown in Figure 3.(2)We choose a set of parameters as follows: 𝑟=1, 𝑘=5.5, 𝛼=4, 𝛽=1,𝛽₁=0.7, 𝜌=0.1, 𝑐₁=0.1, 𝑚=0.5, 𝑝=0.7, 𝛾=0.1, 𝑑=0.1, 𝜏=0.2, and 𝐶=0.1. It easy to compute that 𝑃^∗=0.8,𝑍^∗𝑍=4.10,_∗<𝑚𝐶+𝛾𝑚𝑐₁𝑃𝑝−𝜏=6,_∗<𝑘<𝛼+2𝑃^∗=5.6. (52)So condition (3)(ii) of Theorem 1 holds. The equilibrium 𝑆₂(0.8,4.10,0) is locally asymptotically stable, which is shown in Figure 4(a).(3)Choose 𝑘=7 and other parameters are the same as case (2); then 𝑃^∗=0.8,𝑍^∗𝑍=4.25,_∗<𝑚𝐶+𝛾𝑚𝑐₁𝑝−𝜏=6,𝑘>𝛼+2𝑃^∗=5.6. (53)Therefore, the equilibrium 𝑆₂(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: 𝑟=16, 𝑘=1.6, 𝛼=1, 𝛽=1, 𝛽₁=0.9, 𝜌=0.1, 𝑐₁=0.1, 𝑚=0.5, 𝑝=0.7, 𝛾=0.01, 𝑑=0.1, 𝜏=0.2, and𝐶=0.1. It is easy to verify that 𝑘=1.6>𝑘_𝑡2=0.1615,𝜏=0.2<min𝑝,𝛼𝑚𝑟𝑐₁𝑝−𝛽𝑚𝐶+𝛾𝛼𝑚𝑟𝑐₁=0.625. (54)So, the conditions of Theorem 5 are satisfied; then interior equilibrium 𝑆₃=(1.5,2.4,3.8) is globally asymptotically stable, which is shown in Figure 5.

Figure 3 

The stability of the equilibrium 𝑆1.

(a)

(b)

Figure 4 

The stability of the equilibrium 𝑆2. (a) The equilibrium 𝑆2(0.8,4.10,0) is locally asymptotically stable when 𝑘=5.5<5.6. (b) The equilibrium 𝑆2(0.8,4.25,0) exhibits a Hopf bifurcation when 𝑘=7>5.6.

(a)

(b)

Figure 5 

The globally asymptotical stability of positive equilibrium point. (a) Time series portrait. (b) The phase portrait with different initial values.

From the numerical examples discussed above, we may note the following points: the equilibrium 𝑆₂ exists in absence of harvesting effort but at a lower population level (𝑃^∗,𝑍^∗,0)=(0.8,4.10,0) compared to equilibrium 𝑆₃=(𝑃^∗₁,𝑍^∗₁,𝐸^∗₁)=(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)Let 𝑟=4, 𝑘=1.6, 𝛼=1, 𝛽=1, 𝛽₁=0.9, 𝜌=0.1, 𝑐₁=0.1, 𝑚=0.5, 𝑝=0.7, 𝛾=0.01, 𝑑=0.1, and𝐶=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 𝑆^∗(0.676,3.871,2.227) is found to be stable because condition of Theorem 1 holds in this case.

Figure 6 

The optimal equilibrium 𝑆∗(0.676,3.871,2.227) is stable with different initial levels.

In Figures 7 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)

(b)

Figure 7 

(a) Variation of phytoplankton population with time for different tax levels. (b) Variation of zooplankton population with time for different tax levels.

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 equilibrium 𝑆₂ has been established with carrying capacity as the bifurcation parameter. That is, when carrying capacity 𝑘 crosses a threshold value 𝑘_ℎ, 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors wish to express their gratitude to the editors and the reviewers for the helpful comments. This work is supported in part by NNSF of China (no. 11371112) and by the Heilongjiang Provincial Natural Science Foundation (no. A201208).