ARTICLE IN PRESS Theoretical Population Biology 71 (2007) 367–379 www.elsevier.com/locate/tpb Local facilitation, bistability and transitions in arid ecosystems Sonia Kéfia,, Max Rietkerka, Minus van Baalenb, Michel Loreauc a Department of Environmental Sciences, Utrecht University, Copernicus Institute, P.O. Box 80115, 3508 TC, Utrecht, The Netherlands Laboratoire d’Ecologie—UMR 7625, Université Pierre et Marie Curie, Bât. A, 7ème Etage, Case 237, 7 Quai St.-Bernard, 75252 Paris Cedex 05, France. c Department of Biology, McGill University, 1205 ave Docteur Penfield, Montreal, Qué., Canada b Received 1 June 2006 Available online 29 September 2006 Abstract Arid ecosystems are liable to undergo sudden discontinuous transitions from a vegetated to a desert state as a result of human pressure and climate change. A predictive framework about the conditions under which such transitions occur is lacking. Here, we derive and analyze a general model describing the spatial dynamics of vegetation in arid ecosystems considering local facilitation as an essential process. We investigate the conditions under which continuous or discontinuous transitions from a vegetated to a desert state are likely to occur. We focus on arid ecosystems but our approach is sufficiently general to be applied to other ecosystems with severe environmental conditions. The model exhibits bistability and vegetation patchiness. High local facilitation decreases the risk of discontinuous transitions. Moreover, for arid ecosystems where local facilitation is a driving process, vegetation patchiness indicates proximity to a transition point, but does not allow distinguishing between continuous and discontinuous transitions. r 2006 Elsevier Inc. All rights reserved. Keywords: Bistability; Cellular automaton; Continuous transition; Desertification; Discontinuous transition; Pair approximation; Spatial selforganization; Vegetation patchiness 1. Introduction In the extreme conditions of arid areas, amelioration of the micro-environment underneath the canopy of a ‘nurse plant’ (Niering et al., 1963) significantly improves the conditions for establishment and growth of other plants (e.g., Mc Auliffe, 1988; Pugnaire et al., 1996; Aguiar and Sala, 1999; Holzapfel and Mahall, 1999; Armas and Pugnaire, 2005). Vegetation cover plays the role of physical protection that decreases thermal amplitudes, radiation exposures, wind dessication and soil erosion and prevents soil crust formation (Noy-Meir, 1973; Callaway and Walker, 1997; Greene et al., 2001). Moreover, water and/ or nutrient availability are higher under the canopy of certain woody species and perennial grasses as compared to areas deprived of vegetation, leading to the development of ‘islands of fertility’ (Jackson and Caldwell, 1993; Smith Corresponding author. Fax: +31302532746. E-mail addresses: kefi@geo.uu.nl (S. Kéfi), rietkerk@geo.uu.nl (M. Rietkerk), minus.van.baalen@ens.fr (M. van Baalen), michel.loreau@mcgill.ca (M. Loreau). 0040-5809/$ - see front matter r 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.tpb.2006.09.003 et al., 1994; Schlesinger et al., 1996). Here, we use the term local facilitation to describe this assortment of physical and biological mechanisms that impact the space below and close to the canopy of a plant. In this term, local stresses the spatial component of these mechanisms, and facilitation refers to their positive effects on other plants (Stachowicz, 2001; Bruno et al., 2003). Water is the main limited resource in arid ecosystems. Combined with the harsh conditions for establishment and survival of individuals, this prevents the full extension of the vegetation cover (e.g., Aguiar and Sala, 1999). In line with these mechanisms, vegetation cover in arid ecosystems is characterized by vegetation patchiness, meaning that sparsely populated, or bare, patches coexist with dense vegetation patches (Aguiar and Sala, 1999). Vegetation patches of different sizes, shapes and spatial distributions have been observed in arid ecosystems throughout the world (e.g., Aguiar and Sala, 1999; Rietkerk et al., 2000; Alados et al., 2004; Webster and Maestre, 2004). These observations have raised questions about the origin of this patchiness and the mechanisms that maintain it. ARTICLE IN PRESS 368 S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 Besides its influence on the spatial organization of vegetation, local facilitation constitutes a positive feedback between vegetation abundance and local environmental conditions that can theoretically lead to bistability (Rietkerk and van de Koppel, 1997; Scheffer et al., 2001; van de Koppel et al., 2001, 2002). Bistable systems can be in one of two alternative stable states depending on the initial conditions (Connel and Sousa, 1983). Bistable systems potentially undergo sudden, discontinuous transitions, corresponding to switches from one stable state to the other because of small, continuous changes in parameter values (Scheffer et al., 2001). For example, with decreasing rainfall, such a system may shift from a patched vegetation state to a bare state (von Hardenberg et al., 2001; Rietkerk et al., 2002). An increase in rainfall to values for which this shift occurs does not necessarily lead to a recovery of the vegetation state. This hysteresis effect is specific to discontinuous transitions and does not occur along continuous transitions, where no bistability takes place. Hence, it has been hypothesized that discontinuous transitions are associated with vegetation patchiness, and that certain shapes of patches could indicate imminent discontinuous transition towards the bare state (Rietkerk et al., 2004). Local facilitation by sessile organisms is also ubiquitous in physically harsh habitats other than arid ecosystems, such as salt marshes, alpine areas and intertidal zones (Bertness and Callaway, 1994; Bertness and Leonard, 1997; Stachowicz, 2001; Callaway et al., 2002; van de Koppel et al., 2005). The objective of this work is to gain insight into the link between bistability, vegetation patchiness and local facilitation in harsh environments. More specifically, we investigate the conditions under which local facilitation leads to continuous or discontinuous transitions under environmental and human pressure changes, and how these transitions are associated with vegetation patchiness. This is important for understanding and predicting the effect of global change on physically harsh habitats, including desertification in the case of arid ecosystems. In this paper, we focus on arid ecosystems, although our approach remains sufficiently general to be applied to other ecosystems with severe environmental conditions. Therefore, we constructed a lattice-structured model for arid ecosystems including local facilitation as an essential process. Each cell directly interacts with its nearest neighbors, which corresponds to a plant’s sphere of influence. Besides the full spatial model, two analytic approximations, a mean field (ignoring spatial interactions; Durret and Levin, 1994) and a pair approximation (Sato and Iwasa, 2000; van Baalen, 2000), leading to two systems of ordinary differential equations, were derived and analyzed for conditions leading to bistability. Comparing the results of the mean field and of the pair approximation models allows understanding the role of local interactions on the bistability of the system. Comparing the results of the pair approximation model and of the full spatial model made it possible to assess whether continuous and discontinuous transitions are associated with different vegetation patchiness. 2. The model Consider an arid ecosystem as a lattice-structured habitat, consisting of an infinite number of sites, each roughly the size of a square meter. Each site can be in one of three possible states: vegetated (denoted +), unoccupied (denoted 0) and degraded (denoted ); there is no other sub-structure within a site. The difference between unoccupied and degraded sites is that the latter cannot be colonized whereas in the former seeds can germinate. The basic task is then to formulate a model that allows us to keep track of the areas covered by vegetation and degraded soil, on the basis of the fundamental processes that occur at the site level. The fundamental processes that occur in our model are colonization (events that change {0}-sites into {+}-sites), mortality (events that change {+}-sites into {0}-sites), soil degradation (events that change {0}-sites into {}-sites) and recovery (events that change {}-sites into {0}-sites). At the site level, our model is fully stochastic. This means that each of the possible events that applies to a site in a given state does not occur with certainty but only with a given probability per unit time (i.e., a rate). These rates can be dependent not only on global vegetation cover but also on local configurations (state of the z nearest neighbors). Here we define neighborhood of a site as its four nearest neighbors (von Neumann neighborhood, z ¼ 4). Below, we give mathematical specifications for these rates. These formulae allow us to derive differential equations that describe the changes in expected cover (or mean densities) if we disregard spatial structures. These differential equations, or mean field model, can be easily analyzed to study the conditions for bistability. Not surprisingly, ignoring the spatial structure can lead to incorrect conclusions about the outcome of the population dynamics, as we show using stochastic spatial simulations of a cellular automaton. To better understand the role of spatial structure we formulate and analyze a so-called pair approximation of the cellular automaton that tracks local correlations as well as global averages. 2.1. Basic events and their associated rates Global density of sites in a given state (which is equivalent to a fraction of the total available area) is denoted by ri where i is either +, 0 or . To model colonization and facilitation we also need a way to represent local densities: following Matsuda et al. (1992) we represent these by the probabilities of finding a neighbor of a given site i in state j, in mathematical form qj|i (where i and j are either +, 0 or ). In what follows, all symbols stand for parameters and variables that have nonnegative values (Appendix A). ARTICLE IN PRESS S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 2.1.1. Colonization Plants reproduce by spreading seeds across the lattice, but the recruitment of a new individual has a probability of being successful only if the seeds have arrived at an unoccupied site. Research on seed dispersal indicates that the seed bank is mainly concentrated in the vegetation patches (Mauchamp et al., 2001; Montaña et al., 2001; van Rheede van Oudtshoorn and van Rooyen, 1999). A nonnegligible part of the seed rain, however, is dispersed further away by animals and surface runoff (van Rheede van Oudtshoorn and van Rooyen, 1999). Thus, we suppose that a fraction of the seeds produced by a vegetation site is locally dispersed, while the rest is globally dispersed. We represent the rate of colonization of an unoccupied site by w{0,+}, and we assume that
(1) wf0;þg ¼ b drþ þ ð1  dÞqþj0 GðÞ, where b is a positive constant, which represents the intrinsic seed production rate per vegetation site multiplied by the survival and the germination probabilities; d is the fraction of seeds that disperses globally and G is a function that describes how germination depends on competition for resources. The spatial scale at which competition operates in arid ecosystems remains an area of current research. Existing studies have already shown that the spatial scale of competition is larger than the one of facilitation (e.g., Lejeune et al., 1999; Rietkerk et al., 2002). For simplicity, we made competition occur at a global scale by assuming that G() is a linear function of global vegetation cover Gðrþ Þ ¼   grþ , (2) where e is the maximum germination rate and g represents the strength of density dependence. We chose goe, so that G() is always positive. For brevity, we write be as b and bg as c. Then
(3) wf0;þg ¼ drþ þ ð1  dÞqþj0 ðb  crþ Þ. 2.1.2. Mortality Just as germination turns an unoccupied site into a vegetation site, mortality turns it back into an unoccupied site again. We assume that mortality is density-independent (e.g., Klausmeier, 1999; von Hardenberg et al., 2001; Rietkerk et al., 2002; Meron et al., 2004) wfþ;0g ¼ m, (4) where m is a vegetation characteristic that may be influenced by anthropogenic factors such as grazing by cattle. 2.1.3. Degradation An unoccupied site can not only be recolonized by vegetation, it can also be eroded by rain or wind. As for mortality, we assume that degradation of the soil of an unoccupied site occurs at a density-independent rate wf0;g ¼ d; (5) 369 where d is a positive constant depending on soil type and climatic characteristics. Plant seeds cannot germinate in these sites. 2.1.4. Local facilitation Because of the positive effect of vegetation on its microenvironment, regeneration of a degraded site is assumed to be faster when there are more vegetation sites in the neighborhood. We assume regeneration of a degraded site occurs at rate w{,0} wf;0g ¼ r þ fqþj , (6) where r is a positive constant corresponding to the spontaneous regeneration rate of a site, supposed to be a function of rainfall (the more rainfall, the easier the regeneration) and soil type. Facilitation is represented by f which is the maximum effect of the neighborhood on regeneration realized when all nearest neighbors sites are {+}. Note that another facilitative term could have been included in the degradation rate, where vegetation could prevent the degradation of recolonizable sites in its neighborhood. However, in order to isolate the effect of local positive interactions, and to keep the model as simple as possible, we decided to retain only one type of facilitative effect in the model. Whenever possible we set parameters to realistic values, based on the literature (see Appendix A). For parameters for which no values could be found in the literature, we were careful to assign values with a realistic order of magnitude (Table 1, Appendix A). 2.2. Vegetation dynamics Incorporating the transition rates as defined above, we can derive the exact equations for global vegetation dynamics (changes in densities represent expected changes in cover), drþ ¼ rþ wfþ;0g þ r0 wf0;þg , dt
dr0 ¼ rþ wfþ;0g  r0 wf0;þg þ wf0;g þ r wf;0g , dt (7) (8) dr ¼ r0 wf0;g  r wf;0g . (9) dt These equations are simple bookkeeping equations based on the fact that transition rates transfer sites from one type to another. It can be shown that this set of equations represent the expected dynamics of mean densities provided that the number of sites is sufficiently large (Matsuda et al., 1992; van Baalen and Rand, 1998; van Baalen, 2000). At this stage, however, we cannot calculate these differential equations, as they depend on (at least some) local densities qj|i which depend on the spatial structures that develop. 2.2.1. Mean field analysis The so-called mean field approximation assumes that there is no spatial structure. Hence, local densities are ARTICLE IN PRESS 370 S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 identical to global densities, qi|j ¼ ri. Vegetation dynamics is then given by
drþ ¼ rþ m þ r0 drþ þ ð1  dÞrþ ðb  crþ Þ, dt (10)


