Early Warning Signals of Ecological Transitions:
Methods for Spatial Patterns

Sonia Kéfi, Vishwesha Guttal, William Brock, Stephen Carpenter, Aaron M Ellison, Valerie Livina, David Seekell, Marten Scheffer, Egbert van Nes, Vasilis Dakos

► To cite this version:

Sonia Kéfi, Vishwesha Guttal, William Brock, Stephen Carpenter, Aaron M Ellison, et al. Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns. PLoS ONE, 2014, 9 (3), 10.1371/journal.pone.0092097. [hal-01938665](https://hal.journals.pone.org/doi/10.1371/journal.pone.0092097)

HAL Id: hal-01938665

https://hal.umontpellier.fr/hal-01938665

Submitted on 25 May 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns

Sonia Kéfi¹*⁹, Vishwesha Guttal²⁹, William A. Brock³,⁴, Stephen R. Carpenter⁴,⁵, Aaron M. Ellison⁶, Valerie N. Livina⁷, David A. Seekell⁸, Marten Scheffer⁹, Egbert H. van Nes⁹, Vasilis Dakos¹⁰

¹Institut des Sciences de l'Evolution, CNRS, Université de Montpellier II, Montpellier, France, ²Centre for Ecological Sciences, Indian Institute of Science, Bangalore, India, ³Department of Economics, University of Wisconsin, Madison, Wisconsin, United States of America, ⁴Department of Economics, University of Missouri, Columbia, Missouri, United States of America, ⁵Center for Limnology, University of Wisconsin, Madison, Wisconsin, United States of America, ⁶Harvard Forest, Harvard University, Petersham, Massachusetts, United States of America, ⁷National Physical Laboratory, Hampton Road, Teddington, United Kingdom, ⁸Department of Environmental Sciences, University of Virginia, Charlottesville, Virginia, United States of America, ⁹Department of Aquatic Ecology and Water Quality Management, Wageningen University, Wageningen, The Netherlands, ¹⁰Integrative Ecology Group, Estacion Biologica de Donana, Sevilla, Spain

## Abstract

A number of ecosystems can exhibit abrupt shifts between alternative stable states. Because of their important ecological and economic consequences, recent research has focused on devising early warning signals for anticipating such abrupt ecological transitions. In particular, theoretical studies show that changes in spatial characteristics of the system could provide early warnings of approaching transitions. However, the empirical validation of these indicators lag behind their theoretical developments. Here, we summarize a range of currently available spatial early warning signals, suggest potential null models to interpret their trends, and apply them to three simulated spatial data sets of systems undergoing an abrupt transition. In addition to providing a step-by-step methodology for applying these signals to spatial data sets, we propose a statistical toolbox that may be used to help detect approaching transitions in a wide range of spatial data. We hope that our methodology together with the computer codes will stimulate the application and testing of spatial early warning signals on real spatial data.

**Citation:** Kéfi S, Guttal V, Brock WA, Carpenter SR, Ellison AM, et al. (2014) Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns. PLoS ONE 9(3): e92097. doi:10.1371/journal.pone.0092097

**Editor:** Ricard V. Solé, Universitat Pompeu Fabra, Spain

**Received:** August 2, 2013; **Accepted:** February 19, 2014; **Published:** March 21, 2014

**Copyright:** © 2014 Kéfi et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Funding:** The authors are grateful to the Santa Fe Institute and the Arizona State University for the financial support. SK acknowledges support from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 283068 (CASCADE). VG acknowledges support from a Ramalingaswami Fellowship, Department of Biotechnology, Government of India, and the Ministry of Environment and Forests, Government of India. SRC's work is supported by NSF. AME is supported by NSF award 11-44056. VNL is supported by NERC (NE/F005474/1) postdoctoral fellowship of the AXA Research Fund and grant of the National Measurement Office (2013–2016). VD acknowledges a RUBICON (NWO funded) grant and an EU Marie Curie fellowship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing Interests:** The authors have declared that no competing interests exist.

* E-mail: sonia.kefi@univ-montp2.fr

☞ These authors contributed equally to this work.

## Introduction

A range of ecosystems, from lakes and forests to rangelands and coral reefs, can exhibit multiple stable states [1]. In such ecosystems, abrupt shifts between ecological states may lead to ecological and economic losses. This happens when ecosystems reach a 'tipping point', at which they may rapidly reorganize into an alternative state with contrasting features. Such shifts have been documented not only in ecosystems, but also in a wide spectrum of complex systems including physiological systems, financial markets, and human societies [1]. However, the enormous complexity of such systems and the lack of detailed understanding of their underlying processes make it difficult to identify the points at which these systems may experience major changes. To circumvent this problem, recent research has focused on devising early warning signals of imminent transitions [2].

A number of early warning signals for ecological transitions has been proposed based on a phenomenon called 'critical slowing down' that generally occurs prior to a 'bifurcation' [3,4]. The closer a system is to a bifurcation point, the longer time it takes to

recover to its stable state upon any disturbance. Theoretical studies of ecological models suggest that either a direct measure of slow recovery rate [3–5] or its manifestations in the temporal and spatial dynamics of the system can potentially act as generic early warning signals of an impending transition [2,6–8]. This phenomenon of slowing down is expected to occur before a broad range of transitions, including, but not limited to, the so-called 'catastrophic shifts' [9]. Catastrophic shifts are a particular case of transitions that are especially relevant because of their possible association with hysteresis and their lack of reversibility [1].

Generic early warning signals evaluated on time series have attracted a lot of attention in the literature [2]. However, recent theoretical studies suggest that for ecosystems that are not well-mixed (such as drylands, boreal wetlands, or heterogeneous habitats which host mobile predators), changes in spatial characteristics of the system could provide early warnings of approaching transitions as well [5,10–13]. More generally, the spatial structure of ecosystems can provide information about the ecosystem degradation level [14–22]. Spatial information allows us to devise additional kinds of indicators thus adding to our arsenal
of early warning signals. At the same time, well-resolved spatial data is becoming increasingly available at low cost due to improved technology (such as remote sensing).

Empirical verification of these indicators, however, has not been able to keep pace with the rapid growth in theoretical studies, and a number of recent studies question the ability and practical efficiency of these indicators to anticipate upcoming regime shifts in real systems [23,24]. A few recent laboratory and field experiments [25–29], as well as analysis of climatic paleo-records [30], suggest that generic leading indicators (i.e. variance, skewness and autocorrelation at-lag-1) may indeed be detected in time series of real systems prior to transitions, but the empirical validation of the spatial indicators remains scarce [15, 17, 20, 21, 25, 31].

The discrepancy between theoretical developments and their empirical validation arises from a number of issues, such as the lack of sufficiently resolved and long term data, as well as the lack of a coherent methodological framework that outlines the steps and statistical tools necessary to detect those signals. For indicators based on time series, these issues have been addressed in a recent paper which provides a methodological guide to practitioners and managers to detect early warnings in time series [32] (see also http://www.early-warning-signals.org).

Here, we complement the previous work on detecting early warnings in time series data by providing a step-by-step methodology for detecting early warning signals in spatial data sets. We gather, for the first time, all the early warning signals proposed in the literature so far in a spatial context. We apply these metrics on model-generated data sets along a degradation gradient, and we discuss their interpretation based on a few potential null models. Our analysis mimicks a situation where an ecosystem would be degrading and where we would have access to several snapshots of an ecosystem's spatial structure taken over a period of time, or at different locations along a degradation gradient. We hope that our methodology together with the computer codes will stimulate testing and applications of spatial early warning signals on spatial data (R-code for the spatial analysis can be found at https://github.com/earlywarningtoolbox/spatial_warnings).

# Methods

## Spatial indicators

We first give a brief overview of the spatial early warning signals proposed in the literature so far. More details about the indicators and their precise mathematical formulation are provided in Appendix S1. Table 1 summarizes the spatial indicators and their expected trends along a degradation gradient.

**Slowing-down based indicators: spatial correlation and spectral properties (DFT).** Due to increased recovery time to local equilibrium after a perturbation, neighboring units become more like each other when a system approaches a bifurcation point, i.e. they become increasingly correlated [5]. The increasing spatial coherence can be quantified by the spatial correlation function, or Moran's I, between ecological states separated by a certain distance. The near-neighbor spatial correlation, the analog of autocorrelation at lag 1 for time series, is calculated for the distance between nearest neighboring units of the system.

Spatial spectral properties change as the system approaches a tipping point [13]. To quantify spectral properties, we compute the Discrete Fourier Transform (DFT) that decomposes spatial data into components of sine and cosine waves of different wavenumbers [33]. A wavenumber can be thought of as a 'spatial frequency', or the number of times that a pattern is repeated in a unit of spatial length. In a spatial data set, periodicity is visualized

as wavelength, which is inversely related to wavenumber (i.e., a small wavenumber corresponds to a large wavelength and vice-versa). As DFT is generally a complex number, we often plot the power spectrum (also 2D-periodogram) which is the magnitude of the complex DFT matrix (see Appendix S1 for details). Increased memory manifests itself as spectral reddening, i.e. spatial variation becomes increasingly concentrated at low wavenumbers [13]; in other words, long wavelength fluctuations become dominant prior to a transition [34].

Two metrics that help characterize spatial patterns for periodicity and directionality can be evaluated from the power spectrum. First, the radial-spectrum (r-spectrum) is obtained by summing the power spectrum at constant distances from the origin of the power spectrum, i.e. along concentric circles at different distances from the center. It allows evaluating the periodicity of the patterns. Periodic patterns are characterized by a peak in the power spectrum. The wavenumber at which a peak occurs corresponds to the number of times that a pattern reproduces itself within a unit area of the spatial data, and therefore contains information about the scale of the pattern. Second, the angular-spectrum (θ-spectrum) is obtained by summing the values of the power spectrum using angular sectors. It allows evaluating the isotropy (orientation) of the spatial pattern. For an isotropic data set, the θ-spectrum will show uniform amplitude at all angles, whereas for an anisotropic data set, the amplitude of the spectrum will show strong amplitudes for specific orientations [35].

**Variability based indicators: Spatial variance and spatial skewness.** Increased recovery time enroute to a bifurcation point may lead to stronger fluctuations around the equilibrium state of the system [36]. This can cause spatial variance of the system to increase prior to a transition [10,11]. Spatial variance is formally defined as the second moment around the spatial mean of the state variable. It has also been shown that the fluctuations around the mean can become increasingly asymmetric as the system approaches a bifurcation point. This is because the fluctuations in the direction of the alternative stable state take longer to return back to the equilibrium than those in the opposite direction [11]; this asymmetry can also arise due to local flickering events (i.e. occasional jumps of local units between their current and alternative state) [37]. The spatial asymmetry can be measured by spatial skewness, which is the third central moment scaled by the standard deviation.

**Patch based indicators: shapes and sizes of patches.** Many ecological systems, such as shrublands in semi-arid ecosystems and mussel beds in the intertidal, exhibit striking spatial self-organized patterns [38]. It has been suggested that the nature of local ecological interactions, such as the relative scales of competition and facilitation, can strongly influence the type of emerging spatial patterns, leading to i) regular, periodic patches with a characteristic patch size [22,38–40], or to ii) no characteristic scale of patchiness [15,39–42]. These different types of spatial structures have been observed in a range of ecosystems, however their use as potential indicators of degradation has mostly been developed in the case of drylands, where both types of spatial structures exist. In drylands, it has been shown that the early warnings depend on the type of patchiness exhibited by the ecological system [12].

