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

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.

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 transporting CO₂ 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−𝐾−𝑎𝑋(𝑡)𝑌(𝑡)(𝑡),(𝑡)𝑏+𝑋𝑑𝑌=(𝑡)𝑌(𝑡)𝑑𝑡𝑐𝑎𝑋(𝑡)(𝑡),𝑏+𝑋−𝑑𝑌 (1) where 𝑋(𝑡) and 𝑌(𝑡) represent density of prey population and predator population, respectively. 𝑟 denotes the intrinsic growth rate of the prey population and the carrying capacity of the prey population is 𝐾. 𝑎 is the searching efficiency or the predation rate of the predator population; 𝑏 denotes the half saturation constant and 𝑐<1 is the conversion factor. 𝑑 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 𝑁(𝑡), 𝑃(𝑡), and 𝑍(𝑡) denote nutrient biomass, phytoplankton population density, and zooplankton population density with respect to time 𝑡, respectively. Let 𝐾 represent the carrying capacity of the nutrient biomass within marine ecosystem. 𝑎 is the searching efficiency of phytoplankton population, and ℎ is the predation rate; 𝑏₁ and 𝑏₂ denote the half saturation constant. 𝑐₁<1 is conversion factor from nutrient biomass to phytoplankton population. 𝑐₂<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 scalar 𝐸≥0 denotes the harvesting effort to phytoplankton population, constant 𝑞 is the catchability coefficient of phytoplankton population, and the harvesting term 𝑞𝐸𝑁(𝑡) 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 𝑑₁𝑃(𝑡) and 𝑑₂𝑍(𝑡) while through recycling (nutrient biomass decomposition). 𝑟₁ and 𝑟₂ represent amount of phytoplankton and zooplankton population that converts back into nutrient biomass, respectively. It follows from the biological interpretations of parameters that 𝑑₁>𝑟₁ and 𝑑₂>𝑟₂.(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 biomass 𝑟 is assumed to be a periodically varying function of time and is considered as 𝑟=𝑟₀(1+Γsin(𝜔𝑡)). It should be noted that the average value of 𝑟 is assumed to be 𝑟₀ and Γ represents the degree of seasonality; 𝑟₀Γ is the magnitude of the perturbation and 𝜔 is the angular frequency of the fluctuations caused by seasonality. Since 𝑟 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: 𝑑𝑁(𝑡)𝑑𝑡=𝑟₀((𝑁(𝑡)𝑁(𝑡)1+Γsin𝜔𝑡))1−𝐾−𝑎𝑁(𝑡)𝑃(𝑡)𝑏₁+𝑁(𝑡)+𝑟₁𝑃(𝑡)+𝑟₂𝑍(𝑡),𝑑𝑃(𝑡)=𝑐𝑑𝑡¹𝑎𝑁(𝑡)𝑃(𝑡)𝑏₁(𝑡)−(𝑡)𝑍(𝑡)+𝑁ℎ𝑃𝑏₂+𝑃(𝑡)−𝑑₁𝑃(𝑡)(𝑡),(𝑡)−𝑞𝐸𝑃𝑑𝑍=𝑐𝑑𝑡²ℎ𝑃(𝑡)𝑍(𝑡)𝑏₂+𝑃(𝑡)−𝑑₂𝑍(𝑡), (2) where all parameters share the same interpretations introduced in hypotheses (H1)–(H4); model system (2) is analyzed with the following initial conditions: 𝑁(0)>0,𝑃(0)>0,𝑍(0)>0. (3)

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) where 𝑈(𝑡)=(𝑁(𝑡),𝑃(𝑡),𝑍(𝑡))^𝑇∈ℝ³, and 𝐻(𝑈(𝑡)) are given as follows: 𝐻(𝑈(𝑡=𝐻))₁(𝑈(𝑡𝐻))₂(𝑈(𝑡𝐻))₃(𝑈(𝑡=𝑟))⁰((𝑁(𝑡)𝑁(𝑡)1+Γsin𝜔𝑡))1−𝐾−𝑎𝑁(𝑡)𝑃(𝑡)𝑏₁+𝑁(𝑡)+𝑟₁𝑃(𝑡)+𝑟₂𝑍(𝑡)𝑐𝑎𝑁(𝑡)𝑃(𝑡)𝑏₁(𝑡)−(𝑡)𝑍(𝑡)+𝑁ℎ𝑃𝑏₂+𝑃(𝑡)−𝑑₁𝑃(𝑡)(𝑡)(𝑡)𝑍(𝑡)−𝑞𝐸𝑃ℎ𝑃𝑏₂+𝑃(𝑡)−𝑑₂𝑍(𝑡). (5) Let ℝ³₊=[0,∞)³ be the positive octant in ℝ³; then 𝐿∶ℝ₊³⁺¹→ℝ³ is locally Lipschitz and satisfies the condition 𝐻_𝑖(𝑈(𝑡))|_𝑈∈ℝ³≥0, 𝑖=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,𝐴₀) for some 𝐴₀>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: 𝑁(𝑡)(𝑡)(𝑡(𝑡)(𝑡)(𝑡)≤𝐾Ω=(,𝑃,𝑍))∣𝑁+𝑃+𝑍2𝑟₀+𝑒²8𝑒𝑟₀, (6) where 𝑒=min{𝑞𝐸+𝑑₁−𝑟₁,𝑑₂−𝑟₂}.

