Nonlinear Dynamics of a Nutrient-Phytoplankton Model with Time Delay

We consider a nutrient-phytoplanktonmodel with aHolling type II functional response and a time delay. By selecting the time delay used as a bifurcation parameter, we prove that the system is stable if the delay value is lower than the critical value but unstable when it is above this value. First, we investigate the existence and stability of the equilibria, as well as the existence of Hopf bifurcations. Second, we consider the direction, stability, and period of the periodic solutions from the steady state based on the normal form and the center manifold theory, thereby deriving explicit formulas. Finally, some numerical simulations are given to illustrate the main theoretical results.


Introduction
Ecological systems are characterized by the relationships between species and their natural environment.One of the key factors that affect population dynamics is predator.Due to its universal existence and importance in nature, the dynamic interactions between predators and their prey have remained one of the dominant themes in ecological population dynamics since the origin of this discipline [1][2][3][4].
Phytoplankton has vital effects in aquatic ecosystems where it plays a significant role as the base of the food chain.Phytoplankton controls the global carbon cycle which has an important effect on climate regulation [5].In some situations, marine waters, reservoir, and lakes may exhibit plankton or algal blooms [6,7], where drinking water may be polluted and algal toxins can affect human health through the food chain, while the marine biodiversity may also be decreased.
Nutrient-phytoplankton systems have a very important effect on aquatic ecosystems as part of predator-prey systems.Recently, many researchers have investigated a predator-prey model that involves nutrients, phytoplankton, and zooplankton [8][9][10][11], and various studies have simulated the dynamics of plankton systems.
Huppert et al. [12] proposed a generic bottom-up nutrient-phytoplankton model: where the forced  compartment model comprises two variables: the nutrient concentration  and the phytoplankton biomass .The initial conditions for this model are (0) =  0 > 0 and (0) =  0 > 0.  is an external source of nutrients that flows into the system from the environment,  and  are uptake rates, and  is the loss rate of nutrients.Seasonal environment conditions such as the water temperature, salinity, light, and thermocline depth are often important to different degrees depending on the system under investigation.The modulation of phytoplankton growth  can be represented as a periodic function () = ( + ), where  is the period of forcing, which is assumed to be annual in the present study.Analyses and numerical simulations can obtain insights into the mechanism of bloom recurrence, where modifications to the equations via the inclusion of appropriate functional forms can generate more realistic dynamics.Mäler [13] considered a better formulation of the dynamics of the eutrophication process: where  is the runoff of nutrients into the lake and  is the constant rate of the natural removal of nutrients.After reducing the levels of green plants, the bottom sediments will be more vulnerable to winds, waves, and bottom-feeding fishes.Finally, nutrient sediments are released into the water and they contribute to the further growth of phytoplankton.
The term  2 /1 +  2 represents the feedback from the bottom sediments.Depending on the numerical values for the parameters and specific simulations of the system, different negative externalities can affect the system in various ways.
In this study, based on the idea of Huppert et al. [12] and Mäler [13], we consider a nutrient-phytoplankton model with delay: where two variables  and  represent the concentrations of nutrients and phytoplankton, respectively.All of the parameters are positive constants, where  is a positive delay, that is, the time required to convert nutrients into phytoplankton,  is a constant concentration, that is, the nutrient runoff from the environment,  is the loss rate of the nutrient concentration,  is the predation rate of nutrients by phytoplankton,  denotes the maximum nutrient intake rate of phytoplankton,  is the mortality of phytoplankton,   represents the maximum predation rate of zooplankton on phytoplankton, and   is the half-saturation constant for phytoplankton, and   (/(  + )) represents the predation of zooplankton on phytoplankton.The initial conditions are as follows: () = () ≥ 0, () = () ≥ 0,  ∈ [−, 0), (0) > 0, and (0) > 0, where (), () are continuous bounded functions in the interval [−, 0).Huppert et al. [12] described the dynamics of the system above without considering the effect of delay.Time delay plays an important role in reflecting the real dynamical behaviors of biological systems and many studies have addressed this issue in recent years [14][15][16][17][18][19][20][21][22][23].
In this study, in order to investigate the effects of a time delay on the system, we selected the delay  as a bifurcation parameter.The remainder of this paper is organized as follows.In Section 2, we prove that positive equilibria exist under certain conditions and we then analyze the local stability of the boundary equilibrium.We also analyze the stability of the positive equilibria and the critical conditions when the Hopf bifurcation occurs.In Section 3, we derive an explicit algorithm and sufficient conditions for determining the direction of the Hopf bifurcation and the stability of the periodic solutions.In Section 4, we present some numerical simulations to illustrate our theoretical results.In Section 5, we give some brief conclusions.
In the analysis above, we only prove that positive equilibria exist in the system when certain conditions are satisfied but we do not indicate that positive equilibria do not exist when the conditions are not satisfied.
The Jacobian matrix of the system without time delay at the equilibrium  0 ( 0 , 0) is The index of equilibrium  0 is +1 when the condition 2 0 < (1 +  0 2 ),  0   <   +   is satisfied, which is stable.According to the eigenvalues of  1 , it is locally asymptotically stable when the above conditions are satisfied.In particular, the equilibrium  0 ( 0 , 0) is unstable, when 2 0 > (1 +  0 2 ) or  0   >   +   .From the above discussion, we find a positive equilibrium in the system without time delay under some preconditions, defined by  * ( * ,  * ).The Jacobian matrix of the system at the equilibrium  * ( * ,  * ) is It is locally asymptotically stable using the Routh-Hurwitz criteria when the above conditions are satisfied.In particular, the equilibrium , which is a saddle.

