Time Delayed Equations as Models in Nature and SocietyView this Special Issue
Research Article | Open Access
Nonlinear Dynamics of a Nutrient-Phytoplankton Model with Time Delay
We consider a nutrient-phytoplankton model with a Holling 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.
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–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 . 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–11], and various studies have simulated the dynamics of plankton systems.
Huppert et al.  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 and . 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  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 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.  and Mäler , 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: , , , , and , where , are continuous bounded functions in the interval .
Huppert et al.  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–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.
2. Stability of the Positive Equilibria and Existence of the Hopf Bifurcation
2.1. The Existence of Positive Equilibria
In this section, we consider the existence of the positive equilibria for system (3). We define , . The equations for the equilibria can be depicted as follows:
According to , a vertical isocline can be obtained, : , which can be written as . For , when the value of is close to positive infinite, the value of is close to . By contrast, when the value of approaches zero from the right-hand side, the value of is approximately positive infinite. Due to the vertical isocline, is continuous and derivable when and ; thus we can compute the derivative of ; that is, , and four roots exist when , as follows:From the analysis above, (i) if , we can obtain for an arbitrary real number ; thus decreases monotonically when ; (ii) if , four real number roots exist when . Obviously, . According to the second derivative of , we can find that , are the points with the extreme maximum value and , are the points with the extreme minimum value. We can find that decreases monotonically when and increases monotonically when .
According to , the horizontal isocline : and , which can be written as and . When the horizontal isocline : , that is, , it is easy to find that the line, , is the asymptote of . Obviously, is continuous when , and we can find that decreases monotonically when . For , holds when , whereas holds when .
First, when the system comprises only nutrients without phytoplankton in the system, according to the vertical isocline, : , and according to the horizontal isocline, : . From the analysis above, if , we find that there is always a boundary equilibrium. If , one boundary equilibrium exists at least. The boundary equilibrium can be written as , .
In addition, according to the vertical isocline, : , and according to the horizontal isocline, : . From the discussion above, if , we find that a positive equilibrium exists in the system when the condition that holds. Thus, we can find that there is no positive equilibrium in the system when the condition that holds. If , one positive equilibrium exists at least when holds. The positive equilibrium can be written as .
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 is
The index of equilibrium is when the condition , is satisfied, which is stable. According to the eigenvalues of , it is locally asymptotically stable when the above conditions are satisfied. In particular, the equilibrium is unstable, when or .
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
When , the index of equilibrium is if . It is locally asymptotically stable using the Routh-Hurwitz criteria when the above conditions are satisfied. In particular, the equilibrium is unstable, when , which is a saddle.
2.2. Local Stability of Equilibrium
The linearization of (3) at the equilibrium isTherefore, the characteristic matrix at equilibrium isIts characteristic equation iswhere , , , and .
Theorem 1. If and hold, the boundary equilibrium is locally asymptotically stable.
Proof. There are two eigenvalues where and . When all the roots of the characteristic equation at are negative, the boundary equilibrium is a stable node. Thus, the boundary equilibrium is locally asymptotically stable. The proof is complete.
2.3. Local Stability and the Hopf Bifurcation of Equilibrium
In the following, we focus on the existence of a local Hopf bifurcation at positive equilibrium . If we let , the linearization of (3) at the equilibrium is as follows:Thus, the characteristic matrix at equilibrium isSet , , , and .
Its characteristic equation isorwhere , , and .
Theorem 2. For system (3), the positive equilibrium is locally asymptotically stable for all and the system undergoes a Hopf bifurcation at the positive equilibrium for , , when the conditions that and hold.
Lemma 3. If , , and hold, the positive equilibrium is locally asymptotically stable.
Proof. When , the characteristic equation isIf and , 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. LetIfthen (13) has a pair of purely imaginary roots .
Proof. When , we assume that is a root of (13). By substituting this into (13), we obtainorIf we let and , (19) is equal toand thusIf we let , , we can obtain ; therefore, , and then .
However, we can solve the first formula in (20) to obtainThus, when , (13) has a pair of purely imaginary roots . The proof is complete.
Lemma 5. For arbitrary , , 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 obtainThus, we haveFrom (20) and (21), we obtainbecause ofand thusThis shows that all of the roots cross the imaginary axis at from left to right as increases. The proof is complete.
3. 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 , , using the method introduced by Hassard et al. .
If we let , , , and , system (3) can be written aswhere , , , for , andwhere denotes “higher order terms.” By the Riesz representation theorem, a function of bounded variation exists for such thatIn fact, we can choosewhere
In the following discussion, we omit the term in (31) because determining the direction and stability of the Hopf bifurcation only requires up to the third term. For , we defineThen, (29) can be rewritten asFor , we define the adjoint operator as follows:For and , we define the bilinear form:where . In Section 2, we found that and are eigenvalues of and . Thus, we need to calculate the eigenvectors of and that correspond to and , respectively. Suppose that and are the eigenvectors that correspond to and , respectively. LetUsing the condition that , that is,thus we can obtainSimilarly, suppose that is the eigenvector of that corresponds to ; then , and from the definition of , we haveTherefore, we obtainLet , and due towe can obtain
In the following, we first construct the coordinates to describe the center manifold at (i.e., ). We defineOn the center manifold , we havewhere and are local coordinates for the center manifold at . From (36), (37), (38), (39), and (47), we can obtainwhereBased on (37) and (47), we havewhereNear the origin on the center manifold , by (49)–(51), we haveBy (49), (52), (53), and (54) and by comparing the coefficients, we haveWe considerand because ofwe obtainThus, by comparing the coefficients with (51), we obtain and are still unknown; thus, we need to compute them. For , we haveBy comparing the coefficients, we haveIt follows from (55) thatSolving for , we obtainSimilarly, we obtainwhere and are both two-dimensional vectors, which can be determined by setting in . For , we have ; thusAccording to (55) and (56), we haveBy substituting (68) and (69) into (65) and (66), respectively, we obtainwhere is the identity matrix. From the above, we already know and , and can be expressed in terms of the parameters and delay ; thus, we can compute the values as follows:
Theorem 6. , , and are defined above; thus the following are true:(i)if , then the Hopf bifurcation is supercritical (subcritical);(ii)if , then the periodic solutions are stable (unstable);(iii)if , then the period of the bifurcating periodic solution of system (2) increases (decreases).
4. 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 , , , , , and , 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 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 can be obtained by . We set the parameter and the others are the same as the above. By combining the analysis above, we can obtain the critical value: . 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 . The Hopf bifurcation occurred at the critical value of . Thus, the oscillation will occur if . 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 .