# Complex Dynamics in a Harvested Nutrient-Phytoplankton-Zooplankton Model with Seasonality

Authors: Chao Liu; Peiyong Liu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A nutrient-phytoplankton-zooplankton model is established, where harvest effort on phytoplankton population and seasonal intrinsic growth of nutrient biomass are considered. Positivity and boundedness of solutions of model system are studied. By analyzing the characteristic equation of linearized system corresponding to the proposed model system, local asymptotic stability around interior equilibrium is studied, which reveals that interior equilibrium loses its stability at some critical value of predation rate. By utilizing Poincaré surface of section and computation of Lyapunov exponents, numerical experiments are performed to discuss the complex dynamical behavior of model system. Model system shows rich dynamic behavior including limit cycle and chaotic attractor due to variation of predation rate and seasonal intrinsic growth rate. Further numerical experiments show that reasonable harvesting has stabilizing effect on population dynamics that prevents chaotic behavior. Numerical simulations are carried out to show consistency with theoretical analysis obtained in this paper.


---

## Body

## 1. Introduction

Plankton are microscopic organisms that float freely within marine system. They are made up of tiny plants (called phytoplankton) and tiny animals (called zooplankton). Phytoplankton population is the dominant primary producer in the marine system, which produces approximately forty percent oxygen through photosynthesis. By primary production, death, and sinking, they exert a global scale influence on climate by effectively transportingCO2 from the surface layer to deep marine sediments [1]. Phytoplankton population feeds on inorganic nutrients such as nitrate and phosphate for growth [2]. Furthermore, phytoplankton population comprises a large number of different species and in the bottom of the food chain that support the life of all higher marine organisms such as zooplankton and fish. Consequently, the abundance of phytoplankton population plays a significant role in marine reserves and fishery management [3].The dynamics of marine system have long been one of the dominant themes in mathematical biology. Marine population models are generally based on compartmental models consisting of many coupled nonlinear differential equations. Some of these models seek to contain as many as physical and biological mechanisms. In recent years, there is a growing explicit biological physiological body of evidence [4–13] in marine ecosystem, where nutrient recycling is considered and authors have dealt with a nutrient-plankton model in an aquatic environment in the context of phytoplankton bloom. In [4], infective disease among phytoplankton is considered and dynamical behavior of nutrient-phytoplankton population is investigated. However, interaction mechanism of phytoplankton and zooplankton is not discussed. In [5, 6, 8], authors pay particular attention to the population dynamics of a layer of single phytoplankton species growing over a pool of nutrients. In [7, 9, 10], models of nutrient-plankton interaction with a toxic substance that inhibits either the growth rate of phytoplankton, zooplankton, or both trophic levels are proposed and studied. In [11–13], nutrient-phytoplankton dynamical models with instantaneous and time delayed recycling are proposed, authors investigate the dynamics and examine the responses to model complexities, while effects of predation from zooplankton and nutrient recycling within nutrient-phytoplankton-zooplankton ecosystem are not studied.One of the most interesting topics in mathematical ecology concerns the seasonal intrinsic growth of nutrient biomass in marine ecosystem. Countless organisms live in seasonal marine environment, and many parameters of such system may oscillate simultaneously and not necessarily in plane. It is pertinent to note that intrinsic growth of nutrient biomass is usually periodically adjusted, which may be caused by periodic variation of fecundity of nutrient within marine ecosystem. Recently, the model systems of marine ecosystem with periodic external force have gathered new attention because of their complex dynamical behavior, especially chaotic behavior. A variety of studies have been performed on the interactions between seasonality and internal biological rhythms, whose dynamical effects on nutrient-phytoplankton ecosystem are discussed in [14–17], where complex dynamical behavior such as multiplicity of attractors, catastrophes, and deterministic chaos are found. However, seasonality discussed in the above work [14–17] only concentrates on intrinsic growth of phytoplankton population. As introduced in the above section, phytoplankton population density directly relates to the seasonal variational abundance of nutrient, it is necessary to investigate dynamics of model system due to seasonal variation of nutrient biomass within marine ecosystem. As far as knowledge goes, there is no work focusing on seasonality of nutrient biomass. Furthermore, control strategies for chaotic behavior are not discussed in [14–17]. Generally, ecological system is often deeply perturbed by human exploiting activities [3, 18–20]. In recent years, there has been growing interest in the study of harvested phytoplankton ecosystem, and several dynamical models with harvest effort are investigated in [21] and the references therein. A phytoplankton-zooplankton model with harvesting is proposed and investigated in [21], dynamical behavior and optimal harvesting policy is investigated, while nutrient recycling within marine ecosystem is not studied. To the best of author's knowledge, nobody has explicitly proposed a harvested nutrient-phytoplankton-zooplankton model with seasonality on nutrient biomass and studied its complex dynamical behavior.The organization of this paper is as follows. By incorporating harvest effort on phytoplankton population and seasonal intrinsic growth of nutrient biomass, we extend the model proposed in [11–13, 17, 21]; a harvested nutrient-phytoplankton-zooplankton model with seasonality is established in the second section. Positivity and boundedness of solutions of model system are discussed in the third section. Qualitative analysis of model system is performed in the fourth section; local asymptotic stability analysis of model system around the interior equilibrium is studied, which is utilized to discuss the population dynamics. Transcritical bifurcation and Hopf bifurcation phenomenon is investigated due to variation of intrinsic growth rate and predation rate, respectively. Based on Poincaré surface-of-section technique and computation of Lyapunov exponents, complex dynamical behavior of model system is investigated due to variation of predation rate and seasonal intrinsic growth of nutrient biomass. Further numerical experiments are performed to discuss the stabilizing effect of harvest effort on population dynamics. Numerical simulations are provided to support the analytical findings obtained in this paper. Finally, this paper ends with a conclusion.


## 2. Model Formulation

