# Spatiotemporal Patterns Formed by a Discrete Nutrient-Phytoplankton Model with Time Delay

Authors: Feifan Zhang; Wenjiao Zhou; Lei Yao; Xuanwen Wu; Huayong Zhang
Journal: Complexity
 (2020)

---

## 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.


---

## Body

## 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. [2–4]. 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 [5–10]. 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 [13–19] 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 [25–33]. Time delay always exists and plays an important role in dynamical systems, and some dynamical phenomena could be affected by the time delay [34–37]. 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 [40–43]. 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]:(1)∂N∂t=d1ΔN+I−bNa+NP−qN+εmP,∂P∂t=d2ΔP+αbNτa+NτP−mP,where N and P represent nutrient concentration and the population biomass of phytoplankton, respectively. All parameters are non-negative and Nτ=Nx,t−τ. Let f1N,P=I−bN/a+NP−qN+εmP and g1N,P=αbNτ/a+NτP−mP.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 an×n lattice, and the two variables can be expressed as Ni,jt and Pi,jt (i,j∈1,2,3,…,n and t∈Z+), which represent the nutrient concentration and the phytoplankton density in lattice i,j at time t, 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 as(2)Ni,j,t′=Ni,jt+D1∇d2Ni,jt,Pi,j,t′=Pi,jt+D2∇d2Pi,jt,where D1=d1e/h2, D2=d2e/h2, e, and h are the time interval and space interval. ∇d2 denotes the discrete form of the Laplacian operator:(3)∇d2Ni,jt=Ni+1,jt+Ni−1,jt+Ni,j+1t+Ni,j−1t−4Ni,jt,∇d2Pi,jt=Pi+1,jt+Pi−1,jt+Pi,j+1t+Pi,j−1t−4Pi,jt.Then, we consider the reaction stage:(4)Ni,jt+1=fNi,j,t′,Pi,j,t′=Ni,j,t′+ef1Ni,jt,Pi,jt,Pi,jt+1=gNi,j,t′,Pi,j,t′=Pi,j,t′+eg1Ni,jt,Pi,jt.Boundary conditions are set as periodic conditions as follows:(5)Ni,0t=Ni,kt,Ni,1t=Ni,k+1t,N0,jt=Nk,jt,N1,jt=Nk+1,jt,i,j∈1,2,…,k,t∈Z+.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 i,j, and t, we have(6)∇d2Ni,jt=0,∇d2Pi,jt=0.From (6), together with (2)–(4), then the homogeneous dynamics satisfies the following equation:(7)Ni,jt+1=Ni,jt+eI−bNi,jta+Ni,jtPi,jt−qNi,jt+εmPi,jt,Pi,jt+1=Pi,jt+eαbNτi,jta+Nτi,jtPi,jt−mPi,jt,where Nτi,jt=Ni,j,t−τ, we rewrite equation (7) into the form of maps:(8)NP⟶N+eI−bNa+NP−qN+εmPP+eαbNτa+NτP−mP.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: E0=I/q,0 and E∗=N∗,P∗, where P∗=αIαb−Im−qma/m1−εααb−m and N∗=ma/αb−m. The sufficient and necessary conditions for that N∗,P∗ is positive are as follows:(9)αb−m>0,Iαb−Im−qma>0,εα<1.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 pointN∗,P∗, we just concern the stability of N∗,P∗. The corresponding Jacobian is given by(10)Je=JN∗,P∗=1−eq+αb−mαbI−Im−qmaabmα1−εαeεm−mαeαb−mαbI−Im−qmaabm1−εα1.The two eigenvalues of matrix (10) are as follows:(11)λ1,2=tre±tre2−4Δe2,where(12)tre=2−q+αb−mαbI−Im−qmaabmα1−εαe,Δe=1−q+αb−mαbI−Im−qmaabmα1−εαe+αb−mαbI−Im−qmaαabe2.According to the criterion of stability of the fixed point, ifλ1<1 and λ2<1, the fixed point is locally asymptotically stable. As shown in the previous research [48], this criterion of stability is equivalent to Δe<1,−1+Δe<tre<1+Δe. Through calculations, the conditions for N∗,P∗ to be locally asymptotically stable are determined as follows:(13)Ae2−Be<0,A>0,Ae2−2Be+4>0,where(14)A=αb−mαbI−Im−qmaαab,B=q+αb−mαbI−Im−qmaabmα1−εα.The dynamics of the discrete model without spatial diffusions can be shown in Figure1, 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 τ=0 is considered, as when time delay increases, the following bifurcation will occur.Figure 1
Fixed points of (8). Parameter values: I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, q=0.01, ε=0.001, and τ=0. ∗Asymptotically stable fixed point E∗; “○” represents the boundary fixed point E0.