dr0 ¼ rþ m  r0 drþ þ ð1  dÞrþ ðb  crþ Þ þ d þ r r þ f rþ , dt (11)
dr (12) ¼ r0 d  r  r þ f r þ . dt Using the fact that global densities sum to unity, we can represent one density in terms of the others. Therefore, only two non-linear differential equations are necessary to describe the system. For instance r0 ¼ 1r+r yields drþ ¼ rþ ðb  crþ Þð1  rþ  r Þ  mrþ , dt (13) dr ¼ dð1  rþ  r Þ  ðr þ f rþ Þr . dt (14) 2.2.2. Spatial simulations Simulations of the full spatial model (or cellular automata) were carried out using a stochastic asynchronous update algorithm (discrete time). The initial conditions were generated randomly for given initial values of r+, r0 and r. We assume that the lattice is isotropic, meaning that there is no preferential direction (no slope for example). As previously described, interactions between the cells of the cellular automaton were either with the four nearest neighbors or with the entire lattice. See figure legends for further details about the size of the lattices and the boundary conditions. 2.2.3. Pair approximation The pair approximation method makes it possible to analyze lattice models by constructing a dynamical system of ordinary differential equations (Matsuda et al., 1992; Sato and Iwasa, 2000; van Baalen, 2000) that tracks not only global densities but also local correlations. More specifically, differential equations are derived for the densities of pairs of neighbors: rij represents the probability of finding in a given pair of the lattice one neighbor in state i and the other in state j. We do not assume any asymmetry in the interaction between sites, meaning that rij ¼ rji (Sato and Iwasa, 2000). Since formally qj|i is a conditional probability (given that one neighbor is in state i, what is the probability to find the other in state j?) we can express it in terms of densities of pairs and ‘singletons,’ as rij qjji ¼ . (15) ri Making use again of conservation laws to reduce the number of equations (see Appendix B), and assuming that the dynamics of pairs does not depend on the densities of triplets or higher-order configurations (van Baalen, 2000) we can obtain our pair approximation model from drþþ ¼ 2r0þ w0;þ  2rþþ wþ;0 , dt (16) drþ ¼ r0 w0;þ þ r0þ w0;  rþ ðwþ;0 þ w;0 Þ, dt (17) dr ¼ 2r0 w0;  2r w;0 , dt (18) drþ ¼ r0 w0;þ  rþ wþ;0 , dt (19) dr (20) ¼ r0 w0;  r w;0 dt after substituting the transition rates as defined above (Appendix B). This is a set of five closed equations that can be studied using standard methods for ordinary differential equations. Note that here the transition rates are not affected by z, the number of nearest neighbors; but z intervenes in the equation of pairs (Appendix B; Sato and Iwasa, 2000). 3. Results The mean field model has one trivial equilibrium, corresponding to a system deprived of vegetation that we call desert (r+ ¼ 0, r0 ¼ r/(d+r), r ¼ d/(d+r)) which is stable if bom(d+r)/r. Non-trivial equilibria are solutions of a cubic equation (Appendix B), so that there are in total three possible outcomes for the mean field model: (1) the desert equilibrium (r+ ¼ 0) is the only globally stable equilibrium; (2) a vegetation equilibrium (r+40) is the only globally stable equilibrium; (3) the vegetation and the desert equilibria coexist and their domains of attraction are separated by an unstable equilibrium (Fig. 1). In this latter case, the system is bistable. Analysis of the mean field model suggests that the system is bistable for a certain range of parameter values (Fig. 1). However, this result stems for the assumption that there is no spatial structure. In spatially explicit simulations, vegetation patches are present over a certain range of parameter space (Fig. 2). We define a vegetation patch, or cluster of vegetation, as a group of vegetation sites connected to each other by neighborhood distances, where the neighborhood of a site is identical to that used in the model, i.e. the four nearest neighbors (Pascual et al., 2002). Starting from favorable environmental conditions (high b) and high facilitation (f), low mortality (m) or low degradation (d), the vegetation cover is homogeneous (meaning that there is no conspicuous spatial structure). When environmental conditions deteriorate, clusters of bare cells appear and remain in the vegetation cover, which then becomes more fragmented and is eventually composed of a few isolated patches of vegetation in a bare soil matrix before the system becomes a desert. High facilitation (f), low mortality rate (m) and low degradation rate (d) favor the maintenance of the ARTICLE IN PRESS S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 371 Fig. 1. Effects of the parameters on the bifurcation diagrams of the mean field model. On the ordinate, (r+)eq is the global density of vegetation sites at equilibrium. For clarity, only non-trivial equilibria are traced. A lower b value reflects more severe environmental conditions. In black: stable equilibria, in gray: unstable equilibria. The system is bistable for a certain range of environmental conditions. (a): Effect of facilitation (f). m ¼ 0.1, d ¼ 0.1, c ¼ 0.3, d ¼ 0.2, r ¼ 0.05. (b): Effect of soil type (r and d). m ¼ 0.01, c ¼ 0.3, d ¼ 0.1, f ¼ 0.3. 1: unsusceptible soil r ¼ 0.09, d ¼ 0.1, 2: r ¼ 0.03, d ¼ 0.3, 3: susceptible soil r ¼ 0.01, d ¼ 0.9. (c): Effect of mortality (m). d ¼ 0.1, r ¼ 0.01, c ¼ 0.3, d ¼ 0.2, f ¼ 0.3. (d): Effect of competition intensity (c). m ¼ 0.1, d ¼ 0.1, d ¼ 0.1, r ¼ 0.01, f ¼ 0.9. See Table 1 for the interpretation of the parameters that are not defined here. Fig. 2. Spatial organization of the vegetation predicted by the full spatial model depending on plant (f and m) and soil (d and r) characteristics. The lattices are characterized by a gradient in two parameters: facilitation (f), degradation (d) or regeneration (r) on the vertical axis, and b on the horizontal axis. A lower b value reflects more severe environmental conditions. Each cell of the lattice is submitted to a different combination of the two parameters. The simulations were run on a lattice of 40,000 (200  200) cells, with absorbing boundaries (for this type of analysis, large lattices were required so that the variation step of the parameters from one site to the next one was not too large). In black: sites occupied by vegetation, in gray: recolonizable cells, in white: degraded cells. Vegetation patches are present over a certain range of parameter space. (a) Effect of facilitation (f), and environmental conditions (b). m ¼ 0.1, r ¼ 0, d ¼ 0.2, c ¼ 0.3, d ¼ 0.1. (b) Effect of mortality (m), and environmental conditions (b). m varies from value 0.005 to 0.3, f ¼ 0.9, r ¼ 0, d ¼ 0.2, c ¼ 0.3, d ¼ 0.1. (c) Effect of the degradation rate (d) of a recolonizable site, and environmental conditions (b). m ¼ 0.1, f ¼ 0.9, r ¼ 0, c ¼ 0.3, d ¼ 0.1, d varies from 0.02 to 1. (d) Effect of the regeneration rate (r) of a degraded site, and environmental conditions (b). m ¼ 0.1, f ¼ 0.9, d ¼ 0.2, c ¼ 0.3, d ¼ 0.1. See Table 1 for the interpretation of the parameters that are not defined here. vegetation in the system. An increase in facilitation (f) leads to an increase in recolonizable sites in the neighborhood of vegetation sites. These {o}-sites can then be recolonized by seeds dispersed by their occupied neighbors, thus forming clusters of vegetation sites (Fig. 2(a)). When mortality (m) decreases, vegetation sites are more likely to recolonize {o}sites before they die (Fig. 2(b)). In a similar way, when degradation (d) decreases, {o}-sites have a higher chance of being recolonized by vegetation before being degraded (Fig. 2(c)). The regeneration rate of a degraded site that has no vegetation in its neighborhood (r) seems to hardly affect the spatial organization of the vegetation (Fig. 2(d)). Analysis of the pair approximation model provides analytical insight into the behavior of the full spatial ARTICLE IN PRESS 372 S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 model. Simulation results of the pair approximation model and the full spatial model closely resemble each other, while they both differ from the results of the mean field model (Fig. 3(a) and (b)). The ratio c++ ¼ r++/r2+ was calculated to measure the patchiness of the vegetation sites (Iwasa, 2000; van Baalen, 2000; Hui et al., 2006). In case there is no spatial correlation in the system r++ ¼ r2+ and c++ ¼ 1. If vegetation sites aggregate, the probability that two nearest neighbors are {+} is higher than in a random lattice (i.e., a lattice in which the population is uniformly distributed) and c++41. We then talk about patchiness, or clustering, of the vegetation. Simulations of both the full spatial model and the pair approximation model show that when f6¼0 or d6¼1, vegetation clustering occurs (c++41) (Fig. 3(c)). From random initial conditions, vegetation sites clump together despite the stochasticity of the full spatial model, and form vegetation patches. For stability analyses of the pair approximation model, we investigated the conditions for which an initially rare vegetation state (r+ E0) can increase in the system. This is determined by a so-called invasion condition (Appendix C; Iwasa et al., 1998). If the condition is satisfied, vegetation initially increases and is always maintained in the system at equilibrium (Iwasa et al., 1998). On the contrary, if the condition is not satisfied, the system can either be a desert at equilibrium, or be bistable. In case the state degraded (r E0) is initially rare in the system, the condition for invasion of such sites can also be derived. This condition is always satisfied (Appendix C), meaning that there are always degraded sites in the system at equilibrium. Thus, the same three outcomes as those of the mean field model are possible for the pair approximation model, i.e., a stable desert system, a stable vegetation system or a bistable system (Fig. 4). Before going on with the results, we now explain intuitively how bistability arises in the mean field as well as in the pair approximation models. We call the transition rate from a vegetation to degraded site loss rate, and the transition rate from a degraded to vegetation site recovery rate. If the recovery rate in the absence of vegetation (i.e., such that r+ E0), the basic recovery rate, is greater than the loss rate, then vegetation always invades the system and maintains itself. However, in case the loss rate exceeds the basic recovery rate, the balance between recovery and loss rates depends on both global and local densities of vegetation, and bistability may arise. Now, the system is bistable if the recovery rate can outbalance the loss rate for certain values of global and local densities of vegetation. From this general explanation, the influence of changing parameter values on the occurrence of bistability in both models can be understood. If a region of bistability is crossed when changing the bifurcation parameters, the system undergoes a discontinuous transition. A higher mortality (m) favors discontinuous transitions (Fig. 4(c)). Systems with high degradation rate (d) and low basic regeneration rate (r) are also more likely to undergo discontinuous transitions (Fig. 4(b)). The intensity of competition (c) modifies the slope of the upper branch of the bifurcation diagram (i.e., the resistance of the system, Walker et al., 2002) but does not influence the type of transition (continuous or discontinuous) (Fig. 4(d)). This is because competition intensity hardly influences the system at low vegetation densities. In the pair approximation model the effect of facilitation is reversed compared to the mean field model and facilitation favors continuous transitions (Fig. 4(a)). So, if environmental conditions worsen in a system with high local facilitation, the vegetation density gradually decreases until it reaches zero. After the loss of vegetation in the whole system, the desert system is composed of {}- and {0}-sites. Now, if environmental conditions improve again, the appearance of only one vegetation site in the system is already sufficient Fig. 3. Comparison of the full spatial model with the two approximations. Evolution of site densities and clustering intensity during time. Example of a simulation for m ¼ 0.1, c ¼ 0.3, b ¼ 0.4, d ¼ 0.2, f ¼ 0.8, r ¼ 0, d ¼ 0. For the cellular automaton, the simulation was run on a lattice of 50  50 cells with periodic boundary conditions. See Table 1 for the interpretation of the parameters that are not defined here. Simulations of the full spatial model and of the pair approximation model closely resemble each other while they both differ from the results of the mean field model. (a) Global environment. Global density of vegetated sites (r+) in gray, and global density of degraded sites (r) in black. Normal lines: pair approximation model, thick lines: cellular automaton, dashed lines: mean field model. (b) Local environment of vegetation sites. q+|+ in gray and q|+ in black. Normal lines: pair approximation model, thick lines: cellular automaton, dashed lines: mean field model. (c) Variation of the clustering intensity (c++). Grey: pair approximation model, black: cellular automaton. ARTICLE IN PRESS S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 373 Fig. 4. Effect of the parameters on the bifurcation diagram of the pair approximation model. On the ordinate, (r+)eq is the global density of vegetation sites at equilibrium. For clarity, only non-trivial equilibria are traced. A lower b value reflects more severe environmental conditions. In black: stable equilibria, in gray: unstable equilibria. Parameters have a similar effect on the bifurcation diagrams of the pair approximation model as the mean field model, except for facilitation. Higher facilitation (f) decreases the risk of discontinuous transitions in the pair approximation model. (a): Effect of facilitation (f). m ¼ 0.1, d ¼ 0.1, c ¼ 0.2, d ¼ 0.1, r ¼ 0.01. (b): Effect of soil type (r and d). m ¼ 0.2, c ¼ 0.02, d ¼ 0.1, f ¼ 1. 1: unsusceptible soil r ¼ 0.04, d ¼ 0.1, 2: r ¼ 0.02, d ¼ 0.2, 3: susceptible soil r ¼ 0.01, d ¼ 0.4. (c): Effect of mortality (m). d ¼ 0.1, r ¼ 0.01, c ¼ 0.2, d ¼ 0.1, f ¼ 0.9. The spatial organization observed at the points indicated by the arrows on (c) are displayed Fig. 6. (d): Effect of competition intensity (c): m ¼ 0.1, d ¼ 0.1, d ¼ 0.1, r ¼ 0.01, f ¼ 0.9. See Table 1 for the interpretation of the parameters that are not defined here. Bifurcation diagrams for the pair approximation model were numerically derived with the software package CONTENT (Kuznetsov and Levitin, 1997). Fig. 5. State diagrams: comparison of the mean field (a, b, c) and the pair approximation models (d, e, f). Gray: vegetation system at equilibrium, black: bistability, white: desert system at equilibrium. Decreasing facilitation (f) from left to right. A lower b value reflects more severe environmental conditions. m is the mortality rate of a vegetation site. Spatial interactions lead to a strong reduction of the bistability area in the pair approximation model as compared to the mean field model. r ¼ 0.01, c ¼ 0.2, d ¼ 0.1, d ¼ 0.1. (a) and (d): f ¼ 0.9, (b) and (e): f ¼ 0.3, (c) and (f): f ¼ 0.1. See Table 1 for the interpretation of the parameters that are not defined here. for vegetation recovery. Indeed, the probability is high that empty neighbors of this vegetation site also become occupied by vegetation. Thus, patches of vegetation can easily form and maintain themselves even from very low initial vegetation densities. Comparing the state diagrams of the mean field with the pair approximation model (Fig. 5), the desert and the vegetation domains are larger with spatial structure than without. This leads to a strong reduction of the bistability domain in the pair approximation model, regardless of the ARTICLE IN PRESS 374 S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 Fig. 6. Spatial organization of the vegetation predicted by the full spatial model along discontinuous (top) and continuous transitions (bottom) In black: cells occupied by vegetation, in gray: recolonizable cells, in white: degraded cells. Decreasing b from left to right (a lower b value reflects more severe environmental conditions). The simulations were run on lattices of 50  50 cells with periodic boundary conditions, which was appropriate for the observation of the shapes of the patches. When inspected by eye, the sequence of patch shape along discontinuous transition looks similar to the sequence of patch shape along continuous transition. r ¼ 0, c ¼ 0.2, d ¼ 0.1, d ¼ 0.1, f ¼ 0.9. Top row: Discontinuous transition, m ¼ 0.2, (a): b ¼ 0.8, (b): b ¼ 0.62, (c): b ¼ 0.57. Bottom row: Continuous transition, m ¼ 0.1, (d): b ¼ 0.5, (e): b ¼ 0.35, (f): b ¼ 0.3. See Table 1 for the interpretation of the parameters that are not defined here. intensity of local facilitation. As explained in the previous paragraph, in the area where the pair approximation predicts the maintenance of vegetation and where the mean field approximation predicts bistability, low mortalities leave enough time for patches to arise from vegetation sites because of local facilitation. At high mortalities (m) however, vegetation sites disappear before having had time to form patches: local processes, namely local facilitation (f) and local seed dispersal (1–d) become a constraint and vegetation goes always extinct in the pair approximation model whereas it can maintain itself for high initial vegetation densities in the mean field model. Vegetation patches were followed along both continuous and discontinuous transitions using the full spatial model. Vegetation patches arise along continuous and discontinuous transitions. The sequences of patch shape along continuous as compared to discontinuous transition look similar when inspected by eye (Fig. 6). 4. Discussion Our model exhibits two main properties: bistability under certain environmental conditions and vegetation clustering. These two properties arise in an isotropic environment as a result of soil characteristics (d and r), plant characteristics (m, d and b) and plant interactions (f and g). Vegetation clustering occurs when local processes are included in the model (f6¼0 or d6¼1). When environmental conditions become more severe, the vegetation cover decreases, and vegetation sites clump together, forming patches. This clustering is the result of internal dynamics and arises from random initial conditions, despite the stochasticity of the full spatial model. It may therefore be considered as a form of spatial self-organization (Perry, 1995; Rohani et al., 1997). The pair approximation model is capable of predicting the behavior of the full spatial model even when the mean field model fails. The mean field model is in particular unsuccessful in predicting the bistability area, which is always smaller in the pair approximation than in the mean field model. This illustrates the effect of local interactions on the stability of the system. In the mean field model, a higher positive feedback (f) allows vegetation to survive under harsher environmental conditions and favors discontinuous transitions. This result is in line with the mean field models results of Rietkerk and van de Koppel (1997), Lefever and Lejeune (1997), Lejeune and Tlidi (2002) and Meron et al. (2004). However, in our pair approximation model, which takes local interactions into account, a stronger facilitation (f) diminishes the risks of discontinuous transitions, because patches of vegetation can more easily form, spread and maintain themselves. According to both the mean field and the pair approximation models, systems with more susceptible soils (high d, low r) and with higher mortality (m) are more likely to be bistable. Susceptible soils correspond to clay rather than sandy soils (Walker et al., 2002). Low facilitation and high mortality characterize herbaceous rather than woody vegetation and annual rather than perennial vegetation (e.g., Pugnaire ARTICLE IN PRESS S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 et al., 1996; Holzapfel and Mahall, 1999). Thus, our model suggests that systems with dominant woody or perennial vegetation and sandy soils are less likely to be bistable than systems with dominant herbaceous or annual vegetation and clayey soils. In certain arid ecosystems, so-called regular vegetation patterns are found. Bands (‘‘tiger’’ vegetation) are observed on terrain with slopes above a critical threshold value. Spots, labyrinths or gaps in uniform vegetation cover (‘‘leopard’’ vegetation) are observed on less inclined terrain (d’Herbès et al., 2001). Note that we reserve the term pattern(s) for regular spatial organization of the vegetation (following Levin and Segel, 1985). Various mathematical models succeed in reproducing these regular patterns (Lejeune et al., 1999; von Hardenberg et al., 2001; Lejeune and Tlidi, 2002; Rietkerk et al., 2002; Meron et al., 2004). These models assume that the driving mechanism of regular vegetation pattern formation is the difference in water infiltration and/or evaporation rates between vegetation and bare sites within the system (von Hardenberg et al., 2001; Rietkerk et al., 2002; Meron et al., 2004). Resource concentration underneath the vegetation, where the infiltration rate is high and the evaporation rate is low, causes resource depletion farther away, leading to scaledependent feedback. These models also exhibit bistability of certain regular vegetation patterns and a homogeneous bare state for a certain range of parameter values (Lejeune et al., 1999; Rietkerk et al., 2002; Meron et al., 2004). These models have provided insight into how such striking regular patterns can arise in ecosystems because the patterns’ occurrence can be linked with specific mechanisms leading to scale-dependent feedbacks. In our model, facilitation is considered as a combination of physical and biological processes occurring in the local environment of a plant only, and does not trigger effects farther away. Therefore, local facilitation and competition are two independent, decoupled mechanisms, contrary to models based on scale-dependent feedback, where facilitation’s intensity and spatial scale are closely coupled with competition’s intensity and spatial scale. Local facilitation and scale-dependent feedback both include local positive feedbacks, due to different causes which are well-established empirically. However, they differ in the causes and spatial scale of competition, which has not been well investigated yet in arid ecosystems. Our general model results are in line with those of the models based on scale-dependent feedback in that we obtain bistability (under certain environmental conditions) and spatial self-organization of the vegetation (Lefever and Lejeune, 1997; Lejeune et al., 1999; von Hardenberg et al., 2001; Lejeune and Tlidi, 2002; Rietkerk et al., 2002; Meron et al., 2004). However, the spatial organization is different in our system. The patches are not regularly distributed in space, and do not have easily recognized (by eye) geometric shapes (circles, labyrinths, etc.). So, the two possible mechanisms, namely local facilitation and scale-dependant feedback, lead to two different spatial organizations of the 375 vegetation in models, and these outcomes can be connected to observations in the field. Therefore, we hypothesize that, in the field, systems with vegetation patches that are not regular and which do not have easily recognized geometric shapes are systems were local facilitation rather than scaledependent feedback operates. Conversely, systems with regular vegetation patterns could be systems that are driven by scale-dependent feedback. The link between spatial self-organization and discontinuous transitions, as hypothesized by Rietkerk et al. (2004) for scale-dependent feedback systems, does not occur in our model. In our model, both continuous and discontinuous transitions are associated with vegetation patchiness, and the two types of transitions cannot be distinguished merely by observing changes in the shape of vegetation patches. Vegetation patchiness and the type of transition are disconnected from each other in our model. Thus, in the case of systems driven by local facilitation, clusters of vegetation in a bare matrix indicate proximity to a transition to desert, but further knowledge on the stability of the system is required to assess whether such a transition will be discontinuous. On the other hand, the link between self-organized pattern formation and discontinuous transitions in case of coupled scale-dependent feedback remains an area of future research. Studies remain to be done on which facilitative mechanism occurs where and under which conditions in arid ecosystems. Our study already suggests that observation of spatial organization of vegetation in the field could give information about the type of facilitative mechanism involved. This, in turn, could provide relevant information about a possible link between such organization and the type of transition. Finally, the mechanism of local facilitation takes diverse forms in harsh environments depending on the type of stress that predominates. In subalpine and alpine plant communities, at high elevations surrounding vegetation maintains warmer temperatures, prevents from frost, and protects from strong winds and desiccation (Choler et al., 2001; Callaway et al., 2002). In salt marshes, species tolerant of extreme edaphic conditions prevent salt accumulation in the soil because their cover limits surface evaporation by shading (Bertness and Hacker, 1994). Certain species also aerate the anoxic soils by rhizosphere oxidation (Bertness and Leonard, 1997). In interdidal areas, aquatic and marine vascular plants reduce physical disturbance of the flow and enhance depositions and stabilization of the substrate (Bruno, 2000). In rock shores, shading by conspecifics renders temperature less extreme and variable and dislodgments of sessile invertebrates such as mussels by waves are reduced in dense clumps of individuals (Bertness and Leonard,1997). These ecosystems are also characterized by a certain degree of clustering. The facilitative mechanisms operating in these ecosystems could be modeled in the same way as we have done in our model. So, our model could be further generalized to study the link between local facilitation, spatial organization and type of transition in a wider range of harsh environments. ARTICLE IN PRESS 376 S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 A remaining issue is the ability to identify intrinsically bistable systems as well as systems that can become bistable through external forces. For this purpose, finding indicators of discontinuous transitions would be of general interest in ecology and relevant for ecosystem management and restoration. Acknowledgments The manuscript greatly benefited from the comments of Peter de Ruiter. We would like to thank Maarten Eppinga and Yolanda Pueyo for interesting discussions, and Ryan Chaves, Birka Wicke and Joseph Wilde-Ramsing for helping improving the clarity of the manuscript. The research of S.K. and M.R. is supported by a personal VIDI grant of the Netherlands Organization of Scientific Research/Earth and Life Sciences (NWO-ALW) to M. R. Appendix A. Estimation of the parameter values The parameters and their definitions are displayed in Table 1. The average lifespan of the vegetation of arid ecosystems is very variable (Mauchamp et al., 2001), depending on the species. Precise information on species is difficult to find in the literature. The lifespan of perennial grasses may be a few years while most of the shrubs and trees are long-lived (50–120 years). For example, a mature adult of Acacia Aneura, the tree component dominant in Australian banded landscapes (d’Herbès et al., 2001), is usually more than 100 years old and can live until 200 years. Thus the mortality rate, m, ranges between 0.005 and 1 in our model. Accelerated erosion caused by decreased vegetation cover or increased precipitation occurs at a time scale of 10–50 years. Other degradation processes such as saliniza- tion, acidification, and decomposition of organic matter occur at a shorter time scale (respectively: 0.1, 0.5–5 and 5–20 years) (Walker et al., 2002). Thus the degradation rate, d, ranges between 0.02 and 1 in our simulations. There appear to be no specific studies about seed production in arid ecosystems (Montaña et al., 2001). We have not found even an order of magnitude estimate in the literature. It is however known that, in arid ecosystems, recruitment of new individuals is a very rare event, even more when the vegetation is dominated by long-lived species. This is partly due to the great number of seeds lost in the system, either because they do not arrive at a recolonizable site, or because they leave the system with water runoff or animals. To reflect these phenomena, b (the intrinsic seed production rate per vegetation site  survival probability  germination probability) is kept between 0 and 1. The few studies which have been published on seeds in arid ecosystems all concern seed dispersal. The seed bank is patchily distributed and mainly concentrated in the vegetation patches. It has been shown that 90% of the seeds of F. Cernua, the dominant vegetation cover of bands in Mexico, are beneath the adult crown, while 10% of the seeds are dispersed by runoff water, wind, domestic and wild animals (Montaña et al., 2001). This order of magnitude is confirmed by other studies (Montaña et al., 1990). A study in northern Burkina Faso (Rasmussen et al., 2001) showed that the vegetation took 20–30 years to recover after the drought of the 1970s. In the studied area, the soil is sandy and the mean annual rainfall is around 500 mm. Therefore, this recovery time is an optimal value, rather than an average value. Another study (Walker et al., 2002) suggests indeed the century rather than the decade as a recovery time scale. Thus the recovery rate, r, ranges between 0 and 0.05 in our simulations. Table 1 Table of parameters Class Symbol Interpretation Parameter m f b e b d g c r d z Mortality rate of a vegetated cell Local facilitation, i.e., maximum effect of neighbors on the regeneration of a degraded site realized when all the nearest neighbors are {+}-sites Intrinsic seed production rate per vegetated cell  survival probability  germination probability Establishment probability of seeds on a {0}-site in a system without competition be Fraction of seeds globally dispersed Competitive effect of the global density of {+}-sites on the establishment of new individuals bg Regeneration rate of a {}-site without vegetation in its neighborhood Degradation rate of {0}-sites Number of nearest neighbors of a cell (here, z ¼ 4) r+ r r0 rsl c++ Density of vegetated sites Density of degraded sites Density of recolonizable sites Density of pairs of neighboring sites sl (s and l can be +, 0 or ) Clustering intensity of the vegetated patches Variable ARTICLE IN PRESS S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 The parameters c (competition) and f (facilitation) correspond to two different types of effects of the vegetation on the system. Their values are kept between 0, for no effect, and 1, for a maximal effect. They could not really be extracted from the literature. Intuitively, high values of c and f may however correspond to vegetation dominated by trees and shrubs while low values may represent a system dominated by perennial grasses. The parameter e is the probability of recruitment of a seed on a {0}-site in a system without competition (i.e., a system with no vegetation). This parameter reflects the ‘‘quality of the environment’’. We varied it between 0 and 1, the two possible extreme values, to see the effect of a wide range of different environmental conditions on the system. Appendix B. Systems of equations B.1. The pair approximation model The four conservation equations are: rþ þ r0 þ r ¼ 1, (21) rþþ þ rþ0 þ rþ ¼ rþ , (22) r00 þ r0þ þ r0 ¼ r0 , (23) r þ r0 þ rþ ¼ r . (24) In the system, there are nine variables in total. Three are singletons variables (rs), and six are distinct doublet variables (rss0 ), since rss0 ¼ rs0 s. Because there are also 4 conservation equations, we need five (9–4) equations to solve the system. We chose dr++/dt, dr+/dt, dr/dt, dr+/dt, and dr/dt. The principle of deriving the pair approximation equations from Eqs. (16–20) is explained in detail in Iwasa (2000) and van Baalen (2000). Using the conservation equations Eqs. (21–24) and the transition probabilities Eqs. (3–6) leads to the following closed system of equations: drþ ¼ dðrþ  rþ  rþþ Þ þ ðr  r  rþ Þ dt   ðrþ  rþ  rþþ Þ z1 ð1  dÞ  drþ þ z ð1  rþ  r Þ ðb  crþ Þ   f z  1 rþ f þm ,  rþ r þ þ z z r   dr z  1 rþ f ¼ 2dðr  r  rþ Þ  2r r þ , z r dt (27)   drþ r  rþ  rþþ ¼ drþ þ ð1  dÞ þ dt 1  rþ  r ðb  crþ Þð1  rþ  r Þ  mrþ ,   r dr ¼ dð1  rþ  r Þ  r þ f þ r , r dt ð28Þ (29) B.2. The mean field model The mean field model (Eqs. (13) and (14)) has one trivial equilibrium (r+ ¼ 0, r0 ¼ r/(d+r), r ¼ d/(d+r)) which is stable if bom(d+r)/r. Non trivial equilibria are solutions of the system:     b r m b r br    r3þ  r2þ 1 þ  þ rþ c f c c f fc br  mr  md þ ¼ 0, cf m ð30Þ r ¼ 1  rþ  b  crþ with the constraints:rþ 2 ½0; 1 and r 2 ½0; 1. Appendix C. Invasion analysis for the pair approximation model C.1. Invasion of the degraded state in a system dominated by the vegetation state Let us suppose that the lattice is in a state such that degraded sites are rare in the system (rE0). We ask whether the initially rare state ‘‘degraded’’ can increase in the system. This is determined by the sign of dr/dt. r increases if dr ¼ dð1  rþ  r Þ  ðr þ fqþj Þr 40. (31) dt This condition is always satisfied when r tends to 0, meaning that there are always some degraded sites in the lattice at equilibrium. C.2. Invasion of the vegetation state in a system dominated by the degraded state. ð25Þ drþþ ¼ 2ðrþ  rþ  rþþ Þ dt ! 1d z  1 ðrþ  rþ  rþþ Þ þ ð1  dÞ  drþ þ z z ð1  rþ  rÞ ðb  crþ Þ  2rþþ m, 377 ð26Þ In case the vegetated state is initially rare in the system, its density increases if dr+/dt40 Using the approximation r+E0 and Eq. (19), we obtain the following invasion condition: qojþ 4 m dr  . ð1  dÞb ðd þ rÞð1  dÞ (32) This condition depends on the local density q0|+. To estimate this quantity near the equilibrium, we followed the ARTICLE IN PRESS 378 S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 following procedure (Iwasa et al., 1998):   the equilibrium values of the global and local densities of the ‘‘resident’’ state (r and q|) can be calculated by using Eqs. (18), (20), the Bayes’ formula (i.e., d Eq. (15)) and by neglecting r+. We obtain: r  dþr , d qj  dþr; these results and the approximation r+E0 provide an autonomous system of differential equations of the two variables q|+ and q+|+. With x ¼ q+|+, y ¼ q|+ and w ¼ 1q|+q+|+, this system can be written as follows b r ¼ 0, 2wð1  dÞ  mx  bxwð1  dÞ  bxd z dþr drb z  1 wðd þ rÞ ðd þ ð1  dÞ z r ðd þ rÞ2   f r ¼ 0, y rþ  bywð1  dÞ  ydb z d þr dw þ  x and y can then be expressed as x¼ y¼ 2bwð1  dÞ r zðm þ bwð1  dÞ þ db dþr Þ bdr z1 wðdþrÞ wd þ ðdþrÞ 2 ðd þ ð1  dÞ z r Þ dbr r þ fz þ bwð1  dÞ þ dþr . So that using w ¼ 1xy, we can write 1¼wþ þ 2bwð1  dÞ r zðm þ bwð1  dÞ þ db dþr Þ bdr z1 wðdþrÞ wd þ ðdþrÞ 2 ðd þ ð1  dÞ z r dbr r þ fz þ bwð1  dÞ þ dþr ¼ funcðwÞ. ð33Þ func is continuous and differentiable on [0,1]. Moreover: 8wA[0,1] func’(w)40, func(0)o1 and func(1)41 Thus, there is one unique solution of Eq. (33) between 0 and 1 (Intermediate values theorem). The invasion condition is thus: m dr w4 ð1dÞb  ðdþrÞð1dÞ (using Eq. (32) and the fact that w ¼ 1q|+q+|+ ¼ qo|+)   m dr  3f ðwÞ4f ð1  dÞb ðd þ rÞð1  dÞ   m dr  2b b dþr m 3 þ 2zm ð1  dÞb     d m dr d drb 1d b  dþr þ zðdþrÞ ðz  1Þm þ dþr þ r dr f  þ þ mo1. ð34Þ ð1  dÞðd þ rÞ z When the invasion condition (34) is satisfied, vegetation is always maintained in the system at equilibrium. On the contrary, when the invasion condition (34) is not satisfied, the system is either a desert at equilibrium or bistable. Numerical simulations can then be done to decide between desert or bistable system. References Aguiar, M.R., Sala, O., 1999. Patch structure, dynamics and implications for the functioning of arid ecosystems. Trends Ecol. Evol. 14 (7), 273–277. Alados, C.L., ElAich, A., Papanastasis, V.P., Ozbek, H., Navarro, T., Freitas, H., Vrahnakis, M., Larrosi, D., Cabezudo, B., 2004. Ecol. Model 180, 523–535. Armas, C., Pugnaire, F.I., 2005. Plant interactions govern population dynamics in a semi-arid plant community. J. Ecol. 93, 978–989. Bertness, M.D., Callaway, R., 1994. Positive interactions in communities. Trends Ecol. Evol. 9 (5), 191–193. Bertness, M.D., Hacker, S.D., 1994. Physical stress and positive associations among marsh plants. Am. Nat. 144 (3), 363–372. Bertness, M.D., Leonard, G.H., 1997. The role of positive interactions in communities: lessons from intertidal habitats. Ecology 78 (7), 1976–1989. Bruno, J.F., 2000. Facilitation of cobble beach plant communities through habitat modification by Spartina Alterniflora. Ecology 81, 1179–1192. Bruno, J.F., Stachowicz, J.J., Bertness, M., 2003. Inclusion of facilitation into ecological theory. Trends Ecol. Evol. 18, 119–125. Callaway, R.M., Walker, L.R., 1997. Competition and facilitation: a synthetic approach to interactions in plant communities. Ecology 78 (7), 1958–1965. Callaway, R.M., Brooker, R.W., Choler, P., Kikvidze, Z., Lortie, C.J., Michalet, R., et al., 2002. Positive interactions among alpine plants increase with stress. Nature 417, 844–847. Choler, P., Michalet, R., Callaway, R.M., 2001. Facilitation and competition on gradients in alpine plants communities. Ecology 82 (12), 3295–3308. Connel, J.H., Sousa, W.P., 1983. On the evidence needed to judge ecological stability or persistence. Am. Nat. 121, 789–824. d’Herbès, J.-M., Valentin, C., Tongway, D.J., Leprun, J.C., 2001. Banded vegetation patterns and related structures. In: Tongway, D.J., Valentin, C., Seghieri, J. (Eds.), Banded vegetation patterning in arid and semiarid environments. Ecological processes and consequences for management. Ecological Studies, vol. 149. Springer, Berlin, pp. 1–19. Durret, R., Levin, S., 1994. The importance of being discrete (and spatial). Theor. Popul. Biol. 46, 363–394. Greene, R.S.B., Valentin, C., Esteves, M., 2001. Runoff and erosion processes. In: Tongway, D.J., Valentin, C., Seghieri, J. (Eds.), Banded vegetation patterning in arid and semiarid environments. Ecological processes and consequences for management. Ecological Studies, vol. 149. Springer, Berlin, pp. 52–76. Holzapfel, C., Mahall, B.E., 1999. Bidirectional facilitation and interference between shrubs and annuals in the majove desert. Ecology 80 (5), 1966–1975. Hui, C., McGeoch, M.A., Warren, M., 2006. A spatially explicit approach to estimate species occupancy and spatial correlations. J. Anim. Ecol. 75, 140–147. Iwasa, Y., Nakamaru, M., Levin, S.A., 1998. Allelopathy of bacteria in a lattice population: competition between colicin-sensitive and colicinproducing strains. Evol. Ecol. 12, 785–802. Iwasa, Y., 2000. Lattice models and pair approximation in ecology. In: Diekmann, U., Law, R., Metz, J.A.J. (Eds.), The Geometry of Ecological Interactions. Simplifying Spatial Complexity. Cambridge University Press, Cambridge, pp. 227–251. Jackson, R.B., Caldwell, M.M., 1993. The scale of nutrient heterogeneity around individual plants and its quantification with geostatistics. Ecology 74 (2), 612–614. Klausmeier, C.A., 1999. Regular and irregular patterns in semiarid vegetation. Science 284, 1826–1828. ARTICLE IN PRESS S. Kéfi et al. / Theoretical Population Biology 71 (2007) 367–379 Kuznetsov, Y.A., Levitin, V.V., 1997. CONTENT: a multiplatform environment for continuation and bifurcation analysis of dynamical systems. Centrum voor Wiskunde en Informatica, Amsterdam. Lefever, R., Lejeune, O., 1997. On the origin of tiger bush. B. Math. Bio. Sci. 59 (2), 263–294. Lejeune, O., Couteron, P., Lefever, R., 1999. Short range co-operativity competing with long range inhibition explains vegetation patterns. Acta Oecol. 20 (3), 171–183. Lejeune, O., Tlidi, M., 2002. Localized vegetation patches: a selforganized response to resource scarcity. Phys. Rev. E. 66, 010901. Levin, S.A., Segel, L.A., 1985. Pattern generation in space and aspect. SIAM Rev. 27 (1), 45–67. Matsuda, H., Ogita, N., Sasaki, A., Sato, K., 1992. Statistical mechanics of population.The lattice Lotka–Volterra Model. Prog. Theor. Phys. 88 (6), 1035–1049. Mauchamp, A., Rambal, S., Ludwig, J.A., Tongway, D.J., 2001. Multiscale modelling of vegetation bands. In: Tongway, D.J., Valentin, C., Seghieri, J. (Eds.), Banded vegetation patterning in arid and semiarid environments. Ecological processes and consequences for management. Ecological Studies, vol. 149. Springer, Berlin, pp. 146–165. Mc Auliffe, J.R., 1988. Markovian dynamics of simple and complex desert plant communities. Am. Nat. 131 (4), 459–490. Meron, E., Gilad, E., von Hardenberg, J., Shachak, M., Zoumi, Y., 2004. Vegetation patterns along a rainfall gradient. Chaos Solitons Fract. 19, 367–376. Montaña, C., Lòpez-Portillo, J., Monchamp, A., 1990. The resource of two woody species to the conditions created by a shifting ecotone in an arid environment. J. Ecol. 78, 789–798. Montaña, C., Seghieri, J., Cornet, A., 2001. Vegetation dynamics: recruitment and regeneration in two-phase mosaics. In: Tongway, D.J., Valentin, C., Seghieri, J. (Eds.), Banded vegetation patterning in arid and semiarid environments. Ecological processes and consequences for management. Ecological studies 149. Springer, pp. 133–145. Niering, W.A., Whittaker, R.H., Lowe, C.H., 1963. The Saguaro: a population in relation to environment. Science 142, 15–23. Noy-Meir, I., 1973. Desert ecosystems: environment and producers. Annu. Rev. Ecol. Syst. 4, 25–51. Pascual, M., Roy, M., Guichard, F., Flierl, G., 2002. Cluster size distributions: signatures of self-organization in spatial ecologies. Phil. Trans. R. Soc. Lond. B. 357, 657–666. Perry, D.A., 1995. Self-organizing systems across scales. Trends Ecol. Evol. 10 (6), 241–244. Pugnaire, F.I., Haase, P., Puigdefabregas, J., 1996. Facilitation between higher plant species in a semiarid environment. Ecology 77 (5), 1420–1426. Rasmussen, K., Fog, B., Madsen, J.E., 2001. Desertification in reverse? Observations from northern Burkina Faso. Global Environ. Chang. 11, 271–282. Rietkerk, M., van de Koppel, J., 1997. Alternate stable states and threshold effects in semi-arid grazing systems. Oikos 79, 69–76. Rietkerk, M., Ketner, P., Burger, J., Hoorens, B., Olff, H., 2000. Multiscale soil and vegetation patchiness along a gradient of herbivore impact in a semi-arid grazing system in West Africa. Plant Ecol 148, 207–224. 379 Rietkerk, M., Boerlijst, M.C., van Langevelde, F., HilleRisLambers, R., van de Koppel, J., Kumar, L., Prins, H.H.T., de Roos, A.M., 2002. Self-organization of vegetation in arid ecosystems. Am. Nat. 160 (4), 524–530. Rietkerk, M., Dekker, S.C., de Ruiter, P.C., van de Koppel, J., 2004. Selforganized patchiness and catastrophic shifts in ecosystems. Science 305, 1926–1929. Rohani, P., Lewis, T.J., Grünbaum, D., Ruxton, G.D., 1997. Spatial selforganization in ecology: pretty patterns or robust reality? Trends Ecol. Evol. 12 (2), 70–74. Sato, K., Iwasa, Y., 2000. Pair approximation for lattice-based ecological models. In: Diekmann, U., Law, R., Metz, J.A.J. (Eds.), The Geometry of Ecological Interactions. Simplifying Spatial Complexity. Cambridge University Press, Cambridge, pp. 341–358. Scheffer, M., Carpenter, S., Foley, J.A., Folke, C., Walker, B., 2001. Catastrophic shifts in ecosystems. Nature 413, 591–596. Schlesinger, W.H., Raikes, J.A., Hartley, A.E., Cross, A.F., 1996. On the spatial pattern of soil nutrients in desert ecosystems. Ecology 77 (2), 364–374. Smith, J.L., Halvorson, J.J., Bolton, H., 1994. Spatial relationships of Soil Microbial Biomass and C and N Mineralization in a semiarid shrubsteppe ecosystem. Soil Boil. Biochem. 26 (9), 1151–1159. Stachowicz, J.J., 2001. Mutualism, facilitation, and the structure of ecological communities. Bioscience 51 (3), 235–246. van Baalen, M., 2000. Pair approximation for different spatial geometries. In: Diekmann, U., Law, R., Metz, J.A.J. (Eds.), The Geometry of Ecological Interactions. Simplifying Spatial Complexity. Cambridge University Press, Cambridge, pp. 359–387. van Baalen, M., Rand, D.A., 1998. The unit of selection in viscous populations and the evolution of altruism. J. Theor. Biol. 143, 631–648. van de Koppel, J., Herman, P.M.J., Thoolen, P., Heip, C.H.R., 2001. Do alternative stable states occur in natural ecosystems? Evidence from a tidal flat. Ecology 82 (12), 3449–3461. van de Koppel, J., Rietkerk, M., van Langevelde, F., Kumar, L., Klausmeier, C.A., Fryxell, J.M., Hearne, J.W., van Andel, J., de Ridder, N., Skidmone, A., Stroosnijder, L., Prins, H.H., 2002. Spatial heterogeneity and irreversible vegetation change in semiarid grazing systems. Am. Nat. 159 (2), 209–218. van de Koppel, J., van der Wal, D., Bakker, J.P., Herman, P.M.J., 2005. Self-organization and vegetation collapse in salt marsh ecosystems. Am. Nat. 165 (1), E1–E12. van Rheede van Oudtshoorn, K., van Rooyen, M.W., 1999. Dispersal Biology of Desert Plants. In: Cloudsley-Thompson, J.L. (Ed.). Springer, Berlin Heidelberg. von Hardenberg, J., Meron, E., Shachak, M., Zarmi, Y., 2001. Diversity of vegetation patterns and desertification. Phys. Rev. Lett. 87 (19), 198101. Walker, B.H., Abel, N., Stafford Smith, D.M., Langridge, J.L., 2002. A framework for the determinant of degradation in arid ecosystems. In: Reynolds, J.F., Stafford Smith, D.M. (Eds.), Global Desertification: Do Humans Cause Deserts? Dahlem University Press, Berlin, Germany, pp. 75–94. Webster, R., Maestre, F.T., 2004. Spatial analysis of semi-arid patchy vegetation by the cumulative distribution of match boundary spacings and transition probabilities. Environ. Ecol. Stat. 11, 257–281.