The classical Rosenzweig-MacArthur prey-predator model [3] is as follows:
(1)dX(t)dt=rX(t)(1-X(t)K)-aX(t)Y(t)b+X(t),dY(t)dt=caX(t)Y(t)b+X(t)-dY(t),
where X(t) and Y(t) represent density of prey population and predator population, respectively. r denotes the intrinsic growth rate of the prey population and the carrying capacity of the prey population is K. a is the searching efficiency or the predation rate of the predator population; b denotes the half saturation constant and c<1 is the conversion factor. d represents the death rate of the predator population. All parameters are all positive constants.In the following part, some hypotheses are proposed, which are utilized to construct the mathematical model and discuss population dynamics of a harvested nutrient-phytoplankton-zooplankton ecosystem with seasonality.(H1)
It is assumed that the phytoplankton population feeds on nutrient and the zooplankton population grazes on the phytoplankton for survival. Zooplankton population is solely dependent upon phytoplankton as their most favorable food source and a variation in phytoplankton population density has a great impact on the growth of zooplankton population. In this paper, we will consider a three-tier model of nutrient-phytoplankton-zooplankton with the usage of model system (1). Let N(t), P(t), and Z(t) denote nutrient biomass, phytoplankton population density, and zooplankton population density with respect to time t, respectively. Let K represent the carrying capacity of the nutrient biomass within marine ecosystem. a is the searching efficiency of phytoplankton population, and h is the predation rate; b1 and b2 denote the half saturation constant. c1<1 is conversion factor from nutrient biomass to phytoplankton population. c2<1 is conversion factor from phytoplankton population to zooplankton population.(H2)
Due to its extensive utilization in the field of food and hygienic field, phytoplankton population is extensively exploited in the real world. In this paper, we will consider the harvest effort on phytoplankton population. A scalarE≥0 denotes the harvesting effort to phytoplankton population, constant q is the catchability coefficient of phytoplankton population, and the harvesting term qEN(t) follows the catch per unit effort hypothesis [22].(H3)
It is observed that nutrient biomass is recycled within nutrient-phytoplankton-zooplankton ecosystem [1], due to death, washout, or some natural calamity; let the loss of phytoplankton and zooplankton population be given by d1P(t) and d2Z(t) while through recycling (nutrient biomass decomposition). r1 and r2 represent amount of phytoplankton and zooplankton population that converts back into nutrient biomass, respectively. It follows from the biological interpretations of parameters that d1>r1 and d2>r2.(H4)
It is assumed that the intrinsic growth of nutrient biomass varies periodically, which may be caused by periodic variation of fecundity of nutrient within marine ecosystem. Consequently, the intrinsic growth rate of nutrient biomassr is assumed to be a periodically varying function of time and is considered as r=r0(1+Γsin(ωt)). It should be noted that the average value of r is assumed to be r0 and Γ represents the degree of seasonality; r0Γ is the magnitude of the perturbation and ω is the angular frequency of the fluctuations caused by seasonality. Since r can not be negative, Γ lies in the interval [0,1]. Γ=0 relates to absence of seasonality, and Γ=1 means the maximum value of the parameter is twice its average value.Based on model system (1) and hypotheses (H1)–(H4), a nutrient-phytoplankton-zooplankton model with harvest effort and seasonality takes the following form:
(2)dN(t)dt=r0(1+Γsin(ωt))N(t)(1-N(t)K)-aN(t)P(t)b1+N(t)+r1P(t)+r2Z(t),dP(t)dt=c1aN(t)P(t)b1+N(t)-hP(t)Z(t)b2+P(t)-d1P(t)-qEP(t),dZ(t)dt=c2hP(t)Z(t)b2+P(t)-d2Z(t),
where all parameters share the same interpretations introduced in hypotheses (H1)–(H4); model system (2) is analyzed with the following initial conditions:
(3)N(0)>0,P(0)>0,Z(0)>0.


## 3. Positivity and Boundedness of Solution

Theorem 1.
All solutions of model system (2) with initial conditions are positive.Proof.
The model system (2) can be interpreted as the matrix form
(4)dU(t)dt=H(U(t)),
where U(t)=(N(t),P(t),Z(t))T∈ℝ3, and H(U(t)) are given as follows: (5)H(U(t))=(H1(U(t))H2(U(t))H3(U(t)))=(r0(1+Γsin(ωt))N(t)(1-N(t)K)-aN(t)P(t)b1+N(t)+r1P(t)+r2Z(t)caN(t)P(t)b1+N(t)-hP(t)Z(t)b2+P(t)-d1P(t)-qEP(t)hP(t)Z(t)b2+P(t)-d2Z(t)). Let ℝ+3=[0,∞)3 be the positive octant in ℝ3; then L:ℝ+3+1→ℝ3 is locally Lipschitz and satisfies the condition Hi(U(t))|U∈ℝ3≥0, i=1,2,3.Due to the lemma in [23] and Theorem A.4 in [22], any solution of the model system (2) with positive initial conditions exist uniquely and each component of the solution remains within the interval [0,A0) for some A0>0. Standard and simple arguments show that solutions of model system (2) always exist and stay positive. Hence, this completes the positivity of solutions of model system (2).Theorem 2.
Solutions of model system (2) are bounded in the following region:
(6)Ω={K(2r0+e)28er0(N(t),P(t),Z(t))∣N(t)+P(t)+Z(t)≤K(2r0+e)28er0},
where e=min⁡{qE+d1-r1,d2-r2}.Proof.
LetW(t)=N(t)+P(t)+Z(t); it follows from model system (2) and (H1)–(H4) that
(7)dW(t)dt≤2r0N(t)(1-N(t)K)-a(1-c1)N(t)P(t)b1+N(t)-h(1-c2)P(t)Z(t)b2+P(t)-(d1-r1+qE)P(t)-(d2-r2)Z(t).
According to Theorem1, it derives that
(8)dW(t)dt≤(2r0+e)N(t)-2r0N2(t)K-eW(t),≤K(2r0+e)28r0-eW(t),
where e=min⁡{qE+d1-r1,d2-r2}.
Based on the above analysis, it can be obtained that(9)W(t)≤K(2r0+e)28er0.


## 4. Qualitative Analysis of Model System

In this section, qualitative analysis of model system (2) is performed. By choosing appropriate parameter as a bifurcation parameter, local asymptotic stability including transcritical bifurcation and Hopf bifurcation around interior equilibrium is, respectively, studied in the absence of seasonality. Complex dynamical behavior including limit cycle and chaotic attractor are discussed due to variation of predation rate with the help of numerical experiments. Furthermore, dynamic effect of seasonality and harvest effort on population dynamics of model system is investigated based on Poincaré surface-of-section technique and computation of Lyapunov exponents.


### 4.1. Local Stability Analysis of Model System without Seasonality