### 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, Nτ=Nx,t−τ∂Nx,t/∂t, put it into equation (1) and obtain(15)∂N∂t=d1ΔN+I−bNa+NP−qN+εmP,∂P∂t=d2ΔP+αbPRNx,t−τ∂Nx,t∂t−mP,in which(16)RN=NN+a.Taylor expansion ofRNx,t−τ∂Nx,t/∂t at N,P without considering the higher order nonlinear part and equation (1) can be written as(17)∂Nx,t∂t=d1ΔN+f1N,P,∂Px,t∂t=−αabPτN+a2d1ΔN+d2ΔP+g2N,P−αabPτN+a2f1N,P,where f1N,P=I−bN/a+NP−qN+εmP and g2N,P=αbN/a+NP−mP. Therefore, equation (8) can be expressed as the following map:(18)NP⟶N+eI−ebNN+aP−eqN+eεmPP+eαbNPa+N−emP−αabPτN+a2eI−bNN+aP−qN+εmP.According to [48], the first condition of Neimark–Sacker bifurcation requires the eigenvalues at a fixed point, λ1 and λ2 are conjugate and their modules are 1, shown as(19)λ2=λ1¯,λ1=λ2=1,in which(20)λ1,2=tr0±tr02−4Δ02,(21)tr0=a11+a22=A0+A1+A2τ,Δ0=a11a22−a12a21=A0A1+A2τ−B0B1+B2τ,A0,A1,A2,B0,B1,B2,a11,a12,a21,anda22 have been shown in Appendix A.Then, we can get(22)tr02−4Δ0<0,Δ0=1.By calculation, we can obtain the Neimark–Sacker bifurcation pointτ0:(23)τ0=1+B0B1−A0A1A0A2−B0B2.Under the satisfaction of conditions (22) and (23), the fixed point E∗ of map (18) is translated to the origin by the following translation:(24)x=N−N∗,y=P−P∗.Then, map (18) is transformed into(25)xy⟶a11x+a12y+a13x2+a14y2+a15xy+a16x3+a17x2y+a18xy2+Ox+y4a21x+a22y+a23x2+a24y2+a25xy+a26x3+a27x2y+a28xy2+Ox+y4,where a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25,a26,a27,anda28 have been shown in Appendix A with τ=τ0.The second condition for Neimark–Sacker bifurcation requires(26)d=dλτdττ=τ0=A0A2−B0B22A0A1−B0B1+A0A2−B0B2τ0≠0,λτ0θ≠1,θ=1,2,3,4,in which(27)λτ0,λ¯τ0=tr0τ02±i24Δ0τ0−tr02τ0=β±iγ,i=−1,where tr0τ0 and Δ0τ0 are described in equation (21) with τ=τ0. Then, we can obtain(28)A0+A1+A2τ0≠0,−1.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:(29)xy=a120β−a11−γx˜y˜.To map (25), then the map becomes(30)x˜y˜⟶β−γγβx˜y˜+1a12γFx˜,y˜Gx˜,y˜,where(31)Fx˜,y˜=a12γa12a13+a15β−a11x˜2−a12a15γ2x˜y˜+a122γa12a16+a17β−a11x˜3−a122a17γ2x˜2y˜+Ox˜+y˜4,Gx˜,y˜=a122a13β−a11−a12a23+a14β−a11−a12a24β−a112+a12a15β−a11−a122a25β−a11x˜2−a12a24γ2y˜2+a12γ2a24−a15β−a11+a12a25x˜y˜+a123a16β−a11−a124a26+a122a17β−a11−a123a27β−a11−a122a28β−a112x˜3−a122γa17β−a11−a12a27−2a28β−a11x˜2y˜−a122a28γ2x˜y˜2+Ox˜+y˜4.The second-order and third-order partial derivatives ofFx˜,y˜ and Gx˜,y˜ at 0,0 are calculated as(32)Fx˜x˜=2a122a13γ+2a12a15γβ−a11,Fx˜y˜=−a12a15γ2,Fy˜y˜=0,Fx˜x˜x˜=6a123a16γ+6a122a17γβ−a11,Fx˜x˜y˜=−a122a17γ2,Fx˜y˜y˜=0,Fy˜y˜y˜=0,Gx˜x˜=2a122a13β−a11−a12a23+2a14β−a11−a12a24β−a112+2β−a11a12a15β−a11−a122a25,Gx˜y˜=2γa12a24β−a11−γa12a15β−a11−a122a25,Gy˜y˜=−2a12a24γ2,Gx˜x˜x˜=6a123a16β−a11−a124a26+6a122a17β−a11−a123a27β−a11−6a122a28β−a112,Gx˜x˜y˜=−2γa122a17β−a11−a123a27+4γa122a28β−a11,Gx˜y˜y˜=−2a122a28γ2,Gy˜y˜y˜=0.The third condition of Neimark–Sacker bifurcation requires(33)R=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21≠0,in which(34)ξ11=14a12γFx˜x˜+Fy˜y˜+iGx˜x˜+Gy˜y˜,ξ20=18a12γFx˜x˜−Fy˜y˜+2Gx˜y˜+iGx˜x˜−Gy˜y˜−2Fx˜y˜,ξ02=18a12γFx˜x˜−Fy˜y˜−2Gx˜y˜+iGx˜x˜−Gy˜y˜+2Fx˜y˜,ξ21=116a12γFx˜x˜x˜+Fx˜y˜y˜+Gx˜x˜y˜+Gy˜y˜y˜+iGx˜x˜x˜+Gx˜y˜y˜−Fx˜x˜y˜−Fy˜y˜y˜.After calculation, the third condition for Neimark–Sacker can be expressed as(35)R=−1A∗B∗Fx˜x˜+Fy˜y˜Fx˜x˜−Fy˜y˜+2Gx˜y˜−Gx˜x˜+Gy˜y˜Gx˜x˜−Gy˜y˜−2Fx˜y˜+C∗Gx˜x˜+Gy˜y˜Fx˜x˜−Fy˜y˜+2Gx˜y˜−Fx˜x˜+Fy˜y˜Gx˜x˜−Gy˜y˜−2Fx˜y˜−132a122γ2Fx˜x˜+Fy˜y˜2+Gx˜x˜+Gy˜y˜2+116a12γβFx˜x˜x˜+Fx˜y˜y˜+Gx˜x˜y˜+Gy˜y˜y˜+γGx˜x˜x˜+Gx˜y˜y˜−Fx˜x˜y˜−Fy˜y˜y˜−164a122γ2Fx˜x˜−Fy˜y˜−2Gx˜y˜2+Gx˜x˜−Gy˜y˜+2Fx˜y˜2≠0,in which(36)A∗=32a122γ21−β+γ2,B∗=1−3β+2β2−2γ2β2−γ2+6β−8β2γ2,C∗=2βγ1−β1−2β−2γ2−γ3−4ββ2−γ2.Based on the above calculations, when the conditions (23), (28), and (35) are satisfied, Neimark–Sacker bifurcation occurs at the fixed point E∗. Moreover, when R<0 and d>0, then an attracting invariant closed curve bifurcates from the fixed point E∗ for τ>τ0; otherwise, when R>0 and d>0, a repelling invariant closed curve bifurcates from the fixed point E∗ for 0<τ<τ0.In Figure2(a), the parameter space of time delay and input rate I−τ is shown. Zone I represents that when the input rate of nutrient is low, positive fixed point E∗ does not exist. Zone II and III represent no bifurcation and bifurcation area around the fixed point E∗. And Figure 2(b) shows the bifurcation diagram when I=1.6 and e=0.2. The Hopf bifurcation point is not suitable here with a time step e 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): e cannot be equal to 0.2, but in the discrete model, when e=0.2, Neimark–Sacker bifurcation occurs and drives the system from stable fixed point to quasiperiodic dynamics. When time delay τ=3 satisfy both τ>τN−S and τ>τHopf, the dynamics of the continuous model (with e=0.01 as in [46]) and the discrete model (e=0.2) 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.Figure 2
(a) Parameter space graph of Neimark–Sacker bifurcation and Hopf bifurcation withI and τ; (b) bifurcation diagram of Neimark–Sacker bifurcation with τ and phytoplankton; (c–e) phase portraits of Neimark–Sacker bifurcation with parameters (c) e=0.01, τ=3; (d) e=0.2, τ=2.45; (e) e=0.2, τ=3. Other parameter values: I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, q=0.01, and ε=0.001. ∗Asymptotically stable fixed point E∗.
(a)
(b)
(c)
(d)
(e)


#### 2.3.2. Turing Bifurcation Analysis