Proof. Let 𝑊(𝑡)=𝑁(𝑡)+𝑃(𝑡)+𝑍(𝑡); it follows from model system (2) and (H1)–(H4) that d𝑊(𝑡)d𝑡≤2𝑟₀𝑁(𝑡)𝑁(𝑡)1−𝐾−𝑎1−𝑐₁𝑁(𝑡)𝑃(𝑡)𝑏₁(𝑡)−ℎ+𝑁1−𝑐₂𝑃(𝑡)𝑍(𝑡)𝑏₂(𝑡)−𝑑+𝑃¹−𝑟₁𝑃(𝑡)−𝑑+𝑞𝐸₂−𝑟₂𝑍(𝑡). (7)
According to Theorem 1, it derives that d𝑊(𝑡)≤d𝑡2𝑟₀+𝑒𝑁(𝑡)−2𝑟₀𝑁²(𝑡)𝐾(𝑡),≤𝐾−𝑒𝑊2𝑟₀+𝑒²8𝑟₀−𝑒𝑊(𝑡), (8) where 𝑒=min{𝑞𝐸+𝑑₁−𝑟₁,𝑑₂−𝑟₂}.
Based on the above analysis, it can be obtained that 𝑊(𝑡)≤𝐾2𝑟₀+𝑒²8𝑒𝑟₀. (9)

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: 𝑑𝑁(𝑡)𝑑𝑡=𝑟₀𝑁(𝑡)𝑁(𝑡)1−𝐾−𝑎𝑁(𝑡)𝑃(𝑡)𝑏₁+𝑁(𝑡)+𝑟₁𝑃(𝑡)+𝑟₂𝑍(𝑡),𝑑𝑃(𝑡)=(𝑡)𝑃(𝑡)𝑑𝑡𝑐𝑎𝑁𝑏₁(𝑡)−(𝑡)𝑍(𝑡)+𝑁ℎ𝑃𝑏₂+𝑃(𝑡)−𝑑₁𝑃(𝑡)(𝑡),(𝑡)−𝑞𝐸𝑃𝑑𝑍=(𝑡)𝑍(𝑡)𝑑𝑡ℎ𝑃𝑏₂+𝑃(𝑡)−𝑑₂𝑍(𝑡). (10)

