Abstract
A delay differential system is investigated based on a previously proposed nutrient-phytoplankton model. The time delay is regarded as a bifurcation parameter. Our aim is to determine how the time delay affects the system. First, we study the existence and local stability of two equilibria using the characteristic equation and identify the condition where a Hopf bifurcation can occur. Second, the formulae that determine the direction of the Hopf bifurcation and the stability of periodic solutions are obtained using the normal form and the center manifold theory. Furthermore, our main results are illustrated using numerical simulations.
1. Introduction
Phytoplankton plays a very important role as the first trophic level in aquatic ecosystems. To describe the complex dynamics of phytoplankton populations, the dynamic relationship between phytoplankton and nutrients has been investigated theoretically for a long time, as well as experimentally. Since the pioneering work of Riley et al. [1], various nutrient-phytoplankton models have been proposed and analyzed [2–9]. Several of these models have been shown to predict phytoplankton dynamics successfully in specific situations.
The model proposed by Taylor et al. [10] in 1986 describes the nutrient-dependent growth of a single phytoplankton population by considering sinking and variable vertical mixing: where and are the concentrations of the nutrient and phytoplankton, respectively; the specific growth rate of the phytoplankton is taken to be the product of two parts, a dependence () on incident light () and a function () of the concentration of a single nutrient; is the death rate of phytoplankton; is the regeneration efficiency; is the sinking rate of phytoplankton; is the turbulent diffusion coefficient. The diffusion coefficients for phytoplankton and nutrients are assumed to be the same for simplicity. It is well known that the abundance of the phytoplankton population is affected by many environmental factors, such as the water temperature, salinity, and sunlight intensity [11]. In system (1), is the light intensity and numerical simulation indicates that the light intensity can affect the result.
In this study, we consider an approximated model of system (1) with delay effect as a model of a layer of phytoplankton growing over a pool of nutrients: where is the specific growth rate of phytoplankton; is the thickness of the layer; is the turbulent diffusion coefficient; is a positive delay, that is, the time required to convert nutrients into phytoplankton; and are the concentrations of nutrients and phytoplankton below the layer, respectively. We assume that all parameters are positive, except that is negligible and equal to zero.
Taylor et al. [10] described the dynamics of the above system without considering the effect of delay, mainly using numerical methods. Pardo [12] conducted a mathematical study of the local and global stability of the equilibria in the same system. He proved the positivity and boundedness of the solutions, which made sure that the model is biologically sound. He found that the interior equilibrium point of the system is locally and globally asymptotically stable if it exists and that the boundary equilibrium point is also globally asymptotically stable if the system has only one equilibrium point.
Time delay is known to play important roles in biological dynamical systems, which have been studied by many researchers in recent years [13–27]. Our aim is to investigate how time delay may affect the system (2). Thus, we use the delay as a bifurcation parameter.
The remainder of this paper is organized as follows. In Section 2, we consider the local stability of the equilibria and the condition where Hopf bifurcation can occur based on the characteristic equation. In Section 3, we derive an explicit algorithm to determine the direction of the Hopf bifurcation and the stability of the periodic solutions. The results of numerical simulations are presented to support the theoretical results in Section 4. The paper ends with a brief conclusion.
2. Stability Analysis of the Equilibria
In this section, we mainly consider the existence and stability of the nonnegative equilibria of system (2). The equations for the equilibria are as follows:
There are two solutions: , if is not equal to zero and exists, where
2.1. Local Stability of Equilibrium
We linearize (2) about to obtain the linear system
There are two eigenvalues where is always negative and only if . Therefore, we have the following theorem.
Theorem 1. The equilibrium is stable if and unstable if .
We recall that is equivalent to , which is the existing condition of the unique interior equilibrium. Thus, is stable only if does not exist.
2.2. Local Stability and the Hopf Bifurcation of Equilibrium
Similarly, we linearize (2) about to obtain the linear system
The corresponding characteristic equation is or where
When , the characteristic equation is and the two eigenvalues satisfy which indicates that is locally asymptotically stable when .
When , we assume that is a root of (8). Substituting this in (8) we have or
Let and , and hence (13) is equal to
Hence,
Let , and then .
However, we can solve the first formula in (13) to obtain
Substituting in (8) and taking the derivatives with respect to , we have
For ,
This implies that all of the roots cross the imaginary axis at from left to right as increases. Therefore, the transversality condition holds.
Based on the above, we have the following theorem.
Theorem 2. The interior equilibrium exists only if , and it is stable for but unstable for . Furthermore, system (2) undergoes a Hopf bifurcation at the positive equilibrium for , .
3. Direction and Stability of Hopf Bifurcations
In this section, we consider the direction, stability, and period of the periodic solutions from the steady state using the method introduced by Hassard et al. [28]. Let , , , , and , , where is a real number.
Then (2) can be written as that is, where
According to the Riesz representation theorem, there exists a matrix , which is bounded for such that
In particular, we can select where
For , we define
Then, (19) can be rewritten as
For , we define the adjoint operator as
For and , we define
Suppose that and are the eigenvectors that correspond to and , respectively. Let
With the condition , that is, we can obtain
Similarly, let and with that is, we can obtain
Given the condition , and since we obtain
In the following, we compute the coordinates to describe the center manifold at according to [28]. We define
On the center manifold , we have where
Then, where Based on (26) and (38), we have where
Comparing the coefficients, we have
Note that
Thus, according to the definition of , we can easily obtain
In the final formula, , , , and are still unknown. Thus, we need to compute and accurately, as follows.
Consider
Comparing the coefficients, we have
From (45), we can obtain
Solving the formula, we have
Using the same method, we can obtain where and are both two-dimensional vectors.
For , ; thus
According to (45) and (46), we have
Substituting (56) and (57) into (53) and (54), respectively, we can obtain
From the above, we already know and , and is also expressed using these parameters. Thus, we can compute the following values:
Theorem 3. , , and are defined above, so 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 Pardo [12], if there exists an interior equilibrium in system (2) without a delay effect, then the equilibrium is always globally asymptotically stable. According to our study, however, the equilibrium is not always stable due to the delay effect. That is, an oscillation will occur in system (2) under some conditions. These conditions are obtained using an analytical technique. To verify the validity of the results, a series of numerical results is given in this section. In the following, the parameters are fixed, excluding , , and . The three parameters are taken as control parameters because we want to know how the concentration of nutrients below the layer and the thickness of the layer, but mainly the delay, affect the system.
Consider
In this study, we are interested in the interior equilibrium, though the interior does not always exist. Thus, the field where there exists an interior equilibrium is given when , change, as shown in Figure 1(a). In Figure 1(a), both the interior equilibrium and boundary equilibrium exist in zone I (yellow), whereas there is only a boundary equilibrium in zone II (green), which is separated by the blue line. Based on the theoretical results, an oscillation will occur when the delay is larger than (the critical value). The critical value is determined by the parameters of system (2) except , and the relationship of with respect to , is shown in Figure 1(b). In Figure 1(b), when is fixed, will decrease with the increase of . If is fixed, decreases first with the increase in (), before it decreases. It should be noted that does not always exist when and is much smaller, because of the reason illustrated in Figure 1(a).
(a)
(b)
Based on Figure 1, the parameter is taken as and is set as , so is obtained according to arccos . To investigate the effect of a delay on system (2), we set as and initially. The numerical solutions for phytoplankton are shown in Figure 2(a). In Figure 2(a), the solution eventually converges to the interior equilibrium, but the solution with the delay oscillates initially. The interior equilibrium is globally asymptotically stable in system (2) without delay, so the solution must converge to the interior equilibrium. If the solution with a delay converges to the interior equilibrium because , an oscillation does not appear. In the following, we set as and . The numerical solutions for phytoplankton are shown in Figure 2(b). Hopf bifurcation occurs at , so an oscillation appears because . From Figure 2(b), it is obvious that an oscillation occurs.
(a)
(b)
Furthermore, to investigate the relationship between nutrients and phytoplankton, the numerical solutions for nutrients and phytoplankton are shown in Figure 3, where , , and . Figure 3(a) shows that the solutions oscillate. Thus, the relationship between nutrients and phytoplankton is mutually constrained. The phase of system (2) is shown in Figure 3(b).
(a)
(b)
To determine how , , and affect system (2), the bifurcation of the stable state is considered and Figure 4 shows the results. In Figure 4(a), the parameters , are fixed, and the delay varies. The solid line denotes the stable state, the symbol “” is the Hopf bifurcation point, and the dashed line represents the unstable state. There is only one stable state when is smaller than . However, the stable state becomes unstable when is larger than . Next, a periodic solution arises from the stable state, which is stable when . It is obvious that a delay destroys the stability of the equilibrium. In addition, the delay is fixed, where . We consider the effects of , on system (2), as shown in Figure 4(b). The solid line denotes the stable state, the symbol “” is the Hopf bifurcation point, and the dashed line represents the unstable state. The red line and the blue line both denote the equilibrium (phytoplankton), while the magenta line represents the minimum amount of phytoplankton. There are two Hopf bifurcation points when and , but there is only one when and . When , the equilibrium is stable at first, but the stable state becomes unstable when reaches and a Hopf bifurcation occurs, before a stable periodic solution emerges. When reaches , another Hopf bifurcation occurs and the stable periodic solution disappears, before the equilibrium becomes stable again. In particular, the minimum periodic solution decreases initially with the increase in , before it increases. When and , the minimum periodic solution increases with the increase in and the unstable equilibrium becomes stable via a Hopf bifurcation; that is, the oscillation disappears.
(a)
(b)
5. Conclusion
In this study, we investigate a biological system and consider the effect of time delay. Our results show that the time delay plays a vital role in system (2). Model (1), which was proposed by Taylor et al. [10], includes nutrient recycling without any nutrient loss and vertical variation. Using a numerical simulation, Taylor discussed how the light intensity, diurnal variation in light, and diurnal variation in turbulence affected the system. This is a very realistic model and it is useful for studying the dynamic relationship between nutrients and phytoplankton. However, because it considered too many realistic factors, the model dynamics were very complicated and difficult to discern in the full model. Therefore, an approximated model was proposed, which considered a layer of phytoplankton growing over a pool of nutrients. The model was the same as (2) but it lacked the time delay. Similar to the analysis of system (1), the condition where an oscillation occurred was considered based on a theoretical analysis, including how factors such as the phytoplankton loss rate affected the system based on a numerical simulation. Most of this analysis was based on the results obtained using a numerical method.
Pardo [12] focused on the approximated model in 2000. In a mathematical study, Pardo proved that if the interior equilibrium point exists, it is always globally asymptotically stable. The boundary equilibrium point is also globally asymptotically stable if the system has only one equilibrium point.
In our study, we consider the model with delay as (2). Our theoretical analysis shows that the boundary equilibrium is still stable if , but the unique positive equilibrium will lose its stability if the time delay exceeds a critical value and an oscillation occurs. Furthermore, we also prove that system (2) undergoes a Hopf bifurcation with a specific delay. These results are verified using numerical simulation based on published parameters [10] to ensure that they are realistic. Using a computer simulation, we find that the concentration of phytoplankton will increase to a large value within a short period of time because of the effect of delay, which is very similar to the phytoplankton blooms found in nature. Thus, the time delay may be a reason for phytoplankton blooms, but more research is needed to confirm this hypothesis.
The limitations of our study are that we only consider a discrete delay in the phytoplankton increase and we assumed that the recycling of nutrients is instantaneous. A model with delayed recycling would be more complicated but more similar to reality. Moreover, the model we considered is only an approximate model based on (1) and we consider that is a constant. Thus, we did not consider the influence of the light intensity, although the light intensity is depth independent. Therefore, further research is required to consider these aspects fully.
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. 31370381), 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).