In the case ofΓ=0, model system (2) without seasonality takes the following form:
(10)dN(t)dt=r0N(t)(1-N(t)K)-aN(t)P(t)b1+N(t)+r1P(t)+r2Z(t),dP(t)dt=caN(t)P(t)b1+N(t)-hP(t)Z(t)b2+P(t)-d1P(t)-qEP(t),dZ(t)dt=hP(t)Z(t)b2+P(t)-d2Z(t).In order to reduce number of parameters and facilitate dynamical analysis of model system, model system (10) is dimensionalized with the following scaling:
(11)N¯=NK,P¯=PK,Z¯=ZK,T=r0t.
Then model system (10) takes the following form:
(12)dN¯(T)dT=N¯(T)(1-N¯(T))-αN¯(T)P¯(T)1+β1N¯(T)+f1P¯(T)+f2Z¯(T),dP¯(T)dT=m1N¯(T)P¯(T)1+β1N¯(T)-h1P¯(T)Z¯(T)1+β2P¯(T)-ϵ1P¯(T)-E1P¯(T),dZ¯(T)dT=m2P¯(T)Z¯(T)1+β2P¯(T)-ϵ2Z¯(T),
where α=aK/b1r0, β1=K/b1, m1=ac1K/b1r0, f1=r1/r0, f2=r2/r0, h1=hK/b2r0, m2=c2hK/b2r0, ϵ1=d1/r0, ϵ2=d2/r0, β2=K/b2, and E1=qE/r0.Since the interior equilibrium biologically relates to simultaneous survival of nutrient biomass, phytoplankton population, and zooplankton population, we will mainly concentrate on dynamical analysis of model system around the interior equilibrium in this paper.According to model system (12), the interior equilibrium M*(N¯*,P¯*,Z¯*) is as follows:
(13)P¯*=ϵ2m2-ϵ2β2,Z¯*=m2[m1N¯*-(ϵ1+E1)(1+β1N¯*)]h1(m2-ϵ2β2)(1+β1N¯*),
where N¯* satisfies the following equation:
(14)N¯3+C1N¯2+C2N¯+C3=0,
where C1=(1-β1)/β1, C2=(m2f2[β1(ϵ1+E1)-m1]+h1[ϵ2(β2+α-f1β1)-m2])/h1β1(m1-ϵ2β2), and C3=(f2m2(ϵ1+E1)-ϵ2f1h1)/h1β1(m2-ϵ2β2).Based on Routh-Hurwitz criterion [3], a simple condition that (14) has at least one positive roots is as follows:
(15)C1C2-C3≤0.Furthermore, according to the biological interpretations of interior equilibrium, phytoplankton population and zooplankton population all exist provided thatP¯*>0, Z¯*>0, which follows that
(16)m2-ϵ2β2>0,N¯*>ϵ1+E1m1-β1(ϵ1+E1).In Theorem3, local stability analysis around the interior equilibrium M* is performed.Theorem 3.
Model system (12) is locally asymptotically stable around the interior equilibrium M* if the following conditions hold:
(17)A1>0,A3>0,A1A2-A3>0,
where Ai, i=1,2,3 are defined as follows:
(18)A1=((1+β1N¯*)2[ϵ2β2(ϵ1+E1)+h1(2N¯*-1)](m2-ϵ2β2)×(1+β1N¯*)2-m1ϵ2β2(m2-ϵ2β2)N¯*×(1+β1N¯*)+αh1ϵ2(1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)2)-1,A2=((h1+β2(1-2N¯*-ϵ2))-αβ2ϵ2(1+β1N¯*)2](1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(h1+β2(1-2N¯*-ϵ2))-αβ2ϵ2(1+β1N¯*)2](1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)3)-1,A3=((1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(2N¯*-1)+αϵ2-m1f2](1+β1N¯*)2)×(h1(1+β1N¯*)2)-1.Proof.
Firstly, the Jacobian matrix of the model system (12) around M*(N¯*,P¯*,Z¯*) is derived as follows: (19)J=(1-2N¯*-αN¯*(1+β1N¯*)2-αN¯*1+β1N¯*+f1f2m1P¯*(1+β1N¯*)2m1N¯*1+β1N¯*-h1Z¯*(1+β2P¯*)2-ϵ1-E1-h1P¯*1+β2P¯*0m2Z¯*(1+β2P¯*)2m2P¯*1+β2P¯*-ϵ2).
By virtue of (13) and (19), the characteristic equation of model system (12) around the interior equilibrium M* is as follows:
(20)λ3+A1λ2+A2λ+A3=0,
where Ai, i=1,2,3 are defined as follows:
(21)A1=((1+β1N¯*)2[ϵ2β2(ϵ1+E1)+h1(2N¯*-1)](m2-ϵ2β2)×(1+β1N¯*)2-m1ϵ2β2(m2-ϵ2β2)N¯*×(1+β1N¯*)+αh1ϵ2(1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)2)-1,A2=((1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(h1+β2(1-2N¯*-ϵ2))-αβ2ϵ2(1+β1N¯*)2](1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)3)-1,A3=((1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(2N¯*-1)+αϵ2-m1f2(1+β1N¯*)2](1+β1N¯*)2)×(h1(1+β1N¯*)2)-1.
Based on Routh-Hurwitz criterion [3], model system (12) is locally asymptotically stable around M* provided that A1>0, A3>0, A1A2-A3>0.By denoting the quantitiesAi(h1), i=1,2,3 as smooth functions of the parameter h1 (predation rate), local stability switch and dynamical behavior is discussed due to the variation of predation rate h1.Theorem 4.
When the predation rateh1 crosses a critical value h1*, model system (12) enters into Hopf bifurcation around the interior equilibrium M* if the following conditions hold:(i)
A1(h1*)>0,A3(h1*)>0;(ii)
A1(h1*)A2(h1*)-A3(h1*)=0;(iii)
d[A1(h1)A2(h1)-A3(h1)]/dh1∣h1=h1*<0.Proof.
By choosing predation rateh1 as bifurcation parameter, we will discuss local stability switch due to variation of h1; we can see that if there exists a critical value h1* such that A1(h1*)>0, A3(h1*)>0, A1(h1*)A2(h1*)-A3(h1*)=0, and d[A1(h1)A2(h1)-A3(h1)]/dh1|h1=h1*<0. It is easy to show that A2(h1*)>0 based on conditions (i) and (ii). For the Hopf bifurcation to occur at h1=h1*, the characteristic equation takes the following form:
(22)(ρ2(h1*)+A2(h1*))(ρ(h1*)+A1(h1*))=0,
which has three roots ρ1(h1*)=iA2(h1*), ρ2=-iA2(h1*), and ρ3=-A1(h1*)<0.
Furthermore, in order to show if Hopf bifurcation occurs ath1=h1*, we need to verify the transversality condition [24]:
(23)[dRe(ρ(h1))dh1]|h1=h1*≠0.
Assuming all the roots of (22) are in general of the following form:
(24)ρ1(h1)=μ(h1)+iν(h1),ρ2(h1)=μ(h1)+iν(h1),ρ3(h1)=-A1(h1),
where μ(h1) and ν(h1) are smooth functions with respect to h1.
Now we will verify the transversality condition:(25)[dRe(ρi(h1))dh1]|h1=h1*≠0,i=1,2.
Substitutingρi(h1), i=1,2 into (22) and calculating the derivative, it follows that
(26)B1(h1)μ′(h1)-B2(h1)ν′(h1)+B3(h1)=0,B2(h1)μ′(h1)+B1(h1)ν′(h1)+B4(h1)=0,
where Bi(h1) are defined as follows:
(27)B1(h1)=3μ2(h1)+2A1(h1)μ(h1)+A2(h1)-3ν2(h1),B2(h1)=6μ(h1)ν(h1)+2A1(h1)ν(h1),B3(h1)=μ2(h1)A1′(h1)+A2′(h1)μ(h1)+A3′(h1)-A1′(h1)ν2(h1),B4(h1)=2μ(h1)ν(h1)A1′(h1)+A2′(h1)ν(h1),
where Ai′(h1), i=1,2,3 denotes the derivative of Ai(h1) with respect to h1.
It follows from (22) that μ(h1*)=0, ν(h1*)=A2(h1*). By further computation, it can be obtained that
(28)B1(h1*)=-2A2(h1*),B2(h1*)=2A1(h1*)A2(h1*),B3(h1*)=A3′(h1*)-A1′(h1*)A2(h1*),B4(h1*)=A2′(h1*)A2(h1*).
By solving (26), it can be obtained that
(29)[dRe(ρi(h1))dh1]|h1=h1*=μ′(h1)h1=h1*=-B2(h1*)B4(h1*)+B1(h1*)B3(h1*)B12(h1*)+B22(h1*)=A3′(h1*)-A1′(h1*)A2(h1*)-A1(h1*)A2′(h1*)2[A12(h1*)+A2(h1*)].
If condition (iii) holds, then(30)[dRe(ρi(h1))dh1]|h1=h1*>0,
which implies that transversality condition holds and Hopf bifurcation occurs at h1=h1*.Theorem 5.
By choosingα as a bifurcation parameter, model system (12) undergoes transcritical bifurcation around M* and bifurcation value is α=(m1f2-(m2-ϵ2β2)(1+β1N¯*)2(2N¯*-1))/ϵ2.Proof.
Whenα=(m1f2-(m2-ϵ2β2)(1+β1N¯*)2(2N¯*-1))/ϵ2, it follows from simple computation that (20) has a geometrically simple zero eigenvalue with left eigenvector vl=(h1P¯*/(1+β2P¯*),f2,1) and right eigenvector vr=(1,0,1)T, and
(31)vlDUDαH(U(t))vr|M*=(h1P¯*1+β2P¯*,f2,1)×(β1N¯*(1+2β1N¯*)-1(1+β1N¯*)400000000)(101)=h1P¯*[β1N¯*(1+2β1N¯*)-1](1+β2P¯*)(1+β1N¯*)4,(32)vlDU2H(U(t))(vr,vr,vr)|M*=vl∑i=13‍(eivlTDU(DUHi)Tvr)|M*=α(β1N¯*-1)-2f2m1β1P¯*(1+β1N¯*)3-2,
where H(U(t)), U(t) have been defined in Theorem 1 of this paper and ei, i=1,2,3, are unit vectors. It follows from [25] and the above computation that model system (12) undergoes transcritical bifurcation around M*.


