Emergence of spatiotemporal chaos driven by far-field breakup of spiral waves in the plankton ecological systems

TL;DR

This paper investigates how spiral wave breakup leads to spatiotemporal chaos in a reaction-diffusion model of plankton, revealing complex patterns driven by parameter changes and wave interactions.

Contribution

It provides a detailed analysis of spiral wave breakup and chaos emergence in a plankton model using simulations and bifurcation theory, expanding understanding of ecological pattern formation.

Findings

Spiral waves exist in the mixed Turing-Hopf bifurcation region.

Far-field breakup of spiral waves causes spatiotemporal chaos.

Chaos does not spread uniformly across the entire spatial domain.

Abstract

Alexander B. Medvinsky \emph{et al} [A. B. Medvinsky, I. A. Tikhonova, R. R. Aliev, B.-L. Li, Z.-S. Lin, and H. Malchow, Phys. Rev. E \textbf{64}, 021915 (2001)] and Marcus R. Garvie \emph{et al} [M. R. Garvie and C. Trenchea, SIAM J. Control. Optim. \textbf{46}, 775-791 (2007)] shown that the minimal spatially extended reaction-diffusion model of phytoplankton-zooplankton can exhibit both regular, chaotic behavior, and spatiotemporal patterns in a patchy environment. Based on that, the spatial plankton model is furtherly investigated by means of computer simulations and theoretical analysis in the present paper when its parameters would be expected in the case of mixed Turing-Hopf bifurcation region. Our results show that the spiral waves exist in that region and the spatiotemporal chaos emerge, which arise from the far-field breakup of the spiral waves over large ranges of diffusion coefficients of phytoplankton and zooplankton. Moreover, the spatiotemporal chaos arising from the far-field breakup of spiral waves does not gradually involve the whole space within that region. Our results are confirmed by means of computation spectra and nonlinear bifurcation of wave trains. Finally, we give some explanations about the spatially structured patterns from the community level.

Figures

Equations

$$\frac{\partial p}{\partial t}=rp(1-p)-\frac{ap}{1+bp}h+d_p{\nabla }^{2}p,$$

$$\frac{\partial h}{\partial t}=\frac{ap}{1+bp}h-mh-f\frac{n{h}^2}{n^2+h^2}+d_h{\nabla }^{2}h,$$

$$g_1(p,h)=rp(1-p)-\frac{ap}{1+bp}h,$$

$$g_2(p,h)=\frac{ap}{1+bp}h-mh-f\frac{n{h}^2}{n^2+h^2}.$$

$$\frac{\partial p}{\partial t}=d_p{\nabla }^{2}p+g_1(p,h),$$

$$\frac{\partial h}{\partial t}=d_h{\nabla }^{2}h+g_2(p,h).$$

$$\frac{\partial p}{\partial t}=d_p{\nabla }_{\rho ,\theta }^{2}p+\omega \frac{\partial p}{\partial \theta }+g_1(p^{\ast },h^{\ast }),$$

$$\frac{\partial h}{\partial t}=d_h{\nabla }_{\rho ,\theta }^{2}h+\omega \frac{\partial h}{\partial \theta }+g_2(p^{\ast },h^{\ast }),$$

$$d_p{\nabla }_{\rho ,\theta }^{2}p+\omega \frac{\partial p}{\partial \theta }+g_1^p(p^{\ast },h^{\ast })p+g_1^h(p^{\ast },h^{\ast })h=\lambda p,$$

$$d_h{\nabla }_{\rho ,\theta }^{2}h+\omega \frac{\partial h}{\partial \theta }+g_2^p(p^{\ast },h^{\ast })p+g_2^h(p^{\ast },h^{\ast })h=\lambda h,$$

$$d_p{\nabla }_{\rho ,\rho }^{2}p+\omega \frac{\partial p}{\partial \theta }+g_1^p(p^{\ast },h^{\ast })p+g_1^h(p^{\ast },h^{\ast })h=\lambda p,$$

$$d_h{\nabla }_{\rho ,\rho }^{2}h+\omega \frac{\partial h}{\partial \theta }+g_2^p(p^{\ast },h^{\ast })p+g_2^h(p^{\ast },h^{\ast })h=\lambda h,$$

$$d_p{k}^2{\nabla }_{\xi ,\xi }^{2}p+\omega p_{\xi }+g_1^p(u_{wt}(\xi ))p+g_1^h(u_{wt}(\xi ))h=\lambda p,$$

$$d_h{k}^2{\nabla }_{\xi ,\xi }^{2}h+\omega h_{\xi }+g_2^p(u_{wt}(\xi ))p+g_2^h(u_{wt}(\xi ))h=\lambda h.$$

