## Abstract

Young mussel beds on soft sediments can display large-scale regular spatial patterns. This phenomenon can be explained relatively simply by a reaction–diffusion–advection (RDA) model of the interaction between algae and mussel, which includes the diffusive spread of mussel and the advection of algae. We present a detailed analysis of pattern formation in this RDA model. We derived the conditions for differential-flow instability that cause the formation of spatial patterns, and then systematically investigated how these patterns depend on model parameters. We also present a detailed study of the patterned solutions in the full nonlinear model, using numerical bifurcation analysis of the ordinary differential equations, which were obtained from the RDA model. We show that spatial patterns occur for a wide range of algal concentrations, even when algal concentration is much lower than the prediction by linear analysis in the RDA model. That is to say, spatial patterns result from the interaction of nonlinear terms. Moreover, patterns with different wavelength, amplitude and movement speed may coexist. The results obtained are consistent with the previous observation that self-organization allows mussels to persist with algal concentrations that would not permit survival of mussels in a homogeneous bed.

## 1. Introduction

The study of pattern formation in reaction–diffusion (RD) systems is a very active area of research. Following Turing's (1952) work, many studies have focused on the mechanism behind the formation of regular patterns throughout science. Today Turing's scenario lies at the heart of many applications attempting to explain patterns in fish skin, mammalian coat markings, phyllotaxis, predator–prey systems, terrestrial vegetation, plankton, intertidal communities and so on (Hassell *et al*. 1991; Murray 1993; Rovinsky *et al*. 1997; Gurney *et al*. 1998; van de Koppel & Crain 2006; Rietkerk & van de Koppel 2008). RD models can exhibit a great variety of spatial patterns (see Cross & Hohenberg (1993), Murray (1993) and Maini *et al*. (1997) for reviews). However, their relevance to biological systems has been criticized due to their extreme sensitivity to parameter values. Murray (1993) has illustrated this by linear stability analysis and showed that the parameter regime in which a number of RD models gave rise to diffusion-driven instability is really very small. Hence, there is a clear need for a thorough mathematical analysis of case studies that reveal Turing-like patterns in biological systems.

The role of flow in pattern formation has been studied only recently. Differential flow, which is flow of species at different rates (Rovinsky & Menzinger 1992; Rovinsky *et al*. 1997; Sherratt 2005; Liu *et al*. 2008), has been shown to have a similar destabilizing effect on the homogeneous steady state in reaction–diffusion–advection (RDA) systems as differential diffusion has in the Turing case. The potential result of differential-flow instability (DIFI) is periodic travelling waves. Analysis of a generic RDA model revealed the domain in which spatial patterns occur and clarified the relation between Turing and DIFI patterns. In particular, it was shown that DIFI waves exist in parameters space that lies outside the Turing domain (Satnoianu & Menzinger 2000; Satnoianu *et al*. 2001). There is a broad potential application for DIFI mechanisms in explaining wave-like patterns, as is indicted by the occurrence of travelling wave patterns in animal growth, flow of cells, gene oscillations, early embryo development (Palmeirim *et al*. (1997) and Kaern *et al*. (2002), in which these cases correspond to the special case: the flow rate of one species is zero and the other one is not), skin patterning or gametogenesis. Moreover, DIFI patterns are very robust, with a wide domain of pattern occurrence in parameter space when compared with stationary pattern forming systems such as Turing patterns in RD systems.

Recently, many examples have emerged in the literature pointing at regular self-organized patterns in a broad range of ecosystems, such as arid vegetation, peat lands and mussel beds (Klausmeier 1999; Rietkerk *et al*. 2004; van de Koppel *et al*. 2005; Rietkerk & van de Koppel 2008). In a number of these studies, the mechanisms that were highlighted followed the RDA or DIFI framework. An example is the blue mussel *Mytilus edulis*, which is an important component of many intertidal ecosystems and a main resource for aquaculture (King *et al*. 1990; Dobretsov & Wahl 2008). Beds of young *M. edulis* in the Wadden Sea are typically characterized by banded patterns that are aligned perpendicular to the average incoming tidal flow, and have a wavelength of approximately 6 m. Within these bands, mussel can reach high densities of over 30 kg m^{−2} fresh weight, alternating with almost bare sediment in the interband areas. Van de Koppel *et al*. (2005) provided a potential explanation for the emergence of banded patterns, based on the scale difference in facilitation between neighbouring mussels and competition among the mussels for algal food sources. The scale difference was caused by the flow of the algae in the water with the tidal currents, causing mussels to compete over extensive distances.

