Nonlinear Dynamics of a Nutrient-Phytoplankton Model with Time Delay

Abstract

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.

1. 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–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–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 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 [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 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: , , , , and , where , are continuous bounded functions in the interval .

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–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 is Therefore, the characteristic matrix at equilibrium is Its characteristic equation is where , , , 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 is Set ,  ,  , and .

Its characteristic equation is or where , , 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 is If 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. Let If then (13) has a pair of purely imaginary roots .

Proof. When , we assume that is a root of (13). By substituting this into (13), we obtain or If we let and , (19) is equal to and thus If we let , , we can obtain ; therefore, , and then .
However, we can solve the first formula in (20) to obtain Thus, 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 obtain Thus, we have From (20) and (21), we obtain because of and thus 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.

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. [24].

If we let ,  ,  , and , system (3) can be written as where , , , for , and where denotes “higher order terms.” By the Riesz representation theorem, a function of bounded variation exists for such that In fact, we can choose where

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 define Then, (29) can be rewritten as For , 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. Let Using the condition that , that is, thus we can obtain Similarly, suppose that is the eigenvector of that corresponds to ; then , and from the definition of , we have Therefore, we obtain Let , and due to we can obtain

In the following, we first construct the coordinates to describe the center manifold at (i.e., ). We define On the center manifold , we have where and are local coordinates for the center manifold at . From (36), (37), (38), (39), and (47), we can obtain where Based on (37) and (47), we have where Near the origin on the center manifold , by (49)–(51), we have By (49), (52), (53), and (54) and by comparing the coefficients, we have We consider and because of we obtain Thus, by comparing the coefficients with (51), we obtain and are still unknown; thus, we need to compute them. For , we have By comparing the coefficients, we have It follows from (55) that Solving for , we obtain Similarly, we obtain where and are both two-dimensional vectors, which can be determined by setting in . For , we have ; thus According to (55) and (56), we have By substituting (68) and (69) into (65) and (66), respectively, we obtain where 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.

(a)

(b)

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 .

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 .

(a)

(b)

Figure 2 

The solution for nutrients where : (a) (red dashed line), (blue solid line); (b) (red dashed line), (blue solid line).

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(b).

(a)

(b)

Figure 3 

(a) Solutions for nutrients and phytoplankton where , . (b) Phase diagram.

5. 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 31170338), by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001), and by the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510).