Local Stability of Equilibrium
Therefore, the characteristic matrix at equilibrium  0 ( 0 , 0) is Its characteristic equation is where Proof.There are two eigenvalues where − + (2 0 /(1 + ( 0 ) 2 )) and  0 −−  /  .When all the roots of the characteristic equation at  0 are negative, the boundary equilibrium  0 is a stable node.Thus, the boundary equilibrium  0 is locally asymptotically stable.The proof is complete.
Proof.When  = 0, the characteristic equation is If − < 0 and  +  0 > 0, then all the roots of characteristic equation ( 15) have negative real parts.Thus, the positive equilibrium  * ( * ,  * ) is locally asymptotically stable.The proof is complete.

Lemma 4.
Let If then ( 13) has a pair of purely imaginary roots ± 0 .
Lemma 5.For arbitrary  =   ,  = 0, 1, 2, . .., that satisfy the above conditions, one can find that the derivative of the real number part of the eigenvalue is greater than zero; that is, We can find that the roots of ( 13) cross the imaginary axis in a transverse direction from left to right as  increases.
Proof.Differentiating characteristic equation ( 13) with respect to , we obtain Thus, we have From ( 20) and ( 21), we obtain This shows that all of the roots cross the imaginary axis at  from left to right as  increases.The proof is complete.
From Lemmas 3-5 and the Hopf bifurcation theorem, the proof is complete.

Direction and Stability of the Hopf Bifurcations
In the previous section, we obtained some conditions for the occurrence of Hopf bifurcations.In this section, we consider the direction, stability, and period of the periodic solutions from the steady state, as well as deriving the explicit formulae that determine these factors at the critical value  =   ,  = 0, 1, 2, . .., using the method introduced by Hassard et al. [24].
By substituting (68) and ( 69) into (65) and (66), respectively, we obtain where  is the 2×2 identity matrix.From the above, we already know  20 () and  11 (), and  21 can be expressed in terms of the parameters and delay ; thus, we can compute the values as follows: Theorem 6.  2 ,  2 , and  2 are defined above; thus the following are true: , then the periodic solutions are stable (unstable); (iii) if  2 > 0 ( 2 < 0), then the period of the bifurcating periodic solution of system ( 2) increases (decreases).

