Discrete Dynamics in Nature and Society
Volume 2015 (2015), Article ID 939187, 12 pages
http://dx.doi.org/10.1155/2015/939187

Research Article

Nonlinear Dynamics of a Nutrient-Phytoplankton Model with Time Delay

DeBing Mei,1,2 Min Zhao,2,3 Hengguo Yu,1,2 Chuanjun Dai,2,3 and Yi Wang1,2

¹School of Mathematics and Information Science, Wenzhou University, Wenzhou, Zhejiang 325035, China
²Zhejiang Provincial Key Laboratory for Water Environment and Marine Biological Resources Protection, Wenzhou University, Wenzhou, Zhejiang 325035, China
³School of Life and Environmental Science, Wenzhou University, Wenzhou, Zhejiang 325035, China

Received 1 April 2015; Accepted 9 June 2015

Academic Editor: Luca Guerrini

Copyright © 2015 DeBing Mei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

References

A. A. Berryman, “The origins and evolution of predator-prey theory,” Ecology, vol. 73, no. 5, pp. 1530–1535, 1992.

H.-B. Shi, W.-T. Li, and G. Lin, “Positive steady states of a diffusive predator-prey system with modified Holling-Tanner functional response,” Nonlinear Analysis. Real World Applications, vol. 11, no. 5, pp. 3711–3721, 2010.

M. Baurmann, T. Gross, and U. Feudel, “Instabilities in spatially extended predator-prey systems: spatio-temporal patterns in the neighborhood of Turing-Hopf bifurcations,” Journal of Theoretical Biology, vol. 245, no. 2, pp. 220–229, 2007.

G.-Q. Sun, Z. Jin, Q.-X. Liu, and L. Li, “Dynamical complexity of a spatial predator-prey model with migration,” Ecological Modelling, vol. 219, no. 1-2, pp. 248–255, 2008.

S. Abbas, M. Banerjee, and N. Hungerbühler, “Existence, uniqueness and stability analysis of allelopathic stimulatory phytoplankton model,” Journal of Mathematical Analysis and Applications, vol. 367, no. 1, pp. 249–259, 2010.

Y. P. Wang, M. Zhao, C. J. Dai, and X. H. Pan, “Nonlinear dynamics of a nutrient-plankton model,” Abstract and Applied Analysis, vol. 2014, Article ID 451757, 10 pages, 2014.

M. Sandulescu, C. López, E. Hernández-García, and U. Feudel, “Plankton blooms in vortices: the role of biological and hydrodynamic timescales,” Nonlinear Processes in Geophysics, vol. 14, no. 4, pp. 443–454, 2007.

C. J. Dai, M. Zhao, and L. S. Chen, “Complex dynamic behavior of three-species ecological model with impulse perturbations and seasonal disturbances,” Mathematics and Computers in Simulation, vol. 84, pp. 83–97, 2012.

W. Wang, Q.-X. Liu, and Z. Jin, “Spatiotemporal complexity of a ratio-dependent predator-prey system,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 75, no. 5, Article ID 051913, 9 pages, 2007.

M. R. Garvie, “Finite-difference schemes for reaction-diffusion equations modeling predator-prey interactions in MATLAB,” Bulletin of Mathematical Biology, vol. 69, no. 3, pp. 931–956, 2007.

W. W. Zhang and M. Zhao, “Dynamical complexity of a spatial phytoplankton-zooplankton model with an alternative prey and refuge effect,” Journal of Applied Mathematics, vol. 2013, Article ID 608073, 10 pages, 2013.

A. Huppert, B. Blasius, R. Olinky, and L. Stone, “A model for seasonal phytoplankton blooms,” Journal of Theoretical Biology, vol. 236, no. 3, pp. 276–290, 2005.

K.-G. Mäler, “Development, ecological resources and their management: a study of complex dynamic systems,” European Economic Review, vol. 44, no. 4–6, pp. 645–665, 2000.

M. Zhao, X. Wang, H. Yu, and J. Zhu, “Dynamics of an ecological model with impulsive control strategy and distributed time delay,” Mathematics and Computers in Simulation, vol. 82, no. 8, pp. 1432–1444, 2012.

F. Lian and Y. Xu, “Hopf bifurcation analysis of a predator-prey system with Holling type IV functional response and time delay,” Applied Mathematics and Computation, vol. 215, no. 4, pp. 1484–1495, 2009.

Y. Su, J. Wei, and J. Shi, “Hopf bifurcations in a reaction-diffusion population model with delay effect,” Journal of Differential Equations, vol. 247, no. 4, pp. 1156–1184, 2009.

H. Yu, M. Zhao, and R. P. Agarwal, “Stability and dynamics analysis of time delayed eutrophication ecological model based upon the Zeya reservoir,” Mathematics and Computers in Simulation, vol. 97, pp. 53–67, 2014.

Y. Chen and F. Zhang, “Dynamics of a delayed predator-prey model with predator migration,” Applied Mathematical Modelling, vol. 37, no. 3, pp. 1400–1412, 2013.

S. S. Chen, J. P. Shi, and J. J. Wei, “Time delay-induced instabilities and Hopf bifurcations in general reaction-diffusion systems,” Journal of Nonlinear Science, vol. 23, no. 1, pp. 1–38, 2013.

M. Wang and G. Lv, “Entire solutions of a diffusive and competitive Lotka-Volterra type system with nonlocal delays,” Nonlinearity, vol. 23, no. 7, pp. 1609–1630, 2010.

X. H. Pan, M. Zhao, C. J. Dai, and Y. P. Wang, “Stability and Hopf bifurcation analysis of a nutrient-phytoplankton model with delay effect,” Abstract and Applied Analysis, vol. 2014, Article ID 471507, 10 pages, 2014.

J. Z. Zhang, Z. Jin, J. R. Yan, and G. Q. Sun, “Stability and Hopf bifurcation in a delayed competition system,” Nonlinear Analysis. Theory, Methods & Application, vol. 70, no. 2, pp. 658–670, 2009.

S. Chen, J. Shi, and J. Wei, “The effect of delay on a diffusive predator-prey system with Holling type-II predator functional response,” Communications on Pure and Applied Analysis, vol. 12, no. 1, pp. 481–501, 2013.

B. D. Hassard, N. D. Kazarinoff, and Y. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, UK, 1981.