Abstract

A nutrient-phytoplankton model with multiple delays is studied analytically and numerically. The aim of this paper is to study how the delay factors influence dynamics of interaction between nutrient and phytoplankton. The analytical analysis indicates that the positive equilibrium is always globally asymptotically stable when the delay does not exist. On the contrary, the positive equilibrium loses its stability via Hopf instability induced by delay and then the corresponding periodic solutions emerge. Especially, the stability switches for positive equilibrium occur as the delay is increased. Furthermore, the numerical simulations show that periodic-2 and periodic-3 solutions can appear due to the existence of delays. Numerical results are consistent with the analytical results. Our results demonstrate that the delay has a great impact on the nutrient-phytoplankton dynamics.

1. Introduction

Some phytoplankton, for example, Cyanobacteria, can form dense and sometimes toxic blooms in freshwater and marine environments, which threaten ecological balance, drinking water, fisheries, and even human health [1]. However, the mechanism, by which phytoplankton blooms occur, is currently not very clear, which contribute to the difficulty to prevent or mitigate the proliferation of phytoplankton blooms. These have stimulated lots of researches aiming to understand the growth mechanisms of phytoplankton.

In recent years, dynamics in phytoplankton growth have drawn increasing attention from experimental ecologists, as well as mathematical ecologists. Some results from experiments and field observations imply that many factors affecting the dynamics of phytoplankton growth are bound to exist, such as nutrient [2], light [3], temperature [4], iron supply [5], zooplankton [6]. Especially, due to the effects of limiting factors including temperature, light, and day length, it has been indicated by Rhee and Gotham [7] that the population dynamics of phytoplankton in aquatic environments can change with season, latitude, and depth. Among factors affecting phytoplankton growth, nutrient has been an essential element [810], mainly including nitrogen and phosphate. Results reported by Ryther [11] indicated that phytoplankton indeed consumes lots of nitrogen and phosphate in their growth process, but reducing the nitrogen content in aquatic cannot slow the eutrophication. Using data from 17 lakes, Smith [8] analysed the influence of ratio of total nitrogen to phosphorus on the growth of blue-green algae (Cyanophya) and showed that controlling the ratio can help us improve the quality of aquatic environment very well. Obviously, the production process of phytoplankton is more complex.

However, due to the complexity and nonlinearity of aquatic ecosystem, there are some difficulties in understanding nutrient-phytoplankton dynamics only depending on experiment or field observation, which makes it necessary to use models to provide quantitative insights into dynamic mechanism of phytoplankton growth. For different aquatic environments, we can use various modifications of the classical prey-predator models by introducing functional responses to model nutrient-phytoplankton dynamics [1214]. For example, Huppert et al. [15] describe the dynamics of nutrient-driven phytoplankton blooms by a simple model and identify, using the model analysis, an important threshold effect that a bloom will only be triggered when nutrients exceed a certain defined level. Additionally, most nutrient-phytoplankton models reveal that phytoplankton population and nutrient population can coexist at equilibrium globally under some conditions [16, 17]. However, Sherratt and Smith [18] have reported that a constant population density may not exist in reality because of the existence of some factors, such as noise and physical factors. Actually, experiments and field observations show that the changes of phytoplankton population density usually possess oscillatory behaviour [19, 20].

For the single cell phytoplankton species, in most studies of nutrient-phytoplankton models, it is usually assumed that the processes, such as conversion process of nutrient, in the dynamics of phytoplankton growth are instantaneous [1417, 1923]. It may be doubtful whether there exists the delay in the growth of phytoplankton over the large area or not. Yet, J. Caperon [24] studied time lag in population growth response of Isochrysis galbana, a phytoplankton species, to a variable nitrate environment by both experiments and model, and demonstrated the existence of delay in the growth of Isochrysis galbana. Hence, the delay may indeed exist in the phytoplankton growth, which means that it is necessary to consider delay in nutrient-phytoplankton models. An approach that has been attempted by researchers to model the dynamics of phytoplankton is the role of delay since delay appears as an important component in biosystems and ecosystems [2530].