$$U(\rho ,\cdot )=e^{\lambda t}u(k\rho -\omega t),U(\rho ,T)=e^{\lambda T}u(k\rho -2\pi ).$$

$$\frac{dp}{d\rho }$$

$$\frac{dh}{d\rho }$$

$$\frac{d{p}_1}{d\rho }$$

$$\frac{d{h}_1}{d\rho }$$

$$\mathcal{A}(\lambda ,k):=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ \frac{1}{k^2{d}_p}(\lambda -g_1^p(u_{wt}(\xi )))& -\frac{1}{k^2{d}_p}{g}_1^h(u_{wt}(\xi ))& -\frac{1}{k^2{d}_p}\omega & 0\\ -\frac{1}{k^2{d}_h}{g}_2^p(u_{wt}(\xi ))& \frac{1}{k^2{d}_h}(\lambda -g_2^h(u_{wt}(\xi )))& 0& -\frac{1}{k^2{d}_h}\omega \end{array}\right),$$

$${\Sigma }_{ess}=\{\lambda \in \mathcal{C};\text{det}(\mathcal{A}(\lambda ,k)-\text{i}k)=0\text{for some}k\in \mathcal{R}\}.$$

$$c_p=\frac{\omega }{k}.$$

$$d_p\frac{d^{2}P}{d{z}^2}+c_p\frac{dP}{dz}+g_1(P,H)=0,$$

$$h_p\frac{d^{2}H}{d{z}^2}+c_p\frac{dH}{dz}+g_2(P,H)=0.$$

Full text

Spiral waves; Spatio-temporal pattern; Plankton dynamics; Reaction-diffusion system

I Introduction

There is a growing interest in the spatial pattern dynamics of ecological systems amritkar:258102 ; pekalski:021913 ; Sayama2003 ; gilad:098105 ; mobilia2006 ; mobilia2005 ; Blasius1999 ; PhysRevLett.87.198101 ; PhysRevE.67.056602 ; PhysRevE.64.021915 ; medvinsky:311 ; Gurney1998 ; Murray2002 . However, many mechanisms of the spatio-temporal variability of natural plankton populations are not known yet. Pronounced physical patterns like thermoclines, upwelling, fronts and eddies often set the frame for the biological process. Measurements of the underwater light field are made with state-of-the-art instruments and used to calculate concentrations of phytoplankton biomass (as chlorophyll) as well as other forms of organic matter. Very high diffusion of the marine environment would prevent the formation of any stable patch spatial distribution with much longer life-time than the typical time of biodynamics. Meanwhile, in addition to very changeable transient spatial patterns, there also exist other spatial patterns in marine environment, much more stable spatial structure associated with ocean fronts, spatiotemporal chaos medvinsky:311 ; PhysRevE.64.021915 ; Petrovskii , cyclonic rings, and so called meddies Laurence . In fact, it is significant to create the biological basis for understanding spatial patterns of plankton EsaRanta11281997 . For instance, the impact of space on the persistence of enriched ecological systems was proved in laboratory experiments Luckinbill . Recently, it has been shown both in laboratory experiments Holyoak and theoretically Jansen ; Jansen2000 ; Jansen2001 ; Petrovskii that the existence of a spatial structure makes a predator-prey system less prone to extinction. This is due to the temporal variations of the density of different sub-populations can become asynchronous and the events of local extinction can be compensated due to re-colonization from other sites in the space Allen2001 . During a long period of time, all the spiral waves have been widely observed in diverse physical, chemical, and biological systems RevModPhys.65.851 ; PhysRevLett.87.068101 ; Sawai2005 ; winfree:303 . However, a quite limited number of documents Gurney1998 ; medvinsky:311 ; Biktashev ; Garvie2007 ; garvie:775 concern the spiral wave pattern and its breakup in the ecological systems.

The investigation of transition from regular patterns to spatiotemporally chaotic dynamics in spatially extended systems remains a challenge in nonlinear science Markus2004 ; PhysRevLett.82.1160 ; RevModPhys.65.851 ; Petrovskii . In a nonlinear ecology system, the two most commonly seen patterns are spiral waves and turbulence (spatio-temporal chaos) for the level of the community Craigjohnson2002 . It has been recently shown that spontaneous spatiatemoporal pattern formation is an instrinsic property of a predator-prey system Pascual1993 ; Sherratt1997 ; Petrovskii2001 ; medvinsky:311 ; Petrovskii ; Sherratt2001 and spatiotemporal structures play an important role in ecological systems. For example, spatially induced speciation prevents the extinction of the predator-prey models Savill1998 ; Gurney1998 ; medvinsky:311 . So far, plankton patchiness has been observed on a wide range of spatial temporal scales Abraham1998 ; Folt1999 . There exist various, often heuristic explanations of the spatial patterns phenomenon for these systems. It should be noted that, although conclusive evidence of ecological chaos is still to be found, there is a growing number of indications of chaos in real ecosystems MartenScheffer01011991 ; Hanski1993 ; Ellner1995 ; Dennis2001 .

