Dynamics Induced by Delay in a Nutrient-Phytoplankton Model with Multiple Delays

1Zhejiang Provincial Key Laboratory for Water Environment and Marine Biological Resources Protection, Wenzhou University, Wenzhou, Zhejiang, 325035, China 2College of Life and Environmental Science, Wenzhou University, Wenzhou, Zhejiang, 325035, China 3College of Mathematics, Physics and Electronic Information Engineering, Wenzhou University, Wenzhou, Zhejiang, 325035, China 4Environmental Engineering Program, University of Northern British Columbia, Prince George, British Columbia, V2N 4Z9, Canada


Introduction
Some phytoplankton, for example, Cyanobacteria, can form dense and sometimes toxic blooms in freshwater and marine environments, which threaten ecological balance, drinking water, fisheries, and even human health [1].However, the mechanism, by which phytoplankton blooms occur, is currently not very clear, which contribute to the difficulty to prevent or mitigate the proliferation of phytoplankton blooms.These have stimulated lots of researches aiming to understand the growth mechanisms of phytoplankton.
In recent years, dynamics in phytoplankton growth have drawn increasing attention from experimental ecologists, as well as mathematical ecologists.Some results from experiments and field observations imply that many factors affecting the dynamics of phytoplankton growth are bound to exist, such as nutrient [2], light [3], temperature [4], iron supply [5], zooplankton [6].Especially, due to the effects of limiting factors including temperature, light, and day length, it has been indicated by Rhee and Gotham [7] that the population dynamics of phytoplankton in aquatic environments can change with season, latitude, and depth.Among factors affecting phytoplankton growth, nutrient has been an essential element [8][9][10], mainly including nitrogen and phosphate.Results reported by Ryther [11] indicated that phytoplankton indeed consumes lots of nitrogen and phosphate in their growth process, but reducing the nitrogen content in aquatic cannot slow the eutrophication.Using data from 17 lakes, Smith [8] analysed the influence of ratio of total nitrogen to phosphorus on the growth of blue-green algae (Cyanophya) and showed that controlling the ratio can help us improve the quality of aquatic environment very well.Obviously, the production process of phytoplankton is more complex.
However, due to the complexity and nonlinearity of aquatic ecosystem, there are some difficulties in understanding nutrient-phytoplankton dynamics only depending on experiment or field observation, which makes it necessary 2 Complexity to use models to provide quantitative insights into dynamic mechanism of phytoplankton growth.For different aquatic environments, we can use various modifications of the classical prey-predator models by introducing functional responses to model nutrient-phytoplankton dynamics [12][13][14].For example, Huppert et al. [15] describe the dynamics of nutrient-driven phytoplankton blooms by a simple model and identify, using the model analysis, an important threshold effect that a bloom will only be triggered when nutrients exceed a certain defined level.Additionally, most nutrientphytoplankton models reveal that phytoplankton population and nutrient population can coexist at equilibrium globally under some conditions [16,17].However, Sherratt and Smith [18] have reported that a constant population density may not exist in reality because of the existence of some factors, such as noise and physical factors.Actually, experiments and field observations show that the changes of phytoplankton population density usually possess oscillatory behaviour [19,20].
For the single cell phytoplankton species, in most studies of nutrient-phytoplankton models, it is usually assumed that the processes, such as conversion process of nutrient, in the dynamics of phytoplankton growth are instantaneous [14][15][16][17][19][20][21][22][23].It may be doubtful whether there exists the delay in the growth of phytoplankton over the large area or not.Yet, J. Caperon [24] studied time lag in population growth response of Isochrysis galbana, a phytoplankton species, to a variable nitrate environment by both experiments and model, and demonstrated the existence of delay in the growth of Isochrysis galbana.Hence, the delay may indeed exist in the phytoplankton growth, which means that it is necessary to consider delay in nutrient-phytoplankton models.An approach that has been attempted by researchers to model the dynamics of phytoplankton is the role of delay since delay appears as an important component in biosystems and ecosystems [25][26][27][28][29][30].
Actually, growing evidence shows that there exists time lag in some conversion processes from one state to another in some systems, and delay is an important factor because it can affect the dynamics of these systems.Volterra [31] considered time delay in a prey-predator model first and found oscillatory behaviour for the spatial distribution.For a long time, it has been recognized that delays can give rise to destabilizing effect of the dynamics of systems, where periodic solutions, as well as chaos, may emerge [32][33][34][35].Models incorporating delays in diverse biological and ecological models are extensively studied [36][37][38][39][40][41][42].Especially, the characteristic equation with respect to the linearized system of delay differential equations plays a key role in dynamic analysis, by which we can obtain some information on the stability of equilibrium.In addition, using the normal form theory, one can carry out the bifurcation analysis, such as the direction and stability of periodic solutions arising through Hopf bifurcation [43,44].
The main purpose of this paper is to consider the effects of multiple delays on the nutrient-phytoplankton dynamics.In [15], Huppert et al. presented a simple model to investigate effect of nutrient on phytoplankton blooms, and much better results are obtained.Here, this model is extended into a "two preys-one predator" type to describe nitrogen-phosphorusphytoplankton dynamics, as follows: where , , and  represent nitrogen, phosphorus, and phytoplankton population density at time , respectively;  1 is the nitrogen nutrients input flowing into the system and  2 is the phosphorus nutrients input flowing into the system;  1 is the loss rate of the nitrogen nutrients, and  2 is the loss rate of the phosphorus nutrients;  1 is nitrogen nutrient uptake rate of phytoplankton, and  2 is phosphorus nutrient uptake rate of phytoplankton;  1 and  2 denote the efficiency of nutrient utilization;  1 and  2 are time delay parameters;  is the mortality rate of phytoplankton.Although the function, which describes nutrient uptake dynamics, is not a Michaelis-Menten function, but Lotka-Volterra type, Huppert et al. [15] have indicated that the Lotka-Volterra term is a good first approximation to the Michaelis-Menten type.From biological viewpoint, all parameters are nonnegative.(), (), and () ≥ 0 are continuous on − ≤  < 0, where  = max( 1 ,  2 ) and (0), (0), and (0) > 0.
The paper is organized as follows.In Section 2, we analyze the existence and stability of positive equilibrium in model (1) without delays.In Section 3, we discuss stability of positive equilibrium and Hopf bifurcation under five different cases for delay effect.Subsequently, the direction of bifurcation and the stability of periodic solutions arising through Hopf bifurcation are given in Section 4. In order to analyze further how delay effects influence nutrient-phytoplankton dynamics, a series of numerical simulations are carried out in Section 5. Finally, the paper ends with conclusion in Section 6.