In this part, we consider Turing bifurcation of equations (2)–(4) when τ=0. 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:(37)∇2Xij+λXij=0,with periodic boundary conditions:(38)Xi,0=Xi,k,Xi,1=Xi,k+1,X0,j=Xk,j,X1,j=Xk+1,j,i,j∈1,2,…,k.Whenx≠0 and y≠0, we set Xij=xiyj, put it into equation (38), and obtain xk=yk=1. Then, we can get(39)xr=exp2r−1πki′,yl=exp2l−1πki′,r,l∈1,2,…,k,i′=−1.SubstitutingXij=xiyj into equation (37), we obtain that(40)xi−1yj+xi+1yj+xiyj−1+xiyj+1−4xiyj+λxiyj=0.By calculatingλ=−1−x2/x−1−y2/y, then, using equation (39) we can get(41)λrl=−1−xr2xr−1−yl2yl=4sin2r−1πk+sin2l−1πk=Srl2.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 stateE∗. The equation for spatially nonuniform perturbations can be expressed as(42)N˜i,jt=Ni,jt−N∗,P˜i,jt=Pi,jt−P∗.Noticing∇d2N˜i,jt=∇d2Ni,jt and ∇d2P˜i,jt=∇d2Pi,jt, and the values of the two cannot be constant to 0. Substituting equation (42) into equations (2)–(4) leads to the following equation:(43)N˜i,jt+1=fNN˜i,jt+D1∇d2N˜i,jt+fPP˜i,jt+D2∇d2P˜i,jt+ON˜i,jt+P˜i,jt2,P˜i,jt+1=gNN˜i,jt+D1∇d2N˜i,jt+gPP˜i,jt+D2∇d2P˜i,jt+ON˜i,jt+P˜i,jt2.The high-order terms in the above equations can be ignored when the disturbance is very weak. Using the corresponding characteristic functionXrlij of the eigenvalue λrl to multiply equation (43), we obtain(44)XrlijN˜i,jt+1=fNXrlijN˜i,jt+fPXrlijP˜i,jt+D1fNXrlij∇d2N˜i,jt+D2fPXrlij∇d2P˜i,jt,XrlijP˜i,jt+1=gNXrlijN˜i,jt+gPXrlijP˜i,jt+D1gNXrlij∇d2N˜i,jt+D2gPXrlij∇d2P˜i,jt.Summing equation (44) for all i and j gives the following equation:(45)ΣXrlijN˜i,jt+1=fNΣXrlijN˜i,jt+fPΣXrlijP˜i,jt+D1fNΣXrlij∇d2N˜i,jt+D2fPΣXrlij∇d2P˜i,jt,ΣXrlijP˜i,jt+1=gNΣXrlijN˜i,jt+gPΣXrlijP˜i,jt+D1gNΣXrlij∇d2N˜i,jt+D2gPΣXrlij∇d2P˜i,jt.LetN¯t=∑i,j=1kXrlijN˜ijt+1 and P¯t=∑i,j=1kXrlijPijt+1, combined boundary conditions (38), equation (45) can be transformed into the following form:(46)N¯t+1=fN1−D1λrlN¯t+fP1−D2λrlP¯t,P¯t+1=gN1−D1λrlN¯t+gP1−D2λrlP¯t.At this time, we know that a progressive solution of equations (2)–(4) with boundary conditions (5) is Nijt=N¯tXrlijandPijt=P¯tXrlij. 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:(47)λ±r,l,e=12trr,l,e±12trr,l,e2−4Δr,l,e,in which(48)trr,l,e=a11e+a22e−eh2a11ed1+a22ed2λrl,Δr,l,e=a11ea22e−a12ea21e1−eh2d1λrl1−eh2d2λrl,JN∗,P∗=a11ea12ea21ea22e=1−eq−eabN∗+a2P∗eεm−ebN∗N∗+aeαabN∗+a2P∗1−em+eαbN∗N∗+a.Here,a11e, a12e, a21e, and a22e are denoted for reminding that they are dependent on e. Based on the two eigenvalues, we define(49)λmr,l,e=maxλ+r,l,λ−r,l,λme=maxr=1,l=1nλmr,l,e.The occurrence of Turing bifurcation needs to satisfyλme=1. According to this condition, we can determine the critical value e′ of Turing bifurcation occurring, which can be described as follows:(1)
Iftrr,l,e2>4Δr,l,e is satisfied in a small neighborhood of e=e′, the critical value e′ satisfies maxr=1,l=1ntrr,l,e′−Δr,l,e′=1(2)
Iftrr,l,e2≤4Δr,l,e is satisfied in a small neighborhood of e=e′, the critical value e′ satisfies maxr=1,l=1nΔr,l,e′=1max
r
=
1
,
l
=
1
n
t
r
r
,
l
,
e
′
−
Δ
r
,
l
,
e
′
=
1 and maxr=1,l=1nΔr,l,e′=1 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 λme>1.Then, we discuss the influence of the diffusion coefficientsd1 and d2 on the model. We use the following parameter values: I=1.6, a=0.5, b=0.8, α=0.7, q=0.01, ε=0.001, m=0.45, d1=0.2, d2=0.2, e=0.2, h=1, τ=0, and t=100. 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 λk,l. In Figure 3(a), we can see that the perturbation parameters k and l have symmetrical effects on the eigenvalues. Therefore, for convenience of observation, let k=l, then we can get the change in the eigenvalue to l, 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 d1=1.220, then the eigenvalue of the system will not make the eigenvalue λk,l greater than 1 with the increase of the perturbation parameter; at this time, Turing bifurcation will not occur. When d1=1.389, the curve is tangent to the straight line λk,l=1, that is, only if the diffusion coefficient d1>1.389, the value of the eigenvalue λk,l will exceed 1 as the disturbance coefficient increases; at this time, Turing bifurcation will occur in the system.Figure 3
Variation of eigenvaluesλk,l with perturbation numbers k and l. Parameter values: I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d2=0.20, q=0.01, τ=0, h=1, ε=0.001, and n=100.
(a)
(b)


## 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]:(1)∂N∂t=d1ΔN+I−bNa+NP−qN+εmP,∂P∂t=d2ΔP+αbNτa+NτP−mP,where N and P represent nutrient concentration and the population biomass of phytoplankton, respectively. All parameters are non-negative and Nτ=Nx,t−τ. Let f1N,P=I−bN/a+NP−qN+εmP and g1N,P=αbNτ/a+NτP−mP.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 an×n lattice, and the two variables can be expressed as Ni,jt and Pi,jt (i,j∈1,2,3,…,n and t∈Z+), which represent the nutrient concentration and the phytoplankton density in lattice i,j at time t, 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 as(2)Ni,j,t′=Ni,jt+D1∇d2Ni,jt,Pi,j,t′=Pi,jt+D2∇d2Pi,jt,where D1=d1e/h2, D2=d2e/h2, e, and h are the time interval and space interval. ∇d2 denotes the discrete form of the Laplacian operator:(3)∇d2Ni,jt=Ni+1,jt+Ni−1,jt+Ni,j+1t+Ni,j−1t−4Ni,jt,∇d2Pi,jt=Pi+1,jt+Pi−1,jt+Pi,j+1t+Pi,j−1t−4Pi,jt.Then, we consider the reaction stage:(4)Ni,jt+1=fNi,j,t′,Pi,j,t′=Ni,j,t′+ef1Ni,jt,Pi,jt,Pi,jt+1=gNi,j,t′,Pi,j,t′=Pi,j,t′+eg1Ni,jt,Pi,jt.Boundary conditions are set as periodic conditions as follows:(5)Ni,0t=Ni,kt,Ni,1t=Ni,k+1t,N0,jt=Nk,jt,N1,jt=Nk+1,jt,i,j∈1,2,…,k,t∈Z+.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 i,j, and t, we have(6)∇d2Ni,jt=0,∇d2Pi,jt=0.From (6), together with (2)–(4), then the homogeneous dynamics satisfies the following equation:(7)Ni,jt+1=Ni,jt+eI−bNi,jta+Ni,jtPi,jt−qNi,jt+εmPi,jt,Pi,jt+1=Pi,jt+eαbNτi,jta+Nτi,jtPi,jt−mPi,jt,where Nτi,jt=Ni,j,t−τ, we rewrite equation (7) into the form of maps:(8)NP⟶N+eI−bNa+NP−qN+εmPP+eαbNτa+NτP−mP.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: E0=I/q,0 and E∗=N∗,P∗, where P∗=αIαb−Im−qma/m1−εααb−m and N∗=ma/αb−m. The sufficient and necessary conditions for that N∗,P∗ is positive are as follows:(9)αb−m>0,Iαb−Im−qma>0,εα<1.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 pointN∗,P∗, we just concern the stability of N∗,P∗. The corresponding Jacobian is given by(10)Je=JN∗,P∗=1−eq+αb−mαbI−Im−qmaabmα1−εαeεm−mαeαb−mαbI−Im−qmaabm1−εα1.The two eigenvalues of matrix (10) are as follows:(11)λ1,2=tre±tre2−4Δe2,where(12)tre=2−q+αb−mαbI−Im−qmaabmα1−εαe,Δe=1−q+αb−mαbI−Im−qmaabmα1−εαe+αb−mαbI−Im−qmaαabe2.According to the criterion of stability of the fixed point, ifλ1<1 and λ2<1, the fixed point is locally asymptotically stable. As shown in the previous research [48], this criterion of stability is equivalent to Δe<1,−1+Δe<tre<1+Δe. Through calculations, the conditions for N∗,P∗ to be locally asymptotically stable are determined as follows:(13)Ae2−Be<0,A>0,Ae2−2Be+4>0,where(14)A=αb−mαbI−Im−qmaαab,B=q+αb−mαbI−Im−qmaabmα1−εα.The dynamics of the discrete model without spatial diffusions can be shown in Figure1, 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 τ=0 is considered, as when time delay increases, the following bifurcation will occur.Figure 1
Fixed points of (8). Parameter values: I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, q=0.01, ε=0.001, and τ=0. ∗Asymptotically stable fixed point E∗; “○” represents the boundary fixed point E0.