In ecosystems exhibiting periodic patterns, as the level of external stress increases, a predictable sequence of self-organized patterns based on 'Turing instability' occurs. In isotropic areas (i.e. no preferred orientation of the pattern) the shape of the patterns shifts from gaps to labyrinths and to spots as the system becomes more degraded. Thus, spotted vegetation patterns have been proposed to be an early warning signal of imminent desertification
<table><caption>**Table 1.** Early warning signals of transitions in spatial data.</caption><thead><tr><th rowspan="2">Method/Indicator</th><th colspan="3">Phenomenon</th><th rowspan="2">Expected trend</th><th rowspan="2">Ref.</th></tr><tr><th>Rising memory</th><th>Rising variability</th><th>Patchiness</th></tr></thead><tbody><tr><td>Spatial correlation</td><td>x</td><td></td><td></td><td>increase</td><td>[5]</td></tr><tr><td>Return time</td><td>x</td><td></td><td></td><td>increase</td><td>[12]</td></tr><tr><td>Discrete Fourier Transform</td><td>x</td><td></td><td></td><td>spectral reddening</td><td>[13]</td></tr><tr><td>Spatial variance</td><td></td><td>x</td><td>x</td><td>increase</td><td>[10,11]</td></tr><tr><td>Spatial skewness</td><td></td><td>x</td><td>x</td><td>peaks (see caption)</td><td>[11]</td></tr><tr><td>Patch-size distributions</td><td></td><td></td><td>x</td><td>change in shape of the dist.</td><td>[15]</td></tr><tr><td>Regular spotted patterns</td><td></td><td></td><td>x</td><td>change in patch shape</td><td>[22]</td></tr><tr><td>Power spectrum</td><td></td><td></td><td>x</td><td>spectral reddening</td><td>[43]</td></tr></tbody></table>

Leading spatial indicator, the primary underlying phenomenon, the expected trend along a degradation gradient, and the original references in which these were proposed. The trend of spatial skewness depends on the nature of test data: It can show a nonmonotonic behavior (thus, a peaking) for transition from a low density state to higher density state. For discrete data, it typically shows a monotonic behavior (increasing or decreasing depending on whether it is transitioning from fully covered to bare state or the other way transition).
doi:10.1371/journal.pone.0092097.t001