Recently developed models show that spatial self-structuring in multispecies systems can meet both criteria and provide a rich substrate for community-level section and a major transition in evolution. In present paper, the scenario in the spatially extended plankton ecological system is observed by means of the numerical simulation. The system has been demonstrated to exhibit regular or chaostic, depending on the initial conditions and the parameter values PhysRevE.64.021915 ; garvie:775 . We find that the far-field breakup of the spiral wave leads to complex spatiotemporal chaos (or a turbulentlike state) in the spatially extended plankton model (1). Our results show that regular spiral wave pattern shifts into spatiotemporal chaos pattern by modulating the diffusion coefficients of the species.

II model

In this paper we study the spatially extended nutrient-phytoplankton-zooplankton-fish reaction-diffusion system. Following Scheffer’s minimal approach Scheffer1991a , which was originally formulated as a system of ordinary diffential equation (ODEs) and later developed models PhysRevE.64.021915 ; medvinsky:311 ; Malchow1993 ; MercedesPascual ; garvie:775 , as a further investigation, we study a two-variable phytoplankton and zooplankton model on the level of the community to describe pattern formation with the diffusion. The dimensionless model is written as

[TABLE]

where the parameters are $r$, $a$, $b$, $m$, $n$, $d_{p}$, $d_{h}$, and $f$ which refer to work in Refs. PhysRevE.64.021915 ; medvinsky:311 . The explanation of model (1) relates to the nutrient-phytoplankton-zooplankton-fish ecological system [see Refs. garvie:775 ; Scheffer1991a ; PhysRevE.64.021915 for details]. The local dynamics are given by

[TABLE]

From the earlier results Malchow1993 about non-spatial system of model (1) by means of numerical bifurcation analysis show that the bifurcation and bistability can be found in the system (1) when the parameters are varied within a realistic range. For the fixed parameters (see the caption of Fig. 1 and 2), we can see that the $f$ controls the distance from Hopf bifurcation. For larger $f$, there exists only one stable steady state. As $f$ is decreased further, the homogeneous steady state undergoes a saddle node bifurcation (SN), that is $f_{SN}=0.658$. In this case, a stable and an unstable steady state become existence. Moreover, the bistability will emerge when the parameter $f$ lies the interval $f_{SN}>f>f_{c}=0.445$ (this value is more than the Hopf onset, $f_{H}=0.3397$). There are three steady states: with these kinetics A and C are linearly stable while B is unstable. Outside this interval, the system (1) has unique nontrivial equilibrium. Recent studies medvinsky:311 ; garvie:775 shown that the systems (1) can well-develop the spiral waves in the oscillation regime, but where the authors only consider the special case, i.e., $d_{p}=d_{h}$. A few important issue have not yet been properly addressed such as the spatial pattern if $d_{p}\neq d_{h}$.

Here we report the result that emergence of spatiotemporal chaos due to breakup in the system under the $d_{h}\neq d_{p}$ case. We may now use the $f$ and diffusion ratio, $\nu=d_{h}/d_{p}$, as control parameters to evaluate the region for the spiral wave. Turing instability in reaction-diffusion can be recast in terms of matrix stability Scheel2003 ; Satnoianu . Such with the help of Maple software assistance algebra computing, we obtain the parameters space $(f,\nu)$ bifurcation diagrams of the spiral waves as showing Fig. 2, in which two lines are plotted, Hopf line (solid) and Turing lines (dotted) respectively. In domain I, located above all three bifurcation lines, the homogeneous steady states is the only stable solution of the system. Domain II are regions of homogeneous oscillation in two dimensional spaces liuqx2007 . In domain III, both Hopf and Turing instabilities occur, (i.e., mixed Turing-Hopf modes arise), in which the system generally produces the phase waves. Our results show that the system has spiral wave in this regions. One can see that a Hopf bifurcation can occur at the steady when the parameter $f$ passes through a critical values $f_{H}$ while the diffusion coefficients $d_{p}=d_{h}=0$ and the bifurcation periodic solutions are stable. From our analysis (see Fig. 2), one could also see that the diffusion can induce Turing type instability for the spatial homogeneous stable periodic solutions and the spatially extended model (1) exhibit spatio-temporal chaos patterns. These spatial pattern formation arise from interaction between Hopf and Turing modes, and their subharmonics near hte codimension-two Hopf-Turing bifucation point. Special, it is interesting that spiral wave and travelling wave will appear when the parameters correspond to the Turing-Hopf bifurcation region III in the spatially extended model (1), i.e., the Turing instability and Hopf bifurcation occur simultaneously.