## 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, Nτ=Nx,t−τ∂Nx,t/∂t, put it into equation (1) and obtain(15)∂N∂t=d1ΔN+I−bNa+NP−qN+εmP,∂P∂t=d2ΔP+αbPRNx,t−τ∂Nx,t∂t−mP,in which(16)RN=NN+a.Taylor expansion ofRNx,t−τ∂Nx,t/∂t at N,P without considering the higher order nonlinear part and equation (1) can be written as(17)∂Nx,t∂t=d1ΔN+f1N,P,∂Px,t∂t=−αabPτN+a2d1ΔN+d2ΔP+g2N,P−αabPτN+a2f1N,P,where f1N,P=I−bN/a+NP−qN+εmP and g2N,P=αbN/a+NP−mP. Therefore, equation (8) can be expressed as the following map:(18)NP⟶N+eI−ebNN+aP−eqN+eεmPP+eαbNPa+N−emP−αabPτN+a2eI−bNN+aP−qN+εmP.According to [48], the first condition of Neimark–Sacker bifurcation requires the eigenvalues at a fixed point, λ1 and λ2 are conjugate and their modules are 1, shown as(19)λ2=λ1¯,λ1=λ2=1,in which(20)λ1,2=tr0±tr02−4Δ02,(21)tr0=a11+a22=A0+A1+A2τ,Δ0=a11a22−a12a21=A0A1+A2τ−B0B1+B2τ,A0,A1,A2,B0,B1,B2,a11,a12,a21,anda22 have been shown in Appendix A.Then, we can get(22)tr02−4Δ0<0,Δ0=1.By calculation, we can obtain the Neimark–Sacker bifurcation pointτ0:(23)τ0=1+B0B1−A0A1A0A2−B0B2.Under the satisfaction of conditions (22) and (23), the fixed point E∗ of map (18) is translated to the origin by the following translation:(24)x=N−N∗,y=P−P∗.Then, map (18) is transformed into(25)xy⟶a11x+a12y+a13x2+a14y2+a15xy+a16x3+a17x2y+a18xy2+Ox+y4a21x+a22y+a23x2+a24y2+a25xy+a26x3+a27x2y+a28xy2+Ox+y4,where a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25,a26,a27,anda28 have been shown in Appendix A with τ=τ0.The second condition for Neimark–Sacker bifurcation requires(26)d=dλτdττ=τ0=A0A2−B0B22A0A1−B0B1+A0A2−B0B2τ0≠0,λτ0θ≠1,θ=1,2,3,4,in which(27)λτ0,λ¯τ0=tr0τ02±i24Δ0τ0−tr02τ0=β±iγ,i=−1,where tr0τ0 and Δ0τ0 are described in equation (21) with τ=τ0. Then, we can obtain(28)A0+A1+A2τ0≠0,−1.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:(29)xy=a120β−a11−γx˜y˜.To map (25), then the map becomes(30)x˜y˜⟶β−γγβx˜y˜+1a12γFx˜,y˜Gx˜,y˜,where(31)Fx˜,y˜=a12γa12a13+a15β−a11x˜2−a12a15γ2x˜y˜+a122γa12a16+a17β−a11x˜3−a122a17γ2x˜2y˜+Ox˜+y˜4,Gx˜,y˜=a122a13β−a11−a12a23+a14β−a11−a12a24β−a112+a12a15β−a11−a122a25β−a11x˜2−a12a24γ2y˜2+a12γ2a24−a15β−a11+a12a25x˜y˜+a123a16β−a11−a124a26+a122a17β−a11−a123a27β−a11−a122a28β−a112x˜3−a122γa17β−a11−a12a27−2a28β−a11x˜2y˜−a122a28γ2x˜y˜2+Ox˜+y˜4.The second-order and third-order partial derivatives ofFx˜,y˜ and Gx˜,y˜ at 0,0 are calculated as(32)Fx˜x˜=2a122a13γ+2a12a15γβ−a11,Fx˜y˜=−a12a15γ2,Fy˜y˜=0,Fx˜x˜x˜=6a123a16γ+6a122a17γβ−a11,Fx˜x˜y˜=−a122a17γ2,Fx˜y˜y˜=0,Fy˜y˜y˜=0,Gx˜x˜=2a122a13β−a11−a12a23+2a14β−a11−a12a24β−a112+2β−a11a12a15β−a11−a122a25,Gx˜y˜=2γa12a24β−a11−γa12a15β−a11−a122a25,Gy˜y˜=−2a12a24γ2,Gx˜x˜x˜=6a123a16β−a11−a124a26+6a122a17β−a11−a123a27β−a11−6a122a28β−a112,Gx˜x˜y˜=−2γa122a17β−a11−a123a27+4γa122a28β−a11,Gx˜y˜y˜=−2a122a28γ2,Gy˜y˜y˜=0.The third condition of Neimark–Sacker bifurcation requires(33)R=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21≠0,in which(34)ξ11=14a12γFx˜x˜+Fy˜y˜+iGx˜x˜+Gy˜y˜,ξ20=18a12γFx˜x˜−Fy˜y˜+2Gx˜y˜+iGx˜x˜−Gy˜y˜−2Fx˜y˜,ξ02=18a12γFx˜x˜−Fy˜y˜−2Gx˜y˜+iGx˜x˜−Gy˜y˜+2Fx˜y˜,ξ21=116a12γFx˜x˜x˜+Fx˜y˜y˜+Gx˜x˜y˜+Gy˜y˜y˜+iGx˜x˜x˜+Gx˜y˜y˜−Fx˜x˜y˜−Fy˜y˜y˜.After calculation, the third condition for Neimark–Sacker can be expressed as(35)R=−1A∗B∗Fx˜x˜+Fy˜y˜Fx˜x˜−Fy˜y˜+2Gx˜y˜−Gx˜x˜+Gy˜y˜Gx˜x˜−Gy˜y˜−2Fx˜y˜+C∗Gx˜x˜+Gy˜y˜Fx˜x˜−Fy˜y˜+2Gx˜y˜−Fx˜x˜+Fy˜y˜Gx˜x˜−Gy˜y˜−2Fx˜y˜−132a122γ2Fx˜x˜+Fy˜y˜2+Gx˜x˜+Gy˜y˜2+116a12γβFx˜x˜x˜+Fx˜y˜y˜+Gx˜x˜y˜+Gy˜y˜y˜+γGx˜x˜x˜+Gx˜y˜y˜−Fx˜x˜y˜−Fy˜y˜y˜−164a122γ2Fx˜x˜−Fy˜y˜−2Gx˜y˜2+Gx˜x˜−Gy˜y˜+2Fx˜y˜2≠0,in which(36)A∗=32a122γ21−β+γ2,B∗=1−3β+2β2−2γ2β2−γ2+6β−8β2γ2,C∗=2βγ1−β1−2β−2γ2−γ3−4ββ2−γ2.Based on the above calculations, when the conditions (23), (28), and (35) are satisfied, Neimark–Sacker bifurcation occurs at the fixed point E∗. Moreover, when R<0 and d>0, then an attracting invariant closed curve bifurcates from the fixed point E∗ for τ>τ0; otherwise, when R>0 and d>0, a repelling invariant closed curve bifurcates from the fixed point E∗ for 0<τ<τ0.In Figure2(a), the parameter space of time delay and input rate I−τ is shown. Zone I represents that when the input rate of nutrient is low, positive fixed point E∗ does not exist. Zone II and III represent no bifurcation and bifurcation area around the fixed point E∗. And Figure 2(b) shows the bifurcation diagram when I=1.6 and e=0.2. The Hopf bifurcation point is not suitable here with a time step e 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): e cannot be equal to 0.2, but in the discrete model, when e=0.2, Neimark–Sacker bifurcation occurs and drives the system from stable fixed point to quasiperiodic dynamics. When time delay τ=3 satisfy both τ>τN−S and τ>τHopf, the dynamics of the continuous model (with e=0.01 as in [46]) and the discrete model (e=0.2) 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.Figure 2
(a) Parameter space graph of Neimark–Sacker bifurcation and Hopf bifurcation withI and τ; (b) bifurcation diagram of Neimark–Sacker bifurcation with τ and phytoplankton; (c–e) phase portraits of Neimark–Sacker bifurcation with parameters (c) e=0.01, τ=3; (d) e=0.2, τ=2.45; (e) e=0.2, τ=3. Other parameter values: I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, q=0.01, and ε=0.001. ∗Asymptotically stable fixed point E∗.
(a)
(b)
(c)
(d)
(e)


