Stability and Hopf Bifurcation Analysis of a Nutrient-Phytoplankton Model with Delay Effect

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.


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 nutrientphytoplankton models have been proposed and analyzed [2][3][4][5][6][7][8][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; (0 < ≤ 1) is the regeneration efficiency; V 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; 0 and 0 are the concentrations of nutrients and phytoplankton below the layer, respectively. We assume that all parameters are positive, except that 0 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.
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.

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: 0 ( 0 , 0), * ( * , * ) if is not equal to zero and 0 > * exists, 2.1. Local Stability of Equilibrium 0 . We linearize (2) about 0 to obtain the linear system There are two eigenvalues where − /ℎ is always negative and 0 − − V/ℎ − /ℎ > 0 only if 0 > * . Therefore, we have the following theorem.
The corresponding characteristic equation is When = 0, the characteristic equation is and the two eigenvalues satisfy which indicates that * is locally asymptotically stable when = 0. When ̸ = 0, we assume that = ( ) + ( ) is a root of (8). Substituting this in (8) we have Let ( ) = 0 and ( ) = 0 > 0, and hence (13) is equal to Hence, Let Δ = 4 + 4 2 > 0, and then 0 = √ (− 2 + √ Δ)/2. Substituting ( ) in (8) and taking the derivatives with respect to , we have 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.

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

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 0 , ℎ, 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 (60) 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 0 , ℎ 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 0 (the critical value). The critical value 0 is determined by the parameters of system (2) except , and the relationship of 0 with respect to 0 , ℎ is shown in Figure 1(b). In Figure 1(b), when ℎ is fixed, 0 will decrease with the increase of 0 . If 0 is fixed, 0 decreases first with the increase in ℎ (0 < ℎ < 30), before it decreases. It should be noted that 0 does not always exist when 0 < ℎ < 30 and 0 is much smaller, because of the reason illustrated in Figure 1(a).
Based on Figure 1, the parameter 0 is taken as 10 and ℎ is set as 1, so 0 is obtained according to 0 = (1/ 0 ) arccos ( 2 0 / ). To investigate the effect of a delay on system (2), we set as 0 and 2 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 < 0 = 3.45, an oscillation does not appear. In the following, we set as 0 and 4. The numerical solutions for phytoplankton are shown in Figure 2(b). Hopf bifurcation occurs at = 0 , so an oscillation appears because > 0 = 3.45. From Figure 2(b), it is obvious that an oscillation occurs.
Furthermore, to investigate the relationship between nutrients and phytoplankton, the numerical solutions for nutrients and phytoplankton are shown in Figure 3, where 0 = 10, ℎ = 1, and = 3.6. 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).
To determine how 0 , ℎ, and affect system (2), the bifurcation of the stable state is considered and Figure 4 shows the results. In Figure 4(a), the parameters 0 , ℎ 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 0 . However, the stable state becomes unstable when is larger than 0 . Next, a periodic solution arises from the stable state, which is stable when 0 < < 5. It is obvious that a delay destroys the stability of the equilibrium. In addition, the delay is fixed, where = 6. We consider the effects of 0 , ℎ 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 0 = 7 and 1 < ℎ < 7, but there is only one when 0 = 10 and 0 = 12. When 0 = 7, 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 0 = 10 and 0 = 12, the minimum periodic solution increases with the increase in ℎ and the unstable equilibrium becomes stable via a Hopf bifurcation; that is, the oscillation disappears.

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 0 < * , 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.