Abstract

In this research, a continuous nutrient-phytoplankton model with time delay and Michaelis–Menten functional response is discretized to a spatiotemporal discrete model. Around the homogeneous steady state of the discrete model, Neimark–Sacker bifurcation and Turing bifurcation analysis are investigated. Based on the bifurcation analysis, numerical simulations are carried out on the formation of spatiotemporal patterns. Simulation results show that the diffusion of phytoplankton and nutrients can induce the formation of Turing-like patterns, while time delay can also induce the formation of cloud-like pattern by Neimark–Sacker bifurcation. Compared with the results generated by the continuous model, more types of patterns are obtained and are compared with real observed patterns.

1. Introduction

Phytoplankton is the main producer of aquatic ecosystems and forms the basis of food webs [1]. The growth of phytoplankton can determine the behavioral space of plankton directly and then even affects the development of the whole aquatic ecosystem. As the spatial and temporal distribution/self-organization of phytoplankton is quite visible, it can reflect the “structure” of the aquatic ecosystem, which also reveals the strategy that phytoplankton takes to adapt to environmental changes. Therefore, investigation on the self-organization of phytoplankton has been essential to the understanding of ecosystem evolution and remains a hot topic in ecosystem dynamics and complexity.

One of the most striking manifestations of biological ocean dynamics is the proliferation of particularly dense phytoplankton, which often forms marvelous patterns with unexpectedly clear boundaries, such as banded patterns, chaotic patterns, and solitons. [24]. These patterns are usually produced by the interaction between the behavior of organisms (sinking, floating, or swimming) and the surrounding physical flow [3]. By examining the structure of dense phytoplankton mass reproduction, the local physical dynamics of living organisms, as well as behavioral and physiological responses, can be obtained. Patterns of some accumulations are the reflection of specific flows; observation of these visible patterns may reveal underlying dynamics [3]. However, the formation mechanism of these patterns is still unclear, and further research is urgently needed.

Since the environmental factors in the phytoplankton ecosystem vary with time, simple statistical models cannot accurately describe the dynamic process of phytoplankton growth in time and space. Therefore, the establishment of a dynamic model of plankton populations is an important method to study the physical and biological processes of plankton [510]. In 1949, Riley et al. established a vertical one-dimensional ecological dynamics model to simulate the seasonal changes of plankton in the North Sea of Europe, which symbolizes the study of aquatic ecosystems from qualitative description to the era of quantitative simulation [11]. In the same year, Monod, the founder of modern cell growth, proposed the famous Monod kinetic model, which describes the effect of substrate concentration on the cell growth rate, laying the foundation for the development of cell growth dynamics, thus promoting the development of phytoplankton growth dynamics [12]. After the Monod model was proposed, it was widely used in the study of phytoplankton growth dynamics [1319] and found that the Monod model is very suitable under certain conditions. Ye et al. used the Monod model to quantitatively study the relationship between phytoplankton biomass and nutrients in Taihu Lake [20]. It was found that total nitrogen is a key factor affecting the growth of phytoplankton, and eutrophication of lakes to phytoplankton growth has served as positive feedback. However, some studies have shown that the relationship between phytoplankton growth and nutrients does not follow the Monod model, which means that the Monod model has limitations [21]. Dugdale used the Michaelis–Menten model to describe the absorption of nutrients by planktonic algae under steady-state conditions [22, 23], revealing the effects of limited nutrient salts on the growth of planktonic algae. Studies have shown that it is feasible to describe the quantitative relationship between phytoplankton growth rate and nutrient salt concentration with Michaelis–Menten model under steady-state conditions [24].