### 2.3.2. Turing Bifurcation Analysis

In this part, we consider Turing bifurcation of equations (2)–(4) when τ=0. 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:(37)∇2Xij+λXij=0,with periodic boundary conditions:(38)Xi,0=Xi,k,Xi,1=Xi,k+1,X0,j=Xk,j,X1,j=Xk+1,j,i,j∈1,2,…,k.Whenx≠0 and y≠0, we set Xij=xiyj, put it into equation (38), and obtain xk=yk=1. Then, we can get(39)xr=exp2r−1πki′,yl=exp2l−1πki′,r,l∈1,2,…,k,i′=−1.SubstitutingXij=xiyj into equation (37), we obtain that(40)xi−1yj+xi+1yj+xiyj−1+xiyj+1−4xiyj+λxiyj=0.By calculatingλ=−1−x2/x−1−y2/y, then, using equation (39) we can get(41)λrl=−1−xr2xr−1−yl2yl=4sin2r−1πk+sin2l−1πk=Srl2.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 stateE∗. The equation for spatially nonuniform perturbations can be expressed as(42)N˜i,jt=Ni,jt−N∗,P˜i,jt=Pi,jt−P∗.Noticing∇d2N˜i,jt=∇d2Ni,jt and ∇d2P˜i,jt=∇d2Pi,jt, and the values of the two cannot be constant to 0. Substituting equation (42) into equations (2)–(4) leads to the following equation:(43)N˜i,jt+1=fNN˜i,jt+D1∇d2N˜i,jt+fPP˜i,jt+D2∇d2P˜i,jt+ON˜i,jt+P˜i,jt2,P˜i,jt+1=gNN˜i,jt+D1∇d2N˜i,jt+gPP˜i,jt+D2∇d2P˜i,jt+ON˜i,jt+P˜i,jt2.The high-order terms in the above equations can be ignored when the disturbance is very weak. Using the corresponding characteristic functionXrlij of the eigenvalue λrl to multiply equation (43), we obtain(44)XrlijN˜i,jt+1=fNXrlijN˜i,jt+fPXrlijP˜i,jt+D1fNXrlij∇d2N˜i,jt+D2fPXrlij∇d2P˜i,jt,XrlijP˜i,jt+1=gNXrlijN˜i,jt+gPXrlijP˜i,jt+D1gNXrlij∇d2N˜i,jt+D2gPXrlij∇d2P˜i,jt.Summing equation (44) for all i and j gives the following equation:(45)ΣXrlijN˜i,jt+1=fNΣXrlijN˜i,jt+fPΣXrlijP˜i,jt+D1fNΣXrlij∇d2N˜i,jt+D2fPΣXrlij∇d2P˜i,jt,ΣXrlijP˜i,jt+1=gNΣXrlijN˜i,jt+gPΣXrlijP˜i,jt+D1gNΣXrlij∇d2N˜i,jt+D2gPΣXrlij∇d2P˜i,jt.LetN¯t=∑i,j=1kXrlijN˜ijt+1 and P¯t=∑i,j=1kXrlijPijt+1, combined boundary conditions (38), equation (45) can be transformed into the following form:(46)N¯t+1=fN1−D1λrlN¯t+fP1−D2λrlP¯t,P¯t+1=gN1−D1λrlN¯t+gP1−D2λrlP¯t.At this time, we know that a progressive solution of equations (2)–(4) with boundary conditions (5) is Nijt=N¯tXrlijandPijt=P¯tXrlij. 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:(47)λ±r,l,e=12trr,l,e±12trr,l,e2−4Δr,l,e,in which(48)trr,l,e=a11e+a22e−eh2a11ed1+a22ed2λrl,Δr,l,e=a11ea22e−a12ea21e1−eh2d1λrl1−eh2d2λrl,JN∗,P∗=a11ea12ea21ea22e=1−eq−eabN∗+a2P∗eεm−ebN∗N∗+aeαabN∗+a2P∗1−em+eαbN∗N∗+a.Here,a11e, a12e, a21e, and a22e are denoted for reminding that they are dependent on e. Based on the two eigenvalues, we define(49)λmr,l,e=maxλ+r,l,λ−r,l,λme=maxr=1,l=1nλmr,l,e.The occurrence of Turing bifurcation needs to satisfyλme=1. According to this condition, we can determine the critical value e′ of Turing bifurcation occurring, which can be described as follows:(1)
Iftrr,l,e2>4Δr,l,e is satisfied in a small neighborhood of e=e′, the critical value e′ satisfies maxr=1,l=1ntrr,l,e′−Δr,l,e′=1(2)
Iftrr,l,e2≤4Δr,l,e is satisfied in a small neighborhood of e=e′, the critical value e′ satisfies maxr=1,l=1nΔr,l,e′=1max
r
=
1
,
l
=
1
n
t
r
r
,
l
,
e
′
−
Δ
r
,
l
,
e
′
=
1 and maxr=1,l=1nΔr,l,e′=1 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 λme>1.Then, we discuss the influence of the diffusion coefficientsd1 and d2 on the model. We use the following parameter values: I=1.6, a=0.5, b=0.8, α=0.7, q=0.01, ε=0.001, m=0.45, d1=0.2, d2=0.2, e=0.2, h=1, τ=0, and t=100. 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 λk,l. In Figure 3(a), we can see that the perturbation parameters k and l have symmetrical effects on the eigenvalues. Therefore, for convenience of observation, let k=l, then we can get the change in the eigenvalue to l, 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 d1=1.220, then the eigenvalue of the system will not make the eigenvalue λk,l greater than 1 with the increase of the perturbation parameter; at this time, Turing bifurcation will not occur. When d1=1.389, the curve is tangent to the straight line λk,l=1, that is, only if the diffusion coefficient d1>1.389, the value of the eigenvalue λk,l will exceed 1 as the disturbance coefficient increases; at this time, Turing bifurcation will occur in the system.Figure 3
Variation of eigenvaluesλk,l with perturbation numbers k and l. Parameter values: I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d2=0.20, q=0.01, τ=0, h=1, ε=0.001, and n=100.
(a)
(b)