Actually, growing evidence shows that there exists time lag in some conversion processes from one state to another in some systems, and delay is an important factor because it can affect the dynamics of these systems. Volterra [31] considered time delay in a prey-predator model first and found oscillatory behaviour for the spatial distribution. For a long time, it has been recognized that delays can give rise to destabilizing effect of the dynamics of systems, where periodic solutions, as well as chaos, may emerge [3235]. Models incorporating delays in diverse biological and ecological models are extensively studied [3642]. Especially, the characteristic equation with respect to the linearized system of delay differential equations plays a key role in dynamic analysis, by which we can obtain some information on the stability of equilibrium. In addition, using the normal form theory, one can carry out the bifurcation analysis, such as the direction and stability of periodic solutions arising through Hopf bifurcation [43, 44].

The main purpose of this paper is to consider the effects of multiple delays on the nutrient-phytoplankton dynamics. In [15], Huppert et al. presented a simple model to investigate effect of nutrient on phytoplankton blooms, and much better results are obtained. Here, this model is extended into a “two preys-one predator” type to describe nitrogen- phosphorus-phytoplankton dynamics, as follows:where , , and represent nitrogen, phosphorus, and phytoplankton population density at time , respectively; is the nitrogen nutrients input flowing into the system and is the phosphorus nutrients input flowing into the system; is the loss rate of the nitrogen nutrients, and is the loss rate of the phosphorus nutrients; is nitrogen nutrient uptake rate of phytoplankton, and is phosphorus nutrient uptake rate of phytoplankton; and denote the efficiency of nutrient utilization; and are time delay parameters; is the mortality rate of phytoplankton. Although the function, which describes nutrient uptake dynamics, is not a Michaelis-Menten function, but Lotka-Volterra type, Huppert et al. [15] have indicated that the Lotka-Volterra term is a good first approximation to the Michaelis-Menten type. From biological viewpoint, all parameters are nonnegative. , , and are continuous on , where and , , and .

The paper is organized as follows. In Section 2, we analyze the existence and stability of positive equilibrium in model (1) without delays. In Section 3, we discuss stability of positive equilibrium and Hopf bifurcation under five different cases for delay effect. Subsequently, the direction of bifurcation and the stability of periodic solutions arising through Hopf bifurcation are given in Section 4. In order to analyze further how delay effects influence nutrient-phytoplankton dynamics, a series of numerical simulations are carried out in Section 5. Finally, the paper ends with conclusion in Section 6.

2. Existence and Stability of Positive Equilibrium in Model (1) without Delays

In this section, it is presented first that the first octant is positive invariant in model (1) without delays and the following lemma holds.

Lemma 1. All the solutions of model (1) with initial conditions that initiate in are positive invariant in the absence of delays.

Proof. From the first equation of model (1), we haveHence, under .
Likewise, from the second equation of model (1), we have under .
In the absence of delays in model (1), from the third equation of model (1), if , it can be obtained thatObviously, all the solutions of model (1) without delays are positive invariant if the initial conditions initiate in .
Then, we complete the proof.

For model (1), it is obvious that the extinction equilibrium, , exists. Moreover, in order to discuss the existence of positive equilibrium, the following function is defined:and then we can obtainwhere is the positive root of (4).

For the function , we have when the condition, , holds, and then there is no positive equilibrium in model (1). Obviously, holds if , and then there exists a unique positive root with respect to , which means that there exists a unique positive equilibrium in model (1) under this condition. However, when , it can be verified directly that and , which implies that there is no positive equilibrium in model (1). Thus, summarizing these results, the following theorem can be obtained.

Theorem 2. If holds, then there exists a unique positive equilibrium in model (1); otherwise, there is no positive equilibrium in model (1) if holds.

Letting the unique positive equilibrium be , then the following theorem holds for model (1) without delays.

Theorem 3. If the unique positive equilibrium exists in model (1) without delays, then it is globally asymptotically stable.