In order to reduce number of parameters and facilitate dynamical analysis of model system, model system (10) is dimensionalized with the following scaling: 𝑁𝑁=𝐾,𝑃𝑃=𝐾,𝑍𝑍=𝐾,𝑇=𝑟₀𝑡. (11) Then model system (10) takes the following form: 𝑑𝑁(𝑇)=𝑑𝑇𝑁(𝑇)1−𝑁(𝑇)−𝛼𝑁(𝑇)𝑃(𝑇)1+𝛽₁𝑁(𝑇)+𝑓₁𝑃(𝑇)+𝑓₂𝑍(𝑇),𝑑𝑃(𝑇)=𝑚𝑑𝑇¹𝑁(𝑇)𝑃(𝑇)1+𝛽₁𝑁(𝑇)−ℎ₁𝑃(𝑇)𝑍(𝑇)1+𝛽₂𝑃(𝑇)−𝜖₁𝑃(𝑇)−𝐸₁𝑃(𝑇),𝑑𝑍(𝑇)=𝑚𝑑𝑇²𝑃(𝑇)𝑍(𝑇)1+𝛽₂𝑃(𝑇)−𝜖₂𝑍(𝑇), (12) where 𝛼=𝑎𝐾/𝑏₁𝑟₀, 𝛽₁=𝐾/𝑏₁, 𝑚₁=𝑎𝑐₁𝐾/𝑏₁𝑟₀, 𝑓₁=𝑟₁/𝑟₀, 𝑓₂=𝑟₂/𝑟₀, ℎ₁=ℎ𝐾/𝑏₂𝑟₀, 𝑚₂=𝑐₂ℎ𝐾/𝑏₂𝑟₀, 𝜖₁=𝑑₁/𝑟₀, 𝜖₂=𝑑₂/𝑟₀, 𝛽₂=𝐾/𝑏₂, and 𝐸₁=𝑞𝐸/𝑟₀.

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 𝑀^∗(𝑁^∗,𝑃^∗,𝑍^∗) is as follows: 𝑃^∗=𝜖₂𝑚₂−𝜖₂𝛽₂,𝑍^∗=𝑚₂𝑚₁𝑁^∗−𝜖₁+𝐸₁1+𝛽₁𝑁^∗ℎ₁𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗, (13) where 𝑁^∗ satisfies the following equation: 𝑁³+𝐶₁𝑁²+𝐶₂𝑁+𝐶₃=0, (14) where 𝐶₁=(1−𝛽₁)/𝛽₁, 𝐶₂=(𝑚₂𝑓₂[𝛽₁(𝜖₁+𝐸₁)−𝑚₁]+ℎ₁[𝜖₂(𝛽₂+𝛼−𝑓₁𝛽₁)−𝑚₂])/ℎ₁𝛽₁(𝑚₁−𝜖₂𝛽₂), and 𝐶₃=(𝑓₂𝑚₂(𝜖₁+𝐸₁)−𝜖₂𝑓₁ℎ₁)/ℎ₁𝛽₁(𝑚₂−𝜖₂𝛽₂).

Based on Routh-Hurwitz criterion [3], a simple condition that (14) has at least one positive roots is as follows: 𝐶₁𝐶₂−𝐶₃≤0. (15)

Furthermore, according to the biological interpretations of interior equilibrium, phytoplankton population and zooplankton population all exist provided that 𝑃^∗>0, 𝑍^∗>0, which follows that 𝑚₂−𝜖₂𝛽₂>0,𝑁^∗>𝜖₁+𝐸₁𝑚₁−𝛽₁𝜖₁+𝐸₁. (16)

In Theorem 3, local stability analysis around the interior equilibrium 𝑀^∗ is performed.

Theorem 3. Model system (12) is locally asymptotically stable around the interior equilibrium 𝑀^∗ if the following conditions hold: 𝐴₁>0,𝐴₃>0,𝐴₁𝐴₂−𝐴₃>0, (17) where 𝐴_𝑖, 𝑖=1,2,3 are defined as follows: 𝐴₁=𝜖₂𝛽₂𝜖₁+𝐸₁+ℎ₁2𝑁^∗−1𝑚₂−𝜖₂𝛽₂×1+𝛽₁𝑁^∗2−𝑚₁𝜖₂𝛽₂𝑚₂−𝜖₂𝛽₂𝑁^∗×1+𝛽₁𝑁^∗+𝛼ℎ₁𝜖₂×ℎ₁𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗2−1,𝐴₂=𝜖₂𝑚𝑁^∗−𝜖₁+𝐸₁1+𝛽₁𝑁^∗×𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗2×ℎ₁+𝛽₂1−2𝑁^∗−𝜖₂−𝛼𝛽₂𝜖₂×ℎ₁𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗3−1,𝐴₃=𝜖₂𝑚𝑁^∗−𝜖₁+𝐸₁1+𝛽₁𝑁^∗×𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗2×2𝑁^∗−1+𝛼𝜖₂−𝑚₁𝑓₂×ℎ₁1+𝛽₁𝑁^∗2−1. (18)