in drylands characterized by periodic patterns [16,22]. In anisotropic areas with band-like patterns, the wavelength increases as the system approaches a transition [43].

In contrast to periodic patterns, there are cases where spatial processes give rise to non-periodic (irregular) patterns. In these cases, we can quantify the size of each patch and calculate the frequency of occurrence of different patch sizes. It is common practice to characterize the patchiness of these systems by a function that best describes the distribution of patch sizes. Irregular patterns may be characterized by a scale-free patch-size distribution, which means that there is no typical patch size in the ecosystem. Such a distribution may be well approximated by a pure power law [44] or by other heavy-tailed functions, such as a log-normal, a stretched exponential or a power law with cutoff [44,45]. Scale-free patch-size distributions have been observed in several ecosystems [15,20,41,42]. Computational models of drylands predict that larger vegetation patches become fragmented into smaller ones as aridity or grazing pressure increases and show an increasing deviation from a theoretical power law as the ecosystem approaches the desertification point [15,46]. Therefore, it has been hypothesized that an increasing deviation from power-law distribution of patch sizes can signal increasing degradation (but see [17,18]).

## Statistical significance tests

There are two steps in computing the spatial indicators. The first one is to compute the spatial metrics for a given snapshot. The second is to evaluate the trend of the spatial metrics along a degradation gradient (see Table 1). In doing so, we need to ensure that the spatial metrics for each snapshot and their trends differ from what would be expected by chance. A standard way to produce null models is to generate surrogate data and compare the trends in the indicators obtained from the original data to the trends obtained from the surrogate data [32]. Here, we discuss ways of obtaining null models for spatial early warning signals.

**Null models.** One way of obtaining a null model is to randomly permute or shuffle the elements of the spatial matrix, and this is also called bootstrapping. This is equivalent to a randomization procedure that removes any spatial structure from the original data but conserves the values of spatial variance and spatial skewness since these moments do not depend on the spatial arrangement of the data points. Therefore, such surrogate data

cannot act as a null model for spatial variance or skewness but only for other metrics such as spatial autocorrelation, DFT, and patch-size distribution.

To devise a null model for spatial variance and skewness, Eby, Guttal and others (unpublished data) propose a coarse-graining method which should be applied for both the reshuffled matrix and the real data matrix. In this method, we first divide the full matrix of dimension $n \times n$ into nonoverlapping submatrices of size $s \times s$. We then replace each submatrix by its average to obtain a smaller ‘coarse-grained matrix’ of size $c_g \times c_g$ (note that $c_g=n/s$). The basic intuition behind the method is as follows: consider any two non-overlapping submatrices of dimension (e.g. $5 \times 5$) from the reshuffled matrix. Since the reshuffled matrix is equivalent to a random matrix, the average of the entries of the two sub-matrices chosen would be roughly equal to the average of the full matrix. This exercise of ‘coarse-graining’ necessarily reduces variability across submatrices in the case of a reshuffled (thus, random) matrix. Now, consider two non-overlapping submatrices in the real data. If we expect that the real data contains a non-random spatial pattern, the average of the entire of the two sub-matrices need not be of comparable value to each other nor with the average of the full real data matrix. Therefore, in contrast to the reshuffled matrix case, coarse-graining will not necessarily reduce variability in the real data, especially if it contains a spatial pattern. Since variability determines spatial variance and skewness, the coarse-graining applied to the reshuffled matrix provides a null model for spatial variance and spatial skewness.

An alternative method of building null models is to construct a spatial matrix from a continuous stochastic process. This method is applicable when continuous data (such as biomass density) is available at each spatial point as in data set 1 and 3. More specifically, we construct a null model matrix where each entry is a random number (e.g., from a normal distribution) whose mean and variance are equal to the mean and variance of the original data matrix, respectively. This approach provides a null model for spatial skewness, correlation, and DFT, but not for spatial variance since variance is, by construction, fixed to be the same as the one from the original matrix (see figures in Appendix S3). However, it may be claimed, following the same arguments as above, that the coarse-graining method can help estimate statistical significance of spatial variance.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Spatial patterns along a degradation gradient in the three data sets. In each row, the system approaches the bifurcation point from left to right (see Fig. S1 in Appendix S2 to visualize the location of the four snapshots along the degradation gradient.) First row: local positive feedback model (data set 1). Middle row: local facilitation model (data set 2). Bottom row: scale-dependent feedback model (data set 3). In each panel, darker cells correspond to higher vegetation biomass.
doi:10.1371/journal.pone.0092097.g001</figcaption></figure>

**The particular case of patchy ecosystems.** A system with patchiness may show characteristic scales in the *r*- and *θ*-spectra. These *r*- and *θ*-spectra of the data can be compared to those obtained from a null model. In the case where there is no spatial structure, it is known analytically that the values of the scaled power spectrum observed in each bin when calculating the *r*- or *θ*-spectrum should be distributed as $\frac{\chi^2_k}{2k}$ with *k* the number of values in the bin (see for instance [35,47]; this gives very similar results as the confidence intervals presented in Appendix S3).

Once the patches are identified and their size evaluated, various heavy tailed (e.g. power law, power law with a cut-off, log-normal, etc) and non-heavy tailed distributions (e.g. exponential) can be fitted to the patch-size distribution. One way of fitting a given distribution to data has been to use ordinary least square regression on the log-log transformed probability distribution function of patch sizes. However, this method is known to have substantial bias in estimating the parameter values, especially for small data sets [44,45]. Maximum likelihood methods or least-square fits of the inverse-cumulative distribution (which quantifies the number of patches whose size is larger than a given value *s* for different values of *s*) provide more accurate estimates of the parameters of most heavy-tailed functions [48].

**Trends.** The above methods inform us about whether the spatial indicators for a given spatial data set are significantly different from those of random patterns. However, to anticipate ecological transitions, we also need to know how these spatial indicators are changing along a degradation gradient. Trend statistics like the Kendall's τ [12,30] or Pearson's correlation coefficient [25,49] can be used to quantify the strength of the trend in the indicator along the gradient.

## Simulated spatial data sets

Spatial data in ecology are typically obtained by field studies, data collecting devices placed at various locations of an ecosystem or extracted from spatial imagery. In any of these cases, the nature

of data at a given spatial location can be of two types: (a) a discrete occupancy data, such as presence or absence of vegetation (or species) at each pixel of an image, or (b) a continuous variable, such as NDVI (Normalized Difference Vegetation Index) at each pixel or nutrient concentration at each sampling point.

Here, to serve better our method-illustration purposes, we chose to work on model-generated data rather than real data. We thereby circumvent limitations of missing or noisy data, and avoid issues of misinterpretation arising from potential insufficient knowledge about the underlying degradation gradients in real data. We generated three synthetic data sets using three representative models of tipping points and self-organization in ecological systems. The three models treat ecological variables in a spatially-explicit framework with stochasticity. They all describe vegetation dynamics under resource limitation or grazing pressure, but they differ in the nature of ecological interactions and the emerging spatial vegetation structure. A detailed description of the models can be found in Appendix S2 but a brief description follows.