III Numerical results

The simulation is done in a two-dimensional (2D) Cartesian coordinate system with a grid size of $600\times 600$. The fourth order Runger-Kutta integrating method is applied with a time step $\Delta t=0.005$ time unit and a space step $\Delta x=\Delta y=0.20$ length unit. The results remain the same when the reaction-diffusion equations were solved numerically in one and two spatial dimensions using a finite-difference approximation for the spatial derivatives and an explicit Euler method for the time integration. Neumann (zero-flux) boundary conditions were emmployed in our simulation. The diffusion terms in Eqs. (1a) and (1b) often describe the spatial mixing of species due to self-motion of the organism. The typical diffusion coefficient of plankton patterns $d_{p}$ is about $0.05$, based on the parameters estimatie of Refs Jorgensen ; Okubo1980 using the relationship between turbulent diffusion and the scale of the space in the sea. In the previous studies PhysRevE.64.021915 ; medvinsky:311 ; Malchow1993 ; MercedesPascual ; garvie:775 , the authors provided a valueable insight into the role of spatial pattern for the system (1) if $d_{p}=d_{h}$. From the biological meaning, the diffusion coefficients should satisfy $d_{h}\geq d_{p}$. However, in nature waters it is turbulent diffusion that is supposed to dominate plankton mixing Sugihara , when $d_{h}<d_{p}$ is allowed. The other reason for choosing such parameter is that it is well-known new patterns, such as Turing patterns, can emerge in reaction-diffusion systems in which there is an imbalance between the diffusion coefficients $d_{p}$ and $d_{h}$ Turing1952 ; RevModPhys.65.851 . Therefore, we set $\nu=d_{h}/d_{p}$, and investigated whether a spiral wave would break up into complex spatiotemporal chaos when the diffusion ratio was varied. Throughout this paper, we fix $d_{p}=0.05$ and $d_{h}$ is a control parameter.

In the following, we will show that the dynamic behavior of the spiral wave qualitatively change as the control parameter $d_{h}$ increases from zero, i.e., the diffusion ratio $\nu$ increases from zero, to more than one. For large $\nu$ ($\nu>1$), the outwardly rotating spiral wave is completely stable everywhere, and fills in the space when the proper parameters are chosen, as shown in Fig. 3(A). Figure 3(A) shows a series of snapshots of a well-developed single spiral wave formed spontaneously for the variable $p$ in system (1). The spiral is initiated on a $600\times 600$ grid by the cross-field protocol (the initial distribution chosen in the form of allocated “constant-gradient” perturbation of the co-existence steady state) and zero boundary conditions are employed for simulations in the two dimensions. From Fig. 3(A) we can see that the well-developed spiral waves are formed firstly by the evolution. Inside the domain, new waves emerge, but are evolved by the spiral wave growing from the center. The spiral wave can steadily grow and finally prevail over the whole domain (a movie illustrating the dynamical evolution for this case movie1 [partly $movie_{-}1$, $movie_{-}2$, and $movie_{-}3$ for $d_{h}=0.2$]). Fig. 3(B) shows that the spiral wave first break up far away from the core center and eventually relatively large spiral fragments are surrounded by a ‘turbulent’ bath remain. The size of the surviving part of the spiral does not shrink when $d_{h}$ is further decreasing until finally $d_{h}$ equals to 0, which is different from phenomenon that is observed previous in the two-dimensional space Belousov-Zhabotinsky and FitzHugn-Nagumo oscillatory system Markus2004 ; PhysRevLett.82.1160 ; xie:026107 ; PhysRevE.62.7708 ; PhysRevLett.80.4811 , in which the breakup gradually invaded the stable region near the core center, and finally the spiral wave broke up in the whole medium. Figure 3(C) is the time sequences (arbitrary units) of the variables $p$ and $h$ at an arbitrary spatial point within the spiral wave region, from which we can see that the spiral waves are caused by the accepted as “phase waves” with substantially group velocity, phase velocity and sinusoidal oscillation rather than the relaxational oscillation with large amplitude. This breakup scenario is similar to the breakup of rotating spiral waves observed in numerical simulation in chemical systems Markus2004 ; PhysRevLett.82.1160 ; xie:026107 ; PhysRevE.62.7708 ; PhysRevLett.80.4811 , and experiments in BZ systems Ouyang1996 ; Ouyang2000 , which shows that spiral wave breakup in these systems was related to the Eckhaus instability and more important, the absolute instability.