## 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, Nτ=Nx,t−τ∂Nx,t/∂t, put it into equation (1) and obtain(15)∂N∂t=d1ΔN+I−bNa+NP−qN+εmP,∂P∂t=d2ΔP+αbPRNx,t−τ∂Nx,t∂t−mP,in which(16)RN=NN+a.Taylor expansion ofRNx,t−τ∂Nx,t/∂t at N,P without considering the higher order nonlinear part and equation (1) can be written as(17)∂Nx,t∂t=d1ΔN+f1N,P,∂Px,t∂t=−αabPτN+a2d1ΔN+d2ΔP+g2N,P−αabPτN+a2f1N,P,where f1N,P=I−bN/a+NP−qN+εmP and g2N,P=αbN/a+NP−mP. Therefore, equation (8) can be expressed as the following map:(18)NP⟶N+eI−ebNN+aP−eqN+eεmPP+eαbNPa+N−emP−αabPτN+a2eI−bNN+aP−qN+εmP.According to [48], the first condition of Neimark–Sacker bifurcation requires the eigenvalues at a fixed point, λ1 and λ2 are conjugate and their modules are 1, shown as(19)λ2=λ1¯,λ1=λ2=1,in which(20)λ1,2=tr0±tr02−4Δ02,(21)tr0=a11+a22=A0+A1+A2τ,Δ0=a11a22−a12a21=A0A1+A2τ−B0B1+B2τ,A0,A1,A2,B0,B1,B2,a11,a12,a21,anda22 have been shown in Appendix A.Then, we can get(22)tr02−4Δ0<0,Δ0=1.By calculation, we can obtain the Neimark–Sacker bifurcation pointτ0:(23)τ0=1+B0B1−A0A1A0A2−B0B2.Under the satisfaction of conditions (22) and (23), the fixed point E∗ of map (18) is translated to the origin by the following translation:(24)x=N−N∗,y=P−P∗.Then, map (18) is transformed into(25)xy⟶a11x+a12y+a13x2+a14y2+a15xy+a16x3+a17x2y+a18xy2+Ox+y4a21x+a22y+a23x2+a24y2+a25xy+a26x3+a27x2y+a28xy2+Ox+y4,where a11,a12,a13,a14,a15,a16,a17,a18,a21,a22,a23,a24,a25,a26,a27,anda28 have been shown in Appendix A with τ=τ0.The second condition for Neimark–Sacker bifurcation requires(26)d=dλτdττ=τ0=A0A2−B0B22A0A1−B0B1+A0A2−B0B2τ0≠0,λτ0θ≠1,θ=1,2,3,4,in which(27)λτ0,λ¯τ0=tr0τ02±i24Δ0τ0−tr02τ0=β±iγ,i=−1,where tr0τ0 and Δ0τ0 are described in equation (21) with τ=τ0. Then, we can obtain(28)A0+A1+A2τ0≠0,−1.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:(29)xy=a120β−a11−γx˜y˜.To map (25), then the map becomes(30)x˜y˜⟶β−γγβx˜y˜+1a12γFx˜,y˜Gx˜,y˜,where(31)Fx˜,y˜=a12γa12a13+a15β−a11x˜2−a12a15γ2x˜y˜+a122γa12a16+a17β−a11x˜3−a122a17γ2x˜2y˜+Ox˜+y˜4,Gx˜,y˜=a122a13β−a11−a12a23+a14β−a11−a12a24β−a112+a12a15β−a11−a122a25β−a11x˜2−a12a24γ2y˜2+a12γ2a24−a15β−a11+a12a25x˜y˜+a123a16β−a11−a124a26+a122a17β−a11−a123a27β−a11−a122a28β−a112x˜3−a122γa17β−a11−a12a27−2a28β−a11x˜2y˜−a122a28γ2x˜y˜2+Ox˜+y˜4.The second-order and third-order partial derivatives ofFx˜,y˜ and Gx˜,y˜ at 0,0 are calculated as(32)Fx˜x˜=2a122a13γ+2a12a15γβ−a11,Fx˜y˜=−a12a15γ2,Fy˜y˜=0,Fx˜x˜x˜=6a123a16γ+6a122a17γβ−a11,Fx˜x˜y˜=−a122a17γ2,Fx˜y˜y˜=0,Fy˜y˜y˜=0,Gx˜x˜=2a122a13β−a11−a12a23+2a14β−a11−a12a24β−a112+2β−a11a12a15β−a11−a122a25,Gx˜y˜=2γa12a24β−a11−γa12a15β−a11−a122a25,Gy˜y˜=−2a12a24γ2,Gx˜x˜x˜=6a123a16β−a11−a124a26+6a122a17β−a11−a123a27β−a11−6a122a28β−a112,Gx˜x˜y˜=−2γa122a17β−a11−a123a27+4γa122a28β−a11,Gx˜y˜y˜=−2a122a28γ2,Gy˜y˜y˜=0.The third condition of Neimark–Sacker bifurcation requires(33)R=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21≠0,in which(34)ξ11=14a12γFx˜x˜+Fy˜y˜+iGx˜x˜+Gy˜y˜,ξ20=18a12γFx˜x˜−Fy˜y˜+2Gx˜y˜+iGx˜x˜−Gy˜y˜−2Fx˜y˜,ξ02=18a12γFx˜x˜−Fy˜y˜−2Gx˜y˜+iGx˜x˜−Gy˜y˜+2Fx˜y˜,ξ21=116a12γFx˜x˜x˜+Fx˜y˜y˜+Gx˜x˜y˜+Gy˜y˜y˜+iGx˜x˜x˜+Gx˜y˜y˜−Fx˜x˜y˜−Fy˜y˜y˜.After calculation, the third condition for Neimark–Sacker can be expressed as(35)R=−1A∗B∗Fx˜x˜+Fy˜y˜Fx˜x˜−Fy˜y˜+2Gx˜y˜−Gx˜x˜+Gy˜y˜Gx˜x˜−Gy˜y˜−2Fx˜y˜+C∗Gx˜x˜+Gy˜y˜Fx˜x˜−Fy˜y˜+2Gx˜y˜−Fx˜x˜+Fy˜y˜Gx˜x˜−Gy˜y˜−2Fx˜y˜−132a122γ2Fx˜x˜+Fy˜y˜2+Gx˜x˜+Gy˜y˜2+116a12γβFx˜x˜x˜+Fx˜y˜y˜+Gx˜x˜y˜+Gy˜y˜y˜+γGx˜x˜x˜+Gx˜y˜y˜−Fx˜x˜y˜−Fy˜y˜y˜−164a122γ2Fx˜x˜−Fy˜y˜−2Gx˜y˜2+Gx˜x˜−Gy˜y˜+2Fx˜y˜2≠0,in which(36)A∗=32a122γ21−β+γ2,B∗=1−3β+2β2−2γ2β2−γ2+6β−8β2γ2,C∗=2βγ1−β1−2β−2γ2−γ3−4ββ2−γ2.Based on the above calculations, when the conditions (23), (28), and (35) are satisfied, Neimark–Sacker bifurcation occurs at the fixed point E∗. Moreover, when R<0 and d>0, then an attracting invariant closed curve bifurcates from the fixed point E∗ for τ>τ0; otherwise, when R>0 and d>0, a repelling invariant closed curve bifurcates from the fixed point E∗ for 0<τ<τ0.In Figure2(a), the parameter space of time delay and input rate I−τ is shown. Zone I represents that when the input rate of nutrient is low, positive fixed point E∗ does not exist. Zone II and III represent no bifurcation and bifurcation area around the fixed point E∗. And Figure 2(b) shows the bifurcation diagram when I=1.6 and e=0.2. The Hopf bifurcation point is not suitable here with a time step e 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): e cannot be equal to 0.2, but in the discrete model, when e=0.2, Neimark–Sacker bifurcation occurs and drives the system from stable fixed point to quasiperiodic dynamics. When time delay τ=3 satisfy both τ>τN−S and τ>τHopf, the dynamics of the continuous model (with e=0.01 as in [46]) and the discrete model (e=0.2) 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.Figure 2
(a) Parameter space graph of Neimark–Sacker bifurcation and Hopf bifurcation withI and τ; (b) bifurcation diagram of Neimark–Sacker bifurcation with τ and phytoplankton; (c–e) phase portraits of Neimark–Sacker bifurcation with parameters (c) e=0.01, τ=3; (d) e=0.2, τ=2.45; (e) e=0.2, τ=3. Other parameter values: I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, q=0.01, and ε=0.001. ∗Asymptotically stable fixed point E∗.
(a)
(b)
(c)
(d)
(e)


