4Environmental Engineering Program, University of Northern British Columbia, Prince George, British Columbia, Canada

Academic Editor: Baltazar Aguirre-Hernandez

Received05 Oct 2018

Revised11 Dec 2018

Accepted09 Jan 2019

Published03 Feb 2019

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 [8–10], 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 [12–14]. 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 [14–17, 19–23]. 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 [25–30].

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 [32–35]. Models incorporating delays in diverse biological and ecological models are extensively studied [36–42]. 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:where