* Data set 1 was obtained from a local positive feedback model resulting in a non-patchy vegetation structure [50,51] (Fig. 1 first row). Space is represented as a two-dimensional lattice [52,53]. Locally, vegetation density grows logistically and is lost due to grazing. Biomass and water are exchanged between neighboring sites at a certain rate, such that a site with high biomass (or water) will have the tendency to diffuse biomass (or water) to its neighboring sites. As rainfall falls below a certain threshold, the ecosystem undergoes an abrupt transition from a globally high vegetation density to a bare state due to the nearly synchronous shifts of each of the sites to a desert state [52].

* Data set 2 is based on a local facilitation model that exhibits spatial patterns characterized by a scale-free patch-size distribution [15]. In this stochastic cellular automaton model, an ecosystem is represented by a grid of cells, each of which
<figure><img src="image_3.png" /><figcaption>Figure 2. Flow chart of analysis to perform on a spatial data set.</figcaption></figure>

doi:10.1371/journal.pone.0092097.g002

can be in one of three discrete states: vegetated (+), empty (0) or degraded (-). Empty cells represent fertile soil whereas degraded cells represented eroded soil patches unsuitable for recolonization by vegetation. A key ecological mechanism is the positive effect of vegetation on its local neighborhood through increased regeneration of degraded cells. Because of this local facilitation, vegetated cells tend to form clusters (Fig. 1 second row). When the environmental conditions become harsher, there is a point at which the vegetation dies out and the system becomes a desert resembling a saddle-node (or fold) bifurcation.

* Data set 3 is based on a scale-dependent feedback model that results in periodic (Turing-like) spatial vegetation patterns [54]. This model is based on a three partial differential equations model describing the dynamics of vegetation biomass, soil water and surface water. Plants grow due to soil water availability and die due to natural mortality and/or grazing. The infiltration rate of water in the soil is higher in areas with vegetation than in bare soil, leading to the accumulation of water under vegetation and to its depletion further away, resulting in a scale-dependent feedback responsible for the formation of regular vegetation patterns [54] (Fig. 1 last row). When water availability becomes limited, a homogeneous vegetated state becomes unstable leading to self-organized patterns such as gaps, labyrinths and spots. A further reduction in water availability leads to a transition into a desert state, again mimicking a fold-like bifurcation.

The three models exhibit a bifurcation from one state (e.g., vegetated) to an alternative state (e.g., desert) as an external parameter (such as rainfall, grazing, etc) changes. See Appendix S2 for underlying mathematical equations and parameter values of

these models. In a previous study [12], it has been shown that the three systems take increasingly longer to recover to their equilibrium after perturbation, thus demonstrating that critical slowing down is a generic feature of the transitions observed in these three ecological models, regardless of their different underlying mechanisms and their different types of spatial structures.

For our analysis, we selected ten snapshots (i.e. two-dimensional space discretized into matrices recording the spatial spread of the vegetation at the end of the simulation) for each of these models at different points along a gradient of degrading conditions prior to the transition. We illustrate our analyses using only ten of these points which are not equally spaced along the gradient (their location is shown on Fig. S1 in Appendix S2).

We are interested in quantifying how the spatial characteristics of these matrices change when approaching a tipping, or bifurcation, point. We are considering cases where the whole ecosystem shifts to an alternative state. In our mathematical representation of the ecosystem, this is equivalent to the entire matrix undergoing a shift. The degradation sequence of the matrices might correspond to snapshots in time (e.g. temperature changing through time) or in space (e.g. herbivory pressure changing in space depending on a distance to a water point). Both types of data are relevant to evaluate and test early warning signals. However, shifts of a given spatial ecosystem in time are more commonly the type of phenomena that we are trying to anticipate.

## Results

We suggest a step by step process to decide which spatial indicators should be used (Fig. 2). The spatial statistics that need to be evaluated depend on the type of data set (with discrete or
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Radial-spectrum along a degradation gradient in the three data sets. In a row, each panel corresponds to the radial-spectrum of the system at a different location along the degradation gradient. The system approaches the bifurcation point from left to right column (as in Fig. 1). First row: local positive feedback model. Middle row: local facilitation model (original data transformed using 5×5 submatrices). Bottom row: scale-dependent feedback model. Gray areas correspond to 95% confidence intervals obtained using 200 simulations of a null model (i.e. data sets of same size generated by reshuffling the original data set).
doi:10.1371/journal.pone.0092097.g003</figcaption></figure>

continuous values) at hand. Two of the three datasets used in this paper provide quantitative, continuous data, i.e. vegetation biomass (data set 1 and 3), whereas data set 2 gives qualitative, discrete data indicating the presence or absence of vegetation in each cell. For this latter data set, we transformed the original matrix by the coarse-graining procedure described in “Statistical significance tests”. More specifically, we used submatrices of 5×5 cells in which we counted the number of cells occupied by vegetation [12]. We obtained matrices that were 25 times smaller than the original ones and where values in each of the cells ranged between 0 and 25, indicating the local abundance of vegetation. The size of the submatrix used to transform the original data may affect the behavior of the indicators.

The first question to ask is whether the patterns observed are periodic or not. The *r*-spectrum obtained from the DFT analysis provides information about whether the patterns are periodic, while the *θ*-spectrum indicates whether the patterns are isotropic (i.e. with no specific orientation) or anisotropic (i.e. with a specific orientation, e.g. band-like patterns). If the patterns are not periodic, the generic leading indicators (i.e. spatial variance, spatial skewness and spatial correlation between nearby sites) may be used [12] and the power spectrum should be checked for possible reddening [13]. In addition, if the patterns are not only irregular but also patchy (e.g. can be characterized by two phases,

one vegetated and one bare), the patch-size distribution may be plotted and estimated [15,46]. If the patterns are periodic and anisotropic, the wavelength of the pattern should be evaluated. The wavelength is equivalent to the dominant length scale of the pattern provided by the *r*-spectrum [43]. For periodic isotropic patterns, the skewness of the distribution of values of the data set (e.g. grey pixels in the case of a greyscale image) indicates the type of patterns (i.e. spots, gaps or labyrinths) [43]. Additionally, it is noteworthy that if the data set includes several replicates at each stress level, potential analysis may be performed [55–57] (see more details in Appendix S1).

Next, we present how this methodology can be applied to our three data sets.

## 1. Distinguishing periodic from non-periodic patterns

We used DFT analysis to estimate the *r*-spectra as a function of wavenumbers for all the three data sets (Fig. 3). The first and second data sets (Fig. 3, first and second row) show a noisy pattern indicating that contribution to *r*-spectra is not significant for all wavenumbers. However, the *r*-spectrum for the third data set, the scale-dependent feedback model, shows a clear peak (Fig. 3 last row) even far from the transition. The peak indicates that there is dominant wavelength (corresponding to a characteristic patch size) which is a signature of periodic patterns. Note that periodicity can
## Data set 1