## 2.3.2. Turing Bifurcation Analysis

In this part, we consider Turing bifurcation of equations (2)–(4) when τ=0. 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:(37)∇2Xij+λXij=0,with periodic boundary conditions:(38)Xi,0=Xi,k,Xi,1=Xi,k+1,X0,j=Xk,j,X1,j=Xk+1,j,i,j∈1,2,…,k.Whenx≠0 and y≠0, we set Xij=xiyj, put it into equation (38), and obtain xk=yk=1. Then, we can get(39)xr=exp2r−1πki′,yl=exp2l−1πki′,r,l∈1,2,…,k,i′=−1.SubstitutingXij=xiyj into equation (37), we obtain that(40)xi−1yj+xi+1yj+xiyj−1+xiyj+1−4xiyj+λxiyj=0.By calculatingλ=−1−x2/x−1−y2/y, then, using equation (39) we can get(41)λrl=−1−xr2xr−1−yl2yl=4sin2r−1πk+sin2l−1πk=Srl2.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 stateE∗. The equation for spatially nonuniform perturbations can be expressed as(42)N˜i,jt=Ni,jt−N∗,P˜i,jt=Pi,jt−P∗.Noticing∇d2N˜i,jt=∇d2Ni,jt and ∇d2P˜i,jt=∇d2Pi,jt, and the values of the two cannot be constant to 0. Substituting equation (42) into equations (2)–(4) leads to the following equation:(43)N˜i,jt+1=fNN˜i,jt+D1∇d2N˜i,jt+fPP˜i,jt+D2∇d2P˜i,jt+ON˜i,jt+P˜i,jt2,P˜i,jt+1=gNN˜i,jt+D1∇d2N˜i,jt+gPP˜i,jt+D2∇d2P˜i,jt+ON˜i,jt+P˜i,jt2.The high-order terms in the above equations can be ignored when the disturbance is very weak. Using the corresponding characteristic functionXrlij of the eigenvalue λrl to multiply equation (43), we obtain(44)XrlijN˜i,jt+1=fNXrlijN˜i,jt+fPXrlijP˜i,jt+D1fNXrlij∇d2N˜i,jt+D2fPXrlij∇d2P˜i,jt,XrlijP˜i,jt+1=gNXrlijN˜i,jt+gPXrlijP˜i,jt+D1gNXrlij∇d2N˜i,jt+D2gPXrlij∇d2P˜i,jt.Summing equation (44) for all i and j gives the following equation:(45)ΣXrlijN˜i,jt+1=fNΣXrlijN˜i,jt+fPΣXrlijP˜i,jt+D1fNΣXrlij∇d2N˜i,jt+D2fPΣXrlij∇d2P˜i,jt,ΣXrlijP˜i,jt+1=gNΣXrlijN˜i,jt+gPΣXrlijP˜i,jt+D1gNΣXrlij∇d2N˜i,jt+D2gPΣXrlij∇d2P˜i,jt.LetN¯t=∑i,j=1kXrlijN˜ijt+1 and P¯t=∑i,j=1kXrlijPijt+1, combined boundary conditions (38), equation (45) can be transformed into the following form:(46)N¯t+1=fN1−D1λrlN¯t+fP1−D2λrlP¯t,P¯t+1=gN1−D1λrlN¯t+gP1−D2λrlP¯t.At this time, we know that a progressive solution of equations (2)–(4) with boundary conditions (5) is Nijt=N¯tXrlijandPijt=P¯tXrlij. 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:(47)λ±r,l,e=12trr,l,e±12trr,l,e2−4Δr,l,e,in which(48)trr,l,e=a11e+a22e−eh2a11ed1+a22ed2λrl,Δr,l,e=a11ea22e−a12ea21e1−eh2d1λrl1−eh2d2λrl,JN∗,P∗=a11ea12ea21ea22e=1−eq−eabN∗+a2P∗eεm−ebN∗N∗+aeαabN∗+a2P∗1−em+eαbN∗N∗+a.Here,a11e, a12e, a21e, and a22e are denoted for reminding that they are dependent on e. Based on the two eigenvalues, we define(49)λmr,l,e=maxλ+r,l,λ−r,l,λme=maxr=1,l=1nλmr,l,e.The occurrence of Turing bifurcation needs to satisfyλme=1. According to this condition, we can determine the critical value e′ of Turing bifurcation occurring, which can be described as follows:(1)
Iftrr,l,e2>4Δr,l,e is satisfied in a small neighborhood of e=e′, the critical value e′ satisfies maxr=1,l=1ntrr,l,e′−Δr,l,e′=1(2)
Iftrr,l,e2≤4Δr,l,e is satisfied in a small neighborhood of e=e′, the critical value e′ satisfies maxr=1,l=1nΔr,l,e′=1max
r
=
1
,
l
=
1
n
t
r
r
,
l
,
e
′
−
Δ
r
,
l
,
e
′
=
1 and maxr=1,l=1nΔr,l,e′=1 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 λme>1.Then, we discuss the influence of the diffusion coefficientsd1 and d2 on the model. We use the following parameter values: I=1.6, a=0.5, b=0.8, α=0.7, q=0.01, ε=0.001, m=0.45, d1=0.2, d2=0.2, e=0.2, h=1, τ=0, and t=100. 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 λk,l. In Figure 3(a), we can see that the perturbation parameters k and l have symmetrical effects on the eigenvalues. Therefore, for convenience of observation, let k=l, then we can get the change in the eigenvalue to l, 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 d1=1.220, then the eigenvalue of the system will not make the eigenvalue λk,l greater than 1 with the increase of the perturbation parameter; at this time, Turing bifurcation will not occur. When d1=1.389, the curve is tangent to the straight line λk,l=1, that is, only if the diffusion coefficient d1>1.389, the value of the eigenvalue λk,l will exceed 1 as the disturbance coefficient increases; at this time, Turing bifurcation will occur in the system.Figure 3
Variation of eigenvaluesλk,l with perturbation numbers k and l. Parameter values: I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d2=0.20, q=0.01, τ=0, h=1, ε=0.001, and n=100.
(a)
(b)


## 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 d1 and d2. 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 Figure4, with given parameter values, the critical value of can be calculated as τ0≈2.2703. Thus, we take τ=0, τ=1, and τ=3 and obtain Figures 4(a)–4(c) respectively. Figures 4(a) and 4(b) show that when 0≤τ<τ0, there is no spatial patterns. Figure 4(c) clearly shows the formation of a Neimark–Sacker type pattern when τ>τ0. From the simulated patterns, we can see that simulations are consistent with the Neimark–Sacker analysis.Figure 4
Spatial distribution of phytoplankton population under different time delay parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, d1=0.20, d2=0.20, q=0.01, ε=0.001, and t=800. (a) τ=0; (b) τ=1; (c) τ=3.
(a)
(b)
(c)Figure5 shows a time series of phytoplankton pattern formation process from initial conditions to self-organization when τ=3. 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)).Figure 5
The process of forming cloud-like patterns over time with a time delay of 3. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, d1=0.2, d2=0.2, q=0.01, τ=3, and ε=0.001. (a) t=0; (b) t=500; (c) t=2500; (d) t=5000; (e) t=7500; (f) t=10000; (g) t=12500; (h) t=15000; (i) t=20000.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


### 3.2. Effects of Diffusion Coefficients on the Pattern Formation