Numerical Simulation
According to the analysis above, the positive equilibrium does not always exist.Thus, based on the numerical technology, we can obtain the positive equilibria that exist under certain conditions.When some parameters are set, we select  = 0.3,  = 0.5,  = 0.75,  = 0.2,   = 0.6, and   = 4, and we can obtain the different results shown in Figure 1. Figure 1(a) shows that different nutrient concentrations correspond to various vertical isoclines and we can determine the different intersection points (i.e., equilibria).Figure 1(b) shows the stability with the change in the nutrient concentration when the other parameters are fixed, where the blue solid line shows that the positive equilibrium is unstable, whereas the red solid line shows that the positive equilibrium is stable.
Based on the two previous sections, we know that the stability of the boundary equilibria  0 is relatively simple.However, the stability of the positive equilibria  * is complex.We studied the effect of a time delay  on the positive equilibria of system (3) and our results showed that the time delay has a vital role in system.The critical value  0 can be obtained by  0 = (1/ 0 ) arccos(( 2 0 −  0 )/).We set the parameter  = 1.3 and the others are the same as the above.By combining the analysis above, we can obtain the critical value:  0 = 3.91.The numerical solutions for nutrients are shown in Figures 2(a) and 2(b).In order to study the effect of a time delay, we set  as 0 and 1 initially, and Figure 2(a) shows that the solution converged to a positive equilibrium in the end.The positive equilibrium is locally asymptotically stable in system (3) without delay.The oscillation does not occur when the system had a delay and the parameter  <  0 = 3.91.The Hopf bifurcation occurred at the critical value of  0 .Thus, the oscillation will occur if  >  0 = 3.91.The numerical solutions for nutrients are shown in Figure 2(b), and we set  as 0 and 4. Figure 2(b) shows that the oscillation occurs when  = 4. Furthermore, in order to study the relationship between nutrients and phytoplankton, we set the parameter  as 5 and the numerical solutions for nutrients and phytoplankton are shown in Figure 3. Figure 3(a) shows the relationship between the spatial distributions of nutrients and phytoplankton, which indicates that the density of phytoplankton reaches its maximal value when the concentration of nutrients reaches its minimal value.The density of phytoplankton affects the concentration of nutrients.As the density of phytoplankton increases, the density of nutrients decreases.Therefore, the density of nutrients and phytoplankton are mutually constrained.The phase of system (3) is shown in Figure 3

Conclusion
In this study, we considered a biological system that comprises nutrients and phytoplankton, where we focused on the effects of a delay on the system dynamics.
In Section 2, we showed that the positive equilibrium  * is locally asymptotically stable if there is no time delay.However, as the time delay increases above the critical value, the equilibrium state becomes unstable and a Hopf bifurcation occurs.Using the method introduced by Hassard et al. [24], we derived the direction and stability of the Hopf bifurcation.
Based on the computer simulation, we showed that the relationship between nutrients and phytoplankton is mutually constrained.Thus, an abundance of nutrients leads to a major increase of phytoplankton, which depleted the nutrients via consumption.Due to the effect of time delay, the system shows oscillation.It reveals that time delay has a vital effect on the system.
The main limitations of this study are that the main theoretical results were confirmed by numerical simulations with hypothetical parameter values.Thus, we aim to obtain some real data to confirm the validity of our system.Furthermore, we did not consider the diffusion and advection of nutrients and phytoplankton; therefore, these aspects should be considered comprehensively in our future research.

Figure 1 :Figure 2 :
Figure 1: (a) The intersections between the vertical isocline and horizontal isocline with different nutrient concentrations , where the red solid line represents the horizontal isocline, the grey dashed line represents the asymptote of the horizontal isocline, and the other lines with different colors represent the vertical isocline that corresponds to different nutrient concentrations .(b) The stability of the positive equilibrium with the change in the nutrient concentration . (b).