Proof. We construct a Lyapunov function, as follows:In the model (1) without delays,Obviously, holds under existence of positive equilibrium and holds if and only if and . The largest invariant subset of the set of the point where is . Therefore, according to LaSalle’s theorem, is globally asymptotically stable.
Then, we complete the proof.

Letting the extinction equilibrium be , then we can obtain the following theorem in model (1) in the absence of delay.

Theorem 4. In the absence of delays, let , so that(i)if , then the extinction equilibrium is locally asymptotically stable;(ii)if , then the extinction equilibrium is unstable;(iii)if , then the model (1) undergoes transcritical bifurcation at the extinction equilibrium .

Proof. For simplicity, letThe Jacobian matrix at is The eigenvalues are
Obviously, if , then the extinction equilibrium is locally asymptotically stable.
If , then the extinction equilibrium is unstable.
When , the Jacobian matrix at is Then has a geometrically simple zero eigenvalue with right eigenvector and left eigenvector .
NowandAccording to [45], the model (1) undergoes transcritical bifurcation at the extinction equilibrium in the absence of delays.
Then, we complete the proof.

Actually, when holds, the positive equilibrium does not exist, and the extinction equilibrium is globally asymptotically stable. Then, the following theorem holds.

Theorem 5. If holds, then the extinction equilibrium is globally asymptotically stable.

Proof. We construct a Lyapunov function, as follows:In the model (1) without delays,Obviously, holds if . Therefore, the extinction equilibrium is globally asymptotically stable when .
Then, we complete the proof.

3. Local Stability Analysis and the Hopf Bifurcation

In this section, we first state the following positive invariant theorem.

Lemma 6. All the solutions of model (1) with initial conditions that initiate in are positive invariant.

Proof. We consider a noncontinuable solution of model (1); see [46], defined on , where . Then we can use the method from [47] to prove that, for all , , , and . Suppose that is not true. Then, there exists such that, for all , , , and and either , , or . According to Lemma 1, for all , we haveandAs is defined and continuous on , there is a such that, for all ,andTaking the limit, as , we can getandwhich contradicts the fact that either , , or . Thus, for all , , , and .
Therefore, all the solutions of model (1) are positive invariant if the initial conditions initiate in .
Then, we complete the proof.

Next, we will discuss the stability of the unique positive equilibrium and existence of Hopf bifurcation in model (1) for five different cases: , ; , ; ; , ; , .

According to Theorem 2, let , , ; the linearized form of model (1) can be obtained as follows:where

Then, we obtain the associated characteristic equation of model (21) as follows:where

Case 1. , .
Due to , , (23) becomes Assuming is the pure imaginary root of (25), then the following can be obtained:ThenNow, we define a function as follows:(i)If (): holds, then, (28) has at least one positive root. Without loss of generality, we denote , , and as the roots of (28); hence , , , , if .(ii)If : holds, let and . When , then (28) has no positive roots. However, when , has two real roots, denoted asObviously, . If and holds, then (28) has two positive real roots if and only if and . In addition, we denote two positive roots of (28) as and ; then, (27) has two positive roots, namely, and . Furthermore, we can have the following results.

Proposition 7. (i)If () holds, then (28) has at least one positive root.(ii)If () and holds, then, (28) has no positive root.(iii)If () and holds, then, (28) has two positive roots if and only if , and . Then, according to (26), the critical delay can be obtained as follows:Letting , from [37] we know that also needs to be proved. Differentiating left side of (26) with respect to , we haveso that the following can be obtained:where

If (i) in Proposition 7 and () hold, then we have . However, if (iii) in Proposition 7 holds, assuming , then we obtain and . Hence, we have , , and . Therefore, we have the following results.

Theorem 8. For model (1) with , ,(i)If () and () both hold, then the positive equilibrium is locally asymptotically stable for and Hopf bifurcation occurs at .(ii)If (ii) in Proposition 7 holds, then the positive equilibrium is locally asymptotically stable for all .(iii)If (iii) in Proposition 7 holds, there exists a nonnegative integer , such that the positive equilibrium is locally asymptotically stable whenever and is unstable whenever . Then, model (1) undergoes Hopf bifurcation around at every , .