Reaction-diffusion equations have been widely used to describe temporal and spatial complexity in ecology [2533]. Time delay always exists and plays an important role in dynamical systems, and some dynamical phenomena could be affected by the time delay [3437]. Turing bifurcation helps with the explanation of self-organization. Usually, Turing analysis is carried out on the basis of continuous models. In ecosystems, discrete models can generate more dynamical behaviors, which undoubtedly reflect ecological complexity [38, 39]. It is worth noting that numerical simulation is actually based on discrete forms of continuous systems, which automatically give numerical simulation algorithms. Therefore, the discrete form is the natural link between the actual model and the simulation. At the same time, in many cases, the environment exhibits considerable spatial variability and is best considered using spatial discrete methods [4043]. Wang et al. studying the discrete forms of the Gierer–Meinhardt system found that discrete systems lead to more complex nonlinear behaviors and new dynamic phenomena than continuous forms. In addition, for the discrete time and space G-M system, they found a new instability mechanism, pronunciation flip-Turing instability [44]. In addition, when we apply Turing analysis on discrete models, there are some differences. First, the solution space of the continuous Turing model is infinite dimension, while the corresponding discretization model is in finite dimension. Secondly, in the linearization analysis of Turing bifurcation, elliptic operators have infinite eigenvalues and eigenvalue sequences are unbounded, while discrete elliptic operators have only finite eigenvalues. Therefore, it is of great significance to further reveal the formation mechanism of space-time self-organization by discretizing the reaction-diffusion equation in space and time. In this research, we use the N-S bifurcation theorem described by Guckenheimer and Holmes [45] to determine the conditions under which bifurcation occurs.

Therefore, in this research, we will investigate the self-organization of phytoplankton in a spatiotemporal discrete nutrient-phytoplankton model with time delay. The model is obtained through the discretization of a continuous nutrient-phytoplankton model proposed by Dai et al. [46]. Based on the discrete model, Neimark–Sacker bifurcation analysis and Turing bifurcation analysis are carried out, and parameter space is obtained correspondingly. Numerical simulations are then carried out using MATLAB (2011a) on the formation of phytoplankton patterns with different parameter values. Considering the pattern types and holistic shapes, the simulated phytoplankton patterns are compared with those observed in the real world. In addition, we also illustrate the results that can be obtained by the discrete model in the research are far more than those by the continuous model in [46].

2. Model and Analysis

2.1. A Discrete Spatiotemporal Nutrient-Phytoplankton Model with Delay

Considering the rate of diffusion of nutrients and phytoplankton is different in real-world systems, therefore, in this paper, we obtain the following spatial nutrient-phytoplankton models based on Ref. [46]:where and represent nutrient concentration and the population biomass of phytoplankton, respectively. All parameters are non-negative and . Let and .

In the Turing instability analysis, people often give theoretical analysis to the continuous model, while the Turing pattern simulation undoubtedly adopts the corresponding discretization model. It is important to discretize the reaction-diffusion equation in space-time, which is of great significance for further revealing the formation mechanism of space-time self-organized structure. Therefore, in this paper, we consider the model on a lattice, and the two variables can be expressed as and ( and ), which represent the nutrient concentration and the phytoplankton density in lattice at time , respectively. According to the former research works of [39, 47], there are two stages, reaction stage and diffusion stage, when we discretize the continuous model (equation (1)). Meanwhile, we use a standard second-order central finite difference scheme to discretize (1) in space. The spatial dispersal stage, that is, diffusion, is considered firstly aswhere , , , and are the time interval and space interval. denotes the discrete form of the Laplacian operator:

Then, we consider the reaction stage:

Boundary conditions are set as periodic conditions as follows:

Equations (2) and (4) including both diffusion and reaction stages are defined as our discrete model. According to the spatial homogeneity of the discrete system, with all , and , we have

From (6), together with (2)–(4), then the homogeneous dynamics satisfies the following equation:where , we rewrite equation (7) into the form of maps:

In the sequel, we use (8) to analyze the homogeneous dynamics of equations (2)–(4).

2.2. Analysis of the Stability of the Homogeneous Steady State

The fixed point of map (8) is a homogeneous steady state of equations (2)–(4). So, we carry out the analysis of the stability about the fixed point of map (8). When both the delay and the diffusion are not considered, it is easy to see that map (8) has two fixed points: and , where and . The sufficient and necessary conditions for that is positive are as follows:

The Jacobian matrix method is used to determine the local stability of fixed points. Since pattern formation requires the nontrivial homogeneous stationary state, which is corresponding to the fixed point , we just concern the stability of . The corresponding Jacobian is given by

The two eigenvalues of matrix (10) are as follows:where

