We analyze a nutrient-plankton system with a time delay. We choose the time delay as a bifurcation parameter and investigate the stability of a positive equilibrium and the existence of Hopf bifurcations. By using the center manifold theorem and the normal form theory, the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions are researched. The theoretical results indicate that the time delay can induce a positive equilibrium to switch from a stable to an unstable to a stable state and so on. Numerical simulations show that the theoretical results are correct and feasible, and the system exhibits rich complex dynamics.

1. Introduction

Plankton is the basis of all aquatic food chains and its importance for marine ecosystems is widely recognized [1]. The dynamics of zooplankton-phytoplankton systems have been discussed by many authors [28]. A wealth of studies have shown that a time delay may have a complex effect on the dynamics of such systems, with effects that include stability switches for equilibria, the existence of Hopf bifurcation periodic solutions, and the direction and stability of Hopf bifurcation [918].

Algal blooms are a common feature of marine ecosystems, and a high nutrient concentration has an important influence on algal blooms [19]. Plankton-nutrient interaction models can provide a better understanding of plankton dynamics. The model formulated by Huppert et al. [20] consists of two variables, nutrient level and phytoplankton biomass , in the following system:Considering the biological significance, all the constants, , , , , and , are assumed to be positive. The phytoplankton growth rate can be represented as a periodic function , where is the forcing period affected by seasonal conditions such as light, salinity, and water temperature.

Harmful algal blooms have become a serious environmental problem worldwide and have been widely studied [8, 21]. On the basis of field studies and mathematical modeling, Chattopadhayay et al. [22] proposed the following phytoplankton-zooplankton model:The authors studied the existence and local stability of steady states and the existence of Hopf bifurcation for system (2) by taking various combinations of and [22]. The results showed that toxin-producing plankton is helpful in terminating planktonic blooms by decreasing the grazing pressure of zooplankton.

To reflect maturation time, capture time, and other factors, a time delay is often included in mathematical models of population dynamics. Incorporation of a time delay can provide a better understanding of the dynamics of biological models. In recent years, many researchers have studied the fact that time delay plays important roles in biological dynamical systems [18, 2328]. Following Huppert et al. [20] and Chattopadhayay et al. [22], we consider the following plausible nutrient-phytoplankton-zooplankton system with a time delay:where is the nutrient concentration, and are the density of the phytoplankton and zooplanktons population, respectively, is the external source of nutrients flowing into the system, is the small loss rate to reflect sinking of nutrients from the epilimnion down to the hypolimnion, and are the capture rates for phytoplankton on nutrient and zooplankton on phytoplankton, and and denote the rates of biomass conversion. is the half-saturation constant, and are the death rates for phytoplankton and zooplankton, and is the interspecies competition coefficient for phytoplankton, which reflects phytoplankton self-limitation. The parameter denotes the time delay, which arises because of the time required by phytoplankton to absorb nutrients.

The remainder of the paper is organized as follows. In Section 2, we determine the stability of the positive equilibrium and the existence of Hopf bifurcation. Section 3 illustrates the direction of Hopf bifurcation and the stability of bifurcating periodic solutions. To verify the theoretical analysis, we present numerical simulations in Section 4. Finally, Section 5 provides some conclusions.

2. Stability and Existence of Hopf Bifurcation

In this section, we study the local stability of the positive equilibrium of system (3). The positive equilibrium corresponds to the coexistence of the three species, whereIt is obvious that exists if and only ifFrom (5), we know that the positive equilibrium of system (3) exists if external nutrients flowing into the system exceed a certain critical value, and the ratio of biomass consumed by zooplankton is greater than the sum of the zooplankton mortality rate and the rate of toxic substance production by phytoplankton. Hereafter we assume that the conditions stated in (5) always hold.

Using the translationswe translate the positive equilibrium to the origin. Then the linearized system for system (3) near iswhereThe associated characteristic equation of system (7) iswhereIt is well known that the positive equilibrium of system (3) is locally asymptotically stable if all roots of (9) have negative real parts and is unstable if (9) has a root with positive real parts. Now, we shall discuss the distribution of the roots of (9). First, when , the characteristic equation becomesAccording to the Routh-Hurwitz criterion, is locally asymptotically stable if and only if