Proof. Firstly, the Jacobian matrix of the model system (12) around 𝑀^∗(𝑁^∗,𝑃^∗,𝑍^∗) is derived as follows: 𝐽=1−2𝑁^∗−𝛼𝑁^∗1+𝛽₁𝑁^∗2−𝛼𝑁^∗1+𝛽₁𝑁^∗+𝑓₁𝑓₂𝑚₁𝑃^∗1+𝛽₁𝑁^∗2𝑚₁𝑁^∗1+𝛽₁𝑁^∗−ℎ₁𝑍^∗1+𝛽₂𝑃^∗2−𝜖₁−𝐸₁−ℎ₁𝑃^∗1+𝛽₂𝑃^∗0𝑚₂𝑍^∗1+𝛽₂𝑃^∗2𝑚₂𝑃^∗1+𝛽₂𝑃^∗−𝜖₂. (19)
By virtue of (13) and (19), the characteristic equation of model system (12) around the interior equilibrium 𝑀^∗ is as follows: 𝜆³+𝐴₁𝜆²+𝐴₂𝜆+𝐴₃=0, (20) where 𝐴_𝑖, 𝑖=1,2,3 are defined as follows: 𝐴₁=𝜖₂𝛽₂𝜖₁+𝐸₁+ℎ₁2𝑁^∗−1𝑚₂−𝜖₂𝛽₂×1+𝛽₁𝑁^∗2−𝑚₁𝜖₂𝛽₂𝑚₂−𝜖₂𝛽₂𝑁^∗×1+𝛽₁𝑁^∗+𝛼ℎ₁𝜖₂×ℎ₁(𝑚₂−𝜖₂𝛽₂)1+𝛽₁𝑁^∗2−1,𝐴₂=𝜖₂𝑚𝑁^∗−𝜖₁+𝐸₁1+𝛽₁𝑁^∗×𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗2×ℎ₁+𝛽₂1−2𝑁^∗−𝜖₂−𝛼𝛽₂𝜖₂×ℎ₁𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗3−1,𝐴₃=𝜖₂𝑚𝑁^∗−𝜖₁+𝐸₁1+𝛽₁𝑁^∗×𝑚₂−𝜖₂𝛽₂1+𝛽₁𝑁^∗2×2𝑁^∗−1+𝛼𝜖₂−𝑚₁𝑓₂×ℎ₁1+𝛽₁𝑁^∗2−1. (21)
Based on Routh-Hurwitz criterion [3], model system (12) is locally asymptotically stable around 𝑀^∗ provided that 𝐴₁>0, 𝐴₃>0, 𝐴₁𝐴₂−𝐴₃>0.

By denoting the quantities 𝐴_𝑖(ℎ₁), 𝑖=1,2,3 as smooth functions of the parameter ℎ₁ (predation rate), local stability switch and dynamical behavior is discussed due to the variation of predation rate ℎ₁.

Theorem 4. When the predation rate ℎ₁ crosses a critical value ℎ^∗₁, model system (12) enters into Hopf bifurcation around the interior equilibrium 𝑀^∗ if the following conditions hold:(i)𝐴₁(ℎ^∗₁)>0,𝐴₃(ℎ^∗₁)>0;(ii)𝐴₁(ℎ^∗₁)𝐴₂(ℎ^∗₁)−𝐴₃(ℎ^∗₁)=0;(iii)d[𝐴₁(ℎ₁)𝐴₂(ℎ₁)−𝐴₃(ℎ₁)]/dℎ₁∣_ℎ1=ℎ^∗₁<0.