Existence and Stability of Positive Equilibrium in Model (1) without Delays
In this section, it is presented first that the first octant is positive invariant in model ( 1) without delays and the following lemma holds.

Lemma .
All the solutions of model (1) with initial conditions that initiate in { 3 + } are positive invariant in the absence of delays.
Proof.From the first equation of model (1), we have Hence, Likewise, from the second equation of model ( 1), we have Obviously, all the solutions of model (1) without delays are positive invariant if the initial conditions initiate in { 3 + }.Then, we complete the proof.
For model (1), it is obvious that the extinction equilibrium, ( 1 / 1 ,  2 / 2 , 0), exists.Moreover, in order to discuss the existence of positive equilibrium, the following function is defined: and then we can obtain where  * is the positive root of (4).
then there exists a unique positive equilibrium in model (1); otherwise, there is no positive equilibrium in model (1) Letting the unique positive equilibrium be  * = ( * ,  * ,  * ), then the following theorem holds for model (1) without delays.
eorem .If the unique positive equilibrium exists in model (1) without delays, then it is globally asymptotically stable.
Then, we complete the proof.

Complexity
Obviously, if  >  * , then the extinction equilibrium  0 is locally asymptotically stable.
Then, we complete the proof.
Actually, when  >  * holds, the positive equilibrium does not exist, and the extinction equilibrium  0 is globally asymptotically stable.Then, the following theorem holds.
Proof.We construct a Lyapunov function, as follows: In the model (1) without delays, Obviously, / < 0 Then, we complete the proof.

Local Stability Analysis and the Hopf Bifurcation
In this section, we first state the following positive invariant theorem.