By Corollary 2.4 in the paper of Ruan and Wei [29], we know that if instability occurs for a particular value of time delay, a characteristic root of (9) must intersect the imaginary axis. Hence, suppose that is a purely imaginary root of (9). Inserting this into (9) and separating the real and imaginary parts, we obtainSquaring and adding these two equations, we havewhere and . Let ; then (14) becomesForm (15), we define a function . It is obvious that if , then When , the equation has two real roots:Noticing that and , then we introduce the following results, which have been proved by many authors [30, 31].

Lemma 1. For (15), we obtain the following results:(i)If , then (15) has no positive roots.(ii)If , then (15) has two positive real roots if and only if and .Hence, if the condition,is satiated, then (15) has two positive roots, denoted by , . Then (14) has two positive roots, namely, and We definewhere satisfiesFor convenience, let Let denote the root of (9) such thatThen we have the following lemma.

Lemma 2. If (17) holds, then

Proof. Differentiating both sides of (9) with respect to , we haveThenwhereAs , then
If (17) holds, (15) has two positive roots, denoted by , , and is the local minimum value. Assuming , we obtain and Hence, , , . This completes the proof.

On the basis of the above analysis, we have the following theorem.

Theorem 3. Suppose that (5) and (12) hold; then we obtain the following:(i)If , then the positive equilibrium of system (3) is locally asymptotically stable for all values of .(ii)If (17) holds, there exists a nonnegative integer , such that the positive equilibrium is locally asymptotically stable whenever and is unstable whenever Then system (3) undergoes Hopf bifurcation around for every ,

3. Direction and Stability of Hopf Bifurcation

From Theorem 3 (ii), we have obtained the conditions for the occurrence of Hopf bifurcation when , For convenience, we define , , and and In this section, we consider the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions by using the center manifold and normal form theories presented by Hassard et al. [32]. Let , , , , , and , where . Using Taylor series expansion about , system (3) can be written aswhere Then is Hopf bifurcation value of system (3).

Let According to the Riesz representation theorem, there exists a matrix of the bounded variation for such thatIn fact, we can selectwhere denote the Dirac delta function:For , we defineThen system (25) can be rewritten aswhere and for

For , we define the adjoint operator asand a bilinear inner productwhere . From the results in the previous section, we know that are the eigenvalues of . Thus, are also the eigenvalues of .

Let be the eigenvector of corresponding to the eigenvalue , and let be the eigenvector of corresponding to the eigenvalue . With the condition , we obtainwhere is an identity matrix of order 3; that is,Then we can obtainSimilarly, with or , we haveSolving the above equations, we obtainHence, Then we can choosesuch that and

In the following, we first compute the coordinates to describe the center manifold at according to the approach of Hassard et al. [32]. Let be the solution of (31) when . We define

On the center manifold , we havewhere and are local coordinates for the center manifold in the direction of and , respectively. In this paper, we only consider real solutions as is real if is real. Since , for the solution of (31), we obtainwhereFrom (31) and (41), we havewhereOn the center manifold , from (42)–(44), we obtainCombining (45)–(47) and (42) and comparing the coefficients, we haveSincenoting that , we haveAccording to the definition of , we obtainComparison of the coefficients with (44) yieldswhere , , and , , and are still unknown. Next, we compute and accurately.

For , we haveComparing the coefficients with (46), we obtainFrom the definition of and (48) and (54), we haveFor , solving (55) yieldsSimilarly, using the same method, we obtain and are three-dimensional vectors that can be determined by setting in

For (45), when , . Thus,whereAccording to (48), we obtainFrom (56), (57), and (60), we haveNow, all can be expressed in terms of parameters. Therefore, we can evaluate the following values:

According to the above analysis, we can obtain the following theorem about the properties of Hopf bifurcation.

Theorem 4. For , and defined as above, the properties of Hopf bifurcation at the critical value are as follows:(i)If , Hopf bifurcation is supercritical (subcritical).(ii)If , the Hopf periodic solutions are stable (unstable).(iii)If , the period of the bifurcation periodic solution of system (3) increases (decreases).