Proof. By choosing predation rate ℎ₁ as bifurcation parameter, we will discuss local stability switch due to variation of ℎ₁; we can see that if there exists a critical value ℎ^∗₁ such that 𝐴₁(ℎ^∗₁)>0, 𝐴₃(ℎ^∗₁)>0, 𝐴₁(ℎ^∗₁)𝐴₂(ℎ^∗₁)−𝐴₃(ℎ^∗₁)=0, and 𝑑[𝐴₁(ℎ₁)𝐴₂(ℎ₁)−𝐴₃(ℎ₁)]/𝑑ℎ₁|_ℎ1=ℎ^∗₁<0. It is easy to show that 𝐴₂(ℎ^∗₁)>0 based on conditions (i) and (ii). For the Hopf bifurcation to occur at ℎ₁=ℎ^∗₁, the characteristic equation takes the following form: 𝜌²ℎ^∗₁+𝐴₂ℎ^∗₁𝜌ℎ^∗₁+𝐴₁ℎ^∗₁=0, (22) which has three roots 𝜌₁(ℎ^∗₁)=𝑖𝐴₂(ℎ^∗₁), 𝜌₂=−𝑖𝐴₂(ℎ^∗₁), and 𝜌₃=−𝐴₁(ℎ^∗₁)<0.
Furthermore, in order to show if Hopf bifurcation occurs at ℎ₁=ℎ^∗₁, we need to verify the transversality condition [24]: 𝑑Re𝜌ℎ₁𝑑ℎ_{1ℎ1=ℎ∗1}≠0. (23)
Assuming all the roots of (22) are in general of the following form: 𝜌₁ℎ₁=𝜇ℎ₁+𝑖𝜈ℎ₁,𝜌₂ℎ₁=𝜇ℎ₁+𝑖𝜈ℎ₁,𝜌₃ℎ₁=−𝐴₁ℎ₁, (24) where 𝜇(ℎ₁) and 𝜈(ℎ₁) are smooth functions with respect to ℎ₁.
Now we will verify the transversality condition: 𝑑Re𝜌_𝑖ℎ₁𝑑ℎ_{1ℎ1=ℎ∗1}≠0,𝑖=1,2. (25)
Substituting 𝜌_𝑖(ℎ₁), 𝑖=1,2 into (22) and calculating the derivative, it follows that 𝐵₁ℎ₁𝜇^′ℎ₁−𝐵₂ℎ₁𝜈^′ℎ₁+𝐵₃ℎ₁𝐵=0,₂ℎ₁𝜇^′ℎ₁+𝐵₁ℎ₁𝜈^′ℎ₁+𝐵₄ℎ₁=0, (26) where 𝐵_𝑖(ℎ₁) are defined as follows: 𝐵₁ℎ₁=3𝜇²ℎ₁+2𝐴₁ℎ₁𝜇ℎ₁+𝐴₂ℎ₁−3𝜈²ℎ₁,𝐵₂ℎ₁=6𝜇ℎ₁𝜈ℎ₁+2𝐴₁ℎ₁𝜈ℎ₁,𝐵₃ℎ₁=𝜇²ℎ₁𝐴^′₁ℎ₁+𝐴^′₂ℎ₁𝜇ℎ₁+𝐴^′₃ℎ₁−𝐴^′₁ℎ₁𝜈²ℎ₁,𝐵₄ℎ₁=2𝜇ℎ₁𝜈ℎ₁𝐴^′₁ℎ₁+𝐴^′₂ℎ₁𝜈ℎ₁, (27) where 𝐴^′_𝑖(ℎ₁), 𝑖=1,2,3 denotes the derivative of 𝐴_𝑖(ℎ₁) with respect to ℎ₁.
It follows from (22) that 𝜇(ℎ^∗₁)=0, 𝜈(ℎ^∗₁)=𝐴₂(ℎ^∗₁). By further computation, it can be obtained that 𝐵₁ℎ^∗₁=−2𝐴₂ℎ^∗₁,𝐵₂ℎ^∗₁=2𝐴₁ℎ^∗₁𝐴₂ℎ^∗₁,𝐵₃ℎ^∗₁=𝐴^′₃ℎ^∗₁−𝐴^′₁ℎ^∗₁𝐴₂ℎ^∗₁,𝐵₄ℎ^∗₁=𝐴^′₂ℎ^∗₁𝐴₂ℎ^∗₁. (28)
By solving (26), it can be obtained that 𝑑Re𝜌_𝑖ℎ₁𝑑ℎ_{1ℎ1=ℎ∗1}=𝜇^′ℎ_1ℎ1=ℎ^∗₁𝐵=−²ℎ^∗₁𝐵₄ℎ^∗₁+𝐵₁ℎ^∗₁𝐵₃ℎ^∗₁𝐵²₁ℎ^∗₁+𝐵²₂ℎ^∗₁=𝐴^′₃ℎ^∗₁−𝐴^′₁ℎ^∗₁𝐴₂ℎ^∗₁−𝐴₁ℎ^∗₁𝐴^′₂ℎ^∗₁2𝐴²₁ℎ^∗₁+𝐴₂ℎ^∗₁. (29)
If condition (iii) holds, then 𝑑Re𝜌_𝑖ℎ₁𝑑ℎ_{1ℎ1=ℎ∗1}>0, (30) which implies that transversality condition holds and Hopf bifurcation occurs at ℎ₁=ℎ^∗₁.