No pattern

<figure><img src="image_5.png" /></figure>

## Data set 2

Scale-free pattern

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Figure 4.** Generic leading indicators in data sets 1 and 2 along a degradation gradient. In each panel, the x-values correspond to the rank of the snapshot of the system along the degradation gradient. This mimics a scenario where we would not know the exact value of the driver but where we can order the data set along a degradation gradient and see Fig. S1 in Appendix S2 to visualize the location of the four snapshots along the degradation gradient. For each of the three data sets, 10 snapshots were used. Left: local positive feedback model (data set 1). Right: local facilitation model (data set 2; original data transformed using 5×5 sub-matrices). First row: spatial variance. Second row: spatial skewness. Third row: spatial correlation at lag one. In each panel, Kendall's τ, quantifying the trend of the indicator, is indicated. Gray areas correspond to 95% confidence intervals obtained using 200 simulations of a null model (i.e. data sets of same size generated by reshuffling the original data set).
doi:10.1371/journal.pone.0092097.g004

also be seen by plotting the power spectrum, where the periodicity is visible as a ring corresponding to the dominant wavelength of the data set (see upcoming paragraph “DFT and reddening”). Therefore, we conclude that the two first data sets are not periodic whereas the last one shows some spatial periodicity.

## 2. Probing spatial early warnings for non-periodic patterns

We first focus on the case of the non-periodic patterns, i.e. data set 1 showing no clear spatial structure and data set 2 with vegetation clusters.

**Spatial correlation, variance and skewness.** In data set 1, spatial correlation at lag 1 and spatial variance increase whereas spatial skewness decreases toward the bifurcation point (Fig. 4), as expected from theory. The behavior of these indicators is very similar for data set 2, except that the spatial variance decreases just before the collapse. This difference in the behaviour of the spatial variance is due to the fact that data set 2 measures only presence or absence of vegetation (i.e qualitative data) whereas data set 1 provides biomass density at each location in space (i.e. quantitative information). We note that the spatial variance and skewness for data set 1 are identical to those of null model which was obtained by a random reshuffling; therefore, we do not see error bars. However, we used the coarse-graining method for the discrete data set 2 which provides a null model for spatial variance and skewness. See a later section on ‘Probing statistical significance’ for further comments.

**DFT and reddening.** The spatial power spectrum (or 2D-periodogram, see Eq. 5 in Appendix S1) shows a reddening of the signal, i.e., the amplitude of the power spectra increases at low wavenumbers, as the system approaches the bifurcation point (Fig. 5 first and second rows). The reddening of the power spectra provides advance warning of the transition in all data sets. That trend is even clearer on the r-spectra, which sums the values of the 2D-periodogram for all the wavenumbers and shows that the lower wavenumbers contribute more to the total variance of the data set as the system approaches the bifurcation point (Fig. 3 first and second rows).

**Non-periodic and patchy: patch-size distribution.** Data set 2 was not only characterized by non-periodic patterns, but our visual examination reveals that it also exhibits distinct patches of vegetation and bare ground. In that case, it makes sense to look at the distribution of patch sizes. A way of plotting such data is to calculate the inverse cumulative distribution, i.e. plotting the number of patches whose size is larger than a given value *s* as a function of *s*. The inverse cumulative distribution is nearly scale-free far from the bifurcation point, while its slope decreases and the distribution becomes bent (toward less large patches) as the system approaches the bifurcation point (Fig. 6 top row) [46].

For comparison, we plotted the inverse cumulative patch-size distribution of data set 3 that is also patchy but periodic. Far from the transition, after the onset of pattern formation, the periodic patterns presented a patch-size distribution characterized by a sharp cutoff (Fig. 6 bottom row). As the system approaches the bifurcation point, the value of the cutoff decreases indicating decreasing patch size in the periodic pattern.

## 3. Probing spatial early warnings for periodic patterns

In contrast to the first two data sets, data set 3 exhibits periodic patterns (Fig. 1 and 3 last rows). The θ-spectra does not indicate a strong and clear signal at any specific angle, suggesting that the patterns do not have a clear orientation, i.e. patterns are isotropic (Fig. 7 first row). When the system approaches the bifurcation point (from left to right panel on Fig. 7 second row), the distribution of values of the data set goes from one peak reflecting the absence of patterns, to a two-peak distribution due to the occurrence of both vegetation and bare soil in the system after the emergence of spatial patterns, and finally to a distribution that is skewed toward small values because of the dominance of bare soil in the system. The last distribution observed before the bifurcation point characterizes spot patterns [43] which has been hypothesized to be an indicator of imminent desertification [22].
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Power spectrum along a degradation gradient in the three data sets. In a row, the system approaches the bifurcation point, from left to right column (as in Fig. 1). For a data set of size $M \times N$, the power spectrum is typically plotted for wavenumbers up to $p = \frac{M}{2}$ and $q = \frac{N}{2}$ [47] and is scaled by the spatial variance $\sigma^2$ (i.e. the scaled power spectrum is evaluated as $\frac{I}{\sigma^2}$) [35]. Red color indicates higher values of the scaled power spectrum, $\frac{I}{\sigma^2}$. The x and y-axis correspond to the wavenumbers along these directions. First row: local positive feedback model (data set 1). Middle row: local facilitation model (data set 2; original data transformed using 5×5 sub-matrices). Bottom row: scale-dependent feedback model (data set 3).
doi:10.1371/journal.pone.0092097.g005</figcaption></figure>

## 4. Probing statistical significance

We make further comments on the statistical significance for indicators of spatial data, especially in the context of null models for spatial variance and spatial skewness. As an illustration of the coarse-graining method proposed in the “Methods” section, we compared the trends in the generic leading indicators on coarse-grained matrix obtained by different dimensions of submatrix, $s \times s$ with $s=1$ (in this case the coarse-grained matrix is identical to the original matrix), 2 and 5 but starting from the same original matrix of size 100 × 100 from our data set 2. As shown in the first row of Fig. 8, the values of spatial variance along the degradation gradient differ substantially from the ones of the null model only in the coarse-grained case (second and third column). The same result seems to be true, although the differences are not as pronounced, for spatial skewness (second row). In contrast, the coarse-graining method does not offer a good null model for spatial correlation (last row, Fig. 8). In summary, a random matrix that has the same dimensions and average as the original spatial data can act as a null model for computing spatial correlation only. On the other hand, a comparison between indicators of coarse-grained matrices of both original and random matrix data can act as a null model for spatial variance and spatial skewness.

## Discussion

In this manuscript, we presented a systematic methodology for applying spatial early warning signals of abrupt ecological

transitions. To demonstrate the methods, we employed data generated from simulations of spatially-explicit ecological models with stochasticity that showed abrupt transitions from one state to an alternative state. It is increasingly being recognized that spatial dynamics pose challenges and provide opportunities for both basic science as well as management [58]. Research on spatial dynamics in ecology is continually uncovering new patterns and mechanisms. Some of these processes are likely related to regime shifts. In this context, the main objective of this manuscript was to provide a methodological guide that can stimulate the application of spatial indicators for ecological transitions on empirical data sets from real case studies.