### 4.2.  Numerical Experiments for Model System without Seasonality

In order to facilitate the following numerical experiments, some preliminaries about Poincaré surface-of-section technique and Lyapunov exponents are introduced in Remarks6 and 7.Remark 6.
A traditional approach to gain preliminary insight into the properties of dynamical system is to carry out a one-dimensional bifurcation analysis. One-dimensional bifurcation diagrams of Poincaré maps present information about the dependence of the dynamics on a certain parameter. The analysis is expected to reveal the type of attractor to which the dynamics will ultimately settle down after passing the initial transient phase and within which the trajectory will then remain forever. On a Poincaré surface of section, the dynamical behavior can be described by a discrete map whose phase-space dimension is less than that of the original continuous flow. Chaotic flows can be understood based on concepts that are convenient for maps such as unstable orbits (see [26–28] and references therein).Case 1.
A finite number of points corresponds to a periodic solution; that is, one point corresponds to a solution of period equal to that of the forcing term, namely,2π/ω (period-one); and n points correspond to a solution (subharmonic) of period 2nπ/ω (period n or a subharmonic of order 1/n).Case 2.
A closed curve corresponds to a quasiperiodic solution, that is, a solution consisting of two incommensurate frequencies or, equivalently, having a trajectory that is dense on tours.Case 3.
A collection of points that is fractal corresponds to chaos, namely, a stranger in phase space; a collection of points that form a cloud that is disorganized, partially organized, or fuzzy may (or may not ) correspond to chaotic attractors.Remark 7.
In a given embedding dimension, the Lyapunov exponent is a measure of the speeds at which initially nearby trajectories of the system diverge. There is a Lyapunov exponent for each dimension of the process, which together constitutes the Lyapunov spectrum for the dynamical system. The Lyapunov exponent is related to the predictability of the system, and the largest Lyapunov exponent of a stable system does not exceed zero. However, a chaotic system has at least one positive Lyapunov exponent. A bounded system with a positive Lyapunov exponent is one operational definition of chaotic behaviors, which presents a quantitative measure of the average rate of separation of nearby trajectories on the attractor. Over the years, a number of methods have been introduced for computation of Lyapunov exponents (see [26, 28]). A positive Lyapunov exponent implies that a chaotic process displays long term unpredictability, with the output being sensitively dependent on the initial conditions. Even slightly different initial values can lead to vastly different system outputs. Furthermore, another characteristic for chaos is fractal dimension, which can be computed with the commutated Lyapunov exponents (see Section  9.4 of reference [29]).In the absence of seasonality, we will perform numerical experiments to observe the dynamics of the model system (2) with the following hypothetical set of parameter values, which are given as follows: r0=1, K=1, a=1.63, r1=0.01, r2=0.001, c1=0.95, b1=0.33, d1=0.4, c2=0.9, d2=0.01, b2=0.5, q=0.1, and E=0.15. In the following part, h1 (predation rate) is chosen as bifurcation parameter. By virtue of scaling introduced in model system (12), values of parameters are obtained as follows: α=4.94, β1=3.03, m1=4.69, f1=0.01, f2=0.001, m2=0.9h1, ϵ1=0.4, ϵ2=0.01, β2=2, and E1=0.015. A bifurcation diagram of Poincaré section is plotted to investigate the dynamic effect of predation rate on model system (12), which can be shown in Figure 1. In order to understand Hopf bifurcation occurring at the critical value h1*=0.0711, bifurcation diagram with (0.07,0.075) is plotted in Figure 2.Figure 1
Bifurcation diagram of Poincaré section for the nutrient and predation rateh1 is varied in [0.068,0.1].Figure 2
Bifurcation diagram of Poincaré section for the nutrient and predation rateh1 is varied in [0.07,0.075].According to (15) and (16), it can be obtained that h1>0.068 guarantees the existence of interior equilibrium. Based on Figures 1 and 2, it is observed that model system (12) is stable around the interior equilibrium when h1 varies within (0.068, 0.0711). It should be noted that h1=0.07 is randomly selected within (0.068,0.0711); the dynamical response and phase portrait of model system (12) with h1=0.07 is plotted in Figures 3 and 4, respectively, which shows that model system (12) is locally asymptotically stable around interior equilibrium M*(0.7022,0.1999,13.9114).Figure 3
Dynamical response of model system (12) with predation rate h1=0.07.Figure 4
Phase portrait of model system (12) with predation rate h1=0.07.Ash1 increases across the critical value h1*=0.0711; model system enters into Hopf bifurcation. Dynamical responses of model system (12) with h1*=0.0711 is plotted in Figure 5, and a limit cycle corresponding to periodic solution is observed in Figure 6. As h1 further increases, it leads to a chaotic region at h1∈(0.0711,0.1). The dynamical response and phase portrait of model system with h1=0.08 is plotted in Figures 7 and 8, respectively. It should be noted that the sensitive dependence of future dynamics on the current state, the signature of chaotic behavior, is apparent from the fact that all the trajectories in the handle of teacup are very close together. Based on Theorem 4, model system undergoes Hopf bifurcation at critical bifurcation value h1*=0.0711. With the increase of bifurcation parameter h1, it can observed that model system enters into period doubling from limit cycle oscillations based on Figure 2, and model system finally settles down into chaotic behavior from period doubling. Hence, it can be concluded that model system admits that the transition to chaos is realized through the cascade of period doubling bifurcations.Figure 5
Dynamical response of model system (12) with predation rate h1=0.0711.Figure 6
A limit cycle corresponding to periodic solution in Figure5.Figure 7
Dynamical response of model system (12) with predation rate h1=0.08.Figure 8
Phase portrait of model system (12) with predation rate h1=0.08.Furthermore, chaotic attractor of model system (12) with h1=0.085,0.09,0.095, and 0.1 are plotted in Figure 9, and the corresponding Poincaré section of Figure 9 is plotted in Figure 10, which coincides with Case 3 introduced in Remark 6 implying the occurrence of chaotic attractor. Lyapunov exponents (denoted by Li, i=1,2,3) with time and the corresponding Lyapunov dimensions can be computed, which are utilized to show the existence of chaotic attractor. Table 1 shows Lyapunov exponents and corresponding Lyapunov dimensions of model system (12) with h1=0.08,0.085,0.09,0.095 and 0.1.Table 1
Lyapunove exponents and corresponding Lyapunov dimensions of model system (12) with h1=0.08, 0.085, 0.09, 0.095, and 0.1.
h
1
0.08
0.085
0.09
0.095
0.1
L
1
0.13761
0.1543
0.1818
0.1913
0.2044
L
2
0
0
0
0
0
L
3
−1.1265
−1.5804
−1.5803
−1.5893
−1.5710
D
2.1222
2.0976
2.115
2.1204
2.1301Chaotic attractor of model system (12) with predation rate h1=0.085,0.09,0.095,0.1.
(a)
(b)
(c)
(d)Corresponding Poincaré section of model system (12) with h1=0.085,0.09,0.095,0.1. Poincaré points 5000–15000 are plotted.
(a)
(b)
(c)
(d)It follows from the Table1 that there is a positive Lyapunov exponent for the model system (12) with h1=0.08,0.085,0.09,0.095, and 0.1, respectively. Furthermore, it follows from the Table 1 that Lyapunov dimensions for model system (12) with h1=0.08,0.085,0.09,0.095, and 0.1 are all fractals. According to the characterization of chaotic attractors introduced in Remark 7, a bounded system with a positive Lyapunov exponent is one operational definition of chaotic behavior, and another characteristic for chaos is fractal dimension.Oscillatory population may be driven to extinction in presence of environmental stochasticity when the population density is very low [30]. Based on the above numerical experiments, chaotic dynamics for increasing the predation rate has been discussed. In the following part, role of harvesting in controlling chaotic dynamics will be investigated. Now we fixed the predation rate h1=0.08 in chaotic region and the harvest effort E1 is chosen as bifurcation parameter. Figure 11 shows the bifurcation diagram of Poincaré section for the nutrient and harvest effort E1 is varied in [0,0.15]. It should be noted that harvest effort E1 should be controlled within [0,0.15] based on (15) and (16), which guarantees the existence of interior equilibrium. From Figure 11, it is observed that chaotic behavior remains for E1 that varies from [0,0.1381). With the increase of E1, chaos enter into limit cycle oscillation at E1=0.1381 and finally the model system (12) enters into stable focus from limit cycle.Figure 11
Bifurcation diagram of Poincaré section for the nutrient and harvest effortE1 is varied in [0,0.15].