Case 2. , .
Since , , (23) becomes Similar to Case 1, let be the pure imaginary root of (34); then we obtainThat is,Letting , we define the following function:(i)If : holds, then (37) has at least one positive root. Without loss of generality, we denote , , as the roots of (37); then , , if .(ii)If : holds, let and . When , then (37) has no positive roots. However, when , has two real roots, denoted as Obviously, . If and holds, then (37) has two positive real roots if and only if and . In addition, we denote two positive roots of (37) as and ; then (36) has two positive roots, namely, and . Furthermore, we can obtain the following results.

Proposition 9. (i)If () holds, then (36) has at least one positive root.(ii)If () and holds, then, (36) has no positive root.(iii)If () and holds, then, (36) has two positive roots if and only if and .Then the critical delay can be derived by (35):Let . Differentiating left side of (34) with respect to , we obtainHence, we obtain the following:where

If (i) in Proposition 9 and : both hold, then is obtained. However, if (iii) in Proposition 9 holds, assuming , we obtain and . Then, we have , , and . Therefore, we have the following theorem.

Theorem 10. For model (1) with , , (i)If and hold, then, the positive equilibrium is locally asymptotically stable for and Hopf bifurcation occur at .(ii)If (ii) in Proposition 9 holds, then the positive equilibrium is locally asymptotically stable for all .(iii)If (iii) in Proposition 9 holds, then there exists a nonnegative integer , such that the positive equilibrium is locally asymptotically stable whenever and is unstable whenever . Then, model (1) undergoes Hopf bifurcation around for every , .

Case 3. .
When , (23) becomes Leting be the pure imaginary root of (43), thenthat is,Let and define the following function:From (46), we can clearly see that ; hence (46) has at least one positive root. Without loss of generality, we denote , , and as the roots of (46); then we have , , if holds. Hence, the critical delay can be derived by (44):Let . Differentiating left side of (43) with respect to , we obtainHence, we obtain the following:whereIf (): holds, then is obtained; hence, we have the following result.

Theorem 11. For model (1), when , if () holds, then the positive equilibrium is locally asymptotically stable for and Hopf bifurcation occurs at .

Case 4. , .
Under this case, is considered as a parameter. The same as Case 1, let be the root of (23); then we have the following:whereAccording to (51), the following holds:whereAssuming (): (53) has finite positive root and denoting as , , then the critical value can be represented as follows:Let , . Differentiating left side of (23) with respect to , the following is obtained:and then, we havewhereSupposing (): , then we have the following.

Theorem 12. For model (1), when and , if both and hold, then the positive equilibrium is locally asymptotically stable for and Hopf bifurcation occurs at .

Case 5. , .
Since , , we consider as a parameter. The same as Case 4, letting be the root of (23), we obtain:whereAccording to (59), the following holds:whereSupposing (): (61) has finite positive root , , then we obtain Assuming , and differentiating left side of (23) with respect to , therefore, the following is obtained:Hence, we havewhereSupposing : holds, then we obtain the following theorem.

Theorem 13. For model (1), when and , if both and hold, then the positive equilibrium is locally asymptotically stable for and Hopf bifurcation occurs at .

4. Properties of Periodic Solution

In this section, we will discuss the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions under Case 4 by using normal form method and center manifold theorem [43], and methods of other four cases are similar to Case 4. Assuming , , and Hopf bifurcation occurs at in model (1) when .

Let , , , , and , and we also denote , , and as , , and . Then, model (1) could be rewritten as follows in :where . and are given as follows: where

According to the Riesz representation theorem, we know that there exists a function of bounded variation for such that , for all . Choosing

For , we define Then, model (1) can be rewritten asFor , the adjoint operator of is defined as follows:Associated with a bilinear inner productwhere and we know that are the eigenvalues of and .

Choose to be the eigenvector of corresponding to the eigenvalue and to be the eigenvector of corresponding to the eigenvalue . By computation, we obtain

Besides, from (74) we have such that and .