In recent years, much effort has been devoted to the search for ‘generic’ indicators based on the idea that there are some common behaviors across a range of complex systems as they approach a bifurcation point [59]. Taking into account the spatial organization of natural systems has revealed that many indicators may not behave in spatially-structured systems as they would in other systems. More specifically, recent theoretical studies have suggested that trends of generic leading indicators (spatial variance, spatial skewness, and spatial correlation between nearby sites) could be different in self organized patterned systems [12]. For such ecosystems, system-specific indicators may be more appropriate. In particular, when spatial patterns are periodic or regular (Turing-like), the shape of the patterns may give an idea of the proximity to the threshold where the system may undergo a regime shift [22]; specifically, spots could warn of approaching
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Inverse cumulative patch-size distributions along a degradation gradient in data sets 2 and 3. Along each row, the system approaches the bifurcation point from left to right column (as in Fig. 1). In the first row: local facilitation model (original data not transformed). In the second row: scale-dependent feedback model. Gray areas correspond to 95% confidence interval obtained using 200 simulations of a null model (i.e. data sets of same size generated by reshuffling the original matrix).
doi:10.1371/journal.pone.0092097.g006</figcaption></figure>

desertification [16,22]. If the patterns are non-periodic or irregular (in particular, if the patch-size distribution is described by a heavy-tailed distribution), the patch size distribution may contain information about the degradation level of the ecosystem but more needs to be known about the underlying ecological mechanisms to interpret the changes in the shape of this distribution [15,20,46,60]. In other words, the nature of spatial organization of ecosystems seems to be a key factor in determining what type of indicators may be employed to detect an impending ecological transition in spatially-structured systems. Therefore, a first and essential step when starting to analyze a spatial data set is to get an idea of the type of spatial organization that one is dealing with. A good knowledge of the system and its underlying ecological mechanisms (specifically those responsible for the spatial structure) are required to know which indicator to use and how to interpret the changes. Theoretical studies have started developing methodologies for inferring underlying mechanisms from a limited number of spatial snapshots [61,62]. Such knowledge will facilitate to assess the risks of ecological transitions while accounting for potential false and failed alarms.

Clearly, the usefulness of the indicators presented here depends greatly on the underlying mechanisms driving the change in the ecosystems studied. In all our analyses, we assumed that a single driver is monotonically changing, while the underlying environmental conditions are assumed relatively stable and the environmental stochasticity relatively weak. This may not always be the case. For example, spatial correlation can increase due to changes in the underlying spatial heterogeneity of the environment, or due to alterations in local ‘mixing’ (or diffusion) levels in the landscape [5]. Effects of increasing spatial variance near a critical point could be confounded by various intrinsic factors such as demographic noise arising from changes in population sizes, state-dependent/multiplicative noise. The nature of the dispersal processes between patches in fragmented landscapes may affect expected trends in the indicators. In addition, extrinsic factors such as environmental fluctuations that vary in space and time can complicate our interpretations of early warning signals. In a similar way, system-specific conditions may affect the behavior and thereby interpretation of the indicators. For example, in Mediterranean drylands, degradation is accompanied by a change in the patch-size distribution toward less large patches, whereas the opposite was
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Analysis of the periodic patterns of data set 3. In a row, the system approaches the bifurcation point from left to right column (as in Fig. 1). First row: θ-spectrum. Gray areas correspond to 95% confidence interval obtained using 200 simulations of a null model (i.e. data sets of same size generated by reshuffling the original matrix). Second row: histogram of the values of the data set (or pixels of the image).
doi:10.1371/journal.pone.0092097.g007</figcaption></figure>

observed in salt marshes because the underlying mechanisms
driving the formation of the patterns differ in the two ecosystems
[15,46,60]. Therefore, our interpretations of leading indicators are
prone to both false positives and false negatives.

To prevent such errors, knowledge of the underlying heterogene-
ity, repeated observations of the system [60], or a sufficiently
well-described system so that it can be modelled are necessary to
know what to expect along a degradation trend. As this is not
usually the case, general models have been suggested to be fitted to
time series data[63,64] to simulate surrogate scenarios for
comparing trends obtained from the original data sets. In this
case the model provides an expectation of i) whether and when a
shift is likely, and ii) what trends should look like as the system is
approaching a shift. Such approaches, however, are yet to be
developed for spatially-explicit systems. Although we presented a
couple of ways to develop null models, these are rather simplistic
as they entirely neglect any underlying spatial structuring.
Therefore, the design and selection of null models for spatial
data, which is an area of active research in ecology, is another
important avenue for further research [65]. There are other
promising avenues of research to be pursued. Composite metrics
that combine spatial patterns with their temporal dynamics could
potentially offer new and potentially more reliable indicators of
imminent transitions [66]. Ecological transitions and leading
indicators in the context of metapopulation dynamics, where
factors that stabilize or make species more vulnerable to extinction

have been extensively studied, may be a fertile ground for further
research.

More importantly, despite a few recent studies, we still lack
empirical tests of spatial early warnings. Researchers have studied
spatial warnings in laboratory populations of microbial organisms
such as Daphnia and yeast [25,31]. In these studies, individual
populations are maintained in locally well-mixed small beakers or
petri-dishes and they are ‘connected’ to other populations by
‘controlled dispersal’ where the researcher transfers a fraction of
the local population to its nearest neighbors. In the field,
researchers have employed space-for-time substitution; in this
approach, it is assumed that spatial patterns at locations with
different values of stressors (e.g. grazing or rainfall) are equivalent
to the dynamics of spatial patterns where the stressor is changing
with time. This is a widely used approach in ecology as an
alternative to long term ecological studies, for example to
investigate ecological succession or how ecosystems may respond
to climate change [67]. In the context of alternative stable states
and tipping points, this method has been employed to establish the
existence of (multiple) stable states in savanna ecosystems as a
function of rainfall [57,68], and to forecast how spatial self-
organization of semi-arid vegetation may respond to increasing
stressors such as grazing [15]. However, the application of space-
for-time method is not without limitations [69]. One needs to be
cautious about the possibility of existence of various other sources
of heterogeneity, both biotic and abiotic, along a spatial gradient
of stressor. As we have argued before, the interpretation of the
## Data set 2

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Generic leading indicators in data set 2 along a degradation gradient using different coarsening of the original data. In each panel, the *x*-value indicates the rank of the snapshot along the degradation gradient. This mimics a scenario where we would not know the exact value of the driver but where we can order the data set along a gradient and See Fig. S1 in Appendix S2 to visualize the location of the four snapshots along the degradation gradient. Left: original data. Middle: original data transformed using 2×2 sub-matrices. Middle: original data transformed using 5×5 sub-matrices. First row: spatial variance. Second row: spatial skewness. Third row: spatial correlation at lag one. In each panel, Kendall's τ quantifies the strength of the trend of the indicator. Gray areas correspond to 95% confidence intervals obtained using 200 simulations of a null model (i.e. data sets of same size generated by reshuffling the original matrix). doi:10.1371/journal.pone.0092097.g008</figcaption></figure>