The corresponding trajectories of the spiral core and the spiral arm (far away from the core center) at $y=300$ are shown in Fig. 4, respectively. From Fig. 4, we can see that the spiral core is not completely fixed, but oscillates with a large amplitude. However, as $d_{h}$ decreases to a critical value, an unstable modulation develops in regions which is far away from the spiral core (cf. the middle column of the Fig. 4). These oscillations eventually grow large enough to cause the spiral arm far away from the core to breakup into complex multiple spiral waves, while the core region remains stable (the corresponding movie can be viewed in the online supplemental in Ref. movie1 [partly $movie_{-}1$ and $movie_{-}2$, and for $d_{h}=0.02$]). Figures 3(B) and 4(B) show the dynamic behavior for $d_{h}=0.02$, i.e., $\nu=0.4$. The regular trajectories far away from the core are now the same as the region of the spatial chaos (cf. the middle column of the Fig. 4). It is shown that an decrease in the diffusion ratio $\nu$ which leads to population oscillations of increasing amplitude (cf. the left column of the Fig. 4). In the tradition explain that the minimum value of the population density decreases and population extinction becomes more probable due to stochastic environmental perturbations. However, from the spatial evolution of system (1) (see Fig. 3), the temporal variations of the density of different sub-population can become asynchronous and the events of local extinction can be compensated due to re-colonization (or diffusion) from other sites.

Furthermore, it is well known that the basic arguments in spiral stability analysis can be carried out by reducing the system to one dimensional space Markus2004 ; PhysRevLett.82.1160 ; xie:026107 ; PhysRevE.62.7708 ; PhysRevLett.80.4811 . Here we show some essential properties of the spiral breakup resulting from the numerical simulation. In the next section we will give the theoretical computation by using the eigenvalue spectra. In this model, it is worth noting that we do not neglect the oscillation of the dynamics in the core as shown in Fig. 4 due to the system exhibiting spatial periodic wave trains when the model is simulated in one-dimensional space. Breakup occurs first far away from the core (the source of waves). The spiral wave breaks towards the core until it gets to some constant distance and then the surviving part of the spiral wave stays stable. These minimal stable wavelengths are called $\lambda_{min}$. So the one-parameter family may be described by a dispersion curve $\lambda(d_{h})$ (see Fig. 5). The minimal stable wavelength $\lambda_{min}$ of the spiral wave are shown in Fig. 5 coming from the simulation in two dimensional space. The results of Fig. 5 can be interpreted as follows: the minimal stable wavelengths decrease with respect to the decrease of $d_{h}$ but eventually stay at a relative constant value, which is that the stable spiral waves are always existing for a larger region values of $d_{h}$. Space-time plots at different times are shown in Fig. 6 for two different $d_{h}$, i.e., different $\nu$, which display the time evolution of the spiral wave along the cross section in the two-dimensional images of Fig. 3(A) and (B). As shown in Fig. 6(A) and (B) for $d_{h}=0.2$ and $d_{h}=0.02$ respectively, the waves far away from the core display unstable modulated perturbation due to convective instability Markus2004 ; PhysRevLett.82.1160 ; xie:026107 ; PhysRevE.62.7708 ; PhysRevLett.80.4811 , but this perturbation is gradually advected to the left and right sides, and finally disappears. The instability manifests itself to produce the wave train breakup several waves from the far-field, as shown in Figs. 6(B).

IV Spectra and Nonlinear bifurcation of the spiral wave

In this section, we concentrate on the linear stability analysis of spiral wave by using the spectrum theory PhysRevE.62.7708 ; sandstede:016217 ; Rademacher2007 ; 2006Whell ; SandstedephysicaD . From the results in Refs. PhysRevE.62.7708 ; 2006Whell we know that the absolute spectrum must be computed numerically for any given reaction-diffusion systems. In practice, such computations only require discretization in one-dimensional space and compare with computing eigenvalues of the full stability problem on a large domain due to the spiral wave exhibitting traveling waves in the plane (see Fig. 6 about the space-time graphes). For spiral waves on the unbounded plane, the essential spectrum is also required to compute, since it determined only by the far-field wave trains of the spiral. The linear stability spectrum consists of point eigenvalues and the essential spectrum that is a continuous spectrum for spiral waves.

For sake of simplicity, the Eqs. (1a) and (1b) can been written as following