In this paper, we provide a systematic analysis of the model of mussel stripe formation proposed by van de Koppel *et al*. (2005), extending the analysis in two directions presented in that paper. First, we compare the parameter regimes by linear analysis in which the mussel patterns occur with nonlinear analysis, focusing on their parametric sensitivity. Second, we assess the influence of the form of the nonlinearity in the kinetics on mussel pattern speed, amplitude and wavelength. We investigate the hypothesis that DIFI yields a richer class of patterns which are found in a larger parameter range than those obtained from Turing instability. Our mathematical analysis may provide a putative explanation for the widespread occurrence of banded patterns in mussel beds and other striped ecosystems.

## 2. Model

The spatio-temporal dynamics of mussel growth and survival in soft sediment mussel beds are determined to a large extent by the availability of algal food sources and loss rates due to dislodgment and predation (Cote & Isabelle 1999; van de Koppel *et al*. 2005). Dislodgment and predation strongly depend on local mussel density, because nearby conspecifics are the main substrate for attachment on soft sediments. The effects of competition are spread out over much large distances because water that is depleted in algae by the mussels is carried over the mussel bed by the tidal currents. Based on a simple mathematical model, van de Koppel *et al*. (2005) have recently hypothesized that this scale-dependent interaction between short-range facilitation and long-range competition for algae could explain the observed patterning in mussel beds.

The hypothesis of van de Koppel *et al*. (2005) has been formulated in two coupled differential equations for mussel biomass *M*(*x*, *t*) on the sediment and algae *A*(*x*, *t*) in the water layer overlying the mussel bed,(2.1a)(2.1b)where *A*_{up} describes the concentration of algae in the benthic boundary layer; *f* is the rate of exchange between the boundary layer and the upper water layers; *c* is a consumption constant; and *h* is the height of the lower water layer. Advection by tidal current is represented by the gradient operator *∂**A*/*∂**X* multiplied by the advection constant *V*. Parameter *e* is the conversion constant of ingested algae to mussel production; *d*_{M} is the maximal *per capita* mussel mortality rate; *k*_{M} is the saturation constant; and *D* is the diffusion rate of the mussels. Here, we assume that movement of mussels is random and therefore adopt the classical diffusion approximation, where diffusion is a linear function of the Laplacian operator *∂*^{2}*M*/*∂**X*^{2}. Note that here we simply assume that flow of water is the main source of transport for algae, and that diffusion (in the horizontal dimension) can be neglected. This is a simplification of course, but it shows that diffusion is not essential to explain the patterns. With parameter values as in table 1, van de Koppel *et al*. (2005) showed numerically that system (2.1*a*) and (2.1*b*) can be used to predict mussel bed stripe formation. Hence, at a general level, the model of van de Koppel *et al*. (2005) provided a potential explanation for the patterning observed in natural mussel beds, but a mathematical analysis of spatial instability and the characteristics of the developing spatial patterns is currently lacking. In this paper, we present a mathematical analysis of the conditions for development of spatial patterns, and analyse how the properties of the periodic solution depend on parameter values.

## 3. Conditions for patterns formation

To simplify, we rescale the system (2.1*a*) and (2.1*b*) as follows:(3.1)These rescalings are taken from van de Koppel *et al*. (2005), although our notation is different. The resulting dimensionless equations are(3.2a)(3.2b)In this case, the non-dimensional model has only four essential parameters that determine the dynamics of the model: renewal rate *α*; water flow rate *β*; scaled potential growth rate *δ*; and mussel maximal mortality *γ*. Therefore, in the remainder of the paper we study dimensionless equations (3.2*a*) and (3.2*b*), focusing on the conditions for patterning, and the way in which the patterns vary with these four parameters.

For all parameter values, system (3.2*a*) and (3.2*b*) has a stable trivial steady state *a*=1, *m*=0, corresponding to the existence of algae and absence of mussels. When *δ*>*γ*>*δα*, a non-trivial homogeneous steady state defined by(3.3)By the linear analysis in appendix A and straightforward manipulation of (A 5) one obtains the dispersion relation expressed as(3.4)where *j*=±1, and Re *λ* and Im *λ* denote the real part and imaginary part of *λ*, respectively. Although there is no Turing instability in system (3.2*a*) and (3.2*b*) (see appendix B), note that the dispersion relation (A 4) is similar to the relation for linear dispersion in RD systems (Satnoianu *et al*. 1998, 2001; Satnoianu & Menzinger 2000; Satnoianu 2003; Míguez *et al*. 2006), in which the Turing case is included. This indicates that we might expect a preferred direction for any possible instability mechanism resulting from equation (A 4) or directly from expression (3.4).