trends of early warnings are crucially dependent on the nature of local ecological interactions and the patterns they produce. Therefore, sources of heterogeneity, especially those that alter the ecological processes generating spatial patterns, may therefore compound the complexity of interpretation of results based on space-for-time substitutions.

In conclusion, both the theory and the application of the spatial indicators lag behind the development of the temporal ones. Spatial patterns may however offer advantages for anticipating or detecting ecological transitions of the types studied here [11,32]. Unlike temporal indicators which require long, unbroken time series of frequent observations, spatial indicators may be evaluated even if measurements are irregular or infrequent over time. While spatial pattern measurements require intensive data collection at each time point, in many cases this may be easier than high-frequency time series sampling [5,11,13]. In both space and time, more empirical validation of the indicators currently proposed in the literature is needed. We do not yet have an example where early warning signals were used to avert an upcoming shift (they

have been used in models, experiments or retroactively). We hope that this work will stimulate further development, testing and application of spatial indicators in a broad range of ecosystems.

## Supporting Information

### Appendix S1 Spatial indicators.

(PDF)

### Appendix S2 The three models used to generate the test data sets.

(PDF)

### Appendix S3 Results compared to the null model based on white noise.

(PDF)

## Acknowledgments

This paper was based on the workshop 'Practical Methods for Analysis of Early Warnings for Regime Shifts' organized by VD, SRC and MS in the
Santa Fe Institute, in Santa Fe, New Mexico, USA, from 10 to 12 of October 2011. We thank Vincent Deblauwe for useful comments on previous versions of the manuscript and for sharing his *r*- and *θ*-spectrum code.

## References

1. Scheffer M, Carpenter SR, Foley JA, Folke C, Walker B (2001) Catastrophic shifts in ecosystems. Nature 413: 591-596.

2. Scheffer M, Bascompte J, Brock W, Brovkin V, Carpenter S, et al. (2009) Early-warning signals for critical transitions. Nature 461: 53-59.

3. Wissel C (1984) A universal law of the characteristic return time near thresholds. Oecologia 65: 101-107.

4. van Nes EH, Scheffer M (2007) Slow recovery from perturbations as a generic indicator of a nearby catastrophic regime shift. American Naturalist 169: 738-747.

5. Dakos V, van Nes E, Donangelo R, Fort H, Scheffer M (2010) Spatial correlation as leading indicator of catastrophic shifts. Theoretical Ecology 3: 163-174.

6. Held H, Kleinen T (2004) Detection of climate system bifurcations by degenerate fingerprinting. Geophysical Research Letters 31: L020972.

7. Carpenter SR, Brock WA (2006) Rising variance: a leading indicator of ecological transition. Ecology Letters 9: 308-315.

8. Guttal V, Jayaprakash C (2008) Changing skewness: an early warning signal of regime shifts in ecological systems. Ecology Letters 11: 450-460.

9. Kéfi S, Dakos V, Scheffer M, Van Nes E, Rietkerk M (2013) Early warning signals also precede non-catastrophic transitions. Oikos 122: 641-648.

10. Obórny B, Meszéná G, Szabó G (2005) Dynamics of populations on the verge of extinction. Oikos 109: 291-296.

11. Guttal V, Jayaprakash C (2009) Spatial variance and spatial skewness: leading indicators of regime shifts in spatial ecological systems. Theoretical Ecology 2: 3-12.

12. Dakos V, Kéfi S, Rietkerk M, van Nes E, Scheffer M (2011) Slowing down in spatially patterned ecosystems at the brink of collapse. The American Naturalist 177: E153-E166.

13. Carpenter S, Brock W (2010) Early warnings of regime shifts in spatial dynamics using the discrete fourier transform. Ecosphere 1: art10.

14. Barbier N, Couteron P, Lejoly J, Deblauwe V, Lejeune O (2006) Self-organized vegetation patterning as a fingerprint of climate and human impact on semi-arid ecosystems. Journal of Ecology 94: 537-547.

15. Kéfi S, Rietkerk M, Alados CL, Pueyo Y, Papanastasis VP, et al. (2007) Spatial vegetation patterns and imminent desertification in mediterranean arid ecosystems. Nature 449: 213-217.

16. Kéfi S, Eppinga M, de Ruiter P, Rietkerk M (2010) Bistability and regular spatial patterns in arid ecosystems. Theoretical Ecology 3: 257-269.

17. Maestre F, Escudero A (2010) Is the patch size distribution of vegetation a suitable indicator of desertification processes? reply. Ecology 91: 3742-3745.

18. Kéfi S, Alados C, Chaves R, Pueyo Y, Rietkerk M (2010) Is the patch size distribution of vegetation a suitable indicator of desertification processes? comment. Ecology 91: 3739-3742.

19. Pueyo S (2011) Desertification and power laws. Landscape ecology 26: 305-309.

20. Lin Y, Han G, Zhao M, Chang S (2010) Spatial vegetation patterns as early signs of desertification: a case study of a desert steppe in inner mongolia, china. Landscape ecology 25: 1519-1527.

21. Litzow M, Urban J, Laurel B (2008) Increased spatial variance accompanies reorganization of two continental shelf ecosystems. Ecological applications 18: 1331-1337.

22. Rietkerk M, Dekker SC, de Ruiter PC, van de Koppel J (2004) Self-organized patchiness and catastrophic regime shifts in ecosystems. Science 305: 1926-1929.

23. Hastings A, Wysham D (2010) Regime shifts in ecological systems can occur with no warning. Ecology Letters 13: 464-472.

24. Boettiger C, Hastings A (2013) Tipping points: From patterns to predictions. Nature 493: 157-158.

25. Drake J, Griffen B (2010) Early warning signals of extinction in deteriorating environments. Nature 467: 456-459.

26. Carpenter S, Cole J, Pace M, Batt R, Brock W, et al. (2011) Early warnings of regime shifts: A whole-ecosystem experiment. Science 332: 1079-1082.

27. Veerat Aj, Faassen EJ, Dakos V, van Nes EH, Lurling M, et al. (2012) Recovery rates reflect distance to a tipping point in a living system. Nature 481: 357-359.

28. Dai L, Vorselen D, Korolev K, Gore J (2012) Generic indicators for loss of resilience before a tipping point leading to population collapse. Science 336: 1175-1177.

29. Wang R, Dearing J, Langdon P, Zhang E, Yang X, et al. (2012) Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature 492: 419-422.

30. Dakos V, Scheffer M, Van Nes E, Brovkin V, Petoukhov V, et al. (2008) Slowing down as an early warning signal for abrupt climate change. Proceedings of the National Academy of Sciences 105: 14308-14312.