[TABLE]

Suppose that $(p^{*},h^{*})$ are a solutions and refer to them as steady spirals of Eq. (3) that rotate rigidly with a constant angular velocity $\omega$, and that are asymptotically periodic along rays in the plane. In a coratating coordinate frame, using the standardized analysis method for the spiral waves SandstedephysicaD ; 2006Whell , the Eq. (3) is given by

[TABLE]

where $(\rho,\theta)$ denote polar coordinates, spirals waves are relative equilibria, then the statianry solutions $p^{*}(\rho,\theta)$ and $h^{*}(\rho,\theta)$ both are $2\pi$-periodic functions with $\theta=\varphi-\omega t$. In Eqs. (4a) and (4b) the operator $\nabla^{2}_{\rho,\theta}$ denotes $\partial_{\rho\rho}+\frac{1}{\rho}\partial_{\rho}+\frac{1}{\rho^{2}}\partial_{\theta\theta}$.

IV.1 Computation of spiral spectra

Next, we commpute the leading part of its linear stability spectrum for the system (4). Consider the linearized evolution equation in the rotating frame, the eigenvalue problem of Eqs. (4a) and (4b) associated with the planar spiral solutions $p^{*}(\rho,\theta)$ and $h^{*}(\rho,\theta)$ are given by

[TABLE]

where $g_{1}^{p}$, $\cdots$, $g_{2}^{h}$ denote the derivatives of the nonlinear functions and $g_{1}^{p}(p,h)=r(1-p)-rp-\frac{ah}{1+bp}+\frac{abph}{(1+bp)^{2}}$, $g_{1}^{h}(p,h)=-\frac{ap}{1+bp}$, $g_{2}^{p}(p,h)=\frac{ah}{1+bp}-\frac{abph}{\left(1+bp\right)^{2}}$, and $g_{2}^{h}(p,h)=\frac{ap}{1+bp}-m-\frac{2fnh}{{n}^{2}+{h}^{2}}+\frac{2fn{h}^{3}}{\left({n}^{2}+{h}^{2}\right)^{2}}$. We shall ignore isolated eigenvalues that belong to the point spectrum, instabilities caused by point eigenvalues lead to meanderingor drifting waves, or to an unstable tip motionin in excitable media and oscillation media PhysRevE.62.7708 ; PhysRevLett.68.2090 ; PhysRevLett.72.164 ; PhysRevLett.72.2316 . This phenomenon is not shown in the present paper. Instead, we focus on the continuous spectrum that is responsible for the spiral wave breakup in the far field (see Fig. 3(b)). By the results in Ref. 2006Whell , it turns out that the boundary of the continuous spectrum depends only on the limiting equation for $\rho\rightarrow\infty$. Thus, we have that $\lambda$ is the boundary of the continuous spectrum if, and only if the limiting equation

[TABLE]

have solutions $p(\rho,\theta)$ and $h(\rho,\theta)$ for $(\rho,\theta)\in R^{+}\times[0,2\pi]$, which are bounded but does not decay as $\rho\rightarrow\infty$. Since spiral waves are rotating waves in the plane, the wave train solutions have the form as $u(t,x,y)=u(\rho,\varphi-\omega t)$ for an appropriate wave numbers $k$ and temporal frequency $\omega$, where we assume that $u$ is $2\pi$-periodic in its argument so that $u(\xi)=u(\xi+2\pi)$ for all $\xi$ and $u=(p,h)^{\text{T}}$. Spiral waves converge to wave trains $u(\rho,\varphi-\omega t)\rightarrow u_{wt}(k\rho+\varphi-\omega t)$ for $\rho\rightarrow\infty$, which are corresponding to asymptotically Archimedean in the two-dimensional space. Assume that $k\neq 0$ and $\omega\neq 0$, and in this case, we can pass from the theoretical frame $\rho$ to the comoving frame $\xi=k\rho+\varphi-\omega t$ ($\xi\in\mathcal{R}$) in which the eigenvalue equation (6) becomes

[TABLE]

Indeed, any nontrivial solution $u(\xi)=(p(\xi),h(\xi))^{\text{T}}$ corresponding to the linearization eigenvalue problem (7) give a solution $U(\rho,\cdot)$ of the eigenvalue problem for the temporal period map of (3) in the corotating frame via

[TABLE]

We write the equations (7) as the first-order systems

[TABLE]

in the radial variable $\rho$. Then the spatial eigenvalues or spatial Floquet exponents are deternined as the roots of the Wronskian

[TABLE]

