#### 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 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: 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 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. is conversion factor from nutrient biomass to phytoplankton population. 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 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 . 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 . relates to absence of seasonality, and 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: where all parameters share the same interpretations introduced in hypotheses (H1)–(H4); model system (2) is analyzed with the following initial conditions:

#### 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
where , and are given as follows: Let be the positive octant in ; then is locally Lipschitz and satisfies the condition , .

Due to the lemma in [23] and Theorem 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 for some . 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:
**
where .*

*Proof. *Let ; it follows from model system (2) and (H1)–(H4) that

According to Theorem 1, it derives that
where .

Based on the above analysis, it can be obtained that

#### 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 , model system (2) without seasonality takes the following form:

In order to reduce number of parameters and facilitate dynamical analysis of model system, model system (10) is dimensionalized with the following scaling: Then model system (10) takes the following form: 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: where satisfies the following equation: where , , and .

Based on Routh-Hurwitz criterion [3], a simple condition that (14) has at least one positive roots is as follows:

Furthermore, according to the biological interpretations of interior equilibrium, phytoplankton population and zooplankton population all exist provided that , , which follows that

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:
**
where , are defined as follows:
*

*Proof. *Firstly, the Jacobian matrix of the model system (12) around is derived as follows:

By virtue of (13) and (19), the characteristic equation of model system (12) around the interior equilibrium is as follows:
where , are defined as follows:

Based on Routh-Hurwitz criterion [3], model system (12) is locally asymptotically stable around provided that , , .

By denoting the quantities , 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)*;*(ii)*;*(iii)*. *

*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 , , , and . It is easy to show that based on conditions (i) and (ii). For the Hopf bifurcation to occur at , the characteristic equation takes the following form:
which has three roots , , and .

Furthermore, in order to show if Hopf bifurcation occurs at , we need to verify the transversality condition [24]:

Assuming all the roots of (22) are in general of the following form:
where and are smooth functions with respect to .

Now we will verify the transversality condition:

Substituting , into (22) and calculating the derivative, it follows that
where are defined as follows:
where , denotes the derivative of with respect to .

It follows from (22) that , . By further computation, it can be obtained that

By solving (26), it can be obtained that

If condition (iii) holds, then
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 .*

*Proof. *When , it follows from simple computation that (20) has a geometrically simple zero eigenvalue with left eigenvector and right eigenvector , and
where , have been defined in Theorem 1 of this paper and , , 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, (period-one); and points correspond to a solution (subharmonic) of period (period or a subharmonic of order ).

*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: , , , , , , , , , , , , and . 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: , , , , , , , , , and . 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 , bifurcation diagram with is plotted in Figure 2.

According to (15) and (16), it can be obtained that 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 is randomly selected within ; the dynamical response and phase portrait of model system (12) with is plotted in Figures 3 and 4, respectively, which shows that model system (12) is locally asymptotically stable around interior equilibrium .

As increases across the critical value ; model system enters into Hopf bifurcation. Dynamical responses of model system (12) with 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 . The dynamical response and phase portrait of model system with 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 . 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.

Furthermore, chaotic attractor of model system (12) with , and 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 , ) 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 and .

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

It follows from the Table 1 that there is a positive Lyapunov exponent for the model system (12) with , and , respectively. Furthermore, it follows from the Table 1 that Lyapunov dimensions for model system (12) with , and 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 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 . It should be noted that harvest effort should be controlled within 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 . With the increase of , chaos enter into limit cycle oscillation at and finally the model system (12) enters into stable focus from limit cycle.

##### 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: , , , , , , , , , , , , , , and . By choosing as a bifurcation parameter, bifurcation diagram of Poincaré section for the nutrient and is varied in , 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 . Model system (2) admits periodic attractor within the range .

It should be noted that and are randomly selected from . Dynamical responses of model system with and are plotted in Figures 13 and 14, respectively. Corresponding Poincaré section of model system (2) with 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 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.

**(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 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.