4. Numerical Simulations

In this section, we verify the theoretical results proved in previous sections using numerical simulations for the following parameter values:According to the analysis above, the existence and stability of the interior equilibrium for system (3) both change with the value of the external source of nutrients , as illustrated in Figure 1. In Figure 1(a), (the conversion frequency of stability switches for ) changes with . From (5), the condition for the existence of is . Region I in Figure 1(a) indicates that system (3) has no positive equilibrium if the external source of nutrients flowing into the system is less than the critical value, 0.1416. In region IV, is unstable when condition (12) is not satisfied. This shows that when is greater than a critical value, 0.755, is always unstable and a time delay will no longer affect system (3). When the value of is in region II, condition (17) is not satisfied and (14) has no positive roots. Hence, is stable only when . In region III, Theorem 3 (ii) holds and there exists a nonnegative integer such that is locally asymptotically stable when and it is unstable when To simulate clear stability switches, Figure 1(b) shows interval values for stability switches and Hopf bifurcation about with the values of selected as 0.355, 0.2991, 0.2948, 0.2819, 0.2776, and 0.2733 and . For example, when , Theorem 3 (ii) is satisfied. Then (14) has two positive roots, and. From (18), we obtainIn order to more clearly show the interval values for stability switches, we use the rectangles to represent the stable and unstable interval of positive equilibrium . Therefore, in Figure 1(b), it is easy to see that the stability switches of positive equilibrium occur with the change of time delay .

On the basis of the above analysis, we present examples in Figures 25. The equilibrium point is locally asymptotically stable in Figure 2. When passes through critical values, loses its stability, and bifurcating periodic solutions occur and are stable, as shown in Figures 3-4. In Figure 5, the bifurcating solutions disappear and chaos occurs when passes through critical values.

The numerical results in Figures 25 show that oscillating solutions are strongly affected by . To generalize the dynamical behavior of system (3) influenced by parameter , we make numerical simulation using a wide range of values of parameter to show the bifurcation diagram (see Figure 6(a)). Results for the nutrient density and phytoplankton density as a function of are similar and are not shown here. From Figure 6(a), the order-2 periodic solution appears in the system (3) by increasing parameter . Keeping increasing , chaos occurs. Subsequently, an order-3 periodic solution bifurcates from chaos. Finally, chaos occurs again with increase of . To investigate the influence of parameter on amplitudes of the oscillations observed when the positive equilibrium undergoes Hopf bifurcations, we simulate the stability switches of positive equilibrium in the component (see Figure 6(b)). In Figure 6(b), the green points are , and , respectively. The black solid lines denote that the positive equilibrium component is stable when . The blue and red curves denote the maximum and minimum values, respectively. The black dashed lines indicate that when , the equilibrium loses its stability and oscillations occur. Therefore, the numerical simulations agree with the theoretical predictions.

5. Conclusion

In this paper, we have studied the dynamics of a nutrient-plankton system with a time delay. We have investigated the stability of the positive equilibrium and the existence of Hopf bifurcation. By using center manifold theory and the normal form method, we determined the direction of Hopf bifurcation and the stability of bifurcating periodic solutions.

In detail, we have found that if some conditions are satisfied, the phenomenon of stability switches arises from system (3). When the time delay passes through some critical values, there exists a nonnegative integer such that the positive equilibrium switches times from stability to instability to stability and so on, and Hopf bifurcations occur at the positive equilibrium. The numerical simulations in Figure 1 indicate that the theoretical results are correct, and the number of stability switches changes with the value of the external source of nutrients . This result indicates that the external source of nutrients flowing into the system has important influence on the complex dynamics of system (3).

Numerical simulations also showed that as the time delay further increases, the periodic solutions disappear and chaos appears. All these results not only will help in further investigating the dynamics of pelagic ecosystem in theory but also are very useful to understand the complex phenomena really happening in marine ecosystem.

Conflict of Interests

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


This work was supported by the National Natural Science Foundation of China (Grant nos. 31570364 and 31370381).