### 4.3.  Numerical Experiments for Model System with Seasonality

With the help of MATLAB, numerical experiments with hypothetical set of parameters are performed to understand the dynamical effect of seasonal intrinsic growth rate of nutrient biomass on model system (2). The values of parameters are given as follows: r0=1, K=1, a=1.63, r1=0.01, r2=0.001, c1=0.95, b1=0.33, d1=0.4, c2=0.9, d2=0.01, b2=0.5, h=0.04, q=0.1, E=0.15, and ω=2. By choosing Γ as a bifurcation parameter, bifurcation diagram of Poincaré section for the nutrient and Γ is varied in [0,1], which can be found in Figure 12. It can be observed that chaotic behavior enters into periodic attractor when bifurcation parameter Γ increases through a critical value Γ=0.8. Model system (2) admits periodic attractor Γ within the range [0.8,1].Figure 12
Bifurcation diagram of Poincaré section for the nutrient andΓ is varied in [0,1].It should be noted thatΓ=0.81 and Γ=0.95 are randomly selected from [0.8,1]. Dynamical responses of model system with Γ=0.81 and Γ=0.95 are plotted in Figures 13 and 14, respectively. Corresponding Poincaré section of model system (2) with Γ=0.81 is plotted in Figure 15(a), which coincides with Case 1 introduced in Remark 6 implying the existence of periodic attractor; and corresponding Poincaré section of model system (2) with Γ=0.95 is plotted in Figure 15(b), which coincides with Case I introduced in Remark 6 implying the existence of periodic attractor. Poincaré points 5000–15000 are plotted. Based on the above numerical experiments, chaotic behavior can be controlled by the enhancement of amplitude of seasonal intrinsic growth of nutrient biomass; such type of mechanism reduces chaotic behavior to periodic attractor.Figure 13
Dynamical response of model system (2) with Γ=0.81.Figure 14
Dynamical response of model system (2) with Γ=0.95.Corresponding Poincaré section of model system (2) with Γ=0.81 is plotted in (a), and corresponding Poincaré section of model system (2) with Γ=0.95 is plotted in (b). Poincaré points 5000–15000 are plotted.
(a)
(b)Remark 8.
It can be observed that the dynamical responses fluctuate with complex structures when the model system admits chaotic behavior. Especially, some dynamical response even comes close to zero within some time interval, which biologically reflects harmful phytoplankton bloom. Consequently, there has been considerable scientific attention towards harmful phytoplankton bloom and its associate control strategy. Based on the theoretical analysis, it is found that high nutrient levels and favorable conditions play a key role in rapid or massive growth of phytoplankton bloom; low nutrient concentration, high predation pressure from zooplankton, and other unfavorable conditions limit phytoplankton growth, which leads to oscillations or recurring bloom in the nutrient-phytoplankton-zooplankton marine ecosystem.Based on the numerical experiments shown in Figure11, harvesting has a stabilizing effect on model dynamics. Furthermore, numerical experiments show that chaotic behavior of the population due to enhancement of predation rate can be prevented by adopting a suitable harvesting policy and model system (12) eventually approaches to a stable steady state. If the harvest effort on phytoplankton population increases, the zooplankton population decreases due to lack of the phytoplankton population and nutrient absorption pressure decreases and in turn the nutrient within marine system increases. On the other hand, harvest effort should be constrained within certain interval, which may guarantee the existence of nutrient biomass, phytoplankton population, and zooplankton population within marine system. The theoretical findings obtained in the above section are of inspiration for regulating some constructive policies to control harmful phytoplankton bloom with the introduction of appropriate harvest effort.


## 4.1. Local Stability Analysis of Model System without Seasonality