Then, according to [44], we obtain the following relevant parameter, which helps to determine the direction and stability of Hopf bifurcation:and where and can be determined by the following, respectively:where

Then, we can compute the following values:This determine the properties of bifurcating periodic solutions and the Hopf bifurcation at . That is, (i) determines the direction of the Hopf bifurcation. Specifically, when , the Hopf bifurcation is supercritical (subcritical).(ii) determines the stability of the bifurcating periodic solutions; when , the bifurcating periodic solution is stable (unstable).(iii) determines the period of the bifurcating periodic solutions; when , the period of bifurcating periodic solution increases (decrease).

5. Numerical Simulations

Due to the complexity of model (1), we perform some numerical simulations in this section to investigate further how the delay influences dynamics in model (1). The following parameter values are used , , , , and . Other parameters are chosen as control parameters.

According to the standard linear analysis, when is equal to zero, the analysis reveals that the parameter plane is divided into four parts (see Figure 1(a)). In Figure 1(a), before reaching black dashed line, there exists in model (1) such that the unique positive equilibrium loses its stability when the condition, , holds. When the locus of is between black dashed line and green dashed line, the stability switches for positive equilibrium do not exist although (28) has two positive roots, which means that there exists in model (1) such that the unique positive equilibrium loses its stability when the condition, , holds. However, the stability switches for positive equilibrium emerge when is beyond green dashed line but it does not reach blue zone. When the locus of is in blue zone, the unique positive equilibrium is always stable, which suggests that cannot influence the stability of the positive equilibrium. When equals zero, the similar results for parameter plane are shown in Figure 1(b), but the sequence is reversed. Additionally, according to results in Section 4, we calculate the values of , , and at with , and the corresponding results are shown in Figure 1(c), where we can find that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable; especially, the period of the bifurcating periodic solutions increases as increases. For other cases of and , the same procedures with respect to calculations of , , and can be performed like Figure 1(c).

As examples corresponding to stability of the positive equilibrium with , taken and in Figure 1(a), respectively, the corresponding numerical solutions are shown in Figure 2. Obviously, the positive equilibrium is stable because is below (see Figure 2(a)). In contrast, due to beyond , a periodic solution exists (see Figure 2(b)). Furthermore, set , then we have . Taken , and in Figure 1(a), respectively, the corresponding numerical solutions are shown in Figure 3. Obviously, the positive equilibrium is stable when and (see Figures 3(a) and 3(c)), but the positive equilibrium is unstable when (see Figure 3(b)), which means that the positive equilibrium can gain its stability again for . In Figure 3, the same initial values are applied, and other parameter values except for are also identical. Clearly, the delay is the principal factor giving rise to the difference among (a), (b), and (c).

Numerical solutions in Figure 3 suggest that the stability switches induced by delay may exist. Hence, the bifurcation diagram in parameter plane is given (see Figure 4(a)). For case , there exists a such that the positive equilibrium with respect to is stable. For case , there exists a such that the positive equilibrium with respect to is stable. Additionally, when , Figure 4(a) shows that exists such that the positive equilibrium with respect to is stable. Likewise, when , Figure 4(a) also display that exists such that the positive equilibrium with respect to is stable. Significantly, Figure 4(a) demonstrates that the stability switches for positive equilibrium with respect to emerge when below is fixed.

Figure 4(b) depicts the dependence of stability of positive equilibrium on delay in the fully nonlinear regime for in sequence and , which is consistent with results in Figure 4(a). However, when is fixed, Figure 4(c) shows that the stability switches emerge with increases. Especially, periodic-2 solutions exist for some values of (see magenta zone in Figure 4(c)). As an example of periodic-2 solution existence, taking , a phase and a time-series are given in the inner of Figure 4(c). Moreover, Figure 4(d) shows that there exist periodic-3 solutions in model (1). According to results in Section 2, the positive equilibrium is globally asymptotically stable in the model (1) without delay if it exists. Obviously, the results shown in Figure 4 are induced by delay.

