A Comparison of the “Reduced Losses” and “Increased Production” Models for Mussel Bed Dynamics

Self-organised regular pattern formation is one of the foremost examples of the development of complexity in ecosystems. Despite the wide array of mechanistic models that have been proposed to understand pattern formation, there is limited general understanding of the feedback processes causing pattern formation in ecosystems, and how these affect ecosystem patterning and functioning. Here we propose a generalised model for pattern formation that integrates two types of within-patch feedback: amplification of growth and reduction of losses. Both of these mechanisms have been proposed as causing pattern formation in mussel beds in intertidal regions, where dense clusters of mussels form, separated by regions of bare sediment. We investigate how a relative change from one feedback to the other affects the stability of uniform steady states and the existence of spatial patterns. We conclude that there are important differences between the patterns generated by the two mechanisms, concerning both biomass distribution in the patterns and the resilience of the ecosystems to disturbances.


Introduction
Pattern formation at the landscape scale is an established feature of many ecosystems across the world. Examples include alternating patches of vegetation and bare ground in arid environments (Bastiaansen et al. 2018;Gandhi et al. 2019); "fairy circles" in Namibian grasslands (Zelnik et al. 2015); patterns of open-water pools in peatlands (Belyea 2007;Eppinga et al. 2009); tussock patterns in freshwater marshes (van de Koppel and Crain 2006;Yu 2010); and labyrinthine patterns in mussel beds (van de Koppel et al. 2005Liu et al. 2014). Mathematical modelling has played an important role in the study of this type of pattern formation. Most commonly, models have been used to show that a hypothesised ecological mechanism can indeed generate spatial patterns. A crucial commonality for these underlying processes is that there is a local positive feedback, where organisms improve the conditions in which they grow, but at the same time lower growth conditions at distance, often by means of competition for resources. This scale-dependent alternation between local positive feedback and larger-scale negative feedback has been proposed for many patterned ecosystems, providing a general principle for regular spatial patterns for ecology and beyond. Despite this common ground, a wide variety of specific mechanisms have been proposed as potential generators of regular patterns. Soft-bottomed mussel beds provided an example of a system in which multiple hypothesised mechanisms have been shown to be potential generators of patterns. These large aggregations of mussels form in intertidal regions, most notably in the Dutch Wadden Sea, and are often patterned, with dense clusters of mussels separated by regions of bare sediment. Two possible mechanisms for this patterning have been proposed: (i) "reduced losses"high mussel density reduces the dislodgement and predation of mussels, because of their adherence to one another; (ii) "increased production"-a high density of mussels increases their growth rate, since the faeces and pseudofaeces they produce raise them towards the algae-rich upper water layers. Mathematical modelling confirms that both of these mechanisms can generate spatial patterns (van de Koppel et al. 2005;Liu et al. 2012). How then can one determine which mechanism is the real driver of the observed pattern formation?
A first step in answering this question was taken by Liu and co-workers (Liu et al. 2012) who calculated bifurcation diagrams for models based on each of the two mechanisms, demonstrating significant differences between the patterns in the two models. However, it is unclear to what extent these differences depend on the formulation of the models and the parameter values, rather than being due to the alternative underlying patterning mechanisms. The objective of this paper is to address this issue in a comprehensive way. We will do this by developing a hybrid model that includes both the "reduced losses" and "increased production" mechanisms, with their relative importance controlled by a single tuning parameter. Our hybrid model is based on the "reduced losses" model proposed by van de Koppel et al. (2005) and the "increased production" model proposed by Liu et al. (2012).
We begin by describing the two models separately, starting with the "reduced losses" model (RLM). Denoting the mussel density by m( x, t ) and the density of algae (the main mussel food source Dolmer 2000;Øie et al. 2002) (1b) Here t is time and x is distance away from the shore; A, B, C, D, E, F, K and V are positive parameters. Since its initial formulation in 2005 (van de Koppel et al. 2005), model (1) and minor extensions have been studied in a number of papers, from both applications and mathematical viewpoints (Wang et al. 2009;Liu et al. 2012;Sherratt 2013;Ghazaryan and Manukian 2015;Cangelosi et al. 2014;Sherratt and Mackenzie 2016). The alternative "increased production" model (IPM) centres around the build-up of faeces and pseudofaeces (van Broekhoven et al. 2015) under mussel beds, which contribute to the underlying sediment. Since the mussels' food source (algae) mainly resides in upper water layer, this sediment deposition raises the mussels towards their food source and thus further promotes their growth (Liu et al. 2012(Liu et al. , 2014. In this case, we denote the mussel and algal densities by m( x, t ) and a( x, t ), respectively, and we use s ( x, t ) to represent the amount of accumulated sediment. Here t and x are time and distance away from the shore, respectively. Then, the model of Liu et al. (2012) is: where A, B, η (< 1), S 0 , C, D m , D s , E, F, P, Q and V are positive parameters. The function s + η S 0 s + S 0 has the property of increasing from a nonzero value ( η) when s = 0 towards the saturation level 1 as s → ∞; this reflects the increase in the growth rate of the mussel population with increased proximity to the (algal-rich) upper water layers. Beyond these properties, the precise functional form is arbitrary. Liu et al. (2012) calculated bifurcation diagrams for models (1) and (2) as a way of comparing their predictions for pattern formation. One of their observations was a distinct difference between the mussel densities in the patterns predicted by the two models. They reported that the average mussel density within patterns in the reduced losses model (1) is considerably greater than that in the spatially uniform steady state, but that in the increased production model (2) the two densities are similar. In this paper, we will present a more detailed study of this difference between the models with the aim of clarifying whether it is due to different model formulations and parameters, or whether it can genuinely be attributed to the different underling mechanisms. The difference is important because it provides a potential avenue for testing which model applies in a particular ecological context; mere observation of patterns is not sufficient for this since both models predict pattern formation for a wide range of ecologically plausible parameters. One issue that necessitates a detailed study is that (as we will show) there are multiple patterned solutions for any given set of ecological parameters.
As a first step, we consider a two-equation analogue of (2), which will facilitate direct comparison between the two model frameworks. Our aim is not to obtain a formal two-equation approximation to (2); rather, we seek a prototypical model of the increased production mechanism, and we use (2) as a starting point for formulating this. In this spirit, we will apply two simplifications to (2), neither of which is a good quantitative approximation, but both of which maintain the key qualitative features of the model. Firstly, we make a quasi-steady-state assumption on the kinetics of s, setting s = P Q m in (2a,b). In order for this to be a good quantitative approximation, it would be necessary that P and Q are significantly larger than comparable rate parameters in other equations. This is not actually the case: the timescale of sediment production and erosion is broadly similar to that of mussel birth and death, and for example the values estimated by Liu et al. (2012) for P and Q are the same. However, s does not play a central role in the pattern formation process, and the only effects of the quasi-steady-state assumption are on quantitative details. Secondly, to simplify the two-equation model we replace s + η S 0 s + S 0 by s S 0 . Quantitative validity of this approximation would require that S 0 s η S 0 , which holds for many of the patterns predicted by (2) using the parameter estimates in Liu et al. (2012), but not all. Note that (Liu et al. 2012) does not actually show any plots of the sediment profile in the patterned solutions of (2); we show an example in Fig. 1. In fact, the choice of functional form for the dependence of mussel birth on sediment level s in Liu et al. (2012) was somewhat arbitrary, and the prediction of pattern formation is not sensitive to this functional form. In particular, the nonzero value of s + η S 0 s + S 0 when s = 0 is not important, for the following reason. In Fig. 1, it is clear that the mussel and sediment patterns are approximately in phase; this is a typical feature of simulations, and it has a mathematical basis from linear stability analysis that we outline in "Appendix". It follows that if s is small then m will also be small, implying that the growth rate will necessarily be small.
These two simplifications give the model Numerical studies show that the qualitative features of pattern formation in (3) are broadly similar to those in (2); an example comparison is shown in Fig. 1. However, this similarity is not a key foundation of our study. Rather, our approach is to regard (3) as a prototypical model for the "increased production" mechanism, which can easily compared to the prototypical model (1) for the "reduced losses" mechanism. . Before solving, we nondimensionalised the models. Rescalings (6) are not convenient for (3), and so we used a different set of rescalings taken from Liu et al. (2012); details are given in "Appendix". This transforms (2) and (3) to ( A.1) and (A.2), respectively. We used the domain 0 <x < 50 with periodic boundary conditions, and with parameter valuesα = 50,β = 200,η = 0.1,ν = 360,Ď = 1, δ = 320,θ = 2.5 andδ = 1. (All notations are defined in "Appendix"). We began by solving (A.1), with the variables (ǎ,m,š) set to the spatially uniform steady state (ǎ s ,m s ,š s ) plus the small mixed-mode perturbation 0.01(ǎ s ,m s ,š s ) i=10 i=1 cos(π ix/50). The solutions were plotted at timesť = 600, 602, 604. We then used the solutions forǎ andm atť = 604 as initial conditions for (A.2), again solving fromť = 0 and plotting the solutions at timesť = 600, 602, 604, although we translated the solutions in the negativě x direction by 3.75 space units in order to facilitate comparison. (Since we are using periodic boundary conditions, the solutions are invariant to translation.) Our next step is to nondimensionalise the two two-equation models. Previous papers on (1) and (2) have reduced them to dimensionless forms (van de Koppel et al. 2005;Wang et al. 2009;Liu et al. 2012); however if one applies these rescalings to (1) and (3), respectively, then the resulting pairs of equations are not easily comparable. We use different rescalings that are constructed specifically to facilitate comparison between the two dimensionless models. We give the rescalings in full because our subsequent work depends crucially on the way in which the (dimensional) algal supply rate affects the dimensionless parameters in the two models. For (1), we substitute Here the asterisks denote dimensionless variables and parameters. The dimensionless parameters α * , ν * , ξ * and β * are defined by a combination of dimensional constants, but readers may find it useful to interpret intuitively α * as corresponding to the algal supply from upper water layers, ν * as corresponding to the tide strength, ξ * as corresponding to the ability of intermussel bonds to reduce dislodgement by waves, and β * as corresponding to the per capita dislodgement rate for isolated mussels. For (3), we substitute Again the asterisks denote dimensionless variables and parameters. Equations (7) are the Klausmeier model, which has been well studied as a model of semi-arid vegetation (Klausmeier 1999;Sherratt and Lord 2007;van der Stelt et al. 2013), and whose kinetics are the same as those of the Gray-Scott model, which has important applications to chemistry (Chen and Ward 2009;Doelman et al. 1997;Gray and Scott 1984). The dimensionless parameters α, ν and β are defined by a combination of dimensional constants, but readers may find it useful to interpret intuitively α as corresponding to the algal supply from upper water layers, ν as corresponding to the tide strength, and β as corresponding to the per capita dislodgement rate of mussels by waves.
The basic approach of this paper is to construct a hybrid model, which includes both the reduced losses and increased production feedback mechanisms. The model will include a new (dimensionless) parameter λ, which takes values between 0 and 1 and which controls the relative importance of the two mechanisms, with the equations  (8). We show the evolution of a pattern following a small random perturbation of a uniform state. The parameter values are α = 0.2, β = 0.1, ν = 100, ξ = 0.5, λ = 0.25 reducing to (5) when λ = 1 and to (7) when λ = 0. This approach is made possible by the close similarity between the two pairs of dimensionless equations (5) and (7), and our hybrid model is: Note that a dimensional model can be recovered from (8) by using either (4) or (6) as rescalings. Figure 2 shows a typical example of a spatially patterned solution of (8). A central question that we will address in this paper is how the mussel density in such patterns changes as the supply rate of algae is varied. In view of this, it is important to note that the algal supply parameters A in (1) and A in (2) do not appear in the rescalings for m * or m * and that the only dimensionless parameters that they affect are α * and α * . Therefore, we can address our question by considering how the (dimensionless) mussel density m in patterned solutions of (8) changes with the (dimensionless) parameter α.
Equation (10) is a cubic polynomial whose roots are very complicated algebraically, and we investigated them using numerical continuation. Starting points for this are provided by the extreme cases λ = 0 (increased production model) and λ = 1 (reduced losses model), for which the spatially uniform steady states have a much simpler form; moreover, their stability can be determined quite easily; here and throughout this section, the "stability" that we consider is to spatially uniform perturbations. For λ = 0, (α, 0) is stable for all parameter values, and if α > 2β there are two other steady states: which is unstable and The stability of (12) depends on parameters, but it is always stable if β < 2 which holds for realistic parameter estimates (Liu et al. 2012). For λ = 1, (α, 0) is stable if and only if α < β. When (α/β) lies between 1 and (1/ξ ), there is one other nonnegative steady state which is stable if ξ < 1. With these expressions for the steady states when λ = 0 and λ = 1 as starting points, we used the software package auto (Doedel 1981;Doedel et al. 1991Doedel et al. , 2006 to track the steady states and their stability as λ is varied. The results are significantly different when ξ is above and below 1, and Figs. 3 and 4 illustrate the two cases. We will discuss first the case ξ < 1 (Fig. 3) and then ξ > 1 (Fig. 4), before using our results to construct a diagram (Fig. 5) showing the existence of a positive stable-steady state in a parameter plane. However, a few general comments about the stability of the steady states are a useful preliminary. The form of (10) shows that the steady states depend on the parameters α and β only through the ratio α/β. In general, their stability does depend on the individual values of α and β, but we will show later that when β < 2, the stability also depends only on the ratio α/β. Realistic parameter estimates are consistent with β < 2 (Liu et al. 2012), and therefore, the stability shown in Figs. 3 and 4 applies under this assumption. In the figures, we plot the mussel density m at the steady states against λ, for a series of values of the parameter ratio α/β. Recall that α and β are most usefully interpreted as representing the rates of algal supply and Steady states for ξ < 1. We plot the steady-state solutions of (8) when ξ = 0.15, which is typical of behaviour when ξ < 1. We plot mussel density m against the parameter λ. Solid dashed lines indicate stable unstable steady states, calculated by numerical continuation using auto (Doedel 1981;Doedel et al. 1991Doedel et al. , 2006, with λ as continuation parameter. Dots indicate the steady states for the reduced losses (λ = 1) and increased production (λ = 0) models; these are (a, m) = (α, 0) plus (11), (12) and (13). The five different cases shown in the figure are separated by four critical values of α/β: 1, 1/λ 0 min ≈ 1.55, 2, and 1/ξ ≈ 6.67; a full description is given in the main text. To improve clarity, the plot for α/β = 8.0 is shown in two parts with different horizontal scalings. The stability of the steady states shown in the figure is valid, provided that β < 2, which applies for realistic parameter estimates (Liu et al. 2012) mussel loss, respectively, although of course they are dimensionless parameter ratios that involve a combination of ecological quantities (see (4) and (6)).
Considering first ξ < 1, five different cases are shown in Fig. 3. When α/β is small, the corresponding algal supply is insufficient to maintain a mussel population for any value of λ, and the only steady state is the mussel-free state, which is always stable. As α/β increases through 1, (13) becomes positive, and there is a stable steady state for sufficiently large values of λ. This solution branch actually has a local minimum, at λ = λ min say. The value of m giving this local minimum is negative when α/β is just above 1, but it increases with α/β and passes through zero when (10) has a double root at m = 0. A straightforward calculation shows that this occurs when α/β = 1/λ 0 min , where λ 0 min is the corresponding value of λ and is given by
The four different cases shown in the figure are separated by three critical values of α/β: 1/ξ ≈ 0.67, 1 and 2; a full description is given in the main text. The steady-state solution branches depend on α and β only through their ratio α/β; their stability does depend on the separate values, but the stability of the steady states shown in the figure is valid whenever 0 < β < 2, which applies for realistic parameter estimates (Liu et al. 2012). Figure 7 shows the corresponding bifurcation diagram for a larger value of β For the value of ξ used in Fig. 3, λ 0 min ≈ 0.646 ⇒ 1/λ 0 min ≈ 1.55. The numerical results indicate a change in stability at the local minimum, as one would expect from general bifurcation theory: λ = λ min is a saddle-node bifurcation point (Kuznetsov 2004, Ch. 3). Therefore, when (α/β) is just above 1/λ 0 min , there are two positive steady states for λ > λ min , one stable and one unstable. Figure 5a shows that λ min decreases as (α/β) increases, reaching 0 at α/β = 2. This heralds the appearance of nonzero steady states for λ = 0 (the increased production model case): see (11) and (12). The final change in qualitative behaviour occurs when α/β = 1/ξ , when the value of m at the nonzero steady state (13) for λ = 1 passes through infinity and becomes negative. For α/β > 1/ξ , the solution branch starting at (12) when λ = 0 does not intersect λ = 1; instead, m → ∞ as λ → 1 − . This is illustrated in the right-hand subpanel of Fig. 3 for α/β = 8.
The onset of patterning in (8) occurs when a spatially uniform steady state becomes unstable to inhomogeneous perturbations-this is a Turing or Turing-Hopf bifurcation. A prerequisite for this is a positive steady state that is stable to homogeneous perturbations, and the above results enable us to determine when such a steady state exists (for ξ < 1). Firstly, we note that there is never more than one positive stable  (8) has a positive steady state that is stable to homogeneous perturbations, for a ξ < 1 and b ξ > 1. The blue curve is the locus of λ min , calculated by a two-parameter numerical continuation using auto (Doedel 1981;Doedel et al. 1991Doedel et al. , 2006, with λ and α as continuation parameters. The red curve is α/β = 1/λ: the mussel-free steady state a = α, m = 0 is stable below this curve and unstable above it. The parameter values used for the plots were β = 0.1 and a ξ = 0.15, b ξ = 1.5 steady state. 1 When α/β > 2, there is a stable steady for all λ, with the exception of λ = 1 when α/β > 1/ξ . For 2 > α/β > 1/λ 0 min , the condition for a positive steady state is λ > λ min , while for 1/λ 0 min > α/β > 1 the condition is λ > (β/α), since m = 0 is a solution of (10) when α/β = 1/λ. Finally when α/β < 1, there are no positive stable steady states. Figure 5a illustrates the parameter region giving a positive stable steady state for one value of ξ between 0 and 1.
We now consider the case of ξ > 1, which is a little simpler (Fig. 4) . Again, when α/β is small, the only steady state is the mussel-free state, which is always stable. As α/β increases through 1/ξ , a branch of positive steady states appears. From the viewpoint of bifurcation theory, this solution branch arises from a transcritical bifurcation at (a, m) = (0, ∞) when λ = 1; for fixed λ < 1, the appearance of the steady states corresponds to a saddle-node bifurcation as α/β is increased. As for ξ < 1, the solution branch has a local minimum, at λ = λ min say, which separates stable and unstable parts of the solution branch. At α/β = 1, the mussel density at the nontrivial steady state (13) for λ = 1 changes sign, and for α/β > 1 the musselfree steady state is unstable for values of λ above the intersection point of the two solution branches. In the language of bifurcation theory, there is another transcritical bifurcation, this time at (a, m) = (α, 0) when λ = 1. The value of λ min decreases as α/β increases, and it reaches zero at α/β = 2, heralding the appearance of positive steady states for λ = 0. Figure 5b shows a typical plot of λ min against α/β for ξ > 1. This curve bounds the parameter region giving a positive steady state, which is shaded in the figure. Note that for α/β > 1 there is a positive steady state for an interval of λ values up to 1 (including λ = 0), but not for λ = 1 itself. Figure 5b illustrates the parameter region giving a positive stable steady state for one value of ξ greater than 1.
To conclude our discussion of steady states, we return to the issue of their stability; specifically, we will justify our assertion that provided β < 2, stability depends on the parameters α and β only through their ratio α/β. For a system of two coupled odes the stability of a steady state depends on the trace and determinant of the Jacobian (a.k.a. stability or community) matrix. For a steady state of (8) with m = 0, the determinant of the Jacobian matrix depends on α and β only through their ratio α/β, while the trace has the form T = F 1 β − F 2 . Here F 1 and F 2 are both strictly positive and depend on λ, ξ and the ratio α/β; their algebraic forms are complicated 2 . However, since F 2 > 0 it follows that for any set of values of λ, ξ and α/β, T < 0 for sufficiently small β. Then, the steady state is stable unstable when the determinant of the Jacobian matrix is positive negative, which depends on α and β only through their ratio α/β. In the latter case, the steady state is unstable for all β, but when the determinant is positive, the steady state will only be stable for β < β crit = F 2 F 1 . Thus for β > β crit there are no stable steady states with positive mussel density.
Our method of numerical continuation gives the value of m at the (unique) positive steady state with positive Jacobian determinant at any point in the (α/β)-λ plane for which this steady state exists; these are the shaded regions in Fig. 5. Using this, we calculated β crit , and the results are shown as colour maps in Fig. 6. For both ξ < 1 and ξ > 1, the minimum value of β crit is 2, occuring at λ = 0 and α/β = 2; here, the value of m = 1. Hence whenever β < 2, T < 0 and stability depends only on the ratio α/β. To illustrate the difference in steady-state stability that can occur for larger values of β, Fig. 7 shows the same bifurcation diagrams as in Fig. 3 but for β = 6. The solution branch curves are the same, but for α/β = 1.8 and 2.5 the stability is clearly different.
Our objective is to study pattern formation in (8) as λ varies between 0 and 1. Therefore, we require that there is a positive stable steady state for all values of λ, including the extreme case λ = 1 (reduced losses model). This only occurs when ξ < 1 and 2 < α/β < 1/ξ (fourth panel in Fig. 3), and we restrict attention to parameter values satisfying these constraints, and also β < 2.

Spatial Pattern Formation
When the parameters α, β, ξ and λ are such that (8) has a positive steady state that is stable to spatially uniform perturbations, there is the possibility of pattern formation in the spatial model (8), because the steady state can be destabilised by the advection and diffusion terms. This is sometimes known as "differential flow instability" (Rovinsky and Menzinger 1992). Conditions for patterning can be found via a standard procedure. We linearise equations (8) about the steady state (a s , m s ) and substitute (a − a s , m − m s ) = (a 0 , m 0 ) exp(ikx + μt). The requirement that the constants a 0 and m 0 are not both zero gives a dispersion relation, i.e. an equation for μ in terms of k. The steady The key message from these plots is that in both cases the minimum value of β crit occurs at α/β = 2 and λ = 0, at which point β crit = 2 state is stable if Re μ < 0 for all real k and unstable if Re μ > 0 for some real k. For a more detailed explanation, see Murray (2003), Perumpanani et al. (1995), Sherratt (2005), Siteur et al. (2014). Using this approach, we calculated the curves in the α-ν plane on which the steady state loses stability, leading to patterns. This is illustrated in Fig. 8, which shows that the α-ν region giving patterns shrinks as λ is increased.
The curves plotted in Fig. 8 give an upper limit on the algal supply parameter α for spatial patterning, for given values of the other parameters. Intuitively, when α exceeds this upper limit the food supply to the mussels is large enough to maintain a spatially uniform mussel population. There is also a minimum value of α below which patterns do not occur, because the food supply is insufficient to maintain even a spatially patterned mussel population; this is discussed in more detail below. Between these two threshold values of α there are spatial patterns. Although the patterns have a constant shape, they are not stationary, but rather they move in the positive x direction (away from the shore) because of the directed movement of the algae (towards the shore). Mathematically, such movement is a standard feature of patterns in reactiondiffusion-advection equations (Rovinsky and Menzinger 1992;Perumpanani et al. 1995;Malchow 2000), and it is reflected by the imaginary part of the growth rate μ being nonzero. Intuitively, the model predicts that algal density is higher on the offshore side of a mussel band compared to the on-shore side, because of consumption in the band. This causes a net growth of mussels on the off-shore side and a net loss on the on-shore side, and this causes a gradual net off-shore migration of the band (Song et al. 2017). Although such movement is intrinsic to patterned solutions of (8), it is not observed in real mussel beds, presumably because the (deliberately) simplistic modelling omits some stabilising feature(s) of the real system-one possibility for this Fig. 7 Steady states for ξ < 1 and β > 2. We plot mussel density m against the parameter λ for the steadystate solutions of (8) when ξ = 0.15 and β = 6. Solid dashed lines indicate stable unstable steady states, calculated by numerical continuation using auto (Doedel 1981;Doedel et al. 1991Doedel et al. , 2006, with λ as continuation parameter. Dots indicate the steady states for the reduced losses (λ = 1) and increased production (λ = 0) models; these are (a, m) = (α, 0) plus (11), (12) and (13). The five different cases shown in the figure are separated by four critical values of α/β: 1, 1/λ 0 min ≈ 1.55, 2, and 1/ξ ≈ 6.67; a full description is given in the main text. To improve clarity, the plot for α/β = 8.0 is shown in two parts with different horizontal magnifications. These plots should be compared with Fig. 3, which shows the corresponding plots when β < 2: the steady states are the same but there are differences in stability for α/β = 1.8 and 2.5 is the oscillatory nature of tidal flow, which has recently been incorporated into the reduced losses model (1) and which does lead to patterns without large-scale migration (Sherratt and Mackenzie 2016).
Being periodic solutions moving with constant shape and speed, the patterns are "periodic travelling waves", and general theory for this type of solution implies that for a given value of the algal supply α (and of the other parameters) there is a family of patterns, with the migration speed and wavelength varying within the family (Sherratt and Smith 2008;Sherratt 2012;van der Stelt et al. 2013;Rademacher and Scheel 2007). A convenient way of illustrating the wave family is to plot the patterning region in the α-speed plane, and some examples are shown in Fig. 9 with calculations done using the software package wavetrain (Sherratt 2012). In these plots, the patterning region is shaded; in keeping with Fig. 8 the region shrinks as λ is increased. The right-hand boundary of the patterning region is a curve on which patterns with a particular speed are initiated: mathematically this is a locus of Hopf bifurcations in the travelling wave odes. The left-hand boundary is for the most part a curve along which the pattern Fig. 8 Pattern onset. We plot the patterning onset curves for (8) in the α-ν plane for β = 0.1 and ξ = 0.5 fixed, for 4 values of the parameter λ. As α and ν are varied from lower right to upper left in the parameter plane, the locally stable positive steady state loses stability as the patterning onset curve is crossed. Note that in contrast to much of our work on steady states, the results in this and subsequent figures depend on the values of both α and β, rather than just their ratio. The curves were calculated by determining the dispersion relation as described in the main text and then solving numerically the condition for this dispersion relation to have a double root. The case λ = 0 provides a starting point for this solution, since (8) then reduces to the Klausmeier model for which patterning conditions have been studied in detail in previous work (e.g. Sherratt 2005;Siteur et al. 2014;Sherratt 2015). From this starting point, we gradually increased λ, solving the double root conditions numerically using the solution at the previous value of λ as an initial guess. Note that to enable direct comparison, the values of β and ξ used in this figure are the same as in Figs. 9,10,11,12 and 13 wavelength becomes infinitely long: mathematically this is a locus of homoclinic solutions. For one of the cases shown in Fig. 9 (λ = 0.95) a small part of the left-hand boundary consists instead of a locus of folds; this implies multiple pattern solutions in a small region of parameter space, but this has no practical significance for realistic parameters.

Mussel Density in Spatial Patterns
Our basic objective is to compare the average mussel density in spatial patterns with the density in the spatially uniform steady state from which the patterns have developed. Therefore, we tracked the form of the pattern along solution branches of (8) as α is varied with the other parameters fixed and compared this with the steady state for the same value of α. Again we used the software package wavetrain (Sherratt 2012) for this calculation. Since there is a range of different patterned solutions for any value of α within the patterning range (see Fig. 9), it is necessary to impose some We also show curves of constant migration speed and of constant wavelength (period) passing through the Turing-Hopf point. The various computations and plots were done using the software package wavetrain (Sherratt 2012). The shaded region includes patterns that are both stable and unstable as solutions of (8).
We include unstable patterns in the plots because they can be ecologically relevant in some situations (see, for example, Sherratt et al. 2009). Note that to enable direct comparison, the same parameters are used in Figs. 9, 10, 11, 12 and 13 . There is no special significance to the four values of λ that we have chosen, and in particular our use of λ = 0.05 and 0.95 rather than 0 and 1 is arbitrary form of constraint in order to specify the solution branch to be followed. A number of simulation-based studies on semi-arid vegetation patterns (Sherratt 2013;Siteur et al. 2014;Siero et al. 2015) have found that as rainfall is varied slowly, the biomass within a vegetation pattern changes, but its wavelength remains the same-except for occasional large and abrupt changes in wavelength that occur when patterns of a given wavelength lose stability. Previous studies by one of us (Sherratt 2013(Sherratt , 2016 found a corresponding result in both the reduced losses model (1) and the increased production model (2) for mussel bed patterns. Therefore, in the present work we considered patterns of constant wavelength as α was decreased; we started our calculations at the pattern onset (Turing-Hopf bifurcation) point. However, our results do not depend critically on this assumption of constant wavelength, and to emphasise this we also considered patterns of constant speed (and therefore varying wavelength) as α was decreased. Both of the resulting solution branches are shown in Fig. 9. Figures 10  and 11 show typical results from these calculations. In both the constant speed and constant wavelength cases, there is an increasing separation between the pattern and steady-state solution branches as λ is increased between 0 and 1. Recall that Liu et al. (2012) reported that the average mussel density in the spatial patterns was much greater than that in the steady state for the "reduced losses" model, while there was relatively little difference in the two densities for the "increased production" model. In our hybrid model (8), there is a gradual transition from the increased production to the reduced losses feedback mechanism as λ is increased at either constant speed or constant wavelength. Therefore, the results in Figs. 10 and 11 are consistent with the findings of Liu et al. (2012). The vertical dashed lines in Figs. 10 and 11 are at values of α that are 20% below the pattern-onset values. In Figs. 12 and 13, we show the corresponding patterns, plotting mussel density against space. The steady-state mussel density is also indicated in these figures. This shows that in contrast to the increased production mechanism, feedback of decreased losses type (λ near 1) causes high mussel densities in the peaks of the patterns; it is this that leads to the high level of mean mussel density, compared to the steady state.
To quantify the difference between the pattern and steady-state mussel densities, we calculated a single numerical measure of this difference. The starting point for our calculation is the variation with α of the mean mussel density in patterns and of the steady-state mussel density; typical examples of this are illustrated in Figs. 10 and 11. We then calculated the average over α of the difference between the two densities and divided this by the average of the steady-state density. We restricted our averaging to values of α for which the pattern and the steady state both exist, and for which the steady state is positive. In some cases, there is a fold in the pattern solution branch (for example the λ = 0.95 case in Fig. 11), and then, we restricted attention to the part of the solution branch between the fold and the pattern onset (Turing-Hopf bifurcation) point. This prevents the results from being skewed by patterns with very low densities that are (typically) unstable as solutions of (8). We repeated this procedure for a sequence of values of α, and for three values of ξ , for both the constant speed and constant wavelength solution branches. The results are shown in Fig. 14. This figure illustrates the broad generality of two conclusions. First, the (average) mussel density in the patterns is greater than that in the steady state-this follows from all of the differences plotted in Fig. 14 being positive. And secondly there is a gradual increase in the average difference in mussel density as λ is increased from 0 to 1, i.e. as the feedback mechanism gradually shifts from increased production to decreased losses.

Fig. 10
Density versus α at constant speed. We show comparisons of the average mussel density in spatial patterns with the steady-state mussel density, when the migration speed is held constant at its value at pattern onset (the Turing-Hopf bifurcation point, indicated by • ). Using the software package wavetrain (Sherratt 2012), we tracked the form of the pattern along solution branches of (8) as α is varied with the other parameters fixed. We plot the maximum, minimum and mean (average) mussel density in the patterns and also the steady-state mussel density. The figure shows an increasing separation between the pattern and the steady-state solution as λ is increased between 0 (increased production model) and 1 (decreased losses model). The dashed vertical lines indicate the values of α used in Fig.12. Note that to enable direct comparison, the same parameters are used in Figs. 9,10,11,12 and 13

Recovery Time
Previous modelling predicts that self-organisation into spatial patterns significantly increases the resilience of mussel beds to disturbances, relative to the resilience of spatially uniform mussel populations (van de Koppel et al. 2005;Liu et al. 2012). Moreover, in their study comparing models based on the reduced losses and increased production mechanisms, Liu et al. (2012) reported that this increased resilience occurs to a greater extent in the reduced losses model than in the model based on increased production. However, it is again unclear to what extent this difference depends on the formulation of the models and the parameter values, rather than being due to the alternative underlying patterning mechanisms. To clarify this, we investigated resilience in our hybrid model (8), by removing a fixed proportion of the total mussel population Fig. 11 Density versus α at constant wavelength. We show comparisons of the average mussel density in spatial patterns with the steady-state mussel density, when the wavelength is held constant at its value at pattern onset (the Turing-Hopf bifurcation point, indicated by • ). Using the software package wavetrain (Sherratt 2012), we tracked the form of the pattern along solution branches of (8) as α is varied with the other parameters fixed. We plot the maximum, minimum and mean (average) mussel density in the patterns and also the steady-state mussel density. As in Fig. 10, this figure shows an increasing separation between the pattern and the steady-state solution as λ is increased between 0 (increased production model) and 1 (decreased losses model). The dashed vertical lines indicate the values of α used in Fig. 13. Note that to enable direct comparison, the same parameters are used in Figs. 9,10,11,12 and 13 from a patterned solution, via a spatially varying perturbation, and monitoring the time taken to return to the original pattern. A typical result is shown in Fig. 15, confirming that the recovery time increases gradually as one varies the parameter λ from λ = 0 (increased production) to λ = 1 (reduced losses). This shows that the difference in resilience can indeed be attributed to the underlying pattern mechanism, rather than being a function of parameter values or other model details.

Discussion
The concept of self-organisation has been particularly successful in providing a framework for the emergence of spatial patterns in a wide range of systems. Despite its generality, a wide array of possible self-organisation models have now emerged for an  (8). We plot patterns with the same migration speed as at the pattern onset (Turing-Hopf bifurcation) point, with the value of α set at 20% below that at the pattern onset point. These values of α are indicated by dashed vertical lines in Fig.10. In each panel, we also show the steady-state mussel density (horizontal line). The patterns were calculated by numerical continuation of the travelling wave odes starting at a Hopf bifurcation point, using the software package wavetrain (Sherratt 2012). Note that to enable direct comparison, the same parameters are used in Figs. 9,10,11,12 and 13 equally large array of example ecosystems, and attempts for unification have been limited. This limits our understanding of the functional differences between self-organised ecosystems in a comprehensive way. Here we analysed a generalised model of pattern formation in ecosystems using a model that combines two general mechanisms for pattern formation: "increased production" and "reduced losses". We showed that both mechanisms predict Turing-type regular patterns to develop, but that there are clear differences between the two models.
In reality, the reduced losses and increased production mechanisms are not alternatives-both will apply to some extent, and they act in concert to generate spatial patterns. Our hybrid model (8) provides a means of including both mechanisms, with a single new parameter, and our detailed study of pattern generation provides a comprehensive framework for understanding the solutions of this model. A key takehome message of our work concerns sensitivity to variations in the balance between the reduced losses and increased production mechanisms. Figures 14 and 15 show that  (8). We plot patterns with the wavelength as at the pattern-onset (Turing-Hopf bifurcation) point, with the value of α set at 20% below that at the pattern onset point. These values of α are indicated by dashed vertical lines in Fig. 11. In each panel, we also show the steady-state mussel density (horizontal line). The patterns were calculated by numerical continuation of the travelling wave odes starting at a Hopf bifurcation point, using the software package wavetrain (Sherratt 2012). Note that to enable direct comparison, the same parameters are used in Figs. 9, 10, 11, 12 and 13 a change in this balance has a much greater effect in a system that is dominated by the reduced losses mechanism than for a system dominated by the increased production mechanism. This highlights that when the reduced losses mechanism is dominant, it is particularly important to investigate the possibility and extent of additional feedback mechanisms.
Our work has been set in the specific context of spatially patterned mussels beds. However, one or other of the "decreased losses" and "increased production" mechanisms lie at the heart of many other patterned ecosystems. Vegetation patterns are a characteristic feature of semi-arid ecosystems, consisting of alternating patches of vegetation and bare ground (Bastiaansen et al. 2018;Gandhi et al. 2019). It is widely accepted that a mechanism of increased production type plays a key role in the generation of these patterns; specifically, higher vegetation levels increase the rate at which rain water infiltrates into the soil, leading to higher plant growth rates (Klausmeier 1999;Rietkerk et al. 2002;Meron 2012). Indeed, as we have already remarked, the (a) (b) Fig. 14 Quantitative details of the comparison between the average mussel density in spatial patterns and the steady-state mussel density. In (a, b) the migration speed and wavelength, respectively, are held constant at their value at pattern onset (the Turing-Hopf bifurcation point); examples of the two corresponding solution branches are shown in Fig. 9. For each pair of λ and ξ values, we used the software package wavetrain (Sherratt 2012) to track the form of the pattern along these solution branches as α is varied with the other parameters fixed. We then calculated the average over α of the difference between the mean mussel density and the steady-state mussel density and divided this by the average of the steady-state mussel density. This gives a single number comparing the mussel density in spatial patterns and in the steady state, and we plot this as a function of λ for ξ = 0.2, 0.5 and 0.8. The plots show an increasing separation between the pattern and steady-state solution branches as λ is increased between 0 (increased production model) and 1 (decreased losses model). However, there is no clear trend in the way in which the difference in mussel densities varies with the parameter ξ . The other parameters are β = 0.1 and ν = 100  (8). A spatial pattern is generated in model (8); then, a single pattern period was extracted and a perturbed version of this was used as an initial condition on a spatial domain of length equal to one wavelength, with periodic boundary conditions. The perturbation consisted of a removal of 20% of the total biomass, varying randomly in space. The "recovery time" is a measure of the time taken for the solution to return to the patterned state; it is defined as the time by which solution is within 99% of the pre-perturbation pattern. We found that the details of the perturbation had a negligible effect on the recovery time. There is a clear trend for this recovery time to increase with the parameter λ, that is as the feedback mechanism gradually shifts from increased production to decreased losses. The parameter values are α = 0.2, β = 0.1, ν = 100 λ = 0 limit of (8) gives the widely used Klausmeier model for vegetation patterning (Klausmeier 1999;Sherratt and Lord 2007;van der Stelt et al. 2013). Another ecosystem in which a mechanism of "increased production" type drives pattern formation is the sequence of ridges and hollows found in many peat bogs (Morris et al. 2011). In fact, two different mechanisms combine to generate these patterns, with each mechanism being of "increased production" type. The elevated ridges are drier than neighbouring hollows, which causes an increased production rate of peat and an increase in the thickness of the surface acrotelm layer. This further amplifies the height difference between the ridges and hollows. In addition, the higher evapotranspiration rates by the vascular plants on the ridges lead to increased nutrient accumulation, which generates faster plant growth and thus a thicker acrotelm layer. Mathematical models based on these mechanisms confirm pattern formation (Morris et al. 2011;Eppinga et al. 2009) and have been verified by comparison with field data (Eppinga et al. 2008). An example of an ecosystem in which pattern formation is driven by a mechanism of "reduced losses" type is provided by patterns of hummocks and hollows on intertidal mudflats. Diatoms accumulate on the top of hummocks, forming a biofilm that is strengthened by extracellular secretions, and this promotes sedimentation; in the hollows, water accumulation inhibits the corresponding processes. This was first proposed as a pattern generation mechanism by Rietkerk and van de Koppel (2008) and has subsequently been tested via both modelling and field data (Weerman et al. 2010).
The occurrence of the two alternative mechanisms in a wide range of ecosystems highlights the importance of understanding the similarities and differences in the patterns that they generate. In particular, the combination of modelling studies and field work has the potential to distinguish the two mechanisms, and our results open up a number of different avenues for such future cross-disciplinary research.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. rescalings proposed by Liu et al. (2012): This transforms (2) to the system ∂ǎ/∂ť = 1 −αǎ −βǎm(š +η)/(š + 1) +ν ∂ǎ/∂x (A.1a) ∂m/∂ť =δǎm(š +η)/(š + 1) −m + ∂ 2m /∂x 2 (A.1b) ∂š/∂ť =m −θš +Ď ∂ 2š /∂x 2 . (A.1c) Applying the same rescalings to (3) gives ∂ǎ/∂ť = 1 −αǎ − β θ ǎm 2 +ν ∂ǎ/∂x (A.2a) ∂m/∂ť = δ θ ǎm 2 −m + ∂ 2m /∂x 2 . (A.2b) In Fig. 1 we compare solutions of (A.1) and (A.2) for parameter values taken from Liu et al. (2012). The uniform steady state (ǎ s ,m s ,š s ), which is used in the initial conditions for the simulations in Fig. 1, is given by One notable feature of Fig. 1 is that in the solution of the three-equation model, the patterns form andš are approximately in phase. Motivation for this small phase difference comes from linear stability analysis of the spatially uniform steady state from which patterns bifurcate. Linearising (A.1) about this steady state and substituting (ǎ −ǎ s ,m −m s ,š −š s ) = (ǎ 0 ,m 0 ,š 0 ) exp(ikx + μť) gives λ +θ +Ďk 2 š 0 =m 0 . Therefore, the phase difference between them andš patterns at the Turing-Hopf bifurcation point is arg Re λ +θ +Ďk 2 + i Im λ .
Since |Im λ/k| is the pattern migration speed, which is small, the phase difference is also small. There is no a priori reason for the phase difference not to increase, potentially significantly, as the pattern amplitude increases, but our numerical simulations suggest that this does not happen.