In the case ofΓ=0, model system (2) without seasonality takes the following form:
(10)dN(t)dt=r0N(t)(1-N(t)K)-aN(t)P(t)b1+N(t)+r1P(t)+r2Z(t),dP(t)dt=caN(t)P(t)b1+N(t)-hP(t)Z(t)b2+P(t)-d1P(t)-qEP(t),dZ(t)dt=hP(t)Z(t)b2+P(t)-d2Z(t).In order to reduce number of parameters and facilitate dynamical analysis of model system, model system (10) is dimensionalized with the following scaling:
(11)N¯=NK,P¯=PK,Z¯=ZK,T=r0t.
Then model system (10) takes the following form:
(12)dN¯(T)dT=N¯(T)(1-N¯(T))-αN¯(T)P¯(T)1+β1N¯(T)+f1P¯(T)+f2Z¯(T),dP¯(T)dT=m1N¯(T)P¯(T)1+β1N¯(T)-h1P¯(T)Z¯(T)1+β2P¯(T)-ϵ1P¯(T)-E1P¯(T),dZ¯(T)dT=m2P¯(T)Z¯(T)1+β2P¯(T)-ϵ2Z¯(T),
where α=aK/b1r0, β1=K/b1, m1=ac1K/b1r0, f1=r1/r0, f2=r2/r0, h1=hK/b2r0, m2=c2hK/b2r0, ϵ1=d1/r0, ϵ2=d2/r0, β2=K/b2, and E1=qE/r0.Since the interior equilibrium biologically relates to simultaneous survival of nutrient biomass, phytoplankton population, and zooplankton population, we will mainly concentrate on dynamical analysis of model system around the interior equilibrium in this paper.According to model system (12), the interior equilibrium M*(N¯*,P¯*,Z¯*) is as follows:
(13)P¯*=ϵ2m2-ϵ2β2,Z¯*=m2[m1N¯*-(ϵ1+E1)(1+β1N¯*)]h1(m2-ϵ2β2)(1+β1N¯*),
where N¯* satisfies the following equation:
(14)N¯3+C1N¯2+C2N¯+C3=0,
where C1=(1-β1)/β1, C2=(m2f2[β1(ϵ1+E1)-m1]+h1[ϵ2(β2+α-f1β1)-m2])/h1β1(m1-ϵ2β2), and C3=(f2m2(ϵ1+E1)-ϵ2f1h1)/h1β1(m2-ϵ2β2).Based on Routh-Hurwitz criterion [3], a simple condition that (14) has at least one positive roots is as follows:
(15)C1C2-C3≤0.Furthermore, according to the biological interpretations of interior equilibrium, phytoplankton population and zooplankton population all exist provided thatP¯*>0, Z¯*>0, which follows that
(16)m2-ϵ2β2>0,N¯*>ϵ1+E1m1-β1(ϵ1+E1).In Theorem3, local stability analysis around the interior equilibrium M* is performed.Theorem 3.
Model system (12) is locally asymptotically stable around the interior equilibrium M* if the following conditions hold:
(17)A1>0,A3>0,A1A2-A3>0,
where Ai, i=1,2,3 are defined as follows:
(18)A1=((1+β1N¯*)2[ϵ2β2(ϵ1+E1)+h1(2N¯*-1)](m2-ϵ2β2)×(1+β1N¯*)2-m1ϵ2β2(m2-ϵ2β2)N¯*×(1+β1N¯*)+αh1ϵ2(1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)2)-1,A2=((h1+β2(1-2N¯*-ϵ2))-αβ2ϵ2(1+β1N¯*)2](1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(h1+β2(1-2N¯*-ϵ2))-αβ2ϵ2(1+β1N¯*)2](1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)3)-1,A3=((1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(2N¯*-1)+αϵ2-m1f2](1+β1N¯*)2)×(h1(1+β1N¯*)2)-1.Proof.
Firstly, the Jacobian matrix of the model system (12) around M*(N¯*,P¯*,Z¯*) is derived as follows: (19)J=(1-2N¯*-αN¯*(1+β1N¯*)2-αN¯*1+β1N¯*+f1f2m1P¯*(1+β1N¯*)2m1N¯*1+β1N¯*-h1Z¯*(1+β2P¯*)2-ϵ1-E1-h1P¯*1+β2P¯*0m2Z¯*(1+β2P¯*)2m2P¯*1+β2P¯*-ϵ2).
By virtue of (13) and (19), the characteristic equation of model system (12) around the interior equilibrium M* is as follows:
(20)λ3+A1λ2+A2λ+A3=0,
where Ai, i=1,2,3 are defined as follows:
(21)A1=((1+β1N¯*)2[ϵ2β2(ϵ1+E1)+h1(2N¯*-1)](m2-ϵ2β2)×(1+β1N¯*)2-m1ϵ2β2(m2-ϵ2β2)N¯*×(1+β1N¯*)+αh1ϵ2(1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)2)-1,A2=((1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(h1+β2(1-2N¯*-ϵ2))-αβ2ϵ2(1+β1N¯*)2](1+β1N¯*)2)×(h1(m2-ϵ2β2)(1+β1N¯*)3)-1,A3=((1+β1N¯*)2ϵ2[mN¯*-(ϵ1+E1)(1+β1N¯*)]×[(m2-ϵ2β2)(1+β1N¯*)2×(2N¯*-1)+αϵ2-m1f2(1+β1N¯*)2](1+β1N¯*)2)×(h1(1+β1N¯*)2)-1.
Based on Routh-Hurwitz criterion [3], model system (12) is locally asymptotically stable around M* provided that A1>0, A3>0, A1A2-A3>0.By denoting the quantitiesAi(h1), i=1,2,3 as smooth functions of the parameter h1 (predation rate), local stability switch and dynamical behavior is discussed due to the variation of predation rate h1.Theorem 4.
When the predation rateh1 crosses a critical value h1*, model system (12) enters into Hopf bifurcation around the interior equilibrium M* if the following conditions hold:(i)
A1(h1*)>0,A3(h1*)>0;(ii)
A1(h1*)A2(h1*)-A3(h1*)=0;(iii)
d[A1(h1)A2(h1)-A3(h1)]/dh1∣h1=h1*<0.Proof.
By choosing predation rateh1 as bifurcation parameter, we will discuss local stability switch due to variation of h1; we can see that if there exists a critical value h1* such that A1(h1*)>0, A3(h1*)>0, A1(h1*)A2(h1*)-A3(h1*)=0, and d[A1(h1)A2(h1)-A3(h1)]/dh1|h1=h1*<0. It is easy to show that A2(h1*)>0 based on conditions (i) and (ii). For the Hopf bifurcation to occur at h1=h1*, the characteristic equation takes the following form:
(22)(ρ2(h1*)+A2(h1*))(ρ(h1*)+A1(h1*))=0,
which has three roots ρ1(h1*)=iA2(h1*), ρ2=-iA2(h1*), and ρ3=-A1(h1*)<0.
Furthermore, in order to show if Hopf bifurcation occurs ath1=h1*, we need to verify the transversality condition [24]:
(23)[dRe(ρ(h1))dh1]|h1=h1*≠0.
Assuming all the roots of (22) are in general of the following form:
(24)ρ1(h1)=μ(h1)+iν(h1),ρ2(h1)=μ(h1)+iν(h1),ρ3(h1)=-A1(h1),
where μ(h1) and ν(h1) are smooth functions with respect to h1.
Now we will verify the transversality condition:(25)[dRe(ρi(h1))dh1]|h1=h1*≠0,i=1,2.
Substitutingρi(h1), i=1,2 into (22) and calculating the derivative, it follows that
(26)B1(h1)μ′(h1)-B2(h1)ν′(h1)+B3(h1)=0,B2(h1)μ′(h1)+B1(h1)ν′(h1)+B4(h1)=0,
where Bi(h1) are defined as follows:
(27)B1(h1)=3μ2(h1)+2A1(h1)μ(h1)+A2(h1)-3ν2(h1),B2(h1)=6μ(h1)ν(h1)+2A1(h1)ν(h1),B3(h1)=μ2(h1)A1′(h1)+A2′(h1)μ(h1)+A3′(h1)-A1′(h1)ν2(h1),B4(h1)=2μ(h1)ν(h1)A1′(h1)+A2′(h1)ν(h1),
where Ai′(h1), i=1,2,3 denotes the derivative of Ai(h1) with respect to h1.
It follows from (22) that μ(h1*)=0, ν(h1*)=A2(h1*). By further computation, it can be obtained that
(28)B1(h1*)=-2A2(h1*),B2(h1*)=2A1(h1*)A2(h1*),B3(h1*)=A3′(h1*)-A1′(h1*)A2(h1*),B4(h1*)=A2′(h1*)A2(h1*).
By solving (26), it can be obtained that
(29)[dRe(ρi(h1))dh1]|h1=h1*=μ′(h1)h1=h1*=-B2(h1*)B4(h1*)+B1(h1*)B3(h1*)B12(h1*)+B22(h1*)=A3′(h1*)-A1′(h1*)A2(h1*)-A1(h1*)A2′(h1*)2[A12(h1*)+A2(h1*)].
If condition (iii) holds, then(30)[dRe(ρi(h1))dh1]|h1=h1*>0,
which implies that transversality condition holds and Hopf bifurcation occurs at h1=h1*.Theorem 5.
By choosingα as a bifurcation parameter, model system (12) undergoes transcritical bifurcation around M* and bifurcation value is α=(m1f2-(m2-ϵ2β2)(1+β1N¯*)2(2N¯*-1))/ϵ2.Proof.
Whenα=(m1f2-(m2-ϵ2β2)(1+β1N¯*)2(2N¯*-1))/ϵ2, it follows from simple computation that (20) has a geometrically simple zero eigenvalue with left eigenvector vl=(h1P¯*/(1+β2P¯*),f2,1) and right eigenvector vr=(1,0,1)T, and
(31)vlDUDαH(U(t))vr|M*=(h1P¯*1+β2P¯*,f2,1)×(β1N¯*(1+2β1N¯*)-1(1+β1N¯*)400000000)(101)=h1P¯*[β1N¯*(1+2β1N¯*)-1](1+β2P¯*)(1+β1N¯*)4,(32)vlDU2H(U(t))(vr,vr,vr)|M*=vl∑i=13‍(eivlTDU(DUHi)Tvr)|M*=α(β1N¯*-1)-2f2m1β1P¯*(1+β1N¯*)3-2,
where H(U(t)), U(t) have been defined in Theorem 1 of this paper and ei, i=1,2,3, are unit vectors. It follows from [25] and the above computation that model system (12) undergoes transcritical bifurcation around M*.


