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

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.


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 2 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][5][6][7][8][9][10][11][12][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][12][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 2 Mathematical Problems in Engineering 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][15][16][17], where complex dynamical behavior such as multiplicity of attractors, catastrophes, and deterministic chaos are found.However, seasonality discussed in the above work [14][15][16][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][15][16][17].Generally, ecological system is often deeply perturbed by human exploiting activities [3,[18][19][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 nutrientphytoplankton-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.

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  < 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 nutrientphytoplankton-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 nutrientphytoplankton-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;  1 and  2 denote the half saturation constant. 1 < 1 is conversion factor from nutrient biomass to phytoplankton population. 2 < 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 where all parameters share the same interpretations introduced in hypotheses (H1)-(H4); model system ( 2) is analyzed with the following initial conditions:

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 () = ((), (), ())  ∈ R 3 , and (()) are given as follows: ) . ( 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,  0 ) for some  0 > 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: 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

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-ofsection technique and computation of Lyapunov exponents.

Local Stability Analysis of Model System without Seasonality.
In the case of Γ = 0, 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 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 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  * > 0,  * > 0, 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   ,  = 1, 2, 3 are defined as follows: By virtue of ( 13) and (19), the characteristic equation of model system (12) around the interior equilibrium  * is as follows: where   ,  = 1, 2, 3 are defined as follows: Based on Routh-Hurwitz criterion [3], model system ( 12) is locally asymptotically stable around  * provided that By denoting the quantities   (ℎ 1 ),  = 1, 2, 3 as smooth functions of the parameter ℎ 1 (predation rate), local stability switch and dynamical behavior is discussed due to the variation of predation rate ℎ 1 .

Numerical Experiments for Model System without
Seasonality.In order to facilitate the following numerical experiments, some preliminaries about Poincaré surface-ofsection 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 onedimensional 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][27][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/ (periodone); 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]).
According to (15) and ( 16), it can be obtained that ℎ 1 > 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 ℎ 1 varies within (0.068, 0.0711).It should be noted that ℎ 1 = 0.07 is randomly selected within (0.068, 0.0711); the dynamical response and phase portrait of model system (12) with ℎ 1 = 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).As ℎ 1 increases across the critical value ℎ * 1 = 0.0711; model system enters into Hopf bifurcation.Dynamical responses of model system (12) with ℎ * 1 = 0.0711 is plotted in Figure 5, and a limit cycle corresponding to periodic solution is observed in Figure 6.As ℎ 1 further increases, it leads to a chaotic region at ℎ 1 ∈ (0.0711, 0.1).The dynamical response and phase portrait of model system with ℎ 1 = 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 ℎ * 1 = 0.0711.With the increase of bifurcation parameter ℎ 1 , 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 ℎ 1 = 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 ℎ 1 = 0.08, 0.085, 0.09, 0.095 and 0.1.1 that there is a positive Lyapunov exponent for the model system (12) with ℎ 1 = 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 ℎ 1 = 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.

It follows from the Table
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 ℎ 1 = 0.08 in chaotic region and the harvest effort  1 is chosen as bifurcation parameter.Figure 11 shows the bifurcation diagram of Poincaré section for the nutrient and harvest effort  1 is varied in [0, 0.15].It should be noted that harvest effort  1 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  1 that varies from [0, 0.1381).With the increase of  1 , chaos enter into limit cycle oscillation at  1 = 0.1381 and finally the model system (12) enters into stable focus from limit cycle.
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.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.

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][5][6][7][8][9][10][11][12][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][5][6]8], interaction mechanism of phytoplankton and zooplankton is discussed in this paper, and the effect of predation from zooplankton and nutrient recycling within nutrientphytoplankton-zooplankton ecosystem are also studied.Compared with the previous related work [14][15][16][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.