where $k\in\mathcal{R}$. The function $U(\rho,\cdot)=e^{\lambda t}e^{\text{i}k\rho}u_{0}(k\rho-\omega t)$ satisfies the equation (3) when the spatial and temporal exponents $\text{i}k$ and $\lambda$ satisfy the complex dispersion relation $\text{det}(\mathcal{A}(\lambda,k)-\text{i}k)=0$ for $\lambda\in\mathcal{C}$. We call the $\text{i}k$ in spectrum of $A(\lambda,k)$ as spatial eigenvalues or spatial Floquet exponents.

The stability of the spiral waves state $(p^{*},h^{*})$ on the plane is determined by the essential spectrum given by

[TABLE]

Now, we compute the continuous spectrum with the equation (IV.1) that are parameterized by the wave number $k$. For each $\lambda$, there are infinitely many stable and unstable spatial eigenvalues. We plot $\lambda$ in the complex plane associated spatial spectrum, see Fig. 7. By the explaination of Sandstede et al sandstede:016217 , one would know that if the real part of the essentail spectra is positive, then the associated eigenmodes grow exponentially toward the boundary, i.e., they correspond to a far-field instability. Note that we find the essentail spectra are not sensitive to temporal frequency, $\omega$.

IV.2 Existence and properties of wave trains

Suppose that a reaction-diffusion system on the one-dimensional space such that the variables equal to a homogeneous stationary solution. If the homogeneous steady-state destabilizes, then its linearization accommodates waves of the form $e^{\text{i}(kx-\omega t)}$ for certain values $k$ and $\omega$. Typically, near the transition to instability, small spatially periodic travelling waves arise for any wave number close to $k_{c}$, which is the critical wavenumber. Their wave speed is approximately equal to $\frac{\omega_{c}}{k_{c}}$, where $\omega_{c}$ is corresponding to $k_{c}$. In present paper, we focus exclusively on the situation where $\omega_{c}=0$ and $k_{c}\neq 0$. The bifurcation with $\omega_{c}=0$ and $k_{c}\neq 0$ is known as the Turing bifurcation, and the bifurcating spatially periodic steady patterns are often referred to as Turing patterns. Another class of moved patterns will appear when the instabilities modulated by Hopf-Turing bifurcation, which is resemble a travelling waves. Moreover, the common feature of the spiral waves in one-dimensional space mentioned above is the presence of wave trains which are spatially periodic travelling waves of the form $p_{wt}(kx-\omega t;k)$ and $h_{wt}(kx-\omega t;k)$, where $p_{wt}(\phi;k)$ and $h_{wt}(\phi;k)$ are $2\pi$-periodic about $\phi$. Typically, the spatial wavenumber $k$ and the temporal frequency $\omega$ are related via the nonlinear dispersion relation $\omega=\omega(k)$ so that the phase velocity is given by

[TABLE]

A second quantity related to the nonlinear dispersion relation is the group velocity, $c_{g}=\frac{d\omega}{dk}$, of the wave train which also play a central role in the spiral waves. The group velocity $c_{g}$ gives the speed of propagation of small localized wave-package perturbations of the wave train sandstede:1 . Here, we are only concerned the existence of travelling wave solution. In fact, the spiral waves move at a constant speed outward from the core (see Fig. 6), so that they have the mathematical form $p(x,t)=P(z)$, and $h(x,t)=H(z)$ where $z=x-c_{p}t$. Substituting these solution forms into Eq. (3) gives the ODEs

[TABLE]

Here, we investigate numerically the existence, speed and wavelength of travelling wave patterns. Our approach is to use the bifurcation package Matcont 2.4 Matcont to study the pattern ODEs (13). To do this, the most natural bifurcation parameters are the wave speed $c_{p}$ and $f$, but they give no information about the stability of travelling wave as solutions of the model PDEs (3).

Our starting point is the homogeneous steady state of Eq. (13) with in the domain III of Fig. 2. The typical bifurcation diagrams are illustrated in Fig. 8, which shows that steady spatially peroidic travelling waves exist for the larger values of the speed $c_{p}$, but it is unstable for small values of $c_{p}$. The changes in stability occur via Hopf bifurcation, from which a branch of periodic orbits emanate. Note that here we use the terms “stable” and “unstable” as referring to the ODEs system (13) rather than the model PDEs. Fig. 8(B) illustrates the maximun stable wavelength against the bifucation parameter, speed $c_{p}$, and the small amplitudes have very long wavelength. It is known that $c_{p}=\frac{\omega}{k}$, hence the tavelling wave solution exist when the $c_{p}\neq 0$, i.e., $k\neq 0,\omega\neq 0$. Using Matcont 2.4 package, it is possible to track the locus of the Hopf bifurcation points and the Limit point (fold) bifurcation in a parameter plane, and a typical example of this for the $c_{p}$-$f$ and $c_{p}$-$d_{h}$ plane are illustrated in Fig. 9. The travelling wave solutions exist for values of $c_{p}$ and $f$ lying in left of Hopf bifurcation locus (see Fig. 9(A)). The same structure about the $c_{p}$-$d_{h}$ plane is shown in Fig. 9(B). These reuslts confirm our previous analysis coming from the algebra computation (see Fig. 2) and the numerical results (see Fig. 6).