Lemma .
All the solutions of model ( 1) with initial conditions that initiate in { 3 + } are positive invariant.
Therefore, all the solutions of model ( 1) are positive invariant if the initial conditions initiate in { 3 + }.Then, we complete the proof.
Then, according to [44], we obtain the following relevant parameter, which helps to determine the direction and stability of Hopf bifurcation: and where  1 and  2 can be determined by the following, respectively: ) , where Then, we can compute the following values: This determine the properties of bifurcating periodic solutions and the Hopf bifurcation at  =  * 20 .That is, (i)  determines the direction of the Hopf bifurcation.Specifically, when  > 0(< 0), the Hopf bifurcation is supercritical (subcritical).
(iii)  determines the period of the bifurcating periodic solutions; when  > 0(< 0), the period of bifurcating periodic solution increases (decrease).
According to the standard linear analysis, when  2 is equal to zero, the analysis reveals that the  2 −  1 parameter plane is divided into four parts (see Figure 1(a)).In Figure 1(a), before  2 reaching black dashed line, there exists  10 in model (1) such that the unique positive equilibrium loses its stability when the condition,  1 >  10 , holds.When the locus of  2 is between black dashed line and green dashed line, the stability switches for positive equilibrium do not exist although (28) has two positive roots, which means that there exists  10 in model (1) such that the unique positive equilibrium loses its stability when the condition,  1 >  10 , holds.However, the stability switches for positive equilibrium emerge when  2 is beyond green dashed line but it does not reach blue zone.When the locus of  2 is in blue zone, the unique positive equilibrium is always stable, which suggests that  1 cannot influence the stability of the positive equilibrium.When  1 equals zero, the similar results for  2 −  2 parameter plane are shown in Figure 1(b), but the sequence is reversed.Additionally, according to results in Section 4, we calculate the values of , , and  at  1 =  10 with  2 ∈ (0, 0.3), and the corresponding results are shown in Figure 1(c), where we can find that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable; especially, the period of the bifurcating periodic solutions increases as  2 increases.For other cases of  1 and  2 , the same procedures with respect to calculations of , , and  can be performed like Figure 1(c).
As examples corresponding to stability of the positive equilibrium with  2 = 0.2, taken  1 = 1 and  1 = 3 in Figure 1(a), respectively, the corresponding numerical solutions are shown in Figure 2. Obviously, the positive equilibrium is stable because  1 = 1 is below  10 (see Figure 2(a)).In contrast, due to 3 =  1 beyond  10 , a periodic solution exists (see Figure 2(b)).Furthermore, set  2 = 0.7, then we have  0 1 ≈ 4.9050 <  0 1 ≈ 19.0615 <  1 1 ≈ 38.6683.Taken  1 = 4,  1 = 18 and  1 = 21 in Figure 1(a), respectively, the corresponding numerical solutions are shown in Figure 3. Obviously, the positive equilibrium is stable when  1 = 4 and  1 = 21 (see Figures 3(a) and 3(c)), but the positive equilibrium is unstable when  1 = 18 (see Figure 3(b)), which means that the positive equilibrium can gain its stability again for  1 >  10 .In Figure 3, the same initial values are applied, and other parameter values except for  1 are also identical.Clearly, the delay is the principal factor giving rise to the difference among (a), (b), and (c).
Numerical solutions in Figure 3 suggest that the stability switches induced by delay may exist.Hence, the bifurcation diagram in  1 −  2 parameter plane is given (see Figure 4(a)).For case  2 = 0, there exists a  * 1 such that the positive equilibrium with respect to  1 ∈ (0,  * 1 ) is stable.For case  1 = 0, there exists a  * 2 such that the positive equilibrium with respect to  2 ∈ (0,  * 2 ) is stable.Additionally, when  2 ∈ (0,  * 2 ), Figure 4(a) shows that  10 exists such that the positive equilibrium with respect to  1 ∈ (0,  10 ) is stable.Likewise, when  1 ∈ (0,  * 1 ), Figure 4(a) also display that  20 exists such that the positive equilibrium with respect to  2 ∈    (0,  20 ) is stable.Significantly, Figure 4(a) demonstrates that the stability switches for positive equilibrium with respect to  1 emerge when  2 below  * 2 is fixed.Figure 4(b) depicts the dependence of stability of positive equilibrium on delay  2 in the fully nonlinear regime for  1 in sequence [0, 2, 10, 30] and  1 =  2 , which is consistent with results in Figure 4(a).However, when  2 is fixed, Figure 4(c) shows that the stability switches emerge with  1 increases.Especially, periodic-2 solutions exist for some values of  2 (see magenta zone in Figure 4(c)).As an example of periodic-2 solution existence, taking  1 = 41, a phase and a time-series are given in the inner of Figure 4(c).Moreover, Figure 4(d) shows that there exist periodic-3 solutions in model (1).According to results in Section 2, the positive equilibrium is globally asymptotically stable in the model ( 1) without delay if it exists.Obviously, the results shown in Figure 4 are induced by delay.
In Figure 4(a), we can find that the number of intervals corresponding to stability of positive equilibrium for small values of  2 equals 3.However, when  2 is beyond dashed line ( * * 2 ), the number of intervals is 2.So the number of intervals for stability switches may be different for diverse  2 .Accordingly, we calculate the number of intervals with respect to parameter  2 , as shown in Figure 5(a).Figure 5(b) shows that there exist 3 stable intervals when  2 = 0.7 and  2 = 0, which is an example of Figure 5 predicted by linear analysis and numerical results in the fully nonlinear regime with ( 2 ,  2 ) at (0, 0.7), where top panel represents numerical results for model (1) in the fully nonlinear regime; bottom panel depicts the results predicted by linear analysis, and blue and red area denote stable and unstable, respectively.