## 4.2.  Numerical Experiments for Model System without Seasonality

In order to facilitate the following numerical experiments, some preliminaries about Poincaré surface-of-section technique and Lyapunov exponents are introduced in Remarks6 and 7.Remark 6.
A traditional approach to gain preliminary insight into the properties of dynamical system is to carry out a one-dimensional bifurcation analysis. One-dimensional bifurcation diagrams of Poincaré maps present information about the dependence of the dynamics on a certain parameter. The analysis is expected to reveal the type of attractor to which the dynamics will ultimately settle down after passing the initial transient phase and within which the trajectory will then remain forever. On a Poincaré surface of section, the dynamical behavior can be described by a discrete map whose phase-space dimension is less than that of the original continuous flow. Chaotic flows can be understood based on concepts that are convenient for maps such as unstable orbits (see [26–28] and references therein).Case 1.
A finite number of points corresponds to a periodic solution; that is, one point corresponds to a solution of period equal to that of the forcing term, namely,2π/ω (period-one); and n points correspond to a solution (subharmonic) of period 2nπ/ω (period n or a subharmonic of order 1/n).Case 2.
A closed curve corresponds to a quasiperiodic solution, that is, a solution consisting of two incommensurate frequencies or, equivalently, having a trajectory that is dense on tours.Case 3.
A collection of points that is fractal corresponds to chaos, namely, a stranger in phase space; a collection of points that form a cloud that is disorganized, partially organized, or fuzzy may (or may not ) correspond to chaotic attractors.Remark 7.
In a given embedding dimension, the Lyapunov exponent is a measure of the speeds at which initially nearby trajectories of the system diverge. There is a Lyapunov exponent for each dimension of the process, which together constitutes the Lyapunov spectrum for the dynamical system. The Lyapunov exponent is related to the predictability of the system, and the largest Lyapunov exponent of a stable system does not exceed zero. However, a chaotic system has at least one positive Lyapunov exponent. A bounded system with a positive Lyapunov exponent is one operational definition of chaotic behaviors, which presents a quantitative measure of the average rate of separation of nearby trajectories on the attractor. Over the years, a number of methods have been introduced for computation of Lyapunov exponents (see [26, 28]). A positive Lyapunov exponent implies that a chaotic process displays long term unpredictability, with the output being sensitively dependent on the initial conditions. Even slightly different initial values can lead to vastly different system outputs. Furthermore, another characteristic for chaos is fractal dimension, which can be computed with the commutated Lyapunov exponents (see Section  9.4 of reference [29]).In the absence of seasonality, we will perform numerical experiments to observe the dynamics of the model system (2) with the following hypothetical set of parameter values, which are given as follows: r0=1, K=1, a=1.63, r1=0.01, r2=0.001, c1=0.95, b1=0.33, d1=0.4, c2=0.9, d2=0.01, b2=0.5, q=0.1, and E=0.15. In the following part, h1 (predation rate) is chosen as bifurcation parameter. By virtue of scaling introduced in model system (12), values of parameters are obtained as follows: α=4.94, β1=3.03, m1=4.69, f1=0.01, f2=0.001, m2=0.9h1, ϵ1=0.4, ϵ2=0.01, β2=2, and E1=0.015. A bifurcation diagram of Poincaré section is plotted to investigate the dynamic effect of predation rate on model system (12), which can be shown in Figure 1. In order to understand Hopf bifurcation occurring at the critical value h1*=0.0711, bifurcation diagram with (0.07,0.075) is plotted in Figure 2.Figure 1
Bifurcation diagram of Poincaré section for the nutrient and predation rateh1 is varied in [0.068,0.1].Figure 2
Bifurcation diagram of Poincaré section for the nutrient and predation rateh1 is varied in [0.07,0.075].According to (15) and (16), it can be obtained that h1>0.068 guarantees the existence of interior equilibrium. Based on Figures 1 and 2, it is observed that model system (12) is stable around the interior equilibrium when h1 varies within (0.068, 0.0711). It should be noted that h1=0.07 is randomly selected within (0.068,0.0711); the dynamical response and phase portrait of model system (12) with h1=0.07 is plotted in Figures 3 and 4, respectively, which shows that model system (12) is locally asymptotically stable around interior equilibrium M*(0.7022,0.1999,13.9114).Figure 3
Dynamical response of model system (12) with predation rate h1=0.07.Figure 4
Phase portrait of model system (12) with predation rate h1=0.07.Ash1 increases across the critical value h1*=0.0711; model system enters into Hopf bifurcation. Dynamical responses of model system (12) with h1*=0.0711 is plotted in Figure 5, and a limit cycle corresponding to periodic solution is observed in Figure 6. As h1 further increases, it leads to a chaotic region at h1∈(0.0711,0.1). The dynamical response and phase portrait of model system with h1=0.08 is plotted in Figures 7 and 8, respectively. It should be noted that the sensitive dependence of future dynamics on the current state, the signature of chaotic behavior, is apparent from the fact that all the trajectories in the handle of teacup are very close together. Based on Theorem 4, model system undergoes Hopf bifurcation at critical bifurcation value h1*=0.0711. With the increase of bifurcation parameter h1, it can observed that model system enters into period doubling from limit cycle oscillations based on Figure 2, and model system finally settles down into chaotic behavior from period doubling. Hence, it can be concluded that model system admits that the transition to chaos is realized through the cascade of period doubling bifurcations.Figure 5
Dynamical response of model system (12) with predation rate h1=0.0711.Figure 6
A limit cycle corresponding to periodic solution in Figure5.Figure 7
Dynamical response of model system (12) with predation rate h1=0.08.Figure 8
Phase portrait of model system (12) with predation rate h1=0.08.Furthermore, chaotic attractor of model system (12) with h1=0.085,0.09,0.095, and 0.1 are plotted in Figure 9, and the corresponding Poincaré section of Figure 9 is plotted in Figure 10, which coincides with Case 3 introduced in Remark 6 implying the occurrence of chaotic attractor. Lyapunov exponents (denoted by Li, i=1,2,3) with time and the corresponding Lyapunov dimensions can be computed, which are utilized to show the existence of chaotic attractor. Table 1 shows Lyapunov exponents and corresponding Lyapunov dimensions of model system (12) with h1=0.08,0.085,0.09,0.095 and 0.1.Table 1
Lyapunove exponents and corresponding Lyapunov dimensions of model system (12) with h1=0.08, 0.085, 0.09, 0.095, and 0.1.
h
1
0.08
0.085
0.09
0.095
0.1
L
1
0.13761
0.1543
0.1818
0.1913
0.2044
L
2
0
0
0
0
0
L
3
−1.1265
−1.5804
−1.5803
−1.5893
−1.5710
D
2.1222
2.0976
2.115
2.1204
2.1301Chaotic attractor of model system (12) with predation rate h1=0.085,0.09,0.095,0.1.
(a)
(b)
(c)
(d)Corresponding Poincaré section of model system (12) with h1=0.085,0.09,0.095,0.1. Poincaré points 5000–15000 are plotted.
(a)
(b)
(c)
(d)It follows from the Table1 that there is a positive Lyapunov exponent for the model system (12) with h1=0.08,0.085,0.09,0.095, and 0.1, respectively. Furthermore, it follows from the Table 1 that Lyapunov dimensions for model system (12) with h1=0.08,0.085,0.09,0.095, and 0.1 are all fractals. According to the characterization of chaotic attractors introduced in Remark 7, a bounded system with a positive Lyapunov exponent is one operational definition of chaotic behavior, and another characteristic for chaos is fractal dimension.Oscillatory population may be driven to extinction in presence of environmental stochasticity when the population density is very low [30]. Based on the above numerical experiments, chaotic dynamics for increasing the predation rate has been discussed. In the following part, role of harvesting in controlling chaotic dynamics will be investigated. Now we fixed the predation rate h1=0.08 in chaotic region and the harvest effort E1 is chosen as bifurcation parameter. Figure 11 shows the bifurcation diagram of Poincaré section for the nutrient and harvest effort E1 is varied in [0,0.15]. It should be noted that harvest effort E1 should be controlled within [0,0.15] based on (15) and (16), which guarantees the existence of interior equilibrium. From Figure 11, it is observed that chaotic behavior remains for E1 that varies from [0,0.1381). With the increase of E1, chaos enter into limit cycle oscillation at E1=0.1381 and finally the model system (12) enters into stable focus from limit cycle.Figure 11
Bifurcation diagram of Poincaré section for the nutrient and harvest effortE1 is varied in [0,0.15].