V conclusions and discussion

We have investigated a spatially extended plankton ecological system within two-dimensional space and found that its spatial patterns exhibit spiral waves dynamics and spatial chaos patterns. Specially, the scenario of the spatiotemporal chaos patterns arising from the far-field breakup is observed. Our research is based on numerical analysis of a kinematic mimicking the diffusion in the dynamics of marine organisms, coupled to a two component plankton model on the level of the community. By increasing (decreasing) the diffusion ratio of the two variables, the spiral arm first broke up into a turbulence-like state far away from the core center, but which do not invade the whole space. From the previous studies in the Belousov-Zhabotinsky reaction, we know the reason causing this phenomenon can be illuminated theoretically by the M. Bär and L. Brusch Markus2004 ; PhysRevLett.82.1160 , as well as by using the spectrum theory that poses by B. Sandstede, A. Scheel et al PhysRevE.62.7708 ; sandstede:016217 ; Rademacher2007 ; Wheeler . The far-field breakup can be verified in field observation and is useful to understand the population dynamics of oceanic ecological systems. Such as that under certain conditions the interplay between wake (or ocean) structures and biological growth leads to plankton blooms inside mesoscale hydrodynamic vortices that act as incubators of primary production. From Fig. 3 and corresponding the movies, we see that spatial peridic bloom appear in the phytoplankton populations, and the details of spatial evolution of the distribution of the phytoplankton population during one bloom cycle, respectively.

In Ref. MaartenBoerlijst , the authors study the optimal control of the model (1) from the spatiotemporal chaos to spiral waves by the parameters for fish predation treated as a multiplicative control variable. Spatial order emerges in a range of spatial models of multispecies interactions. Unsurprisingly, spatial models of multispecies systems often manifests very different behaviors from their mean-field counterparts. Two important general features of spatial models of multispecies systems are that they allow the possibility of global persistence in spite of local extinctions and so are usually more stable than their mean-field equivalents, and have a tendency to self-organzie spatially or regular spatiotemporal patterns Dieckmann ; MaartenBoerlijst . The spatial structures produces nonrandom spatial patterns such as spiral waves and spatiotemporal chaos at scales much larger than the scale of interaction among individuals level. These structures are not explicitly coded but emerge from local interaction among individuals and local diffusion.

As we know that plankton plays an important role in the marine ecosystem and the climate, because of their participation in the global carbon and nitrogen cycle at the base of the food chain Duinker1994 . From the review Pascual2005 , a recently developed ecosystem model incorporates different phytoplankton functional groups and their competition for light and multiple nutrients. Simulations of these models at specific sites to explore future scenarios suggest that global environmental change, including global-warming-induced changes, will alter phytoplankton community structure and hence alter global biogeochemical cycles Litchman . The coupling of spatial ecosystem model to global climate raises again a series of open questions on the complexity of model and relevant spatial scales. So the study of spatial model with large-scale is more important in the ecological system. Basing on numerical simulation on the spatial model, we can draft that the oceanic ecological systems show permanent spiral waves and spatiotemporal chaos in large-scale over a range of parameter values $d_{h}$, which indicates that periodically sustained plankton blooms in the local area. As with all areas of evolutionary biology, theoretical development advances more quickly than does empiraical evidence. The most powerful empirical approach is to conduct experiments in which the spatial pattern can be measured directly, but this is difficulties in the design. However, we can indirectly measured these phenomenona by the simulation and compared with the satellite pictures. For example, the spatiotemporal chaos patterns agree with the perspective observation of the Fig. 3 in Ref. Pascual2005 . Also, some satellite imageries [http://oceancolor.gsfc.nasa.gov] have displayed spiral patterns that represent the phytoplankton [the chlorophyll] biomass and thus demonstrated that plankton patterns in the ocean occur on much broader scales and therefore mechanisms thought diffusion should be considered.

Acknowledgements.

This work is supported by the National Natural Science Foundation of China under Grant No. 10471040 and the Natural Science Foundation of Shan’xi Province Grant No. 2006011009.