In Figure 4(a), we can find that the number of intervals corresponding to stability of positive equilibrium for small values of equals 3. However, when is beyond dashed line (), the number of intervals is 2. So the number of intervals for stability switches may be different for diverse . Accordingly, we calculate the number of intervals with respect to parameter , as shown in Figure 5(a). Figure 5(b) shows that there exist 3 stable intervals when and , which is an example of Figure 5(a). Evidently, parameter can remarkably influence the number of intervals for stability of positive equilibrium.

6. Conclusions

In this paper we proposed a nitrogen-phosphorus-phytoplankton model with multiple delays. The analysis focused on the effect of delay on nutrient-phytoplankton dynamics. In the absence of delay, theoretical analysis indicated that the unique positive equilibrium is globally asymptotically stable in model (1) if it exists. Deng et al. [48] also studied a nitrogen-phosphorus-phytoplankton model without delay, where Holling II function was employed to describe the nutrient uptake dynamics of phytoplankton. Although the function modelling the nutrient uptake dynamics of phytoplankton is different, they get the same results. These results mean that the nutrient-phytoplankton ecosystem will approach the stable equilibrium. However, it has been reported [18] that a constant population density may not exist because of the existence of some factors including noise, interval factors, and physical factors. And ecological studies [49, 50] also criticized this idea of “the balance of nature.” Actually, the existence of nutrient-plankton oscillations has been detected by laboratory experiments and field observation [19, 20]. Additionally, Benincà et al. [49] present the first experimental demonstration of chaos in a long-term experiment with a complex food web, where the food web was consisted of bacteria, several phytoplankton species, herbivorous and predatory zooplankton species, and detritivores. And they also find that the community moved back and forth between stabilizing and chaotic dynamics during the cyclic succession, and their findings provide a field demonstration of nonequilibrium coexistence of competing species through a cyclic succession at the edge of chaos [50]. These reports support that the nonequilibrium dynamics, such as oscillations and chaos, can exist in reality.

In the present paper, we find that the unique positive equilibrium may lose its stability via Hopf bifurcation when delay appears, and then a periodic solution emerges, which means that nutrient-phytoplankton oscillation occurs. Obviously, the factor giving rise to nutrient- phytoplankton oscillation is delay in our studies. And the period and the stability of the bifurcating periodic solutions with respect to delay are discussed by using center manifold argument and normal form theory. In fact, instability induced by delay in nutrient-plankton model has been studied widely, and many studies indicate that the equilibrium is always unstable when delay is beyond a critical value [25, 29, 51, 52]. Yet, it should be emphasized in the present paper that the stability switches induced by delay can occur under some conditions.

Moreover, numerical simulations showed how the delay influences nutrient-phytoplankton dynamics. Numerical results for model (1) in the fully nonlinear regime are consistent with the linear analysis. In numerical simulations, we found that delay indeed gives rise to the emergence of stability switches for the positive equilibrium. Yet, the numerical results show that the parameter intervals for stability switches may depend on other parameters as well, e.g., . Additionally, numerical results also indicated that periodic-2 solutions and periodic-3 solutions can emerge under some conditions for delay, which means that complex dynamics induced by delay exist in model (1). From biological viewpoint, the existence of periodic solutions implies that the fluctuations exist in density of phytoplankton population; that is, nutrient-phytoplankton oscillation emerges. Especially, by Li and York’s theory, periodic-3 solution implies chaos, which means that chaotic density fluctuations can display a variety of different periodicities and the long-term prediction of phytoplankton density can be fundamentally impossible. The chaotic density fluctuations do not contribute to the control of phytoplankton bloom. Consequently, the importance of the present paper is not the precision with which it predicts specific events for phytoplankton blooms but its contribution to the studies on how the delay influences nutrient-phytoplankton dynamics.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The author declares that there are no conflicts of interest.

Acknowledgments

This work was supported by the Zhejiang Provincial Natural Science Foundation of China under Grant no. LQ18C030002, the National Natural Science Foundation of China (Grant nos. 61871293, 31570364, and 41876124), and the Zhejiang Provincial Natural Science Foundation of China under Grant no. LY16B070008.