Theorem 5. By choosing 𝛼 as a bifurcation parameter, model system (12) undergoes transcritical bifurcation around 𝑀^∗ and bifurcation value is 𝛼=(𝑚₁𝑓₂−(𝑚₂−𝜖₂𝛽₂)(1+𝛽₁𝑁^∗)²(2𝑁^∗−1))/𝜖₂.

Proof. When 𝛼=(𝑚₁𝑓₂−(𝑚₂−𝜖₂𝛽₂)(1+𝛽₁𝑁^∗)²(2𝑁^∗−1))/𝜖₂, it follows from simple computation that (20) has a geometrically simple zero eigenvalue with left eigenvector 𝑣_𝑙=(ℎ₁𝑃^∗/(1+𝛽₂𝑃^∗),𝑓₂,1) and right eigenvector 𝑣_𝑟=(1,0,1)^𝑇, and 𝑣_𝑙𝐷_𝑈𝐷_𝛼𝐻(𝑈(𝑡))𝑣_𝑟𝑀∗=ℎ₁𝑃^∗1+𝛽₂𝑃^∗,𝑓₂×𝛽,1₁𝑁^∗1+2𝛽₁𝑁^∗−11+𝛽₁𝑁^∗4101=ℎ00000000₁𝑃^∗𝛽₁𝑁^∗1+2𝛽₁𝑁^∗−11+𝛽₂𝑃^∗1+𝛽₁𝑁^∗4,𝑣_𝑙𝐷²_𝑈𝐻(𝑈(𝑡))𝑣_𝑟,𝑣_𝑟,𝑣_𝑟𝑀∗=𝑣_𝑙³_𝑖=1𝑒_𝑖𝑣^𝑇_𝑙𝐷_𝑈𝐷_𝑈𝐻_𝑖^𝑇𝑣_𝑟𝑀∗=𝛼𝛽₁𝑁^∗−1−2𝑓₂𝑚₁𝛽₁𝑃^∗1+𝛽₁𝑁^∗3−2, (31)(32) where 𝐻(𝑈(𝑡)), 𝑈(𝑡) have been defined in Theorem 1 of this paper and 𝑒_𝑖, 𝑖=1,2,3, are unit vectors. It follows from [25] and the above computation that model system (12) undergoes transcritical bifurcation around 𝑀^∗.

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 Remarks 6 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 𝑛 points correspond to a solution (subharmonic) of period 2𝑛𝜋/𝜔 (period 𝑛 or a subharmonic of order 1/𝑛).

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: 𝑟₀=1, 𝐾=1, 𝑎=1.63, 𝑟₁=0.01, 𝑟₂=0.001, 𝑐₁=0.95, 𝑏₁=0.33, 𝑑₁=0.4, 𝑐₂=0.9, 𝑑₂=0.01, 𝑏₂=0.5, 𝑞=0.1, and 𝐸=0.15. In the following part, ℎ₁ (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, 𝛽₁=3.03, 𝑚₁=4.69, 𝑓₁=0.01, 𝑓₂=0.001, 𝑚₂=0.9ℎ₁, 𝜖₁=0.4, 𝜖₂=0.01, 𝛽₂=2, and 𝐸₁=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 ℎ^∗₁=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 rate ℎ1 is varied in [0.068,0.1].

Figure 2 

Bifurcation diagram of Poincaré section for the nutrient and predation rate ℎ1 is varied in [0.07,0.075].

According to (15) and (16), it can be obtained that ℎ₁>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 ℎ₁ varies within (0.068, 0.0711). It should be noted that ℎ₁=0.07 is randomly selected within (0.068,0.0711); the dynamical response and phase portrait of model system (12) with ℎ₁=0.07 is plotted in Figures 3 and 4, respectively, which shows that model system (12) is locally asymptotically stable around interior equilibrium 𝑀^∗(0.7022,0.1999,13.9114).

Figure 3 

Dynamical response of model system (12) with predation rate ℎ1=0.07.

Figure 4 

Phase portrait of model system (12) with predation rate ℎ1=0.07.

As ℎ₁ increases across the critical value ℎ^∗₁=0.0711; model system enters into Hopf bifurcation. Dynamical responses of model system (12) with ℎ^∗₁=0.0711 is plotted in Figure 5, and a limit cycle corresponding to periodic solution is observed in Figure 6. As ℎ₁ further increases, it leads to a chaotic region at ℎ₁∈(0.0711,0.1). The dynamical response and phase portrait of model system with ℎ₁=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 ℎ^∗₁=0.0711. With the increase of bifurcation parameter ℎ₁, 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 ℎ1=0.0711.

Figure 6 

A limit cycle corresponding to periodic solution in Figure 5.

Figure 7 

Dynamical response of model system (12) with predation rate ℎ1=0.08.

Figure 8 

Phase portrait of model system (12) with predation rate ℎ1=0.08.

Furthermore, chaotic attractor of model system (12) with ℎ₁=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 𝐿_𝑖, 𝑖=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 ℎ₁=0.08,0.085,0.09,0.095 and 0.1.

Table 1 

Lyapunove exponents and corresponding Lyapunov dimensions of model system (12) with ℎ₁=0.08, 0.085, 0.09, 0.095, and 0.1.

(a)

(b)

(c)

(d)

Figure 9 

Chaotic attractor of model system (12) with predation rate ℎ1=0.085,0.09,0.095,0.1.

(a)

(b)

(c)

(d)

Figure 10 

Corresponding Poincaré section of model system (12) with ℎ1=0.085,0.09,0.095,0.1. Poincaré points 5000–15000 are plotted.

It follows from the Table 1 that there is a positive Lyapunov exponent for the model system (12) with ℎ₁=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 ℎ₁=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 ℎ₁=0.08 in chaotic region and the harvest effort 𝐸₁ is chosen as bifurcation parameter. Figure 11 shows the bifurcation diagram of Poincaré section for the nutrient and harvest effort 𝐸₁ is varied in [0,0.15]. It should be noted that harvest effort 𝐸₁ 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 𝐸₁ that varies from [0,0.1381). With the increase of 𝐸₁, chaos enter into limit cycle oscillation at 𝐸₁=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 effort 𝐸1 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: 𝑟₀=1, 𝐾=1, 𝑎=1.63, 𝑟₁=0.01, 𝑟₂=0.001, 𝑐₁=0.95, 𝑏₁=0.33, 𝑑₁=0.4, 𝑐₂=0.9, 𝑑₂=0.01, 𝑏₂=0.5, ℎ=0.04, 𝑞=0.1, 𝐸=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.

(a)

(b)

Figure 15 

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.

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 Figure 11, 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.

Conflict of Interests

All authors of this paper declare that there is no conflict of interests regarding the publication of this paper. The authors have no proprietary, financial, professional, or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, this paper.

Acknowledgments

The authors gratefully acknowledge anonymous reviewers, editors and Professor Lev Ryashko's comments. This work is supported by National Natural Science Foundation of China, Grant no. 61104003, Grant no. 61273008, and Grant no. 61104093; Research Foundation for Doctoral Program of Higher Education of Education Ministry, Grant no. 20110042120016; Hebei Province Natural Science Foundation, Grant no. F2011501023; Fundamental Research Funds for the Central Universities, Grant no. N120423009; Research Foundation for Science and Technology Pillar Program of Northeastern University at Qinhuangdao, Grant no. XNK201301. This work is supported by State Key Laboratory of Integrated Automation of Process Industry, Northeastern University, supported by Hong Kong Admission Scheme for Mainland Talents and Professionals, Hong Kong, Special Administrative Region.