Abstract
We consider the mathematical formulation, analysis, and numerical solution of a nonlinear system of nutrient-phytoplankton, which consists of a series of reaction-advection-diffusion equations. We derive the critical conditions for Turing instability without an advection term and define the range of Turing instability with the change of nutrient concentration. We show that horizontal movement of phytoplankton could influence the system and that it is unstable when the horizontal velocity exceeds a critical value. We also compare reaction-diffusion equations with reaction-advection-diffusion equations through simulations, with spotted, banded, and crenulate patterns produced from our model. We found that different spatial constructions could occur, impacted by the diffusion and sinking of nutrients and phytoplankton. The new model may help us better understand the dynamics of an aquatic community.
1. Introduction
Environmental resources include clean air and water, forests, soil, climate, and many others. Ecological services exist within and are sustained by the ecosystem. They are generated by the continuous connection between organisms, populations, communities, and their physical and chemical environment [1]. Modeling the dynamics of phytoplankton is of great importance to many aspects of human interest, since phytoplankton provides the basis of the food chain in lakes, seas, and oceans [2]. In some cases lakes, reservoirs, and marine waters may experience plankton or algal blooms [3, 4]. Plankton blooms often lead to disorder of aquatic ecosystems, for example, the decrease of aquatic species, and then damage to the biological diversity. Therefore, research into plankton dynamics has become a topic of wide interest.
The historical work of Turing (1952), showed that complex spatial patterns can be produced from the reaction and diffusion of a number of chemicals [5], a phenomenon now known as diffusion-driven instability. A great deal of research has investigated spatial patterns in predator-prey systems. These arise through diffusion-driven instability and depend on significant differences between predator and prey diffusion coefficients [6, 7]. The effect of equal diffusion coefficients of spatial models has also been well studied with reaction-diffusion equations [8, 9].
Many researchers have investigated the predator-prey model among nutrients, phytoplankton, and zooplankton [10–20]. A series of studies simulated the dynamics of the plankton system. Zhang and Wang [21] studied a nutrient-phytoplankton-zooplankton model in an aquatic environment and showed that phytoplankton blooms may occur even if the nutrient input rate is low. Luo [22] investigated phytoplankton-zooplankton dynamics in periodic environments and showed that the eruption of algal blooms largely results from increasing of nutrient concentration. Sae-jie et al. [23] used a general forced dynamical model of nutrient-phytoplankton. They found that sinusoidal and periodic step function input of nutrients caused phytoplankton blooms with periodic behaviour and that changing the frequency of inputs produced blooms with a wide range of different dynamical behavior.
In recent years, numerous ecologists have paid increasing attention to spatial processes in practical contexts [24]. Klausmeier [25], for example, demonstrated a simple model of plant and water dynamics based on ecologically realistic assumptions, which showed that spatial patterns in ecological systems can result from both self-organization and amplification of an underlying heterogeneity. Other researchers have probed ecosystems, including vegetation [26–28] and phytoplankton systems. Asaeda et al. [29] developed a numerical model incorporating phytoplankton, submerged macrophytes, Potamogeton pectinatus L, which belongs to the Monocotyledon, their decomposition process, and the resulting nutrient dynamics in the overlying water, which investigated how external nutrient loading, water temperature, water depth, and retention time of water changed submerged macrophytes and phytoplankton biomass, phosphorous concentrations, and light penetration.
A large number of studies have simulated the processes of patchiness pattern formation [30, 31]. The research may be classified into two groups depending on what is considered as a major cause of these phenomena. One group concentrates on biological factors such as nonlinear growth and grazing interplay between phytoplankton and zooplankton [32]. The other emphasizes the influences of physical factors such as currents, eddies, and the turbulent stirring induced by them [33].
Medvinsky et al. [30] followed the first approach. They developed a model of phytoplankton and zooplankton with equal diffusivities and neglected turbulent stirring. This succeeded in reproducing patchiness-like patterns developed from regular spirals. In their simulation, the initial spiral pattern formation is followed by its collapse, which starts near the center, generating spatiotemporal chaos over the region. They also showed that the processes are almost independent of the initial distributions of the constituents. They found that biological factors play a very important role in the production of plankton patchiness [30].
Most researches following the second contention [31–33] adopt the three-component model comprising nutrients or their substitutes, phytoplankton, and zooplankton [34]. Abraham [31] developed a model consisting of the carrying capacities and population densities of phytoplankton and zooplankton. The maturation time of zooplankton enabled fine structures to emerge in the spatial distributions of the zooplankton population whereas the distribution of the phytoplankton population did not show such fine structures because of their rapid growth. He used the seeded-eddy model [35] as turbulent flows.
Our present research builds on the second contention. We propose a model with a series of reaction-advection-diffusion equations in regard to nutrients and phyto-plankton, which allow spatial phenomena, such as sinking and turbulence, to be described directly, therefore enabling research on spatial structures. The sinking, turbulence, and mixing of phytoplankton have an important influence on the dynamic behavior of the system. Wang et al. [36] developed a relatively simple reaction-diffusion-advection model of the interaction between algae and mussels, which included the diffusive spread of mussels and the advection of algae. Theoretical results also showed the significance of sinking, turbulence, and diffusion [37, 38].
In Section 2 we present the mathematical formulation of the model incorporating the biology in the partial differential equations. In Section 3, we analyze the positive equilibrium state, the necessary conditions of Turing instability are deduced, and we analyze the stable behavior of the nonspatial and spatial systems. We derive the condition where the stationary homogeneous solution is linearly unstable. In Section 4, a number of simulations are discussed. From this we show the key factors and the critical conditions in the system. In Section 5, we discuss the outcomes and draw our final conclusion.
2. Mathematical Formulation of the Model
Fresh water wetland ecosystems depend on the available nutrients. Most often, man can contribute to the nutrient through sewage and agricultural fertilizer runoffs. Phytoplankton growth, in particular, is dependent on light and nutrients [1]. Our model is described aswhere the variables and represent the concentration of nutrients and phytoplankton, respectively, which are functions of the actual time, ; all the parameters are positive constants, with regard to the parameters; is the (constant) input of nutrients from the environment; is the rate of the nutrient concentration loss; is the consumption rate of nutrient by phytoplankton; is the maximum nutrient intake rate for phytoplankton; is the phytoplankton mortality; is the maximum predation rate of zooplankton on phytoplankton; is the half-saturation constant for phytoplankton; is the sinking velocity of phytoplankton, which is caused by turbulent stirring; is the width in the water column; is the depth of water, and are the diffusion rates of nutrients and phytoplankton, respectively, which arise from turbulent flow; and is the usual Laplacian operator in two-dimensional space.
The consumption of nutrients by phytoplankton is represented by the Holling I Functional response function. The implicitly assumed zooplankton grazing on phytoplankton is combined as a nondynamic term, consisting of a Holling II Functional response function.
The term addressing removal of nutrients, , is added to estimate the losses except for the intake by phytoplankton, such as the sedimentation to a deep bottom or a runoff from the system. This term has another meaning in the case of mathematical modelling; if the phytoplankton is finally extinct, the nutrient concentration continues to increase infinitely. To avoid this situation, we add the removal term. However, if the term was omitted, the simulation results for spatial pattern formation are unchanged.
The term of represents the feedback of the bottom sediments [1]. When the bottom phytoplankton is decreased, the bottom sediments will be more easily affected by winds, waves, and bottom zooplankton. Finally, nutrients will be more easily released into the water by turbulent stirring and turbulent flow. It is obvious that the term is nondecreasing in . Furthermore, when increases, it will reach an upper bound due to when all green plants in the bottom disappearing, increasing will not release any more nutrients from the bottom sediments. Generally, ecologists make use of the sigmoid function, and hence is the feedback.
As a response to the interaction between nutrients and phytoplankton, we utilize a Holling type I Functional response function in our model. Holling explained the functional response of predators to prey density and the impact on mimicry and population regulation [39].
It is also important to use a Holling type II Functional response to describe zooplankton grazing on phytoplankton, because this is a nondynamic term. This type of functional response has been broadly applied to represent zooplankton predation in many theoretical studies [30, 39–41] and shows good agreement with experimental data [42, 43].
3. Stability Analysis of the System
3.1. Analysis of the Spatially Uniform System
The system of (1a) and (1b), without spatial derivatives, can be described as
The vertical isocline, where , is , which is continuous when . When is close to zero, is close to the positive infinite. On the other side, when is close to positive infinite, is close to . Because is continuous when , there must exist a point at least that satisfies . If , we can obtain , then it is monotonically decreasing when the . Since when , it follows that when .
The horizontal isocline, where , is , and . When , the line, , is the asymptote of . The horizontal isocline is continuous when . Also, is satisfied when ; is satisfied when .
When the system consists of bare nutrients without phytoplankton in the nonspatial model, , and . From the discussion above, if , there is always a boundary equilibrium state, which can be denoted as .
On the other hand, where and , if , there exists a positive equilibrium state in the nonspatial system when ; whereas there is no positive equilibrium in the nonspatial system when . The positive equilibrium can be written as .
Thus, we confirm that positive equilibrium states exist in the nonspatial system. Of course, it does not mean that there cannot be positive equilibrium states in the nonspatial system when these conditions are not satisfied.
The Jacobian matrix of the nonspatial system at the equilibrium is
The index of equilibrium is when the condition, , , is satisfied, which is stable. According to the eigenvalues of , it is locally asymptotically stable when the above conditions are satisfied. In particular, the equilibrium is unstable, when or .
From the above discussion, we find a positive equilibrium in the nonspatial system under some preconditions, defined by . The Jacobian matrix of the nonspatial system at the equilibrium, is
When , the index of equilibrium is if . It is locally asymptotically stable using the Routh-Hurwitz criteria when the above conditions are satisfied. In particular, the equilibrium is unstable, when , which is a saddle. Figure 1 shows the phase diagram of system which is locally asymptotically stable.
The specific expression of equilibrium can be difficult to obtain. The stability of the nonspatial system is defined only when some parameters are given.
3.2. The Analysis of the Spatial System without Advection
The spatial system of (1a) and (1b) without the advection term isLet us consider the significance of the diffusion rate of nutrients and phytoplankton, represented by and , respectively, we employ a linear analysis technique to focus on the parameters essential for the system behavior. We want to derive the necessary conditions, where and affect the system.
We have considered the stationary homogeneous to be linearly stable in the absence of the diffusion terms. Symmetry breaking occurred when the stationary homogeneous solution was unstable to small spatial perturbations. To analyze the spatial system of (5a) and (5b) and confirm how a small heterogeneous perturbation of the homogeneous steady state changes over a large time period, consider the perturbation [35]where is the wave-number, is the perturbation growth rate, and is the usual imaginary unit (). Substituting (6) into (5a) and (5b) and neglecting all nonlinear terms in and , the characteristic equation for the eigenvalues, iswhere the elements of the Jacobian determinant of the nonspatial system are taken at the positive equilibrium solution , as below:and it is clear that , , and when is at positive equilibrium. The characteristic equation for the system of (5a) and (5b) can be simplified to where and .
Symmetry breaking occurred when the equilibrium is unstable to small spatial perturbations. However it is linearly stable to perturbations in the absence the diffusion terms. Turing bifurcation corresponds to a nonequilibrium phase transition, which transforms the system from a homogeneous steady state to a heterogeneous steady state and produces spatial periodic oscillations. The system is stable to spatial homogeneous perturbations, which requires , and so the system must satisfy
Furthermore, one of the eigenvalues is positive, so the only mechanism for the system to lose stabilization is to produce a saddle node bifurcation. The system is unstable, which corresponds to the small perturbations of modulus. Thus it must satisfy for some values of
And, . Hence we can derive the necessary conditions of Turing instability for the system of (5a) and (5b).
When , Turing instability occurs, if the conditions (11), (12), (13) and (14) are satisfied. According to (11), (12), we can obtain , , but , the condition (13) is not satisfied. Thus, the system does not satisfied the conditions of Turing instability when .
When , and the other parameter values are known, we can derive the stability of the boundary equilibrium and positive equilibrium , which are dependent on , as shown in Figure 2. The curves and straight line represent the stability of the equilibrium in the nonspatial system. For boundary equilibrium , the yellow straight line shows that it is unstable for arbitrary . For positive equilibrium , the black curve shows the region where it is unstable, and the red curve shows where it is stable. For positive equilibrium, the vertical green () and blue () lines divide the plane into three zones of different system behavior: (I) the system is not stable for the nonspatial system; (II) the system is stable for the nonspatial system, and satisfies the conditions for Turing instability; (III) the system is stable for the nonspatial system, but does not satisfy the conditions for Turing instability.
3.3. Analysis of the Spatial System with Advection
Let us now consider the spatial system including the advection and diffusion terms and analyze the interaction with each other. Similar to the previous section, symmetry breaking occurs when the stationary homogeneous solution is unstable to small spatial perturbations with the mutual effect of diffusion and advection, and it is stable to perturbations without those terms. The perturbation considered is that of (6).
Substituting (6) into (1a) and (1b), and neglecting all nonlinear terms in and , the characteristic equation for the eigenvalues, , iswhere the elements of the determinant are the same as in the previous discussion. This may be rewritten aswhere , and .
From (15), it follows thatwhere , and .
The real and imaginary parts of the eigenvalues arewhere . If the solution is stable, the real parts of all eigenvalues are less than zero, that is, ; it is unstable when one or more of the real parts is greater than zero. Therefore, the critical point is , but the analytical expression for the critical point is not easy to be acquire. However, when , (17) can be simplified as and the equilibrium is unstable where the right-hand side is less than zero. Using the parameters given in Figure 2, in zone II, system satisfies the conditions for Turing instability with , , and (i.e., ). Hence, an essential condition for (17) is . If , for some value of wave-number , the right-hand side of (19) is positive for ; if , for some value of wave-number , then the right-hand side of (19) is positive, for . In zone III, the system does not satisfy the conditions for Turing instability. However in this zone, , if , so the right-hand side of (17) is positive, for .
In Figure 2, we described the effect of the nutrient concentration, , on the system of (5a) and (5b) (i.e., without the advection term). In Figure 3, we show the effect of sinking velocity () for the case of fixed nutrient concentration in zone III, . The system is stable, but does not satisfy the conditions for Turing instability. The critical values are when , , and the maximal real part of is equal to zero. The equilibrium is unstable when the sinking velocity exceeds the critical value.
4. Simulation
In general, systems including reaction-diffusion or reaction-advection-diffusion terms (Equations (5a), (5b) and (1a), (1b), resp.), describes the time evolution of the spatial distribution of species density. The nonuniform stationary states of the model, the analytical expression for the spatial models, are difficult to obtain. Analytical methods are not sufficient to understand the system completely, so it is essential to use simulations. If the parameters are chosen appropriately, the simulations will produce meaningful two-dimensional patterns.
To solve the differential equations, we discretize the time and space of the problem, transforming from an infinite-dimensional to a finite-dimensional form. The continuous problem in two dimensions is solved on a rectangular spatial grid of points. The spacing between the spatial grid points is defined by the spatial grid constant . In the discrete system the Laplacian describing diffusion is calculated using finite differences. That is, the derivatives are approached by differences over , and as the differences approximate the derivatives. The time evolution is also discretized, in steps of . The time evolution can be solved using the Euler method, approaching the value of the concentration at the next time step using the exchange rate of the concentration from the former time step.
In our numerical simulations we utilized zero-flux boundary conditions with a system size of space units and parameters , , , , , , , and . The spatiotemporal dynamics of a reaction-diffusion system or reaction-advection-diffusion system depends on the initial conditions, which some authors have taken into account in relation to the problem of biological invasion [30, 44, 45]. The initial conditions are naturally described by some specific functions, describing a randomly perturbed homogeneous state. Systems for (5a), (5b) and (1a), (1b) are solved numerically with a centered finite difference approach for the spatial derivatives with a spatial step size of and the Euler method for time integration with a time step size of . In the simulation, we assume that the diffusion rate of nutrients is greater than that of the phytoplankton.
Figure 4 shows the evolution of the spatial pattern at 0, 125,000, 250,000, 375,000, and 500,000 iterations, with a small random perturbation of the stationary homogeneous solution . The parameters satisfied the conditions for Turing instability in the system of (5a) and (5b), as discussed in Section 3.2. Figure 4 shows the spontaneous emergence of spotted spatial patterns. After 375,000 interation the regular spotted pattern prevails over the whole domain and the progress of the system undergoes no further change over time.
(a)
(b)
(c)
(d)
(e)
Figure 5 shows the evolution of spatial pattern with a small random perturbation of the stationary solution of the spatially homogeneous system. Band patterns develop in the domain, as the iterations increase. In Section 3.3, we found that if the conditions are satisfied, the system will become unstable to small spatial perturbations in the presence of diffusion and advection. Because the sinking velocity is linked to nutrient concentration, the parameters given in Figure 5, are satisfy linearly unstable and oscillation will occure, forming the clear band pattern of nutrients shown.
(a)
(b)
(c)
(d)
(e)
We know from Sections 3.2 and 3.3, that the system of (5a) and (5b) with the parameters as shown in Figure 6 does not satisfy the conditions for Turing instability, whereas the system of (1a) and (1b) would be unstable to spatial perturbations of nutrient concentration and horizontal velocity of phytoplankton. Figure 6 shows that including diffusion and advection terms, (1a) and (1b) produces periodic oscillation, with the crenulate patterns becoming more regular over time.
(a)
(b)
(c)
(d)
(e)
Figures 4 and 5 show the relationship between diffusion-driven instability and diffusion-advection-driven instability, respectively. Figure 4 shows a regular spotted pattern, whereas Figure 5 shows a regular band pattern. These are all unstable to small spatial perturbations, in the case of Figure 4 this causes Turing instability and produces the regular Turing pattern, whereas in the case of Figure 5 they cause linear instability and produce regular parallel patterns from the effect of horizontal velocity.
Figures 5 and 6 show the connection between different nutrient constant and horizontal velocity. Figure 5 shows the band pattern of reaction-diffusion-advection equations with , and . This system satisfies the conditions for Turing instability without the advection term (i.e., ). Figure 6, shows crenulated patterns of reaction-advection-diffusion with , and . However, the system does not satisfy the conditions for Turing instability without the advection term.
Thus, the horizontal flux of phytoplankton has a significant role in the system.
5. Discussion and Conclusion
In the real world, there are many examples of biological patterning, such as, the snow leopard, in which the body and tail are all spots; our national treasure: the panda, whose body is covered with striped skin-hair; banded vegetation patterns in some ecological systems; spiral waves; and so on. These phenomena are caused by internal and physical factors.
We developed a reaction-advection-diffusion model for nutrients and phytoplankton, and studied the relationship between nutrients and phytoplankton. Although our model was only an abstraction of the reality, it was shown to be useful to explain some phenomena. Theoretical studies have researched the conditions of the instability of the systems of (5a), (5b), and (1a), (1b). Our research concentrated on nutrient source and phytoplankton, and probed how the nutrient concentration and diffusion of nutrients and phytoplankton affect the system of (5a) and (5b). We also investigated the effects of horizontal movement of phytoplankton and nutrient inputs on the system of (1a) and (1b).
We have analyzed the conditions for Turing instability, which can be induced by the diffusion of nutrients and phytoplankton in the system of (5a) and (5b). For the system of (1a) and (1b), we showed that the horizontal movement of phytoplankton causes the homogeneous steady state to become unstable. The critical value of the velocity of horizontal phytoplankton relied on the input of nutrients.
Simulations confirmed the occurrences of these instabilities, caused by small spatial perturbations.
In our study, we researched the effect of predation of zooplankton on phytoplankton. Turing instability would not occur when lacking of the predation of zooplankton on phytoplankton. Thus, the predation of zooplankton on phytoplankton could affect the occurrence of Turing instability.
We also investigated the effect of nutrients and small spatial heterogeneous perturbations. The abundant nutrients could lead to the massive growth of phytoplankton. Abundant nutrients and the diffusion of nutrients and phytoplankton could induce the occurrence of Turing instability. The interaction among the diffusion of nutrients and phytoplankton, and the horizontal movement of phytoplankton could also cause the emergence of instability and oscillation.
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), the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510), and the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001).