31. Dai L, Korolev K, Gore J (2013) Slower recovery in space before collapse of connected populations. Nature 496: 355-358.

## Author Contributions

Conceived and designed the experiments: SK VG VD SRC DAS VNL. Performed the experiments: SK VG. Analyzed the data: SK VG. Contributed reagents/materials/analysis tools: SK VG VD. Wrote the paper: SK VG WAB SRC AME VNL DAS MS EHvN VD.

32. Dakos V, Carpenter SR, Ellison AM, Guttal V, Ives AR, et al. (2012) Early warning signals of critical transitions: Methods for time series. PLOS ONE 7: e41010.

33. Legendre P, Legendre L (1998) Numerical ecology. Amsterdam: Elsevier Science.

34. Guttal V (2008) Applications of nonequilibrium statistical physics to ecological systems (phd thesis). The Ohio State University 3: 39-41.

35. Couteron P (2002) Quantifying change in patterned semi-arid vegetation by fourier analysis of digitized aerial photographs. International Journal of Remote Sensing 23: 3407-3425.

36. Sethna JP (2007) Statistical Mechanics: Entropy, Order Parameters, and Complexity. Oxford: Oxford University Press.

37. Dakos V, van Nes EH, Scheffer M (2013) Flickering as an early warning signal. Theoretical Ecology : DOI:10.1007/s12080-013-0186-4.

38. Rietkerk M, van de Koppel J (2008) Regular pattern formation in real ecosystems. Trends in Ecology & Evolution 23: 169-175.

39. Manor A, Shnerb M (2008) Facilitation, competition, and vegetation patchiness: From scale free distribution to patterns. Journal of Theoretical Biology 253: 838-842.

40. von Hardenberg J, Kletter A, Yizhaq H, Nathan J, Meron E (2010) Periodic versus scale-free patterns in dryland vegetation. Proceedings of the Royal Society B: Biological Sciences 2077: 1771-1776.

41. Scanlon T, Caylor K, Levin S, Rodriguez-Iturbe I (2007) Positive feedbacks promote power-law clustering of Kalahari vegetation. Nature 449: 209-212.

42. Foti R, del Jesus M, Rinaldo A, Rodriguez-Iturbe I (2013) Signs of critical transition in the everglades wetlands in response to climate and anthropogenic changes. Proceedings of the National Academy of Sciences 110: 6296-6300.

43. Deblauwe V, Couteron P, Lejeune O, Bogaert J, Barbier N (2011) Environmental modulation of self-organized periodic vegetation patterns in sudan. Ecography 34: 990-1001.

44. Newman M (2005) Power laws, pareto distributions and zipf's law. Contemporary physics 46: 323-351.

45. White E, Enquist B, Green JL (2008) On estimating the exponent of power-law frequency distributions. Ecology 89: 905-912.

46. Kéfi S, Rietkerk M, Roy M, Franc A, De Ruiter P, et al. (2011) Robust scaling in ecosystems and the meltdown of patch size distributions before extinction. Ecology letters 14: 29-35.

47. Mugglesone MA, Renshaw E (1998) Detection of geological lineations on aerial photographs using two-dimensional spectral analysis. Computers and geosciences 24: 771-784.

48. Clauset A, Shalizi CR, Newman ME (2009) Power-law distributions in empirical data. SIAM review 51: 661-703.

49. Seekell D, Carpenter S, Cline T, Pace M (2012) Conditional heteroskedasticity forecasts regime shift in a whole-ecosystem experiment. Ecosystems 15: 741-747.

50. Shnerb NM, Sarah P, Lavee H, Solomon S (2003) Reactive glass and vegetation patterns. Physical Review Letters 90: 038101.

51. Guttal V, Jayaprakash C (2007) Impact of noise on bistable ecological systems. Ecological Mod-elling 201: 420-428.

52. van Nes E, Scheffer M (2005) Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems. Ecology 86: 1797-1807.

53. Keitt T, Lewis M, Holt R (2001) Allee effects, invasion pinning, and species borders. The American Naturalist 157: 203-216.

54. Rietkerk M, Boerlijst MC, van Langevelde F, HilleRisLambers R, van de Koppel J, et al. (2002) Self-organization of vegetation in arid ecosystems. American Naturalist 160: 524-530.

55. Livina V, Kwasniok F, Lenton T (2010) Potential analysis reveals changing number of climate states during the last 60 ky. Climate of the Past 6: 77-82.

56. Livina V, Kwasniok F, Lohmann G, Kantelhardt J, Lenton T (2011) Changing climate states and stability: from plicene to present. Climate dynamics 37: 2437-2453.

57. Hirota M, Holmgren M, Nes EV, Scheffer M (2011) Global resilience of tropical forest and savanna to critical transitions. Science 334: 232-235.

58. Walters C, Pauly D, Christensen V (1999) Ecospace: Prediction of mesoscale spatial patterns in trophic relationships of exploited ecosystems, with emphasis on the impacts of marine protected areas. Ecosystems 2: 539-554.

59. Gilmore R (1981) Catastrophe Theory for Scientists and Engineers. New York: John Wiley and Sons.

60. Weerman E, Van Belzen J, Rietkerk M, Temmerman S, Kéfi S, et al. (2012) Changes in diatom patch-size distribution and degradation in spatially self-organized intertidal mudflat ecosystem. Ecology 93: 608-618.

61. Marcos-Nikolaus P, Martin-González J, Solé R (2002) Spatial forecasting: Detecting determinism from single snapshots. International Journal of Bifurcation and Chaos 12: 369-376.
62. Mocenni C, Facchini A, Vicino A (2010) Identifying the dynamics of complex spatio-temporal systems by spatial recurrence properties. Proceedings of the National Academy of Sciences of the United States of America 107: 8097–8102.

63. Boettiger C, Hastings A (2012) Quantifying limits to detection of early warning for critical transitions. Journal of the Royal Society Interface 9: 1–29.

64. Lade S, Gross T (2012) Early warning signals for critical transitions: a generalized modeling approach. Plos computational biology 8: e1002360.

65. Fortin M, Dale M (2005) Spatial analysis: a guide for ecologists. Cambridge, UK: Cambridge University Press.

66. Corrado R, Cherubini AM, Pennetta C (2013) Signals of critical transitions in ecosystems associated with fluctuations of spatial patterns. In: Noise and Fluctuations (ICNF), 2013 22nd International Conference on. IEEE, pp. 1–4.

67. Blois JL, Williams JW, Fitzpatrick MC, Jackson ST, Ferrier S (2013) Space can substitute for time in predicting climate-change effects on biodiversity. Proceedings of the National Academy of Sciences 110: 9374–9379.

68. Staver AC, Archibald S, Levin SA (2011) The global extent and determinants of savanna and forest as alternative biome states. Science 334: 230–232.

69. Johnson EA, Miyanishi K (2008) Testing the assumptions of chronosequences in succession. Ecology Letters 11: 419–431.