Based on the Turing bifurcation analysis in Section2.3.2, we can see that the diffusion coefficient d1 can induce the occurrence of Turing bifurcations, and the critical value is d1=1.389. 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 d2 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.Figure 6
Spatial distribution of phytoplankton population under different diffusion parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, e=0.2, h=1, q=0.01, ε=0.001, τ=0, and t=500. (a) d1=1.220, d2=0.20; (b) d1=1.390, d2=0.20; (c) d1=1.389, d2=0.21.
(a)
(b)
(c)In order to show the formation process of Turing-type patterns, two sets of values are selected:d1=1.390, d2=0.20 and d1=1.389, d2=0.21. 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.Figure 7
The formation of phytoplankton population spatial distribution pattern under certain diffusion parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=0, h=1, and ε=0.001. (a) t=5; (b) t=25; (c) t=30; (d) t=35; (e) t=40; (f) t=50; (g) t=100; (h) t=400; (i) t=500.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 8
The formation of phytoplankton population spatial distribution pattern under certain diffusion parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.389, d2=0.21, q=0.01, τ=0, h=1, and ε=0.001. (a) t=5; (b) t=25; (c) t=30; (d) t=35; (e) t=50; (f) t=75; (g) t=100; (h) t=300; (i) t=500.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Note that when time delay satisfies the Neimark-Sacker bifurcation condition and coefficientsd1 and d2 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 τ=0 (see Figure 7). Note that as shown in Figure 9(d) with t=50, the same formation stage of patterns can be obtained in Figure 7(e) with t=40. 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].Figure 9
The formation of phytoplankton population spatial distribution pattern under certain parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, h=1, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=3, n=200, and ε=0.001. (a) t=5; (b) t=25; (c) t=40; (d) t=50; (e) t=60; (f) t=75; (g) t=250; (h) t=400; (i) t=500.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


### 3.3. Comparison between Simulated Patterns and Real Observed Patterns

From numerical simulations, mainly two types of patterns are obtained: cloud-like patterns (Figure10(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)).Figure 10
Cloud-like patterns. Parameter values: (a)I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, d1=0.2, d2=0.2, q=0.01, τ=3, ε=0.001, t=20000, n=100, e=0.2, and h=1. (b) On June 24, 2010, in the North Atlantic, phytoplankton was photographed on a medium resolution imaging spectrometer on NASA’s “Aqua” satellite.
(a)
(b)Figure 11
Banded patterns. Parameter values: (a)I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.389, d2=0.21, q=0.01, τ=0, h=1, ε=0.001, t=100, and n=100. (b) On November 25, 2014, the phytoplankton bloom photographed in the sea near Shenzhen, China; (c) I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=0, h=1, ε=0.001, t=100, and n=100. (d) On May 7, 2012, noctilucent algae bloomed in the sea near Rizhao, Shandong, China; (e) I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=3, h=1, ε=0.001, t=500, and n=300. (f) On November 25, 2014, the phytoplankton bloom photographed in the sea near Shenzhen, China.
(a)
(b)
(c)
(d)
(e)
(f)Figures11(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.


## 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 Figure4, with given parameter values, the critical value of can be calculated as τ0≈2.2703. Thus, we take τ=0, τ=1, and τ=3 and obtain Figures 4(a)–4(c) respectively. Figures 4(a) and 4(b) show that when 0≤τ<τ0, there is no spatial patterns. Figure 4(c) clearly shows the formation of a Neimark–Sacker type pattern when τ>τ0. From the simulated patterns, we can see that simulations are consistent with the Neimark–Sacker analysis.Figure 4
Spatial distribution of phytoplankton population under different time delay parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, d1=0.20, d2=0.20, q=0.01, ε=0.001, and t=800. (a) τ=0; (b) τ=1; (c) τ=3.
(a)
(b)
(c)Figure5 shows a time series of phytoplankton pattern formation process from initial conditions to self-organization when τ=3. 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)).Figure 5
The process of forming cloud-like patterns over time with a time delay of 3. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, d1=0.2, d2=0.2, q=0.01, τ=3, and ε=0.001. (a) t=0; (b) t=500; (c) t=2500; (d) t=5000; (e) t=7500; (f) t=10000; (g) t=12500; (h) t=15000; (i) t=20000.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 3.2. Effects of Diffusion Coefficients on the Pattern Formation

Based on the Turing bifurcation analysis in Section2.3.2, we can see that the diffusion coefficient d1 can induce the occurrence of Turing bifurcations, and the critical value is d1=1.389. 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 d2 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.Figure 6
Spatial distribution of phytoplankton population under different diffusion parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, e=0.2, h=1, q=0.01, ε=0.001, τ=0, and t=500. (a) d1=1.220, d2=0.20; (b) d1=1.390, d2=0.20; (c) d1=1.389, d2=0.21.
(a)
(b)
(c)In order to show the formation process of Turing-type patterns, two sets of values are selected:d1=1.390, d2=0.20 and d1=1.389, d2=0.21. 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.Figure 7
The formation of phytoplankton population spatial distribution pattern under certain diffusion parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=0, h=1, and ε=0.001. (a) t=5; (b) t=25; (c) t=30; (d) t=35; (e) t=40; (f) t=50; (g) t=100; (h) t=400; (i) t=500.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 8
The formation of phytoplankton population spatial distribution pattern under certain diffusion parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.389, d2=0.21, q=0.01, τ=0, h=1, and ε=0.001. (a) t=5; (b) t=25; (c) t=30; (d) t=35; (e) t=50; (f) t=75; (g) t=100; (h) t=300; (i) t=500.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Note that when time delay satisfies the Neimark-Sacker bifurcation condition and coefficientsd1 and d2 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 τ=0 (see Figure 7). Note that as shown in Figure 9(d) with t=50, the same formation stage of patterns can be obtained in Figure 7(e) with t=40. 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].Figure 9
The formation of phytoplankton population spatial distribution pattern under certain parameters. Parameter values:I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, h=1, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=3, n=200, and ε=0.001. (a) t=5; (b) t=25; (c) t=40; (d) t=50; (e) t=60; (f) t=75; (g) t=250; (h) t=400; (i) t=500.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 3.3. Comparison between Simulated Patterns and Real Observed Patterns

From numerical simulations, mainly two types of patterns are obtained: cloud-like patterns (Figure10(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)).Figure 10
Cloud-like patterns. Parameter values: (a)I=1.6, a=0.5, b=0.8, α=0.7, m=0.45, d1=0.2, d2=0.2, q=0.01, τ=3, ε=0.001, t=20000, n=100, e=0.2, and h=1. (b) On June 24, 2010, in the North Atlantic, phytoplankton was photographed on a medium resolution imaging spectrometer on NASA’s “Aqua” satellite.
(a)
(b)Figure 11
Banded patterns. Parameter values: (a)I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.389, d2=0.21, q=0.01, τ=0, h=1, ε=0.001, t=100, and n=100. (b) On November 25, 2014, the phytoplankton bloom photographed in the sea near Shenzhen, China; (c) I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=0, h=1, ε=0.001, t=100, and n=100. (d) On May 7, 2012, noctilucent algae bloomed in the sea near Rizhao, Shandong, China; (e) I=1.6, a=0.5, b=0.8, α=0.7, e=0.2, m=0.45, d1=1.390, d2=0.20, q=0.01, τ=3, h=1, ε=0.001, t=500, and n=300. (f) On November 25, 2014, the phytoplankton bloom photographed in the sea near Shenzhen, China.
(a)
(b)
(c)
(d)
(e)
(f)Figures11(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 E0=I/q,0, but the time delay may affect the whole process; (2) when time delay exists and is greater than a certain critical value τ0, the time delay can not only lead to the instability of the stable fixed point E∗, but also form a cloud Neimark–Sacker pattern through Neimark–Sacker instability; (3) when d1≠d2, by changing the value of d1 or d2 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, Figure10(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.



---

*Source: 8541432-2020-01-09.xml*