According to the criterion of stability of the fixed point, if and , the fixed point is locally asymptotically stable. As shown in the previous research [48], this criterion of stability is equivalent to . Through calculations, the conditions for to be locally asymptotically stable are determined as follows:where

The dynamics of the discrete model without spatial diffusions can be shown in Figure 1, and the following bifurcations and simulations are all carried out in this situation. Note that, in the discrete system, the arrows in the figure represent the direction of iteration rather than dynamic flow. And is considered, as when time delay increases, the following bifurcation will occur.

2.3. Bifurcation Analysis of the Homogeneous Steady State
2.3.1. Neimark–Sacker Bifurcation Analysis

In this part, we consider Neimark–Sacker bifurcation of equations (2)–(4) when is small. We first make a corresponding change to equation (1) and obtain equation (17) and then discretize equation (17) to obtain map (18).

If is small, , put it into equation (1) and obtainin which

Taylor expansion of at without considering the higher order nonlinear part and equation (1) can be written aswhere and . Therefore, equation (8) can be expressed as the following map:

According to [48], the first condition of Neimark–Sacker bifurcation requires the eigenvalues at a fixed point, and are conjugate and their modules are 1, shown asin which have been shown in Appendix A.

Then, we can get

By calculation, we can obtain the Neimark–Sacker bifurcation point :

Under the satisfaction of conditions (22) and (23), the fixed point of map (18) is translated to the origin by the following translation:

Then, map (18) is transformed intowhere have been shown in Appendix A with .

The second condition for Neimark–Sacker bifurcation requiresin whichwhere and are described in equation (21) with . Then, we can obtain

On the basis of map (25), the canonical form is studied to obtain the last condition for the occurrence of Neimark–Sacker bifurcation. The invertible transformation is applied:

To map (25), then the map becomeswhere

The second-order and third-order partial derivatives of and at are calculated as

The third condition of Neimark–Sacker bifurcation requiresin which

After calculation, the third condition for Neimark–Sacker can be expressed asin which

Based on the above calculations, when the conditions (23), (28), and (35) are satisfied, Neimark–Sacker bifurcation occurs at the fixed point . Moreover, when and , then an attracting invariant closed curve bifurcates from the fixed point for ; otherwise, when and , a repelling invariant closed curve bifurcates from the fixed point for .

In Figure 2(a), the parameter space of time delay and input rate is shown. Zone I represents that when the input rate of nutrient is low, positive fixed point does not exist. Zone II and III represent no bifurcation and bifurcation area around the fixed point . And Figure 2(b) shows the bifurcation diagram when and . The Hopf bifurcation point is not suitable here with a time step not small enough, while the calculated Neimark–Sacker bifurcation point is more correct compared with the simulated bifurcation point. This difference can be shown in Figure 2(d): cannot be equal to 0.2, but in the discrete model, when , Neimark–Sacker bifurcation occurs and drives the system from stable fixed point to quasiperiodic dynamics. When time delay satisfy both and , the dynamics of the continuous model (with as in [46]) and the discrete model () show differences: Hopf bifurcation induces the continuous system into a limit cycle (Figure 2(c)); Neimark–Sacker bifurcation induces the discrete system into quasiperiodic iteration and the amplitude of variables is larger.

2.3.2. Turing Bifurcation Analysis

In this part, we consider Turing bifurcation of equations (2)–(4) when . Turing bifurcation requires two conditions. Firstly, a nontrivial homogeneous stationary state exists and is stable to spatially homogeneous perturbations, which has been obtained in the above section. Secondly, the stable stationary state is unstable to at least one type of spatially heterogeneous perturbations. This paper is based on the method of Bai and Zhang to do Turing bifurcation analysis of the discrete model [49]. We first consider the following eigenvalue equation:with periodic boundary conditions:

When and , we set , put it into equation (38), and obtain . Then, we can get

Substituting into equation (37), we obtain that

By calculating , then, using equation (39) we can get

Turing bifurcation is caused by the lack of spatial symmetry. To analyze the Turing bifurcation, a spatially nonuniform perturbation is applied at the spatially uniform state . The equation for spatially nonuniform perturbations can be expressed as

Noticing and , and the values of the two cannot be constant to 0. Substituting equation (42) into equations (2)–(4) leads to the following equation:

The high-order terms in the above equations can be ignored when the disturbance is very weak. Using the corresponding characteristic function of the eigenvalue to multiply equation (43), we obtain

Summing equation (44) for all and gives the following equation:

Let and , combined boundary conditions (38), equation (45) can be transformed into the following form:

At this time, we know that a progressive solution of equations (2)–(4) with boundary conditions (5) is . Therefore, we know that if the solution of the equation (46) is unstable, it will cause the equations (2)–(4) to be unstable. The eigenvalues of equations (2)–(4) are as follows:in which

Here, , , , and are denoted for reminding that they are dependent on . Based on the two eigenvalues, we define

The occurrence of Turing bifurcation needs to satisfy . According to this condition, we can determine the critical value of Turing bifurcation occurring, which can be described as follows:(1)If is satisfied in a small neighborhood of , the critical value satisfies (2)If is satisfied in a small neighborhood of , the critical value satisfies

and determine the critical conditions for the occurrence of Turing bifurcation in discrete systems. From this critical condition, the conditions for the formation of the Turing pattern can be obtained. By simplifying the two cases described above, it can be known that Turing instability occurs and leads to the formation of spatial heterogeneity patterns when .

Then, we discuss the influence of the diffusion coefficients and on the model. We use the following parameter values: , , , , , , , , , , , , and . The discrete model is solved numerically in a rectangular spatial grid consisting of 100 × 100 units, and periodic boundary conditions are adopted. As shown in Figure 3, Turing bifurcation can be displayed via the change of eigenvalue . In Figure 3(a), we can see that the perturbation parameters and have symmetrical effects on the eigenvalues. Therefore, for convenience of observation, let , then we can get the change in the eigenvalue to , as shown in Figure 3(b). If there is no disturbance, the system will be at a stable point. But, when there is disturbance, Turing bifurcation can happen. According to Figure 3(b), we can know that if the diffusion coefficient , then the eigenvalue of the system will not make the eigenvalue greater than 1 with the increase of the perturbation parameter; at this time, Turing bifurcation will not occur. When , the curve is tangent to the straight line , that is, only if the diffusion coefficient , the value of the eigenvalue will exceed 1 as the disturbance coefficient increases; at this time, Turing bifurcation will occur in the system.

3. Numerical Simulations

In this section, numerical simulations on the formation of phytoplankton are carried out focusing on three aspects including time delay and diffusion coefficients and . Note that initial conditions are set as homogeneous steady states on 100 × 100 (mostly) lattice with ±5% random perturbations. Then simulated patterns are compared with real observed patterns.

3.1. Effects of Time Delay on the Pattern Formation

According to the stability analysis and Neimark–Sacker bifurcation analysis, time delay can cause the instability of the system. As shown in Figure 4, with given parameter values, the critical value of can be calculated as . Thus, we take , , and and obtain Figures 4(a)4(c) respectively. Figures 4(a) and 4(b) show that when , there is no spatial patterns. Figure 4(c) clearly shows the formation of a Neimark–Sacker type pattern when . From the simulated patterns, we can see that simulations are consistent with the Neimark–Sacker analysis.

Figure 5 shows a time series of phytoplankton pattern formation process from initial conditions to self-organization when . Small patches can be formed at the beginning (Figures 5(b) and 5(c)), and then gradually forming spirals (Figure 5(d)). The spirals will break with time, and reform afterwards (Figures 5(e)5(g)). Then cloud-like patterns will finally be formed (Figures 5(h) and 5(i)).

3.2. Effects of Diffusion Coefficients on the Pattern Formation

Based on the Turing bifurcation analysis in Section 2.3.2, we can see that the diffusion coefficient can induce the occurrence of Turing bifurcations, and the critical value is . Figure 6(a) shows that there is no pattern, while in Figure 6(b) we can clearly see the appearance of a typical banded Turing-type pattern. In addition, diffusion coefficient can also induce the occurrence of Turing bifurcation as shown in Figure 6(c). The simulations are quite consistent with the above analysis of Turing bifurcation.