## 4.3.  Numerical Experiments for Model System with Seasonality

With the help of MATLAB, numerical experiments with hypothetical set of parameters are performed to understand the dynamical effect of seasonal intrinsic growth rate of nutrient biomass on model system (2). The values of parameters are given as follows: r0=1, K=1, a=1.63, r1=0.01, r2=0.001, c1=0.95, b1=0.33, d1=0.4, c2=0.9, d2=0.01, b2=0.5, h=0.04, q=0.1, E=0.15, and ω=2. By choosing Γ as a bifurcation parameter, bifurcation diagram of Poincaré section for the nutrient and Γ is varied in [0,1], which can be found in Figure 12. It can be observed that chaotic behavior enters into periodic attractor when bifurcation parameter Γ increases through a critical value Γ=0.8. Model system (2) admits periodic attractor Γ within the range [0.8,1].Figure 12
Bifurcation diagram of Poincaré section for the nutrient andΓ is varied in [0,1].It should be noted thatΓ=0.81 and Γ=0.95 are randomly selected from [0.8,1]. Dynamical responses of model system with Γ=0.81 and Γ=0.95 are plotted in Figures 13 and 14, respectively. Corresponding Poincaré section of model system (2) with Γ=0.81 is plotted in Figure 15(a), which coincides with Case 1 introduced in Remark 6 implying the existence of periodic attractor; and corresponding Poincaré section of model system (2) with Γ=0.95 is plotted in Figure 15(b), which coincides with Case I introduced in Remark 6 implying the existence of periodic attractor. Poincaré points 5000–15000 are plotted. Based on the above numerical experiments, chaotic behavior can be controlled by the enhancement of amplitude of seasonal intrinsic growth of nutrient biomass; such type of mechanism reduces chaotic behavior to periodic attractor.Figure 13
Dynamical response of model system (2) with Γ=0.81.Figure 14
Dynamical response of model system (2) with Γ=0.95.Corresponding Poincaré section of model system (2) with Γ=0.81 is plotted in (a), and corresponding Poincaré section of model system (2) with Γ=0.95 is plotted in (b). Poincaré points 5000–15000 are plotted.
(a)
(b)Remark 8.
It can be observed that the dynamical responses fluctuate with complex structures when the model system admits chaotic behavior. Especially, some dynamical response even comes close to zero within some time interval, which biologically reflects harmful phytoplankton bloom. Consequently, there has been considerable scientific attention towards harmful phytoplankton bloom and its associate control strategy. Based on the theoretical analysis, it is found that high nutrient levels and favorable conditions play a key role in rapid or massive growth of phytoplankton bloom; low nutrient concentration, high predation pressure from zooplankton, and other unfavorable conditions limit phytoplankton growth, which leads to oscillations or recurring bloom in the nutrient-phytoplankton-zooplankton marine ecosystem.Based on the numerical experiments shown in Figure11, harvesting has a stabilizing effect on model dynamics. Furthermore, numerical experiments show that chaotic behavior of the population due to enhancement of predation rate can be prevented by adopting a suitable harvesting policy and model system (12) eventually approaches to a stable steady state. If the harvest effort on phytoplankton population increases, the zooplankton population decreases due to lack of the phytoplankton population and nutrient absorption pressure decreases and in turn the nutrient within marine system increases. On the other hand, harvest effort should be constrained within certain interval, which may guarantee the existence of nutrient biomass, phytoplankton population, and zooplankton population within marine system. The theoretical findings obtained in the above section are of inspiration for regulating some constructive policies to control harmful phytoplankton bloom with the introduction of appropriate harvest effort.


## 5. Conclusion

It is well known that phytoplankton population and zooplankton population play a key role in large-scale global processes such as ocean-atmosphere dynamics and climate change. Plankton population constitutes the bottom level of the marine and terrestrial life. Recently, harmful algal blooms (HAB) are widely reported and have become a serious environmental problem worldwide as its serious social and economic consequence. Therefore, a better understanding of mechanisms that determine the plankton dynamics is of considerable interest in recent decades [4–13]. Considering the economic interest of exploitation of phytoplankton and seasonal variation of nutrient biomass, a dynamical model is proposed in this paper. It is utilized to discuss the interaction mechanism of nutrient-phytoplankton-zooplankton marine system, where harvest effort on phytoplankton population is considered as well as seasonal intrinsic rate of nutrient biomass. Conditions for positivity and boundedness of solutions are obtained in Theorems 1 and 2, respectively. By choosing predation rate as bifurcation parameter, local stability analysis of model system (12) is performed around the interior equilibrium and sufficient conditions for occurrence of Hopf bifurcation are obtained, which can be found in Theorems 3 and 4, respectively. Further attempts are made to understand the effect of harvest effort and seasonality on population dynamics of the model system. Numerical experiments reveal that harvest effort on phytoplankton population has a stabilizing effect on population dynamics, and the increase of seasonal intrinsic growth of nutrient biomass is responsible for controlling chaotic behavior.Compared with the previous related work [4–6, 8], interaction mechanism of phytoplankton and zooplankton is discussed in this paper, and the effect of predation from zooplankton and nutrient recycling within nutrient-phytoplankton-zooplankton ecosystem are also studied. Compared with the previous related work [14–17, 21], nutrient recycling within marine ecosystem and seasonality of nutrient biomass are considered and investigated in this paper. Furthermore, the stabilizing effect on model system that concentrates on control strategies for chaotic behavior is also discussed. Generally, nutrient-phytoplankton-zooplankton ecosystem is often deeply perturbed by human exploiting activities; some biological interpretations about the theoretical analysis obtained in this paper are also provided that can be found in Remark 8, which are beneficial to discussing interaction and coexistence mechanism of population within the harvested marine ecosystem. It makes the work studied in this paper have some new and positive features.



---

*Source: 521917-2014-03-30.xml*