The condition for a spatial mode (referred as *k*) to be unstable and thus to grow into a pattern is that Re *λ*>0. Therefore, pattern formation will occur if Re *λ*>0 for any values of *k*≠0. In appendix C, we derive that the trivial state (1, 0) is linearly stable to non-uniform perturbations (Re *λ*>0). On the other hand, the non-trivial state (*a*_{s}, *m*_{s}) exhibits a stationary finite wavenumber instability also known as the Turing instability (Turing 1952; Cross & Hohenberg 1993; Murray 1993). For general advection–diffusion equations, this was first investigated by Jorne (1974), with a number of subsequent studies. In particular, previous analysis (Rovinsky & Menzinger 1992; Perumpanani *et al*. 1995; Satnoianu *et al*. 1998; Sherratt 2005) shows that spatial patterns occur if the advection rate (*β* in equations (3.2*a*) and (3.2*b*)) is above critical values, but the authors are unable to obtain an analytical expression for these critical values. So it is important to identify this parameter for the particular problem of mussel stripes. In what follows, we consider the equations (3.2*a*) and (3.2*b*) for the parameter values in table 1 with *β* and *δ* allowed to vary and take *j*=1 in equation (3.4), since we need only consider the maximum value of Re *λ*. At *β*=*β*_{c}≃13.0, there exists *λ*(*k*=*k*_{c})=0, and , where *k*_{c}=0.065 is the critical wavenumber, as is shown in figure 1, as well as in movie-growth 1 (see the electronic supplementary material and http://7y.nuc.edu.cn/jinzhen/English1.aspx). For *β*>*β*_{c}, the uniform state (*a*_{s}, *m*_{s}) becomes unstable to a continuous band of wavenumbers. Note that Re *λ*<0 at *k*=0 reflecting the purely temporal stability of this state in the non-spatial model.

From equation (3.4) we can derive the DIFI neutral curves, i.e. those curves in *β*–*k* space on which Re *λ*=0. After a simple calculation we find that these are given by(3.5)A necessary condition for equation (3.5) to hold is that *k*^{2}<*b*_{22} since *b*_{11}<0 and *Δ*>0. Generically, equation (3.5) has a strictly positive minimum *β*_{c}, and DIFI is predicted for all *β*>max{0, *β*_{c}}. The expression on the right-hand side of equation (3.5) is positive for all *k* within and the neutral curve *β*=*β*(*k*, *δ*) has a minimum at a non-zero value of *β*=*β*_{c}>0. A graph of a typical neutral curve as given by equation (3.5) is shown in figure 2. This figure shows that the neutral curve is convex with a unique minimum in the range having vertical asymptotes at *k*=0 and and 0.1582, respectively. Note that in fact *β*_{c}=*β*_{min}=*β*(*k*_{c}, *δ*) and the most unstable mode *k*_{c} is obtained by solving ∂*β*/∂*k*=0.

From the above analysis, we now deduce two important conclusions based on dispersion relation (A 5). First, Im *λ*≠0 for all (figure 3) and *β*>*β*_{c}, which leads naturally to a pattern classification that is periodic in both space and time. Hence, this is a type *I*_{0} system, following Cross & Hohenberg (1993), also referred to as spatial periodic travelling waves (see §4). Second, the emergence of the spatially periodic patterns (as observed in field given by van de Koppel *et al*. 2005) is found for a wide range of flow rates *β*, with a lower bound of *β*_{c}. This conclusion can be confirmed by using the data in table 1, which derives to the non-dimensional *β*=41.5692, which is more than the critical value of 13.0. In addition, it is instructive to note that the critical DIFI wavenumber does not equal the critical wavenumber provided by Turing theory (Murray 1993), since no stationary space-periodic DIFI waves are found in system (3.2*a*) and (3.2*b*) by stationary bifurcation analysis (see Satnoianu *et al*. (1998, 2001) and Satnoianu & Menzinger (2000) for full details).

We now define an extended DIFI bifurcation as a bifurcation to spatial-periodic solutions in a coupled RDA system from a spatially homogeneous reference state (which may be temporally stable or unstable). Near the bifurcation onset, wavenumber selection can be addressed using a weakly nonlinear analysis (Cross & Hohenberg 1993; Hoyle 2006). However, numerical integration of system (3.2*a*) and (3.2*b*) on a large domain (*L*≫2*π*/*k*_{c}≡*L*_{c}≃96.6644) with either periodic or fixed boundary conditions, slightly above *β*_{c} yields large amplitude periodic patterns whose amplitudes approach (1, 0) and (*a*_{s}, *m*_{s}). This behaviour indicates the presence of a subcritical bifurcation, i.e. the branch of periodic states should bifurcate first, in direction *β*<*β*_{c} (Kuznetsov 2004). Thus, in our case a weakly nonlinear analysis will not provide substantial insights into the nonlinear problem of pattern selection. Note that we also find that there exists a similar relation for parameter *δ*, but the difference is that the large *δ* leads again to a stable state for the homogeneous state. These can be seen from figures 4 and 9, as well as in movie-growth 2 in the electronic supplementary material.

Figure 5 illustrates the results of a numerical analysis of stability within the (*β*, *δ*) space, obtained by computing the maximal Re *λ* directly from formula (A 5). Recall that *δ* is linearly related to the algal concentration *A*_{up}, while *β* is a measure of advection constant *V* in table 1. Intuitively, figure 5 can be explained as follows: if *δ* (e.g. algal concentration) is too small, then mussels will simply die out. On the other hand, if *δ* is sufficiently large, then the competition among mussels for algae will not be very strong, resulting in homogeneous mussel. For intermediate levels of algal concentration, mussels can survive but with strong competition among mussels for algae, leading to the development of a striped mussel bed. The dependence of the lower threshold on *β* occurs because large advection will increase the competitive advantage of mussels in areas of high mussel density over those in low-density areas. So advection must be limited by a lower bound to form spatial patterns.

To understand the pattern selection mechanism, we first use the method of spatial dynamics (Champneys 1998; Yochelis *et al*. 2008) presented below, which was found to be an effective method to explore the properties of the steady-state solutions (Burke & Knobloch 2006; Yochelis *et al*. 2006; Gomila *et al*. 2007; Sherratt & Lord 2007; Knobloch 2008) of the spatial system.

## 4. Existence and properties of pattern solutions

A key issue for mussel bed stripes is to determine their wavelength and pattern speed, and the way in which these vary with model parameters. Note that some insight into the wavelength and speed can be obtained from linear analysis (Sherratt 2005). When patterns occur, one expects that these will be dominated by the fastest growing mode (*k*_{c}) although this only holds when the parameters are sufficiently close to the DIFI bifurcation curve defined by equation (3.4). However, these are properties of the nonlinear system (3.2*a*) and (3.2*b*), and a detailed investigation requires nonlinear analysis (Perumpanani *et al*. 1995; Satnoianu *et al*. 2001), which is the scope of §5.

When Re *λ*≠0 and Im *λ*≠0 in dispersion relation (3.4) (i.e. *λ* is a complex number), a Hopf bifurcation occurs which gives rise to space- and time-periodic travelling DIFI waves. From an ecological viewpoint, some field studies report such movement. Although van de Koppel *et al*. (2005) have showed that the mature mussel bands do not move much in the real world because the mussels fix themselves when they get older, young mussels move around extensively (or move around by the waves), and the wavelength of the pattern is, on average, approximately 6 m (i.e. approx. 103.9230 units in the dimensionless system (3.2*a*) and (3.2*b*)). We use numerical analysis to investigate pattern form, wavelength selection, and their dependence on model parameters in the fully nonlinear regime by using numerical analysis. The numerical analysis suggests that system (2.1*a*) and (2.1*b*) does not have any stationary patterns, with the homogeneous steady states being the only stationary solution. Instead, the patterns move at a constant speed, so that they have the mathematical form *a*(*x*, *τ*)=*u*(*z*) and *m*(*x*, *τ*)=*v*(*z*), where *z*=*x*+*cτ* and *c* is the migration speed; the patterns are periodic travelling wave solutions. Substituting these solution forms into (3.2*a*) and (3.2*b*) gives the ordinary differential equations (ODEs; Murray 1993; Ai & Huang 2005; Yochelis *et al*. 2008)(4.1)where *z* is now treated as a variable. This second-order system of ODEs predicts spatial periodic wave forms in the co-moving frame. Equation (4.1) can then be written as a system of three first-order ODEs in the standard way(4.2)Steady-state solutions of this system are of the form (*u*, *v*, *w*)=(*a*_{s}, *m*_{s}, 0), with (*a*_{s}, *m*_{s}) being the spatially uniform steady-state solutions of equations (3.2*a*) and (3.2*b*). There exist a monotone or oscillatory travelling waves solutions with respect to *z* from −∞ to +∞ in equation (4.2) when the parameters lie outside the Hopf instability, which correspond to a heteroclinic orbit that connects the trivial uniform steady-state solution (1, 0, 0) and non-trivial uniform steady-state solutions (*u*, *v*, *w*)=(*a*_{s}, *m*_{s}, 0). By contrast, the periodic travelling wave solutions will emerge when Hopf instability occurs.

Standard linear analysis of equations (4.2) reveals that a Hopf bifurcation occurs when the wave speed drops below a critical value *c*_{Hopf}. Steady-state solutions are generally stable for speeds above *c*_{Hopf} but unstable for speeds below *c*_{Hopf}. This Hopf bifurcation point is the point of origin for travelling wave families and can be calculated analytically (we omit the calculation here for brevity). Using the software package Auto (Doedel *et al*. 2001) we track the family of travelling wave solutions that arise from *c*_{Hopf} in equations (4.2) and study the changes in the wave family shape caused by varying parameter values. We use Auto to locate a particular periodic solution (a periodic travelling wave) to equations (4.2), and then track how the properties of that periodic travelling wave vary as we gradually change the equation parameters. We first use Auto to compute the eigenvalues for equations (4.2) with a given set of parameters and (*u*, *v*, *w*)=(*a*_{s}, *m*_{s}, 0), for *c* increasing through *c*_{Hopf} identifying the position of the Hopf bifurcations. We then continue along the bifurcation curve by increasing *c*.

We first located the wave speed at which a Hopf bifurcation occurs in equation (4.2), using the homogeneous steady state of system (4.2) as a starting point. When *c*=0 this steady state is stable, but provided *β* is sufficiently large (*β*>*β*_{c}), the system undergoes a Hopf bifurcation as *c* is increased. The periodic solutions that emanate from this Hopf bifurcation are the slowly moving patterns that correspond to mussel stripes, where the spatial wavelength equals the period. As *c* is increased beyond the Hopf point, the steady state is at first unstable, but regains stability through a second Hopf bifurcation. Here, we use the terms ‘stable’ and ‘unstable’ as referring to the ODEs (4.2) rather than the model partial differential equations (PDEs) (3.2*a*) and (3.2*b*). There is no direct relationship between the stability of solutions in one system and the stability in the other. Determination of periodic-travelling wave stability, even numerically, is a notoriously difficult problem (Sandstede & Scheel 1999, 2000; Sandstede 2002).

A typical bifurcation diagram and associated spatial solutions are illustrated in figure 6, which shows that for large values of *β*, changes in stability occur via Hopf bifurcations with respect to speed *c*, from which a single branch extends that links the two Hopf points. When the speed is above the critical values *c*_{c}, equilibrium stability is lost and a periodic orbit emanates. We present the bifurcation diagrams for solution amplitudes in figure 6, using the standard definition for the *L*_{2}-norm(4.3)to distinguish among solutions, where *L* and *U*_{k}(*x*) (*k*=1, 2, 3) are the spatial period and three variables of equation (4.2), respectively. Note that the standard interval [0,1] is used in the interval of integration by Auto (Doedel *et al*. 2001). Careful investigation shows that the periodic orbits are formed by bifurcation with an intermediate range of speeds; for large speed no pattern solutions occur (figure 6*a*). This is also illustrated by the plot of pattern wavelength against speed *c* in (ii) of figure 6*a*, which uses the same parameter values as (i) of figure 6*a*. The wavelength initially increases with *c* being increased, reaches a maximum, and then decreases in the bifurcation diagram.

From a practical viewpoint, one may also be concerned about the effect of other parameters such as *α*, *γ* and *δ* on the solution amplitude and wavelength. Using Auto, we track the locus of both the Hopf bifurcation points and the fixed wavelength in the two-dimensional parameter planes. Figure 7*a* presents a bifurcation analysis of the *c*–*α* parameter plane. As *c* is increased, the two Hopf bifurcation points move towards one another and eventually meet, so that the locus of Hopf points becomes a single curve. The maximum value of *c* along this curve corresponds to the maximal velocity of mussel bed patterning. Patterned solutions exist for all combinations of *c* and *α* lying in between these two loci. Using Auto, we also track the loci of patterns for particular wavelengths (figure 7, thin curves). Each curve belonging to a particular wavelength emanates from the Hopf bifurcation, and ends at the other branch (figure 7*a*). Note that for any given value of the migration speed *c* (below the fold point), large *α* is required for different pattern coexistence. In ecological terms, this suggests that the enhanced maximal rate of exchange between the lower- and upper water layers *f* also will lead to coexistence of different patterns. The corresponding loci are shown for the *c*–*γ* and *c*–*δ* planes in figure 7*b*,*c*, respectively. The basic structure is the same; in both cases, all the loci of patterns of fixed wavelength are close to the Hopf bifurcation curve. The minimum of *δ* and the maximal values of *γ* required for fixed wavelength are exactly predicted by nonlinear analysis, which correspond to the steady state (*a*_{s}, *m*_{s}, 0) becoming unstable to inhomogeneous perturbations.

Next, we track the loci of Hopf bifurcation points for particular wave speed *c* at 0.6433. Figure 8*a*,*b* presents a bifurcation analysis of the *δ*–*γ* and *δ*–*β* parameter planes, respectively. As *δ* is decreased, the two Hopf bifurcation points move towards one another and eventually meet together, at a bifurcation of higher co-dimension (seen as the cusps in figure 8*a*). Patterned solutions exist for values of *δ* and *γ* lying in between these two loci. Figure 8*b* presents the nonlinear bifurcation analysis of predicted *δ*–*β* parameter plane for the fixed wave speed. Note that the ‘U’-shape is the same as the linear prediction results for appropriate speed value due to the equilibrium independent of *β* (cf. figure 5). Hence, the limited value of *β* can be predicted by both linear and nonlinear analysis.

The above results give a detailed understanding of the existence of patterned solutions for a mussel bed, as a function of the model parameters *α*, *γ* and *δ* and the speed *c* of spatial pattern, as well as *β* (see §3). Moreover, we also find that for these equations (4.2) *c*_{Hopf} only exists when *β*>*β*_{c}. If a Hopf bifurcation exists for a chosen value of *β* then the wave family can be calculated by using an identical procedure as described above. We provide an Auto analysis given in the electronic supplementary material (mussel-auto1.tar.gz).

In the following, we concentrate on the effects of algal concentration on the spatial patterns. Figure 9 illustrates the comparison of bifurcation diagrams of pattern solution with different fixed wavelength. For a sufficiently large value of *δ*, no pattern formation occurs, because the homogeneous equilibrium is resistant to small spatially heterogeneous perturbations. As *δ* is decreased, the steady state becomes unstable via a Hopf bifurcation at *δ*≃0.188 (figure 9), which is consistent with the bifurcation analysis of system (3.2*a*) and (3.2*b*) (van de Koppel *et al*. 2005) and our prediction in §3 (figure 5). The branch of periodic solutions emanating from this Hopf bifurcation is spatially as well as temporally periodic, with a wavelength of 100. This is a mode 3 pattern. If we continue this branch, both *δ* and the solution amplitude gradually decrease until reaching a limited point, after which *δ* then increases. The speed approaches zero, while the amplitude tends to a non-zero finite limit.

Note that the upper limit of *δ* with respect to pattern formation is exactly as predicted by the linear analysis in figure 5, and corresponds to steady state (3.3) becoming unstable to inhomogeneous perturbations, owing to DIFI. However, the minimum value of *δ* where DIFI is predicted by the linear analysis (where *δ*=*γ*) exceeds the level where stable patterns are predicted by nonlinear analysis in figure 9, namely the minimal value of *δ* below which a positive homogeneous equilibrium (3.3) does not exist. Hence, spatial patterns occur in a range that well exceeds the range for which DIFI is predicted, as a result of the nonlinear nature of model (3.2*a*) and (3.2*b*). In addition, the patterns with a variety of different wavelengths, speeds and amplitudes may coexist for the intermediate *δ*, namely algal concentration (figure 9). These have been confirmed by numerical simulation with equation (3.2*a*) and (3.2*b*) (figures 10 and 11), as well as the field observation (figure 12). In our understanding, this is caused by the nonlinear interaction in model (3.2*a*) and (3.2*b*).

## 5. Numerical results

In the previous sections, we provided a detailed understanding of the effects of parameters *β*, *δ* and *c* (wave speed) on spatial patterns. Note that the wave speed is not a parameter in system (3.2*a*) and (3.2*b*), but rather is a property of the PDEs (3.2*a*) and (3.2*b*). We find that for system (3.2*a*) and (3.2*b*), the central features of a travelling wave, being wavelength and speed, as well as spatial pattern can be captured by analysing a model where space is represented in a one-dimensional model instead of in two dimensions. Therefore, we focus on the results obtained in one-dimensional space in the previous sections. In this section, we analyse the numerical solution to the PDEs (3.2*a*) and (3.2*b*) in both one- and two-dimensional space. We have taken the approach of discretizing the PDEs in space, applying upwinding to the convective term and Euler integration of the finite-difference equations. We used a rectangular spatial grid consisting of 200×200 cells (50×50 m), with a unidirectional water flow in the *y*-axis direction. We assume the simulated domain to be a section of a large mussel bed. For this reason, we adopted periodic boundary condition (in *y*-axis direction). To start off, we filled a grid of points with algal and mussel levels corresponding the homogeneous equilibrium, for a certain level of *δ* (here 0.14). Next, we perturbed this uniform solution by increasing or decreasing the mussel at some point in a random way. From those initial values, we let equation (3.2*a*) and (3.2*b*) evolve with time, where mussel and algal growth at each grid point slowly evolved towards a fixed pattern. Mussel patterns emerged, depending on the level of *δ* and *β* (figure 10), where the (i) *a*–*c* shows the two-dimensional mussel patterns and (ii) *a*–*c* shows space-temporal pattern, respectively. As we predicted in §4, our numerical results confirm these predictions by linear analysis that the mussel striped pattern occurs when the parameter is within the shaded domain (figure 10*a*). However, the spatial pattern also can survive when *δ*<*γ*, confirming our nonlinear analysis in §3 (figure 9).

Note that from these simulations, we find that long-term solutions presented in figure 10 are travelling waves. In figure 10*a*, we observe the emergence of both parallel crossed patterns, which indicate that patterns of different speed are coexisting in system (3.2*a*) and (3.2*b*). Figure 11 shows a typical spatial pattern after 50 days (approx. 7500 units in the dimensionless system (3.2*a*) and (3.2*b*)) in mussel density on the sediment for equations (2.1*a*) and (2.1*b*). There are periodic patterns of peaks in mussel density *m*, separated by regions in which mussel is almost absent. One can see from figures 10*a* and 11*a* that multiple patterns with different wavelengths, speed and amplitudes coexist for the intermediate algal concentration. In contrast to the lower algal concentration, only a single spatial pattern appears in the space (figures 10*c* and 11*c*).

## 6. Field observation of spatial patterns in mussel beds

Our previous analysis predicts that spatial patterns with different wavelengths may coexist in mussel beds. Here, we study aerial photographs of patterned mussel beds in the Wadden Sea using digitized images to estimate spatial wavelengths. We subsampled a single digital image to show that different spatial wavelengths can occur within a single mussel bed (figure 12). We tested for regular spacing in extracted panels A, B and C in figure 12*a* by applying two-dimensional spectral analysis. Spectral analysis provides a measure of the amount of variance within the image, referred to as *I*_{l,d}, that is explained by a simple cosine with a specific wavelength *l* and direction *d*. To investigate the dependence of *I* on wavelength, we plotted the radial spectrum of *I*, which reveals the dominant wavelength that fits within the analysed image, irrespective of direction. Detailed treatments of the method can be found in articles by Renshaw & Ford (1983, 1984), David (1988), Mugglestone & Renshaw (1998), Couteron & Lejeune (2001) and van de Koppel & Crain (2006).

The spectral analysis reveals that different wavelengths can be observed within a single mussel bed. Specific wavelengths are found to be 5.3, 9.2 and 17.6 m for A, B and C, respectively. Note that although both intrinsic (i.e. depletion of algae) and external processes (high differences) can cause variation in conditions that could explain variation in wavelength, it is unlikely that on the scale of this mussel bed (approx. 200 m across), abiotic conditions such as flow rates vary considerably. Hence, our spectral analysis provides evidence for the co-occurrence of wave patterns with different wavelengths within a single bed. Note that there also exist some wide ranges of different wavelengths within each panel.

## 7. Discussion

Regular pattern formation is an intriguing natural phenomenon found in a broad range of ecosystems, such as mussel beds, arid ecosystems, wetland ecosystems, savannahs, coral reefs, ribbon forests, intertidal mudflats and marsh tussocks (see Rietkerk & van de Koppel (2008) for a review). Recent studies point at an important use of spatial patterns as indicators for eminent catastrophic shifts in response to changing conditions, such as climate change (Rietkerk *et al*. 2004; Kéfi *et al*. 2007; Solé 2007). However, a systematic mathematical analysis of pattern characteristics is lacking for most ecosystems with regular patterns. We provide a mathematical analysis of the formation of spatial, stripe-like regular patterns in mussel beds occurring on soft sediment (Gascoigne *et al*. 2005; van de Koppel *et al*. 2005), and assess their potential use as indicators for the resilience of mussel beds. Using the model proposed by van de Koppel *et al*. (2005), we present the first systematic analysis of the dependence of these striped patterns on key ecological parameters. Our results show that at low mortality or high algal input rates where the resilience of mussel beds is high, many wavelengths can potentially be found. At high mortality or low algal input rates, only high wavelengths are predicted. In this region, disturbances can induce a collapse of the mussel bed, as the boundary equilibrium (*a*_{*}, 0) is stable and single wavelength dominance indicates proximity of a collapse threshold. Furthermore, our results show that for a spatial pattern of a particular wavelength, the amplitude of spatial oscillation is an increasing function of algal concentration. These results point at the potential use of pattern characteristics for assessing the resilience of mussel beds to disturbances, for instance winter storms.

Our combined results suggest that a multitude of wavelengths can coexist for a wide range of algal concentrations. These predictions are in agreement with the results of spectral analysis of aerial photographs obtained from the Wadden Sea. Although we cannot exclude the possibility that underlying gradients in abiotic conditions, such as flow rates, cause the observed variation in wavelength, hydrodynamic conditions within mussel beds are unlikely to vary considerably within the scale of the studied mussel bed and hence wavelength coexistence is a potential explanation. Our mathematical analysis points at the potential usefulness of combined measurement of variation in stripe wavelength, amplitude, at potentially stripe movement speed in assessing the resilience of mussel beds. Moreover, it indicates that integration of numerical hydrodynamic models with ecological benthic–pelagic exchange models provides a feasible approach to the development of an indicator system for assessing the resilience and vulnerability of mussel bed ecosystems.

The idea put forward by Turing that diffusive instability can lead to the formation of patterns is a simple but profound one. His work puts forward that if, in the absence of diffusion, a linearly stable homogeneous steady-state exists, spatially inhomogeneous patterns can evolve through a spatial symmetry-breaking instability. This instability is caused by differential diffusion rates of the system's components, and is therefore called *diffusion driven instability*. We analyse an alternative mechanism that causes a similar instability, arising from DIFI. While diffusion driven instability solely generates stationary spatial patterns, DIFI can generate both stationary and travelling patterns. Linear analysis predicts that spatial patterns may form only within the bifurcation curve, where patterns will be of low amplitude around the equilibrium. However, exponentially growing modes become bounded by nonlinear terms and a spatially inhomogeneous steady-state emerges. In this case, it is well known that linear analysis fails to predict the potential for large amplitude and strongly nonlinear patterns, as it does not take into account the nonlinear interaction among the variables (Satnoianu *et al*. 1998, 2001; Mei 2000). However, for large amplitude patterns, nonlinear effects will be important and may alter the parameter dependencies (see our nonlinear continuum based on package Auto in §4). In addition, these properties extend to other ecosystems that follow the same differential-flow structure, such as Klausmeier's (1999) model by Sherratt (2005) and Sherratt & Lord (2007).

There is a pressing need for suitable and accessible indicators for degradation in systems under anthropogenic stress (Solé 2007). From the pattern characteristics that we have analysed in this paper, wavelength is the most accessible indicator. Wavelength can be determined from aerial observation at a single instance using spectral or wavelet analysis, although georeferencing and hence inclusion of reference points is required for an accurate determination of wavelength. Determination of wave speed is more difficult to assess, as it requires time series of georeferenced aerial images. Amplitude can only be derived from biomass measurements, thus requiring a ground-based study. When combined with mathematical analysis of the bifurcation structure of the mussel bed system, empirical determination of wavelength, wave speed and amplitude provides a potential indication of the state, resilience and vulnerability of mussel bed ecosystems.

## Acknowledgments

We thank Norbert Dankers and Kees Kersting for providing and permitting us to use the satellite images. This work is supported by the National Natural Science Foundation of China under grant no. 60771026, the Program for New Century Excellent Talents in University (NCET050271), the Special Scientific Research Foundation for the Subjects of Doctors in University no. 20060110005 and the Natural Science Foundation of Shan'xi Province grant no. 2006011009.

## Footnotes

- Received September 27, 2008.
- Accepted October 9, 2008.

- © 2008 The Royal Society