In order to show the formation process of Turing-type patterns, two sets of values are selected: , and , . Figure 7 shows that, at the beginning, random distribution of phytoplankton starts to form a labyrinth pattern (Figures 7(b)7(d)), and then some patches die out and form some thin stripes (Figures 7(f) and 7(g)). Finally, the number of stripes decreases over time (Figures 7(h) and 7(i)). Figure 8 shows a similar process of pattern formation with slight differences.

Note that when time delay satisfies the Neimark-Sacker bifurcation condition and coefficients and satisfy Turing bifurcation condition, simulations on pattern formation will be different. As shown in Figure 9, the general process of pattern formation is similar to that when time delay (see Figure 7). Note that as shown in Figure 9(d) with , the same formation stage of patterns can be obtained in Figure 7(e) with . This reveals that the time delay will not change the pattern type in this situation, but will cause the delay of pattern formation stages, which is reasonable and consistent with the research of [46].

3.3. Comparison between Simulated Patterns and Real Observed Patterns

From numerical simulations, mainly two types of patterns are obtained: cloud-like patterns (Figure 10(a)) and banded patterns (Figures 11(a), 11(c), and 11(e)). Some observed patterns of phytoplankton are found and can be briefly compared with the simulated patterns from the holistic shapes. Figure 10(b) is the blooming phenomenon in the Barents Sea and shows similarity to the simulated cloud pattern (Figure 10(a)).

Figures 11(b) and 11(f) are phytoplankton bloom images taken in the waters near Shenzhen, China, in 2014, and Figure 11(d) is a large area phytoplankton bloom image of Noctiluca scintillans taken in the waters near Rizhao in Shandong, China in 2012. The simulated banded patterns are similar to the observed patterns. But note that we can only compare the holistic shapes or types of the patterns due to the unknown of real pattern scales, and further research should be done to explore more detailed comparison on spatial scales of the patterns.

4. Discussion and Conclusions

Through theoretical analysis and numerical simulation, the principle of space-time dynamics of phytoplankton growth and the formation process of spatial pattern distribution of phytoplankton population are explored. We mainly study the effects of time delay and, diffusion coefficients on the spatiotemporal dynamics of phytoplankton growth. The conclusions of theoretical studies and the results of numerical simulations can show that: (1) time delay does not affect the stability of the stable fixed point , but the time delay may affect the whole process; (2) when time delay exists and is greater than a certain critical value , the time delay can not only lead to the instability of the stable fixed point , but also form a cloud Neimark–Sacker pattern through Neimark–Sacker instability; (3) when , by changing the value of or to make all the eigenvalues of the diffusion term are greater than 1, then Turing bifurcation occurs. A band-shaped Turing pattern can be formed.

Both the continuous model (original model in [46]) and the discrete model in this research are reaction-diffusion models, and they have the same functional responses. The differences between the two models are as follows: in the continuous model, reaction process and diffusion process occurs at the same time, while in the discrete model, the occurring order of reaction process and diffusion process can be varied, to be exactly in this research, diffusion process occur first and then reaction process. The main reason for us to calculate the diffusion process first is that spatial movements are quicker than reaction in the nutrient-phytoplankton system in the real world.

From the simulated patterns in this research, we can see that more complex and more types of patterns can be obtained with the discrete model than the continuous model. For example, Figure 10(a) is the simulated pattern with a similar type to the only simulated pattern in [46]. Besides, spotted, banded, and some mixed patterns are obtained in this research. This is reasonable, as Dai et al. [46] only did the Turing bifurcation. However, even with the same bifurcations, discrete models can still generate more complex patterns, due to the flexible of both time and space steps. Although continuous PDEs are still the main type of model that are used in pattern formation, many research studies on discrete dynamic systems have shown that spatiotemporal discrete models are not only more related to realistic processes but also can generate more complex dynamics [38, 39, 43, 44]. And more work could be done to explore the bifurcations induced by time step variation.

Appendix

A. Symbols Not Explained in Section 2

Data Availability

Codes related to the simulations can be made available on request (china907a@163.com).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to acknowledge with great gratitude for the supports of the Major Science and Technology Program for Water Pollution Control and Treatment (no. 2017ZX07101002) and the National Natural Science Foundation of China (Grant nos. 11875126 and 11771140).