Conclusions
In this paper we proposed a nitrogen-phosphorusphytoplankton model with multiple delays.The analysis focused on the effect of delay on nutrient-phytoplankton dynamics.In the absence of delay, theoretical analysis indicated that the unique positive equilibrium is globally asymptotically stable in model (1) if it exists.Deng et al. [48] also studied a nitrogen-phosphorus-phytoplankton model without delay, where Holling II function was employed to describe the nutrient uptake dynamics of phytoplankton.
Although the function modelling the nutrient uptake dynamics of phytoplankton is different, they get the same results.These results mean that the nutrient-phytoplankton ecosystem will approach the stable equilibrium.However, it has been reported [18] that a constant population density may not exist because of the existence of some factors including noise, interval factors, and physical factors.And ecological studies [49,50] also criticized this idea of "the balance of nature."Actually, the existence of nutrient-plankton oscillations has been detected by laboratory experiments and field observation [19,20].Additionally, Benincà et al. [49] present the first experimental demonstration of chaos in a long-term experiment with a complex food web, where the food web was consisted of bacteria, several phytoplankton species, herbivorous and predatory zooplankton species, and detritivores.And they also find that the community moved back and forth between stabilizing and chaotic dynamics during the cyclic succession, and their findings provide a field demonstration of nonequilibrium coexistence of competing species through a cyclic succession at the edge of chaos [50].These reports support that the nonequilibrium dynamics, such as oscillations and chaos, can exist in reality.
In the present paper, we find that the unique positive equilibrium may lose its stability via Hopf bifurcation when delay appears, and then a periodic solution emerges, which means that nutrient-phytoplankton oscillation occurs.Obviously, the factor giving rise to nutrient-phytoplankton oscillation is delay in our studies.And the period and the stability of the bifurcating periodic solutions with respect to delay are discussed by using center manifold argument and normal form theory.In fact, instability induced by delay in nutrientplankton model has been studied widely, and many studies indicate that the equilibrium is always unstable when delay is beyond a critical value [25,29,51,52].Yet, it should be emphasized in the present paper that the stability switches induced by delay can occur under some conditions.Moreover, numerical simulations showed how the delay influences nutrient-phytoplankton dynamics.Numerical results for model (1) in the fully nonlinear regime are consistent with the linear analysis.In numerical simulations, we found that delay indeed gives rise to the emergence of stability switches for the positive equilibrium.Yet, the numerical results show that the parameter intervals for stability switches may depend on other parameters as well, e.g.,  2 .Additionally, numerical results also indicated that periodic-2 solutions and periodic-3 solutions can emerge under some conditions for delay, which means that complex dynamics induced by delay exist in model (1).From biological viewpoint, the existence of periodic solutions implies that the fluctuations exist in density of phytoplankton population; that is, nutrient-phytoplankton oscillation emerges.Especially, by Li and York's theory, periodic-3 solution implies chaos, which means that chaotic density fluctuations can display a variety of different periodicities and the long-term prediction of phytoplankton density can be fundamentally impossible.The chaotic density fluctuations do not contribute to the control of phytoplankton bloom.Consequently, the importance of the present paper is not the precision with which it predicts specific events for phytoplankton blooms but its contribution to the studies on how the delay influences nutrient-phytoplankton dynamics.

Figure 4 :
Figure 4: (a) Bifurcation diagram with  2 = 0.7 for  1 V.S.  2 , where the symbol "S" denotes stable and the symbol "US" denotes unstable; (b) bifurcation diagram with  2 = 0.7 for the effect of  2 on nutrient-phytoplankton dynamics, where dashed line denotes unstable; solid line denotes stable; the green solid square is Hopf bifurcation point, dot-dashed line corresponding to  * 2 in (a), blue line represents equilibrium, and red line represents amplitude of periodic solutions.(c) Bifurcation diagram with  2 = 0.7 and  2 = 1, where the yellow solid square is Hopf bifurcation point, the blue solid circle is bifurcation point for periodic-2 solution, the magenta zone indicates the existence of periodic-2 solutions, and the cyan solid diamond denotes the value of  1 for a phase and a time-series in the inner of (c).(d) A periodic-3 solution for phytoplankton population, where  2 = 0.7,  1 = 100, and  2 = 2.

Figure 5 :
Figure5: (a) The number of intervals for stability switches of positive equilibrium with respect to  2 , where  2 = 0; (b) based on (a), an example for comparison between results predicted by linear analysis and numerical results in the fully nonlinear regime with ( 2 ,  2 ) at (0, 0.7), where top panel represents numerical results for model(1) in the fully nonlinear regime; bottom panel depicts the results predicted by linear analysis, and blue and red area denote stable and unstable, respectively.