Abstract
In this paper, we study a reaction-diffusion nutrient-phytoplankton model with nutrient recycling delay and toxin-producing phytoplankton. For the delay-free case, we prove that the release rate of toxic chemicals by the toxin-producing phytoplankton can destabilize the system and cause periodic oscillation. For the time-delayed case, we observe that nutrient recycling delay may bring about stability switches and Hopf bifurcation. Moreover, we derive the formula for determining the direction and stability of the bifurcating periodic solutions. Finally, we give some numerical simulations to support the theoretical analyses.
1. Introduction
Aquatic ecosystem has a significant effect on human’s living environment. It not only plays an important role in maintaining the global biological matter and water cycles but also undertakes the functions of water source, power source, transportation pollution purification place, and so on. Obviously, it is of great importance to investigate the dynamic behaviors of planktonic ecosystem. In the last few decades, lots of mathematical models have been constructed in terms of differential equations to describe the dynamical interactions between the dissolved nutrient and plankton population, see [1–5] and references cited therein. These models can deepen our knowledge of aquatic ecosystem cycle. Some interesting dynamical behaviors have been investigated, such as uniform persistence, asymptotic stability, Hopf bifurcation, Turing instability, and travelling waves.
In nature, the dead phytoplankton can be partly transformed into the nutrient. By taking into account the recycling factor, Chakraborty et al. [6] proposed the following nutrient-phytoplankton model:where and denote the concentrations of nutrient and phytoplankton population, respectively. All the coefficients are positive constants, is the constant nutrient input flow, is the nutrient uptake rate of phytoplankton, is the nutrient conversion rate of phytoplankton, is the per capita-loss rate of nutrient, is the per capita-mortality rate of phytoplankton, is the portion of the phytoplankton recycled back to the nutrient concentration, is the release rate of toxic chemicals by the toxin-producing phytoplankton population, and is the half-saturation constant. In view of the conversion efficiency in ecological system, we always assume that and throughout the paper.
For system (1), the limit cycle behavior is obtained by numerical simulations in [6]. It is found that inclusion of nutrient recycling enhances the chances of recurring bloom dynamics. After that, some expanded models were investigated by incorporating zooplankton population in [7–11].
As indicated in [12], no life is possible if spatial and temporal processes do not interact. Actually, the individuals can go anywhere in the spatial domain since a well-known fact about Fickian diffusion leading to a Laplacian term. In the lakes or oceans, dissolved nutrient and phytoplankton population can move with the tide. The spatial diffusion has been considered in modelling the biological and ecological systems, such as [13–15]. Reaction-diffusion equations are widely used for modelling chemical reactions, biological systems, population dynamics, and nuclear reactor physics. This kind of models can better describe the spatial and temporal interactions. So, it is more realistic to use the reaction-diffusion equations to model the aquatic system. Then, we have the following model with spatial diffusion:where and are the diffusion coefficients of nutrient and phytoplankton, respectively. The initial functions and are nonnegative and nonidentical to zero. The homogeneous Neumann boundary conditions mean that the nutrient and phytoplankton do not leave the living region. For simplicity, we specify the spatial location as an interval , where denotes the depth of the water column and can be changed into any depth through a suitable transformation. Moreover, some other nutrient-plankton models incorporating spatial effects are investigated in [16–21].
In models (1) and (2), it is assumed that the decomposition process of dead phytoplankton is instantaneous. However, it is not accurate enough. Usually, the interaction between nutrient and phytoplankton is not instantaneous and it does take time to transform dead phytoplankton into nutrient. Time delay can also bring rich dynamical behaviors, especially the stability switches, and hence it cannot be ignored [22–24]. Therefore, we further build the following model with nutrient recycling delay motivated by [25, 26]:where the nutrient recycling delay is the decomposition time of dead phytoplankton population, . The initial functions and are nonnegative and nonidentical to zero. The aim of this paper is to investigate the joint effects of spatial diffusion and time delay on the spatiotemporal dynamics of the nutrient-phytoplankton system.
The highlights of this paper mainly contain two aspects. First, we synthetically consider the nutrient recycling delay and spatial diffusion in modelling the interactions between nutrient concentration and phytoplankton population, and the model is more closer to the real world. Second, by regarding nutrient recycling delay as the bifurcation parameter, we find that time delay may induce stability switches and sustained oscillations, which can help to understand the cyclical burst of algae bloom.
The rest of this paper is organized as follows. In Section 2, we establish the existence of nonnegative constant equilibria. In Section 3, we analyze the stability of positive equilibrium in the absence of time delay, the stability switches, and Hopf bifurcation in the presence of time delay, respectively. In Section 4, the algorithm for determining the Hopf bifurcating periodic solutions is derived. In Section 5, some numerical simulations are presented to verify the theoretical results. Finally, a brief conclusion is given in Section 6.
2. Existence of Equilibria
In this section, we give some preliminary results. For convenience, we make the following assumptions: () . () . () and . () . () . () and for some .
Spatial diffusion and time delay do not affect the number and location of differential system, so systems (1)–(3) have the same nonnegative constant equilibria and the equilibria are the roots of the following algebraic equation:
Obviously, there is always the boundary equilibrium . Next, we establish the existence of positive constant equilibrium. From the second equation of (4), we have
Substituting it into the first equation of (4) leads tothat is,where , , , . Based on the Descartes’ rule of signs, equation (7) has one or three positive roots and assumption can guarantee the existence of unique positive real root of equation (7).
Thus, we can establish the existence of positive equilibrium.
Lemma 1. If () holds, then systems (1)–(3) have the unique positive equilibrium .
Remark 1. From Lemma 1, if the nutrient uptake rate of phytoplankton and release rate of toxic chemicals by toxin-producing phytoplankton population are appropriately large, and the net growth rates of nutrient and phytoplankton are appropriately small, then the unique positive equilibrium of system may exist.
3. Stability Analysis and Occurrence of Hopf Bifurcation
In the following, we analyze the stability of . Linearizing system (1) at , we havewhere
Then, the characteristic equation of (8) iswhere , , . Equation (10) can be turned intowhere
3.1. The Case of
If and , then characteristic equation (11) can be simplified towhich corresponds to ordinary differential system (1) and . By Routh–Hurwitz criterion, the positive equilibrium of (1) is asymptotically stable if the condition () is satisfied.
Next, we explore the influence of parameter on the stability of . Let , then is the potential Hopf bifurcation value. In addition, we have
From above, we can conclude the stability results of system (1).
Theorem 1. (i)If () and () are satisfied, then the positive equilibrium of (1) is asymptotically stable.(ii)If () and () are satisfied, then the positive equilibrium of (1) is asymptotically stable when and unstable when for . Conversely, of (1) is asymptotically stable when and unstable when for . Moreover, is the Hopf bifurcation value.
Remark 2. From Theorem 1, if the equilibrium state of phytoplankton is smaller than the half-saturation constant, then condition () can easily be satisfied and the nutrient and phytoplankton can achieve a balanced coexistence in the absence of diffusion. Otherwise, the release rate of toxic chemicals may induce the periodic variations of nutrient and phytoplankton.
Here, we mainly investigate the Turing instability induced by spatial diffusion. The linearization of system (2) at the positive equilibrium isand the corresponding characteristic equation iswhere
It is obvious that when . Then, the positive equilibrium of diffusive system (2) is stable if for any and unstable if for some . It is not difficult to determine the sign of by using the properties of parabola.
Following Turing’s idea in [27], different diffusion rates of two chemical reactants could sometimes destabilize an otherwise stable equilibrium of the reaction-diffusion system and lead to nonuniform spatial patterns. This kind of diffusion-driven instability is usually called Turing instability.
Theorem 2. Suppose () and () hold, we have the following statements:(i)If () is satisfied, then the positive equilibrium of (2) is locally asymptotically stable.(ii)If () is satisfied, then the positive equilibrium of (2) is unstable. It implies that Turing instability occurs.
Remark 3. From Theorem 2, the diffusion coefficients have significant effect on the stability of positive equilibrium in the presence of spatial diffusion. If the diffusion rate of phytoplankton is larger, then the aquatic system is stable. On the contrary, the system may be unstable.
3.2. The Case of
For , let be a root of (11). We have
Separating the real and imaginary parts can lead to
Squaring both sides of the above equations, we get
From and , it can be found that implies . Moreover, equation (20) has two positive real roots under the assumption (). Thus, characteristic equation (11) has two pairs of purely imaginary roots , where
By direct computation, we have and . Thus, .
In what follows, we shall verify the transversality condition.
Lemma 2. and .
Proof. Differentiating both sides of equation (11) with respect to yieldsThen,
Through the above analysis and motivated by Corollary 2.4 of Ruan and Wei in [28], we can establish the distribution of roots of characteristic equation (11) in the case of .
Lemma 3. If (), (), (), and () are satisfied, then we have the following:(i)If , then all roots of (11) have negative real parts for ; equation (11) has a pair of purely imaginary roots and other roots have negative real parts for ; equation (11) has at least a root with positive real part for .(ii)If , then there exists a positive integer such that all roots of (11) have negative real parts for , and equation (11) has at least one root with positive real part for .
As time delay changes, the stability of the equilibrium solution may also change. The stability may transit from stable to unstable at some critical value and may also transit from unstable to stable at some other critical value. Such phenomena are often referred to as stability switches, i.e., the state can be switched from stable to unstable and back to stable and so on just by progressive increase of the time delay. Then, we can obtain the stability switches and existence of Hopf bifurcation for in (3).
Theorem 3. If (), (), (), and () are satisfied, then the following conclusions are true:(i)If , then positive equilibrium of (3) is asymptotically stable for and is unstable for . System (3) undergoes Hopf bifurcations near when .(ii)If , then there exists a positive integer such that switches times from stability to instability to stability and so on such that is asymptotically stable for and is unstable for . System (3) undergoes Hopf bifurcation around for every .
4. Properties of Hopf Bifurcation
Here, we derive the Hopf bifurcation properties by using the center manifold theorem and normal form theory for partial functional differential equations in [29, 30]. For simplicity, we denote the Hopf bifurcation value by and the corresponding simply purely imaginary characteristic roots by , respectively.
Let , , , , and with , then is the Hopf bifurcation value of system (3). For simplicity, we drop the tilde and rewrite system (3) in the formwhere , , and , are given byand
From the previous section, we know that is a pair of simple purely imaginary eigenvalues of the following linear differential equation:
Due to the Riesz representation theorem, there exists a matrix function , whose elements are of bounded variation functions such thatwhere
For , , we define
Then, and are adjoint operators under the bilinear form
It is calculated that and are eigenvectors of and corresponding to and , respectively, where
Then, is the center subspace of system (3), and is the adjoint subspace.
Let , , and , . The complex-valued inner product on Hilbert space isfor , and , , , ,where . Then, the center subspace of system (27) at is
If , then we can obtain
The flow of system (24) on the center manifold can be written as
Moreover, for of (24) at , we have , where
Following the calculation procedures in [29, 30], we can getwhere
Thus, we can compute
Then, by the Hopf bifurcation theory in [29, 30], we have the following conclusions about the properties of the Hopf bifurcation.
Theorem 4. For system (3),(i) determines the direction of Hopf bifurcation: if , then the supercritical bifurcation bifurcates stable limit cycles; if , then the subcritical bifurcation bifurcates unstable limit cycles.(ii) determines the stability of bifurcating periodic solutions: the periodic solutions are orbitally asymptotically stable if and are unstable if .(iii) determines the period of the bifurcating periodic solutions: the period is monotonically increasing at the time delay when and is monotonically decreasing at the time delay when .
5. Numerical Simulations
In this section, we mainly verify the theoretical results obtained in the previous sections by some numerical examples. For system (3), we set , , , , , , , , , and .
Example 1. If and with the initial conditions and , then the positive constant equilibrium of system (2) is asymptotically stable, as shown in Figure 1.
Example 2. If and with the initial conditions and , then the positive constant equilibrium of system (2) is unstable and Hopf bifurcation occurs, see Figure 2.
For , we fix . Then, system (3) has the unique positive equilibrium solution . Moreover, equation (20) has two positive roots and , . According to Theorem 3, if time delay varies, then the stability of positive equilibrium of system (3) would change.
Example 3. For , we can compute , then and . From Theorem 4, we can know that the Hopf bifurcation is supercritical and the spatially homogeneous periodic solutions are stable. The numerical solutions with initial conditions and are presented in Figures 3 and 4.
Example 4. For , we can compute , then and . The Hopf bifurcation is supercritical and the spatially homogeneous periodic solutions are stable. The numerical solutions with initial conditions and are presented in Figures 5 and 6.
It can be found that nutrient recycling delay can bring about the stability switches in the nutrient-phytoplankton system.
6. Conclusion
In this paper, we have expanded the model proposed in [6] and investigated a novel reaction-diffusion nutrient-phytoplankton ecological model incorporating nutrient recycling delay. We proved that the release rate of toxic chemicals by the toxin-producing phytoplankton population can destabilize the positive constant equilibrium and may generate time-periodic solutions. We also revealed the effect of nutrient recycling delay on our model. By taking time delay as the bifurcation parameter, multiple stability switches of the positive constant equilibrium and Hopf bifurcation phenomena can be observed under some conditions, which contrast sharply with the case in the absence of time delay. The results can help to understand the periodic outbreak of algae.
It is observed that the model is more comprehensive by taking account of the zooplankton population, and rich dynamical behaviors have been presented in [3, 8, 10, 26], such as uniform persistence, Turing bifurcation, and global Hopf bifurcation. However, these models do not investigate the joint effect of spatial diffusion and time delay. In addition, in a reaction-diffusion system, when the Hopf bifurcation curve and Turing bifurcation curve intersect, the Turing-Hopf bifurcation will occur and generate spatiotemporal patterns, see [31–33]. We will do further research on spatiotemporal patters via Turing-Hopf bifurcation in the delayed reaction-diffusion model incorporating three components in the future following the aforementioned references.
Data Availability
No data were used in this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Research Foundation Project (12061033), Natural Science Research Key Project of Universities in Anhui Province (KJ2020A0002 and KJ2020A0016), Anhui Provincial Natural Science Funds (2108085MA04), and First-Class Course of AUFE (acylkc202008).