Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns

Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns 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: Ke´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. Sole´, Universitat Pompeu Fabra, Spain Received August 2, 2013; Accepted February 19, 2014; Published March 21, 2014 Copyright:  2014 Ke´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. recover to its stable state upon any disturbance. Theoretical studies Introduction of ecological models suggest that either a direct measure of slow A range of ecosystems, from lakes and forests to rangelands and recovery rate [3–5] or its manifestations in the temporal and coral reefs, can exhibit multiple stable states [1]. In such spatial dynamics of the system can potentially act as generic early ecosystems, abrupt shifts between ecological states may lead to warning signals of an impending transition [2,6–8]. This ecological and economic losses. This happens when ecosystems phenomenon of slowing down is expected to occur before a broad reach a ‘tipping point’, at which they may rapidly reorganize into range of transitions, including, but not limited to, the so-called an alternative state with contrasting features. Such shifts have been ‘catastrophic shifts’ [9]. Catastrophic shifts are a particular case of documented not only in ecosystems, but also in a wide spectrum of transitions that are especially relevant because of their possible complex systems including physiological systems, financial mar- association with hysteresis and their lack of reversibility [1]. kets, and human societies [1]. However, the enormous complexity Generic early warning signals evaluated on time series have of such systems and the lack of detailed understanding of their attracted a lot of attention in the literature [2]. However, recent underlying processes make it difficult to identify the points at theoretical studies suggest that for ecosystems that are not well- which these systems may experience major changes. To circum- mixed (such as drylands, boreal wetlands, or heterogeneous vent this problem, recent research has focused on devising early habitats which host mobile predators), changes in spatial warning signals of imminent transitions [2]. characteristics of the system could provide early warnings of A number of early warning signals for ecological transitions has approaching transitions as well [5,10–13]. More generally, the been proposed based on a phenomenon called ‘critical slowing spatial structure of ecosystems can provide information about the down’ that generally occurs prior to a ‘bifurcation’ [3,4]. The ecosystem degradation level [14–22]. Spatial information allows us closer a system is to a bifurcation point, the longer time it takes to to devise additional kinds of indicators thus adding to our arsenal PLOS ONE | www.plosone.org 1 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings of early warning signals. At the same time, well-resolved spatial as wavelength, which is inversely related to wavenumber (i.e., a data is becoming increasingly available at low cost due to small wavenumber corresponds to a large wavelength and vice- improved technology (such as remote sensing). versa). As DFT is generally a complex number, we often plot the Empirical verification of these indicators, however, has not been power spectrum (also 2D-periodogram) which is the magnitude of the complex DFT matrix (see Appendix S1 for details). Increased able to keep pace with the rapid growth in theoretical studies, and a number of recent studies question the ability and practical memory manifests itself as spectral reddening, i.e. spatial variation becomes increasingly concentrated at low wavenumbers [13]; in efficiency of these indicators to anticipate upcoming regime shifts in real systems [23,24]. A few recent laboratory and field other words, long wavelength fluctuations become dominant prior to a transition [34]. experiments [25–29], as well as analysis of climatic paleo-records [30], suggest that generic leading indicators (i.e. variance, skewness Two metrics that help characterize spatial patterns for and autocorrelation at-lag-1) may indeed be detected in time series periodicity and directionality can be evaluated from the power of real systems prior to transitions, but the empirical validation of spectrum. First, the radial-spectrum (r-spectrum) is obtained by the spatial indicators remains scarce [15,17,20,21,25,31]. summing the power spectrum at constant distances from the origin of the power spectrum, i.e. along concentric circles at different The discrepancy between theoretical developments and their empirical validation arises from a number of issues, such as the distances from the center. It allows evaluating the periodicity of the patterns. Periodic patterns are characterized by a peak in the lack of sufficiently resolved and long term data, as well as the lack of a coherent methodological framework that outlines the steps power spectrum. The wavenumber at which a peak occurs corresponds to the number of times that a pattern reproduces itself and statistical tools necessary to detect those signals. For indicators based on time series, these issues have been addressed in a recent within a unit area of the spatial data, and therefore contains paper which provides a methodological guide to practitioners and information about the scale of the pattern. Second, the angular- managers to detect early warnings in time series [32] (see also spectrum (h-spectrum) is obtained by summing the values of the http://www.early-warning-signals.org). power spectrum using angular sectors. It allows evaluating the isotropy (orientation) of the spatial pattern. For an isotropic data Here, we complement the previous work on detecting early warnings in time series data by providing a step-by-step set, the h-spectrum will show uniform amplitude at all angles, whereas for an anisotropic data set, the amplitude of the spectrum methodology for detecting early warning signals in spatial data sets. We gather, for the first time, all the early warning signals will show strong amplitudes for specific orientations [35]. Variability based indicators: Spatial variance and spatial proposed in the literature so far in a spatial context. We apply these metrics on model-generated data sets along a degradation skewness. Increased recovery time enroute to a bifurcation point may lead to stronger fluctuations around the equilibrium gradient, and we discuss their interpretation based on a few potential null models. Our analysis mimicks a situation where an state of the system [36]. This can cause spatial variance of the system to increase prior to a transition [10,11]. Spatial variance is ecosystem would be degrading and where we would have access to formally defined as the second moment around the spatial mean of several snapshots of an ecosystem’s spatial structure taken over a the state variable. It has also been shown that the fluctuations period of time, or at different locations along a degradation around the mean can become increasingly asymmetric as the gradient. We hope that our methodology together with the system approaches a bifurcation point. This is because the computer codes will stimulate testing and applications of spatial fluctuations in the direction of the alternative stable state take early warning signals on spatial data (R-code for the spatial longer to return back to the equilibrium than those in the opposite analysis can be found at https://github.com/ direction [11]; this asymmetry can also arise due to local flickering earlywarningtoolbox/spatial_warnings). events (i.e. occasional jumps of local units between their current and alternative state) [37]. The spatial asymmetry can be Methods measured by spatial skewness, which is the third central moment Spatial indicators scaled by the standard deviation. We first give a brief overview of the spatial early warning signals Patch based indicators: shapes and sizes of proposed in the literature so far. More details about the indicators patches. Many ecological systems, such as shrublands in semi- arid ecosystems and mussel beds in the intertidal, exhibit striking and their precise mathematical formulation are provided in Appendix S1. Table 1 summarizes the spatial indicators and their spatial self-organized patterns [38]. It has been suggested that the nature of local ecological interactions, such as the relative scales of expected trends along a degradation gradient. Slowing-down based indicators: spatial correlation and competition and facilitation, can strongly influence the type of spectral properties (DFT). Due to increased recovery time to emerging spatial patterns, leading to i) regular, periodic patches local equilibrium after a perturbation, neighboring units become with a characteristic patch size [22,38–40], or to ii) no more like each other when a system approaches a bifurcation characteristic scale of patchiness [15,39–42]. These different types point, i.e. they become increasingly correlated [5]. The increasing of spatial structures have been observed in a range of ecosystems, spatial coherence can be quantified by the spatial correlation however their use as potential indicators of degradation has mostly function, or Moran’s I, between ecological states separated by a been developed in the case of drylands, where both types of spatial certain distance. The near-neighbor spatial correlation, the analog structures exist. In drylands, it has been shown that the early of autocorrelation at lag 1 for time series, is calculated for the warnings depend on the type of patchiness exhibited by the distance between nearest neighboring units of the system. ecological system [12]. Spatial spectral properties change as the system approaches a In ecosystems exhibiting periodic patterns, as the level of tipping point [13]. To quantify spectral properties, we compute external stress increases, a predictable sequence of self-organized the Discrete Fourier Transform (DFT) that decomposes spatial patterns based on ‘Turing instability’ occurs. In isotropic areas (i.e. data into components of sine and cosine waves of different no preferred orientation of the pattern) the shape of the patterns wavenumbers [33]. A wavenumber can be thought of as a ‘spatial shifts from gaps to labyrinths and to spots as the system becomes frequency’, or the number of times that a pattern is repeated in a more degraded. Thus, spotted vegetation patterns have been unit of spatial length. In a spatial data set, periodicity is visualized proposed to be an early warning signal of imminent desertification PLOS ONE | www.plosone.org 2 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Table 1. Early warning signals of transitions in spatial data. Phenomenon Method/Indicator Expected trend Ref. Rising memory Rising variability Patchiness Spatial correlation x increase [5] Return time x increase [12] Discrete Fourier Transform x spectral reddening [13] Spatial variance x x increase [10,11] Spatial skewness x x peaks (see caption) [11] Patch-size distributions x change in shape of the dist. [15] Regular spotted patterns x change in patch shape [22] Power spectrum x spectral reddening [43] 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 cannot act as a null model for spatial variance or skewness but only anisotropic areas with band-like patterns, the wavelength increases for other metrics such as spatial autocorrelation, DFT, and patch- as the system approaches a transition [43]. size distribution. In contrast to periodic patterns, there are cases where spatial To devise a null model for spatial variance and skewness, Eby, processes give rise to non-periodic (irregular) patterns. In these Guttal and others (unpublished data) propose a coarse-graining cases, we can quantify the size of each patch and calculate the method which should be applied for both the reshuffled matrix frequency of occurrence of different patch sizes. It is common and the real data matrix. In this method, we first divide the full practice to characterize the patchiness of these systems by a matrix of dimension n|n into nonoverlapping submatrices of size function that best describes the distribution of patch sizes. s|s. We then replace each submatrix by its average to obtain a Irregular patterns may be characterized by a scale-free patch-size smaller ‘coarse-grained matrix’ of size c |c (note that c ~n=s). g g g distribution, which means that there is no typical patch size in the The basic intuition behind the method is as follows: consider any ecosystem. Such a distribution may be well approximated by a two non-overlapping submatrices of dimension (e.g. 5|5) from pure power law [44] or by other heavy-tailed functions, such as a the reshuffled matrix. Since the reshuffled matrix is equivalent to a log-normal, a stretched exponential or a power law with cutoff random matrix, the average of the entries of the two sub-matrices [44,45]. Scale-free patch-size distributions have been observed in chosen would be roughly equal to the average of the full matrix. several ecosystems [15,20,41,42]. Computational models of dry- This exercise of ‘coarse-graining’ necessarily reduces variability lands predict that larger vegetation patches become fragmented across submatrices in the case of a reshuffled (thus, random) into smaller ones as aridity or grazing pressure increases and show matrix. Now, consider two non-overlapping submatrices in the an increasing deviation from a theoretical power law as the real data. If we expect that the real data contains a non-random ecosystem approaches the desertification point [15,46]. Therefore, spatial pattern, the average of the entires of the two submatrices it has been hypothesized that an increasing deviation from power- need not be of comparable value to each other nor with the law distribution of patch sizes can signal increasing degradation average of the full real data matrix. Therefore, in contrast to the (but see [17,18]). 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 Statistical significance tests coarse-graining applied to the reshuffled matrix provides a null There are two steps in computing the spatial indicators. The model for spatial variance and spatial skewness. first one is to compute the spatial metrics for a given snapshot. The An alternative method of building null models is to construct a second is to evaluate the trend of the spatial metrics along a spatial matrix from a continuous stochastic process. This method is degradation gradient (see Table 1). In doing so, we need to ensure applicable when continuous data (such as biomass density) is that the spatial metrics for each snapshot and their trends differ available at each spatial point as in data set 1 and 3. More from what would be expected by chance. A standard way to specifically, we construct a null model matrix where each entry is a produce null models is to generate surrogate data and compare the random number (e.g., from a normal distribution) whose mean trends in the indicators obtained from the original data to the and variance are equal to the mean and variance of the original trends obtained from the surrogate data [32]. Here, we discuss data matrix, respectively. This approach provides a null model for ways of obtaining null models for spatial early warning signals. spatial skewness, correlation, and DFT, but not for spatial variance Null models. One way of obtaining a null model is to since variance is, by construction, fixed to be the same as the one randomly permute or shuffle the elements of the spatial matrix, from the original matrix (see figures in Appendix S3). However, it and this is also called bootstrapping. This is equivalent to a may be claimed, following the same arguments as above, that the randomization procedure that removes any spatial structure from coarse-graining method can help estimate statistical significance of the original data but conserves the values of spatial variance and spatial variance. spatial skewness since these moments do not depend on the spatial arrangement of the data points. Therefore, such surrogate data PLOS ONE | www.plosone.org 3 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 The particular case of patchy ecosystems. A system with of data at a given spatial location can be of two types: (a) a discrete patchiness may show characteristic scales in the r- and h-spectra. occupancy data, such as presence or absence of vegetation (or These r- and h-spectra of the data can be compared to those species) at each pixel of an image, or (b) a continuous variable, obtained from a null model. In the case where there is no spatial such as NDVI (Normalized Difference Vegetation Index) at each structure, it is known analytically that the values of the scaled pixel or nutrient concentration at each sampling point. power spectrum observed in each bin when calculating the r-or h- Here, to serve better our method-illustration purposes, we chose to work on model-generated data rather than real data. We 2k spectrum should be distributed as with k the number of values thereby circumvent limitations of missing or noisy data, and avoid 2k in the bin (see for instance [35,47]; this gives very similar results as issues of misinterpretation arising from potential insufficient the confidence intervals presented in Appendix S3). knowledge about the underlying degradation gradients in real Once the patches are identified and their size evaluated, various data. We generated three synthetic data sets using three representative models of tipping points and self-organization in 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 ecological systems. The three models treat ecological variables in a fitted to the patch-size distribution. One way of fitting a given spatially-explicit framework with stochasticity. They all describe vegetation dynamics under resource limitation or grazing pressure, distribution to data has been to use ordinary least square regression on the log-log transformed probability distribution but they differ in the nature of ecological interactions and the function of patch sizes. However, this method is known to have emerging spatial vegetation structure. A detailed description of the substantial bias in estimating the parameter values, especially for models can be found in Appendix S2 but a brief description small data sets [44,45]. Maximum likelihood methods or least- follows. square fits of the inverse-cumulative distribution (which quantifies Data set 1 was obtained from a local positive feedback model the number of patches whose size is larger than a given value s for N resulting in a non-patchy vegetation structure [50,51] (Fig. 1 different values of s) provide more accurate estimates of the first row). Space is represented as a two-dimensional lattice parameters of most heavy-tailed functions [48]. [52,53]. Locally, vegetation density grows logistically and is Trends. The above methods inform us about whether the lost due to grazing. Biomass and water are exchanged between spatial indicators for a given spatial data set are significantly neighboring sites at a certain rate, such that a site with high different from those of random patterns. However, to anticipate biomass (or water) will have the tendency to diffuse biomass (or ecological transitions, we also need to know how these spatial water) to its neighboring sites. As rainfall falls below a certain indicators are changing along a degradation gradient. Trend threshold, the ecosystem undergoes an abrupt transition from a statistics like the Kendall’s t [12,30] or Pearson’s correlation globally high vegetation density to a bare state due to the coefficient [25,49] can be used to quantify the strength of the trend nearly synchronous shifts of each of the sites to a desert state in the indicator along the gradient. [52]. Data set 2 is based on a local facilitation model that exhibits Simulated spatial data sets N spatial patterns characterized by a scale-free patch-size Spatial data in ecology are typically obtained by field studies, distribution [15]. In this stochastic cellular automaton model, data collecting devices placed at various locations of an ecosystem an ecosystem is represented by a grid of cells, each of which or extracted from spatial imagery. In any of these cases, the nature PLOS ONE | www.plosone.org 4 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Figure 2. Flow chart of analysis to perform on a spatial data set. doi:10.1371/journal.pone.0092097.g002 can be in one of three discrete states: vegetated (+), empty (o) these models. In a previous study [12], it has been shown that the three systems take increasingly longer to recover to their or degraded (2). Empty cells represent fertile soil whereas degraded cells represented eroded soil patches unsuitable for equilibrium after perturbation, thus demonstrating that critical slowing down is a generic feature of the transitions observed in recolonization by vegetation. A key ecological mechanism is the positive effect of vegetation on its local neighborhood these three ecological models, regardless of their different underlying mechanisms and their different types of spatial through increased regeneration of degraded cells. Because of structures. this local facilitation, vegetated cells tend to form clusters (Fig. 1 For our analysis, we selected ten snapshots (i.e. two-dimentional second row). When the environmental conditions become space discretized into matrices recording the spatial spread of the harsher, there is a point at which the vegetation dies out and vegetation at the end of the simulation) for each of these models at the system becomes a desert resembling a saddle-node (or fold) different points along a gradient of degrading conditions prior to bifurcation. the transition. We illustrate our analyses using only ten of these Data set 3 is based on a scale-dependent feedback model that points which are not equally spaced along the gradient (their results in periodic (Turing-like) spatial vegetation patterns [54]. location is shown on Fig. S1 in Appendix S2). This model is based on a three partial differential equations We are interested in quantifying how the spatial characteristics model describing the dynamics of vegetation biomass, soil of these matrices change when approaching a tipping, or water and surface water. Plants grow due to soil water bifurcation, point. We are considering cases where the whole availability and die due to natural mortality and/or grazing. ecosystem shifts to an alternative state. In our mathematical The infiltration rate of water in the soil is higher in areas with representation of the ecosystem, this is equivalent to the entire vegetation than in bare soil, leading to the accumulation of matrix undergoing a shift. The degradation sequence of the water under vegetation and to its depletion further away, matrices might correspond to snapshots in time (e.g. temperature resulting in a scale-dependent feedback responsible for the changing through time) or in space (e.g. herbivory pressure formation of regular vegetation patterns [54] (Fig. 1 last row). changing in space depending on a distance to a water point). Both When water availability becomes limited, a homogeneous types of data are relevant to evaluate and test early warning vegetated state becomes unstable leading to self-organized signals. However, shifts of a given spatial ecosystem in time are patterns such as gaps, labyrinths and spots. A further reduction more commonly the type of phenomena that we are trying to in water availability leads to a transition into a desert state, anticipate. again mimicking a fold-like bifurcation. Results The three models exhibit a bifurcation from one state (e.g., vegetated) to an alternative state (e.g., desert) as an external We suggest a step by step process to decide which spatial parameter (such as rainfall, grazing, etc) changes. See Appendix S2 indicators should be used (Fig. 2). The spatial statistics that need to for underlying mathematical equations and parameter values of be evaluated depend on the type of data set (with discrete or PLOS ONE | www.plosone.org 5 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 565 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 continuous values) at hand. Two of the three datasets used in this one vegetated and one bare), the patch-size distribution may be paper provide quantitative, continuous data, i.e. vegetation plotted and estimated [15,46]. If the patterns are periodic and anisotropic, the wavelength of the pattern should be evaluated. biomass (data set 1 and 3), whereas data set 2 gives qualitative, discrete data indicating the presence or absence of vegetation in The wavelength is equivalent to the dominant length scale of the pattern provided by the r-spectrum [43]. For periodic isotropic each cell. For this latter data set, we transformed the original matrix by the coarse-graining procedure described in ‘‘Stastistical 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 significance tests’’. More specifically, we used submatrices of 565 of patterns (i.e. spots, gaps or labyrinths) [43]. Additionally, it is cells in which we counted the number of cells occupied by noteworthy that if the data set includes several replicates at each vegetation [12]. We obtained matrices that were 25 times smaller stress level, potential analysis may be performed [55–57] (see more than the original ones and where values in each of the cells ranged details in Appendix S1). between 0 and 25, indicating the local abundance of vegetation. Next, we present how this methodology can be applied to our The size of the submatrix used to transform the original data may affect the behavior of the indicators. three data sets. The first question to ask is whether the patterns observed are periodic or not. The r-spectrum obtained from the DFT analysis 1. Distinguishing periodic from non-periodic patterns provides information about whether the patterns are periodic, We used DFT analysis to estimate the r-spectra as a function of while the h-spectrum indicates whether the patterns are isotropic wavenumbers for all the three data sets (Fig. 3). The first and (i.e. with no specific orientation) or anisotropic (i.e. with a specific second data sets (Fig. 3, first and second row) show a noisy pattern orientation, e.g. band-like patterns). If the patterns are not indicating that contribution to r-spectra is not significant for all periodic, the generic leading indicators (i.e. spatial variance, wavenumbers. However, the r-spectrum for the third data set, the spatial skewness and spatial correlation between nearby sites) may scale-dependent feedback model, shows a clear peak (Fig. 3 last be used [12] and the power spectrum should be checked for row) even far from the transition. The peak indicates that there is possible reddening [13]. In addition, if the patterns are not only dominant wavelength (corresponding to a characteristic patch size) irregular but also patchy (e.g. can be characterized by two phases, which is a signature of periodic patterns. Note that periodicity can PLOS ONE | www.plosone.org 6 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 Figure 4. Generic leading indicators in data sets 1 and 2 along system approaches the bifurcation point (Fig. 6 top row) [46]. a degradation gradient. In each panel, the x-values correspond to the rank of the snapshot of the system along the degradation gradient. For comparison, we plotted the inverse cumulative patch-size This mimicks a scenario where we would not know the exact value of distribution of data set 3 that is also patchy but periodic. Far from the driver but where we can order the data set along a degradation the transition, after the onset of pattern formation, the periodic gradient and see Fig. S1 in Appendix S2 to visualize the location of the patterns presented a patch-size distribution characterized by a four snapshots along the degradation gradient. For each of the three data sets, 10 snapshots were used. Left: local positive feedback model sharp cutoff (Fig. 6 bottom row). As the system approaches the (data set 1). Right: local facilitation model (data set 2; original data bifurcation point, the value of the cutoff decreases indicating transformed using 565 sub-matrices). First row: spatial variance. Second decreasing patch size in the periodic pattern. row: spatial skewness. Third row: spatial correlation at lag one. In each panel, Kendall’s t, quantifying the trend of the indicator, is indicated. 3. Probing spatial early warnings for periodic patterns Gray areas correspond to 95% confidence intervals obtained using 200 simulations of a null model (i.e. data sets of same size generated by In contrast to the first two data sets, data set 3 exhibits periodic reshuffling the original data set). patterns (Fig. 1 and 3 last rows). The h-spectra does not indicate a doi:10.1371/journal.pone.0092097.g004 strong and clear signal at any specific angle, suggesting that the patterns do not have a clear orientation, i.e. patterns are isotropic also be seen by plotting the power spectrum, where the periodicity (Fig. 7 first row). When the system approaches the bifurcation is visible as a ring corresponding to the dominant wavelength of point (from left to right panel on Fig. 7 second row), the the data set (see upcoming paragraph ‘‘DFT and reddening’’). distribution of values of the data set goes from one peak reflecting Therefore, we conclude that the two first data sets are not periodic the absence of patterns, to a two-peak distribution due to the whereas the last one shows some spatial periodicity. occurrence of both vegetation and bare soil in the system after the emergence of spatial patterns, and finally to a distribution that is 2. Probing spatial early warnings for non-periodic skewed toward small values because of the dominance of bare soil in the system. The last distribution observed before the bifurcation patterns point characterizes spot patterns [43] which has been hypothe- We first focus on the case of the non-periodic patterns, i.e. data sized to be an indicator of imminent desertification [22]. set 1 showing no clear spatial structure and data set 2 with vegetation clusters. PLOS ONE | www.plosone.org 7 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Figure 5. Power spectrum along a degradation gradient in the three data sets. In a row, the system approaches the bifurcation point, from M N left to right column (as in Fig. 1). For a data set of size M|N, the power spectrum is typically plotted for wavenumbers up to p~ and q~ [47] 2 2 and is scaled by the spatial variance s (i.e. the scaled power spectrum is evaluated as ) [35]. Red color indicates higher values of the scaled power spectrum, . 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 565 sub-matrices). Bottom row: scale-dependent feedback model (data set 3). doi:10.1371/journal.pone.0092097.g005 transitions. To demonstrate the methods, we employed data 4. Probing statistical significance generated from simulations of spatially-explicit ecological models We make further comments on the statistical significance for with stochasticity that showed abrupt transitions from one state to indicators of spatial data, especially in the context of null models an alternative state. It is increasingly being recognized that spatial for spatial variance and spatial skewness. As an illustration of the dynamics pose challenges and provide opportunities for both basic coarse-graining method proposed in the ‘‘Methods’’ section, we science as well as management [58]. Research on spatial dynamics compared the trends in the generic leading indicators on coarse- in ecology is continually uncovering new patterns and mecha- grained matrix obtained by different dimensions of submatrix, nisms. Some of these processes are likely related to regime shifts. In s|s with s~1 (in this case the coarse-grained matrix is identical to this context, the main objective of this manuscript was to provide a the original matrix), 2 and 5 but starting from the same original methodological guide that can stimulate the application of spatial matrix of size 100|100 from our data set 2. As shown in the first indicators for ecological transitions on empirical data sets from real row of Fig. 8, the values of spatial variance along the degradation case studies. gradient differ substantially from the ones of the null model only in In recent years, much effort has been devoted to the search for the coarse-grained case (second and third column). The same ‘generic’ indicators based on the idea that there are some common result seems to be true, although the differences are not as behaviors across a range of complex systems as they approach a pronounced, for spatial skewness (second row). In contrast, the bifurcation point [59]. Taking into account the spatial organiza- coarse-graining method does not offer a good null model for tion of natural systems has revealed that many indicators may not spatial correlation (last row, Fig. 8). In summary, a random matrix behave in spatially-structured systems as they would in other that has the same dimensions and average as the original spatial systems. More specifically, recent theoretical studies have suggest- data can act as a null model for computing spatial correlation only. ed that trends of generic leading indicators (spatial variance, On the other hand, a comparison between indicators of coarse- spatial skewness, and spatial correlation between nearby sites) grained matrices of both original and random matrix data can act could be different in self organized patterned systems [12]. For as a null model for spatial variance and spatial skewness. such ecosystems, system-specific indicators may be more appro- priate. In particular, when spatial patterns are periodic or regular Discussion (Turing-like), the shape of the patterns may give an idea of the In this manuscript, we presented a systematic methodology for proximity to the threshold where the system may undergo a applying spatial early warning signals of abrupt ecological regime shift [22]; specifically, spots could warn of approaching PLOS ONE | www.plosone.org 8 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 colum (as in Fig. 1)n. First row: local facilitation model (original data not transformed). Bottom row: scale-dependent feedback model. Gray areas corresponds 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 desertification [16,22]. If the patterns are non-periodic or irregular Clearly, the usefulness of the indicators presented here depends (in particular, if the patch-size distribution is described by a heavy- greatly on the underlying mechanisms driving the change in the ecosystems studied. In all our analyses, we assumed that a single tailed distribution), the patch size distribution may contain information about the degradation level of the ecosystem but driver is monotonically changing, while the underlying environ- mental conditions are assumed relatively stable and the environ- more needs to be known about the underlying ecological mental stochasticity relatively weak. This may not always be the mechanisms to interpret the changes in the shape of this distribution [15,20,46,60]. In other words, the nature of spatial case. For example, spatial correlation can increase due to changes in the underlying spatial heterogeneity of the environment, or due organization of ecosystems seems to be a key factor in determining to alterations in local ‘mixing’ (or diffusion) levels in the landscape what type of indicators may be employed to detect an impending [5]. Effects of increasing spatial variance near a critical point could ecological transition in spatially-structured systems. Therefore, a be confounded by various intrinsic factors such as demographic first and essential step when starting to analyze a spatial data set is noise arising from changes in population sizes, state-dependent/ to get an idea of the type of spatial organization that one is dealing multiplicative noise. The nature of the dispersal processes between with. A good knowledge of the system and its underlying ecological patches in fragmented landscapes may affect expected trends in mechanisms (specifically those responsible for the spatial structure) the indicators. In addition, extrinsic factors such as environmental are required to know which indicator to use and how to interpret fluctuations that vary in space and time can complicate our the changes. Theoretical studies have started developing method- interpretations of early warning signals. In a similar way, system- ologies for inferring underlying mechanisms from a limited specific conditions may affect the behavior and thereby interpre- number of spatial snapshots [61,62]. Such knowledge will facilitate tation of the indicators. For example, in Mediterranean drylands, to assess the risks of ecological transitions while accounting for degradation is accompanied by a change in the patch-size potential false and failed alarms. distribution toward less large patches, whereas the opposite was PLOS ONE | www.plosone.org 9 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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: h-spectrum. Gray areas corresponds 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 observed in salt marshes because the underlying mechanisms have been extensively studied, may be a fertile ground for further driving the formation of the patterns differ in the two ecosystems research. [15,46,60]. Therefore, our interpretations of leading indicators are More importantly, despite a few recent studies, we still lack prone to both false positives and false negatives. empirical tests of spatial early warnings. Researchers have studied To prevent such errors, knowledge of the underlying heteroge- spatial warnings in laboratory populations of microbial organisms neity, repeated observations of the system [60], or a sufficiently such as Daphnia and yeast [25,31]. In these studies, individual well-described system so that it can be modelled are necessary to populations are maintained in locally well-mixed small beakers or know what to expect along a degradation trend. As this is not petri-dishes and they are ‘connected’ to other populations by usually the case, general models have been suggested to be fitted to ‘controlled dispersal’ where the researcher transfers a fraction of time series data[63,64] to simulate surrogate scenarios for the local population to its nearest neighbors. In the field, comparing trends obtained from the original data sets. In this researchers have employed space-for-time substitution; in this case the model provides an expectation of i) whether and when a approach, it is assumed that spatial patterns at locations with shift is likely, and ii) what trends should look like as the system is different values of stressors (e.g. grazing or rainfall) are equivalent approaching a shift. Such approaches, however, are yet to be to the dynamics of spatial patterns where the stressor is changing developed for spatially-explicit systems. Although we presented a with time. This is a widely used approach in ecology as an couple of ways to develop null models, these are rather simplistic alternative to long term ecological studies, for example to as they entirely neglect any underlying spatial structuring. investigate ecological succession or how ecosystems may respond Therefore, the design and selection of null models for spatial to climate change [67]. In the context of alternative stable states data, which is an area of active research in ecology, is another and tipping points, this method has been employed to establish the important avenue for further research [65]. There are other existence of (multiple) stable states in savanna ecosystems as a promising avenues of research to be pursued. Composite metrics function of rainfall [57,68], and to forecast how spatial self- that combine spatial patterns with their temporal dynamics could organization of semi-arid vegetation may respond to increasing potentially offer new and potentially more reliable indicators of stressors such as grazing [15]. However, the application of space- imminent transitions [66]. Ecological transitions and leading for-time method is not without limitations [69]. One needs to be indicators in the context of metapopulation dynamics, where cautious about the possibility of existence of various other sources factors that stabilize or make species more vulnerable to extinction of heterogeneity, both biotic and abiotic, along a spatial gradient of stressor. As we have argued before, the interpretation of the PLOS ONE | www.plosone.org 10 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 mimicks 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 262 sub-matrices. Middle: original data transformed using 565 sub-matrices. First row: spatial variance. Second row: spatial skewness. Third row: spatial correlation at lag one. In each panel, Kendall’s t 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 trends of early warnings are crucially dependent on the nature of have been used in models, experiments or retroactively). We hope local ecological interactions and the patterns they produce. that this work will stimulate further development, testing and application of spatial indicators in a broad range of ecosystems. 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 Supporting Information space-for-time substitutions. Appendix S1 Spatial indicators. In conclusion, both the theory and the application of the spatial (PDF) indicators lag behind the development of the temporal ones. Spatial patterns may however offer advantages for anticipating or Appendix S2 The three models used to generate the test detecting ecological transitions of the types studied here [11,32]. data sets. Unlike temporal indicators which require long, unbroken time (PDF) series of frequent observations, spatial indicators may be evaluated Appendix S3 Results compared to the null model based even if measurements are irregular or infrequent over time. While on white noise. spatial pattern measurements require intensive data collection at (PDF) each time point, in many cases this may be easier than high- frequency time series sampling [5,11,13]. In both space and time, Acknowledgments more empirical validation of the indicators currently proposed in the literature is needed. We do not yet have an example where This paper was based on the workshop ‘Practical Methods for Analysis of early warning signals were used to avert an upcoming shift (they Early Warnings for Regime Shifts’ organized by VD, SRC and MS in the PLOS ONE | www.plosone.org 11 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Santa Fe Institute, in Santa Fe, New Mexico, USA, from 10 to 12 of Author Contributions October 2011. We thank Vincent Deblauwe for useful comments on Conceived and designed the experiments: SK VG VD SRC DAS VNL. previous versions of the manuscript and for sharing his r - and h - spectrum Performed the experiments: SK VG. Analyzed the data: SK VG. code. Contributed reagents/materials/analysis tools: SK VG VD. Wrote the paper: SK VG WAB SRC AME VNL DAS MS EHvN VD. References 1. Scheffer M, Carpenter SR, Foley JA, Folke C, Walker B (2001) Catastrophic 32. Dakos V, Carpenter SR, Ellison AM, Guttal V, Ives AR, et al. (2012) Early shifts in ecosystems. Nature 413: 591–596. warning signals of critical transitions: Methods for time series. PLOS ONE 7: 2. Scheffer M, Bascompte J, Brock W, Brovkin V, Carpenter S, et al. (2009) Early- e41010. warning signals for critical transitions. Nature 461: 53–59. 33. Legendre P, Legendre L (1998) Numerical ecology. Amsterdam: Elsevier 3. Wissel C (1984) A universal law of the characteristic return time near thresholds. Science. Oecologia 65: 101–107. 34. Guttal V (2008) Applications of nonequilibrium statistical physics to ecological 4. van Nes EH, Scheffer M (2007) Slow recovery from perturbations as a generic systems (phd thesis). The Ohio State University 3: 39–41. indicator of a nearby catastrophic regime shift. American Naturalist 169: 738– 35. Couteron P (2002) Quantifying change in patterned semi-arid vegetation by 747. fourier analysis of digitized aerial photographs. International Journal of Remote 5. Dakos V, van Nes E, Donangelo R, Fort H, Scheffer M (2010) Spatial Sensing 23: 3407–3425. correlation as leading indicator of catastrophic shifts. Theoretical Ecology 3: 36. Sethna JP (2007) Statistical Mechanics: Entropy, Order Parameters, and 163–174. Complexity. Oxford: Oxford University Press. 6. Held H, Kleinen T (2004) Detection of climate system bifurcations by 37. Dakos V, van Nes EH, Scheffer M (2013) Flickering as an early warning signal. degenerate fingerprinting. Geophysical Research Letters 31: L020972. Theoretical Ecology : DOI:10.1007/s12080-013-0186-4. 7. Carpenter SR, Brock WA (2006) Rising variance: a leading indicator of 38. Rietkerk M, van de Koppel J (2008) Regular pattern formation in real ecological transition. Ecology Letters 9: 308–315. ecosystems. Trends in Ecology & Evolution 23: 169–175. 8. Guttal V, Jayaprakash C (2008) Changing skewness: an early warning signal of 39. Manor A, Shnerb M (2008) Facilitation, competition, and vegetation patchiness: regime shifts in ecological systems. Ecology Letters 11: 450–460. From scale free distribution to patterns. Journal of Theoretical Biology 253: 9. Ke´fi S, Dakos V, Scheffer M, Van Nes E, Rietkerk M (2013) Early warning 838–842. signals also precede non-catastrophic transitions. Oikos 122: 641–648. 40. von Hardenberg J, Kletter A, Yizhaq H, Nathan J, Meron E (2010) Periodic 10. Obo´rny B, Mesze´na G, Szabo´ G (2005) Dynamics of populations on the verge of versus scale-free patterns in dryland vegetation. Proceedings of the Royal Society extinction. Oikos 109: 291–296. B: Biological Sciences 2077: 1771–1776. 11. Guttal V, Jayaprakash C (2009) Spatial variance and spatial skewness: leading 41. Scanlon T, Caylor K, Levin S, Rodriguez-Iturbe I (2007) Positive feedbacks indicators of regime shifts in spatial ecological systems. Theoretical Ecology 2: 3– promote power-law clustering of Kalahari vegetation. Nature 449: 209–212. 12. 42. Foti R, del Jesus M, Rinaldo A, Rodriguez-Iturbe I (2013) Signs of critical 12. Dakos V, Ke´fi S, Rietkerk M, van Nes E, Scheffer M (2011) Slowing down in transition in the everglades wetlands in response to climate and anthropogenic spatially patterned ecosystems at the brink of collapse. The American Naturalist changes. Proceedings of the National Academy of Sciences 110: 6296–6300. 177: E153–E166. 43. Deblauwe V, Couteron P, Lejeune O, Bogaert J, Barbier N (2011) 13. Carpenter S, Brock W (2010) Early warnings of regime shifts in spatial dynamics Environmental modulation of self-organized periodic vegetation patterns in using the discrete fourier transform. Ecosphere 1: art10. sudan. Ecography 34: 990–1001. 44. Newman M (2005) Power laws, pareto distributions and zipf’s law. Contempo- 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 rary physics 46: 323–351. ecosystems. Journal of Ecology 94: 537–547. 45. White E, Enquist B, Green JL (2008) On estimating the exponent of power-law 15. Ke´fi S, Rietkerk M, Alados CL, Pueyo Y, Papanastasis VP, et al. (2007) Spatial frequency distributions. Ecology 89: 905–912. vegetation patterns and imminent desertification in mediterranean arid 46. Ke´fi S, Rietkerk M, Roy M, Franc A, De Ruiter P, et al. (2011) Robust scaling in ecosystems. Nature 449: 213–217. ecosystems and the meltdown of patch size distributions before extinction. Ecology letters 14: 29–35. 16. Ke´fi S, Eppinga M, de Ruiter P, Rietkerk M (2010) Bistability and regular spatial patterns in arid ecosystems. Theoretical Ecology 3: 257–269. 47. Mugglestone MA, Renshaw E (1998) Detection of geological lineations on aerial 17. Maestre F, Escudero A (2010) Is the patch size distribution of vegetation a photographs using two-dimensional spectral analysis. Computers and geosci- suitable indicator of desertification processes? reply. Ecology 91: 3742–3745. ences 24: 771–784. 18. Ke´fi S, Alados C, Chaves R, Pueyo Y, Rietkerk M (2010) Is the patch size 48. Clauset A, Shalizi CR, Newman ME (2009) Power-law distributions in empirical distribution of vegetation a suitable indicator of desertification processes? data. SIAM review 51: 661–703. comment. Ecology 91: 3739–3742. 49. Seekell D, Carpenter S, Cline T, Pace M (2012) Conditional heteroskedasticity 19. Pueyo S (2011) Desertification and power laws. Landscape ecology 26: 305–309. forecasts regime shift in a whole-ecosystem experiment. Ecosystems 15: 741– 20. Lin Y, Han G, Zhao M, Chang S (2010) Spatial vegetation patterns as early 747. signs of desertification: a case study of a desert steppe in inner mongolia, china. 50. Shnerb NM, Sarah P, Lavee H, Solomon S (2003) Reactive glass and vegetation Landscape ecology 25: 1519–1527. patterns. Physical Review Letters 90: 038101. 21. Litzow M, Urban J, Laurel B (2008) Increased spatial variance accompanies 51. Guttal V, Jayaprakash C (2007) Impact of noise on bistable ecological systems. reorganization of two continental shelf ecosystems. Ecological applications 18: Ecological Mod-elling 201: 420–428. 1331–1337. 52. van Nes E, Scheffer M (2005) Implications of spatial heterogeneity for 22. Rietkerk M, Dekker SC, de Ruiter PC, van de Koppel J (2004) Self-organized catastrophic regime shifts in ecosystems. Ecology 86: 1797–1807. patchiness and catastrophic regime shifts in ecosystems. Science 305: 1926– 53. Keitt T, Lewis M, Holt R (2001) Allee effects, invasion pinning, and species 1929. borders. the American Naturalist 157: 203–216. 23. Hastings A, Wysham D (2010) Regime shifts in ecological systems can occur with 54. Rietkerk M, Boerlijst MC, van Langevelde F, HilleRisLambers R, van de no warning. Ecology Letters 13: 464–472. Koppel J, et al. (2002) Self-organization of vegetation in arid ecosystems. 24. Boettiger C, Hastings A (2013) Tipping points: From patterns to predictions. American Naturalist 160: 524–530. Nature 493: 157–158. 55. Livina V, Kwasniok F, Lenton T (2010) Potential analysis reveals changing 25. Drake J, Griffen B (2010) Early warning signals of extinction in deteriorating number of climate states during the last 60 kyr. Climate of the Past 6: 77–82. environments. Nature 467: 456–459. 56. Livina V, Kwasniok F, Lohmann G, Kantelhardt J, Lenton T (2011) Changing 26. Carpenter S, Cole J, Pace M, Batt R, Brock W, et al. (2011) Early warnings of climate states and stability: from pliocene to present. Climate dynamics 37: regime shifts: A whole-ecosystem experiment. Science 332: 1079–1082. 2437–2453. 27. Veerat AJ, Faassen EJ, Dakos V, van Nes EH, Lurling M, et al. (2012) Recovery 57. Hirota M, Holmgren M, Nes EV, Scheffer M (2011) Global resilience of tropical rates reflect distance to a tipping point in a living system. Nature 481: 357–359. forest and savanna to critical transitions. Science 334: 232–235. 28. Dai L, Vorselen D, Korolev K, Gore J (2012) Generic indicators for loss of 58. Walters C, Pauly D, Christensen V (1999) Ecospace: Prediction of mesoscale resilience before a tipping point leading to population collapse. Science 336: spatial patterns in trophic relationships of exploited ecosystems, with emphasis 1175–1177. on the impacts of marine protected areas. Ecosystems 2: 539–554. 29. Wang R, Dearing J, Langdon P, Zhang E, Yang X, et al. (2012) Flickering gives 59. Gilmore R (1981) Catastrophe Theory for Scientists and Engineers. New York: early warning signals of a critical transition to a eutrophic lake state. Nature 492: John Wiley and Sons. 419–422. 60. Weerman E, Van Belzen J, Rietkerk M, Temmerman S, Ke´fi S, et al. (2012) 30. Dakos V, Scheffer M, Van Nes E, Brovkin V, Petoukhov V, et al. (2008) Slowing Changes in diatom patch-size distribution and degradation in spatially self- down as an early warning signal for abrupt climate change. Proceedings of the organized intertidal mudflat ecosystem. Ecology 93: 608–618. National Academy of Sciences 105: 14308–14312. 61. Marcos-Nikolaus P, Martin-Gonza´lez J, Sole´ R (2002) Spatial forecasting: 31. Dai L, Korolev K, Gore J (2013) Slower recovery in space before collapse of Detecting determinism from single snapshots. International Journal of connected populations. Nature 496: 355–358. Bifurcation and Chaos 12: 369–376. PLOS ONE | www.plosone.org 12 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 62. Mocenni C, Facchini A, Vicino A (2010) Identifying the dynamics of complex 66. Corrado R, Cherubini AM, Pennetta C (2013) Signals of critical transitions in ecosystems associated with fluctuations of spatial patterns. In: Noise and spatio-temporal systems by spatial recurrence properties. Proceedings of the Fluctuations (ICNF), 2013 22nd International Conference on. IEEE, pp. 1–4. National Academy of Sciences of the United States of America 107: 8097–8102. 67. Blois JL, Williams JW, Fitzpatrick MC, Jackson ST, Ferrier S (2013) Space can 63. Boettiger C, Hastings A (2012) Quantifying limits to detection of early warning substitute for time in predicting climate-change effects on biodiversity. for critical transitions. Journal of the Royal Society Interface 9: 1–29. Proceedings of the National Academy of Sciences 110: 9374–9379. 64. Lade S, Gross T (2012) Early warning signals for critical transitions: a 68. Staver AC, Archibald S, Levin SA (2011) The global extent and determinants of generalized modeling approach. Plos computational biology 8: e1002360. savanna and forest as alternative biome states. Science 334: 230–232. 65. Fortin M, Dale M (2005) Spatial analysis: a guide for ecologists. Cambridge, 69. Johnson EA, Miyanishi K (2008) Testing the assumptions of chronosequences in UK: Cambridge University Press. succession. Ecology Letters 11: 419–431. PLOS ONE | www.plosone.org 13 March 2014 | Volume 9 | Issue 3 | e92097 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png PLoS ONE Pubmed Central

Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns

Loading next page...
 
/lp/pubmed-central/early-warning-signals-of-ecological-transitions-methods-for-spatial-4nCnzLJlJs

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Pubmed Central
Copyright
© 2014 Kéfi et al
ISSN
1932-6203
eISSN
1932-6203
DOI
10.1371/journal.pone.0092097
Publisher site
See Article on Publisher Site

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: Ke´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. Sole´, Universitat Pompeu Fabra, Spain Received August 2, 2013; Accepted February 19, 2014; Published March 21, 2014 Copyright:  2014 Ke´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. recover to its stable state upon any disturbance. Theoretical studies Introduction of ecological models suggest that either a direct measure of slow A range of ecosystems, from lakes and forests to rangelands and recovery rate [3–5] or its manifestations in the temporal and coral reefs, can exhibit multiple stable states [1]. In such spatial dynamics of the system can potentially act as generic early ecosystems, abrupt shifts between ecological states may lead to warning signals of an impending transition [2,6–8]. This ecological and economic losses. This happens when ecosystems phenomenon of slowing down is expected to occur before a broad reach a ‘tipping point’, at which they may rapidly reorganize into range of transitions, including, but not limited to, the so-called an alternative state with contrasting features. Such shifts have been ‘catastrophic shifts’ [9]. Catastrophic shifts are a particular case of documented not only in ecosystems, but also in a wide spectrum of transitions that are especially relevant because of their possible complex systems including physiological systems, financial mar- association with hysteresis and their lack of reversibility [1]. kets, and human societies [1]. However, the enormous complexity Generic early warning signals evaluated on time series have of such systems and the lack of detailed understanding of their attracted a lot of attention in the literature [2]. However, recent underlying processes make it difficult to identify the points at theoretical studies suggest that for ecosystems that are not well- which these systems may experience major changes. To circum- mixed (such as drylands, boreal wetlands, or heterogeneous vent this problem, recent research has focused on devising early habitats which host mobile predators), changes in spatial warning signals of imminent transitions [2]. characteristics of the system could provide early warnings of A number of early warning signals for ecological transitions has approaching transitions as well [5,10–13]. More generally, the been proposed based on a phenomenon called ‘critical slowing spatial structure of ecosystems can provide information about the down’ that generally occurs prior to a ‘bifurcation’ [3,4]. The ecosystem degradation level [14–22]. Spatial information allows us closer a system is to a bifurcation point, the longer time it takes to to devise additional kinds of indicators thus adding to our arsenal PLOS ONE | www.plosone.org 1 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings of early warning signals. At the same time, well-resolved spatial as wavelength, which is inversely related to wavenumber (i.e., a data is becoming increasingly available at low cost due to small wavenumber corresponds to a large wavelength and vice- improved technology (such as remote sensing). versa). As DFT is generally a complex number, we often plot the Empirical verification of these indicators, however, has not been power spectrum (also 2D-periodogram) which is the magnitude of the complex DFT matrix (see Appendix S1 for details). Increased able to keep pace with the rapid growth in theoretical studies, and a number of recent studies question the ability and practical memory manifests itself as spectral reddening, i.e. spatial variation becomes increasingly concentrated at low wavenumbers [13]; in efficiency of these indicators to anticipate upcoming regime shifts in real systems [23,24]. A few recent laboratory and field other words, long wavelength fluctuations become dominant prior to a transition [34]. experiments [25–29], as well as analysis of climatic paleo-records [30], suggest that generic leading indicators (i.e. variance, skewness Two metrics that help characterize spatial patterns for and autocorrelation at-lag-1) may indeed be detected in time series periodicity and directionality can be evaluated from the power of real systems prior to transitions, but the empirical validation of spectrum. First, the radial-spectrum (r-spectrum) is obtained by the spatial indicators remains scarce [15,17,20,21,25,31]. summing the power spectrum at constant distances from the origin of the power spectrum, i.e. along concentric circles at different The discrepancy between theoretical developments and their empirical validation arises from a number of issues, such as the distances from the center. It allows evaluating the periodicity of the patterns. Periodic patterns are characterized by a peak in the lack of sufficiently resolved and long term data, as well as the lack of a coherent methodological framework that outlines the steps power spectrum. The wavenumber at which a peak occurs corresponds to the number of times that a pattern reproduces itself and statistical tools necessary to detect those signals. For indicators based on time series, these issues have been addressed in a recent within a unit area of the spatial data, and therefore contains paper which provides a methodological guide to practitioners and information about the scale of the pattern. Second, the angular- managers to detect early warnings in time series [32] (see also spectrum (h-spectrum) is obtained by summing the values of the http://www.early-warning-signals.org). power spectrum using angular sectors. It allows evaluating the isotropy (orientation) of the spatial pattern. For an isotropic data Here, we complement the previous work on detecting early warnings in time series data by providing a step-by-step set, the h-spectrum will show uniform amplitude at all angles, whereas for an anisotropic data set, the amplitude of the spectrum methodology for detecting early warning signals in spatial data sets. We gather, for the first time, all the early warning signals will show strong amplitudes for specific orientations [35]. Variability based indicators: Spatial variance and spatial proposed in the literature so far in a spatial context. We apply these metrics on model-generated data sets along a degradation skewness. Increased recovery time enroute to a bifurcation point may lead to stronger fluctuations around the equilibrium gradient, and we discuss their interpretation based on a few potential null models. Our analysis mimicks a situation where an state of the system [36]. This can cause spatial variance of the system to increase prior to a transition [10,11]. Spatial variance is ecosystem would be degrading and where we would have access to formally defined as the second moment around the spatial mean of several snapshots of an ecosystem’s spatial structure taken over a the state variable. It has also been shown that the fluctuations period of time, or at different locations along a degradation around the mean can become increasingly asymmetric as the gradient. We hope that our methodology together with the system approaches a bifurcation point. This is because the computer codes will stimulate testing and applications of spatial fluctuations in the direction of the alternative stable state take early warning signals on spatial data (R-code for the spatial longer to return back to the equilibrium than those in the opposite analysis can be found at https://github.com/ direction [11]; this asymmetry can also arise due to local flickering earlywarningtoolbox/spatial_warnings). events (i.e. occasional jumps of local units between their current and alternative state) [37]. The spatial asymmetry can be Methods measured by spatial skewness, which is the third central moment Spatial indicators scaled by the standard deviation. We first give a brief overview of the spatial early warning signals Patch based indicators: shapes and sizes of proposed in the literature so far. More details about the indicators patches. Many ecological systems, such as shrublands in semi- arid ecosystems and mussel beds in the intertidal, exhibit striking and their precise mathematical formulation are provided in Appendix S1. Table 1 summarizes the spatial indicators and their spatial self-organized patterns [38]. It has been suggested that the nature of local ecological interactions, such as the relative scales of expected trends along a degradation gradient. Slowing-down based indicators: spatial correlation and competition and facilitation, can strongly influence the type of spectral properties (DFT). Due to increased recovery time to emerging spatial patterns, leading to i) regular, periodic patches local equilibrium after a perturbation, neighboring units become with a characteristic patch size [22,38–40], or to ii) no more like each other when a system approaches a bifurcation characteristic scale of patchiness [15,39–42]. These different types point, i.e. they become increasingly correlated [5]. The increasing of spatial structures have been observed in a range of ecosystems, spatial coherence can be quantified by the spatial correlation however their use as potential indicators of degradation has mostly function, or Moran’s I, between ecological states separated by a been developed in the case of drylands, where both types of spatial certain distance. The near-neighbor spatial correlation, the analog structures exist. In drylands, it has been shown that the early of autocorrelation at lag 1 for time series, is calculated for the warnings depend on the type of patchiness exhibited by the distance between nearest neighboring units of the system. ecological system [12]. Spatial spectral properties change as the system approaches a In ecosystems exhibiting periodic patterns, as the level of tipping point [13]. To quantify spectral properties, we compute external stress increases, a predictable sequence of self-organized the Discrete Fourier Transform (DFT) that decomposes spatial patterns based on ‘Turing instability’ occurs. In isotropic areas (i.e. data into components of sine and cosine waves of different no preferred orientation of the pattern) the shape of the patterns wavenumbers [33]. A wavenumber can be thought of as a ‘spatial shifts from gaps to labyrinths and to spots as the system becomes frequency’, or the number of times that a pattern is repeated in a more degraded. Thus, spotted vegetation patterns have been unit of spatial length. In a spatial data set, periodicity is visualized proposed to be an early warning signal of imminent desertification PLOS ONE | www.plosone.org 2 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Table 1. Early warning signals of transitions in spatial data. Phenomenon Method/Indicator Expected trend Ref. Rising memory Rising variability Patchiness Spatial correlation x increase [5] Return time x increase [12] Discrete Fourier Transform x spectral reddening [13] Spatial variance x x increase [10,11] Spatial skewness x x peaks (see caption) [11] Patch-size distributions x change in shape of the dist. [15] Regular spotted patterns x change in patch shape [22] Power spectrum x spectral reddening [43] 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 cannot act as a null model for spatial variance or skewness but only anisotropic areas with band-like patterns, the wavelength increases for other metrics such as spatial autocorrelation, DFT, and patch- as the system approaches a transition [43]. size distribution. In contrast to periodic patterns, there are cases where spatial To devise a null model for spatial variance and skewness, Eby, processes give rise to non-periodic (irregular) patterns. In these Guttal and others (unpublished data) propose a coarse-graining cases, we can quantify the size of each patch and calculate the method which should be applied for both the reshuffled matrix frequency of occurrence of different patch sizes. It is common and the real data matrix. In this method, we first divide the full practice to characterize the patchiness of these systems by a matrix of dimension n|n into nonoverlapping submatrices of size function that best describes the distribution of patch sizes. s|s. We then replace each submatrix by its average to obtain a Irregular patterns may be characterized by a scale-free patch-size smaller ‘coarse-grained matrix’ of size c |c (note that c ~n=s). g g g distribution, which means that there is no typical patch size in the The basic intuition behind the method is as follows: consider any ecosystem. Such a distribution may be well approximated by a two non-overlapping submatrices of dimension (e.g. 5|5) from pure power law [44] or by other heavy-tailed functions, such as a the reshuffled matrix. Since the reshuffled matrix is equivalent to a log-normal, a stretched exponential or a power law with cutoff random matrix, the average of the entries of the two sub-matrices [44,45]. Scale-free patch-size distributions have been observed in chosen would be roughly equal to the average of the full matrix. several ecosystems [15,20,41,42]. Computational models of dry- This exercise of ‘coarse-graining’ necessarily reduces variability lands predict that larger vegetation patches become fragmented across submatrices in the case of a reshuffled (thus, random) into smaller ones as aridity or grazing pressure increases and show matrix. Now, consider two non-overlapping submatrices in the an increasing deviation from a theoretical power law as the real data. If we expect that the real data contains a non-random ecosystem approaches the desertification point [15,46]. Therefore, spatial pattern, the average of the entires of the two submatrices it has been hypothesized that an increasing deviation from power- need not be of comparable value to each other nor with the law distribution of patch sizes can signal increasing degradation average of the full real data matrix. Therefore, in contrast to the (but see [17,18]). 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 Statistical significance tests coarse-graining applied to the reshuffled matrix provides a null There are two steps in computing the spatial indicators. The model for spatial variance and spatial skewness. first one is to compute the spatial metrics for a given snapshot. The An alternative method of building null models is to construct a second is to evaluate the trend of the spatial metrics along a spatial matrix from a continuous stochastic process. This method is degradation gradient (see Table 1). In doing so, we need to ensure applicable when continuous data (such as biomass density) is that the spatial metrics for each snapshot and their trends differ available at each spatial point as in data set 1 and 3. More from what would be expected by chance. A standard way to specifically, we construct a null model matrix where each entry is a produce null models is to generate surrogate data and compare the random number (e.g., from a normal distribution) whose mean trends in the indicators obtained from the original data to the and variance are equal to the mean and variance of the original trends obtained from the surrogate data [32]. Here, we discuss data matrix, respectively. This approach provides a null model for ways of obtaining null models for spatial early warning signals. spatial skewness, correlation, and DFT, but not for spatial variance Null models. One way of obtaining a null model is to since variance is, by construction, fixed to be the same as the one randomly permute or shuffle the elements of the spatial matrix, from the original matrix (see figures in Appendix S3). However, it and this is also called bootstrapping. This is equivalent to a may be claimed, following the same arguments as above, that the randomization procedure that removes any spatial structure from coarse-graining method can help estimate statistical significance of the original data but conserves the values of spatial variance and spatial variance. spatial skewness since these moments do not depend on the spatial arrangement of the data points. Therefore, such surrogate data PLOS ONE | www.plosone.org 3 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 The particular case of patchy ecosystems. A system with of data at a given spatial location can be of two types: (a) a discrete patchiness may show characteristic scales in the r- and h-spectra. occupancy data, such as presence or absence of vegetation (or These r- and h-spectra of the data can be compared to those species) at each pixel of an image, or (b) a continuous variable, obtained from a null model. In the case where there is no spatial such as NDVI (Normalized Difference Vegetation Index) at each structure, it is known analytically that the values of the scaled pixel or nutrient concentration at each sampling point. power spectrum observed in each bin when calculating the r-or h- Here, to serve better our method-illustration purposes, we chose to work on model-generated data rather than real data. We 2k spectrum should be distributed as with k the number of values thereby circumvent limitations of missing or noisy data, and avoid 2k in the bin (see for instance [35,47]; this gives very similar results as issues of misinterpretation arising from potential insufficient the confidence intervals presented in Appendix S3). knowledge about the underlying degradation gradients in real Once the patches are identified and their size evaluated, various data. We generated three synthetic data sets using three representative models of tipping points and self-organization in 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 ecological systems. The three models treat ecological variables in a fitted to the patch-size distribution. One way of fitting a given spatially-explicit framework with stochasticity. They all describe vegetation dynamics under resource limitation or grazing pressure, distribution to data has been to use ordinary least square regression on the log-log transformed probability distribution but they differ in the nature of ecological interactions and the function of patch sizes. However, this method is known to have emerging spatial vegetation structure. A detailed description of the substantial bias in estimating the parameter values, especially for models can be found in Appendix S2 but a brief description small data sets [44,45]. Maximum likelihood methods or least- follows. square fits of the inverse-cumulative distribution (which quantifies Data set 1 was obtained from a local positive feedback model the number of patches whose size is larger than a given value s for N resulting in a non-patchy vegetation structure [50,51] (Fig. 1 different values of s) provide more accurate estimates of the first row). Space is represented as a two-dimensional lattice parameters of most heavy-tailed functions [48]. [52,53]. Locally, vegetation density grows logistically and is Trends. The above methods inform us about whether the lost due to grazing. Biomass and water are exchanged between spatial indicators for a given spatial data set are significantly neighboring sites at a certain rate, such that a site with high different from those of random patterns. However, to anticipate biomass (or water) will have the tendency to diffuse biomass (or ecological transitions, we also need to know how these spatial water) to its neighboring sites. As rainfall falls below a certain indicators are changing along a degradation gradient. Trend threshold, the ecosystem undergoes an abrupt transition from a statistics like the Kendall’s t [12,30] or Pearson’s correlation globally high vegetation density to a bare state due to the coefficient [25,49] can be used to quantify the strength of the trend nearly synchronous shifts of each of the sites to a desert state in the indicator along the gradient. [52]. Data set 2 is based on a local facilitation model that exhibits Simulated spatial data sets N spatial patterns characterized by a scale-free patch-size Spatial data in ecology are typically obtained by field studies, distribution [15]. In this stochastic cellular automaton model, data collecting devices placed at various locations of an ecosystem an ecosystem is represented by a grid of cells, each of which or extracted from spatial imagery. In any of these cases, the nature PLOS ONE | www.plosone.org 4 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Figure 2. Flow chart of analysis to perform on a spatial data set. doi:10.1371/journal.pone.0092097.g002 can be in one of three discrete states: vegetated (+), empty (o) these models. In a previous study [12], it has been shown that the three systems take increasingly longer to recover to their or degraded (2). Empty cells represent fertile soil whereas degraded cells represented eroded soil patches unsuitable for equilibrium after perturbation, thus demonstrating that critical slowing down is a generic feature of the transitions observed in recolonization by vegetation. A key ecological mechanism is the positive effect of vegetation on its local neighborhood these three ecological models, regardless of their different underlying mechanisms and their different types of spatial through increased regeneration of degraded cells. Because of structures. this local facilitation, vegetated cells tend to form clusters (Fig. 1 For our analysis, we selected ten snapshots (i.e. two-dimentional second row). When the environmental conditions become space discretized into matrices recording the spatial spread of the harsher, there is a point at which the vegetation dies out and vegetation at the end of the simulation) for each of these models at the system becomes a desert resembling a saddle-node (or fold) different points along a gradient of degrading conditions prior to bifurcation. the transition. We illustrate our analyses using only ten of these Data set 3 is based on a scale-dependent feedback model that points which are not equally spaced along the gradient (their results in periodic (Turing-like) spatial vegetation patterns [54]. location is shown on Fig. S1 in Appendix S2). This model is based on a three partial differential equations We are interested in quantifying how the spatial characteristics model describing the dynamics of vegetation biomass, soil of these matrices change when approaching a tipping, or water and surface water. Plants grow due to soil water bifurcation, point. We are considering cases where the whole availability and die due to natural mortality and/or grazing. ecosystem shifts to an alternative state. In our mathematical The infiltration rate of water in the soil is higher in areas with representation of the ecosystem, this is equivalent to the entire vegetation than in bare soil, leading to the accumulation of matrix undergoing a shift. The degradation sequence of the water under vegetation and to its depletion further away, matrices might correspond to snapshots in time (e.g. temperature resulting in a scale-dependent feedback responsible for the changing through time) or in space (e.g. herbivory pressure formation of regular vegetation patterns [54] (Fig. 1 last row). changing in space depending on a distance to a water point). Both When water availability becomes limited, a homogeneous types of data are relevant to evaluate and test early warning vegetated state becomes unstable leading to self-organized signals. However, shifts of a given spatial ecosystem in time are patterns such as gaps, labyrinths and spots. A further reduction more commonly the type of phenomena that we are trying to in water availability leads to a transition into a desert state, anticipate. again mimicking a fold-like bifurcation. Results The three models exhibit a bifurcation from one state (e.g., vegetated) to an alternative state (e.g., desert) as an external We suggest a step by step process to decide which spatial parameter (such as rainfall, grazing, etc) changes. See Appendix S2 indicators should be used (Fig. 2). The spatial statistics that need to for underlying mathematical equations and parameter values of be evaluated depend on the type of data set (with discrete or PLOS ONE | www.plosone.org 5 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 565 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 continuous values) at hand. Two of the three datasets used in this one vegetated and one bare), the patch-size distribution may be paper provide quantitative, continuous data, i.e. vegetation plotted and estimated [15,46]. If the patterns are periodic and anisotropic, the wavelength of the pattern should be evaluated. biomass (data set 1 and 3), whereas data set 2 gives qualitative, discrete data indicating the presence or absence of vegetation in The wavelength is equivalent to the dominant length scale of the pattern provided by the r-spectrum [43]. For periodic isotropic each cell. For this latter data set, we transformed the original matrix by the coarse-graining procedure described in ‘‘Stastistical 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 significance tests’’. More specifically, we used submatrices of 565 of patterns (i.e. spots, gaps or labyrinths) [43]. Additionally, it is cells in which we counted the number of cells occupied by noteworthy that if the data set includes several replicates at each vegetation [12]. We obtained matrices that were 25 times smaller stress level, potential analysis may be performed [55–57] (see more than the original ones and where values in each of the cells ranged details in Appendix S1). between 0 and 25, indicating the local abundance of vegetation. Next, we present how this methodology can be applied to our The size of the submatrix used to transform the original data may affect the behavior of the indicators. three data sets. The first question to ask is whether the patterns observed are periodic or not. The r-spectrum obtained from the DFT analysis 1. Distinguishing periodic from non-periodic patterns provides information about whether the patterns are periodic, We used DFT analysis to estimate the r-spectra as a function of while the h-spectrum indicates whether the patterns are isotropic wavenumbers for all the three data sets (Fig. 3). The first and (i.e. with no specific orientation) or anisotropic (i.e. with a specific second data sets (Fig. 3, first and second row) show a noisy pattern orientation, e.g. band-like patterns). If the patterns are not indicating that contribution to r-spectra is not significant for all periodic, the generic leading indicators (i.e. spatial variance, wavenumbers. However, the r-spectrum for the third data set, the spatial skewness and spatial correlation between nearby sites) may scale-dependent feedback model, shows a clear peak (Fig. 3 last be used [12] and the power spectrum should be checked for row) even far from the transition. The peak indicates that there is possible reddening [13]. In addition, if the patterns are not only dominant wavelength (corresponding to a characteristic patch size) irregular but also patchy (e.g. can be characterized by two phases, which is a signature of periodic patterns. Note that periodicity can PLOS ONE | www.plosone.org 6 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 Figure 4. Generic leading indicators in data sets 1 and 2 along system approaches the bifurcation point (Fig. 6 top row) [46]. a degradation gradient. In each panel, the x-values correspond to the rank of the snapshot of the system along the degradation gradient. For comparison, we plotted the inverse cumulative patch-size This mimicks a scenario where we would not know the exact value of distribution of data set 3 that is also patchy but periodic. Far from the driver but where we can order the data set along a degradation the transition, after the onset of pattern formation, the periodic gradient and see Fig. S1 in Appendix S2 to visualize the location of the patterns presented a patch-size distribution characterized by a four snapshots along the degradation gradient. For each of the three data sets, 10 snapshots were used. Left: local positive feedback model sharp cutoff (Fig. 6 bottom row). As the system approaches the (data set 1). Right: local facilitation model (data set 2; original data bifurcation point, the value of the cutoff decreases indicating transformed using 565 sub-matrices). First row: spatial variance. Second decreasing patch size in the periodic pattern. row: spatial skewness. Third row: spatial correlation at lag one. In each panel, Kendall’s t, quantifying the trend of the indicator, is indicated. 3. Probing spatial early warnings for periodic patterns Gray areas correspond to 95% confidence intervals obtained using 200 simulations of a null model (i.e. data sets of same size generated by In contrast to the first two data sets, data set 3 exhibits periodic reshuffling the original data set). patterns (Fig. 1 and 3 last rows). The h-spectra does not indicate a doi:10.1371/journal.pone.0092097.g004 strong and clear signal at any specific angle, suggesting that the patterns do not have a clear orientation, i.e. patterns are isotropic also be seen by plotting the power spectrum, where the periodicity (Fig. 7 first row). When the system approaches the bifurcation is visible as a ring corresponding to the dominant wavelength of point (from left to right panel on Fig. 7 second row), the the data set (see upcoming paragraph ‘‘DFT and reddening’’). distribution of values of the data set goes from one peak reflecting Therefore, we conclude that the two first data sets are not periodic the absence of patterns, to a two-peak distribution due to the whereas the last one shows some spatial periodicity. occurrence of both vegetation and bare soil in the system after the emergence of spatial patterns, and finally to a distribution that is 2. Probing spatial early warnings for non-periodic skewed toward small values because of the dominance of bare soil in the system. The last distribution observed before the bifurcation patterns point characterizes spot patterns [43] which has been hypothe- We first focus on the case of the non-periodic patterns, i.e. data sized to be an indicator of imminent desertification [22]. set 1 showing no clear spatial structure and data set 2 with vegetation clusters. PLOS ONE | www.plosone.org 7 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Figure 5. Power spectrum along a degradation gradient in the three data sets. In a row, the system approaches the bifurcation point, from M N left to right column (as in Fig. 1). For a data set of size M|N, the power spectrum is typically plotted for wavenumbers up to p~ and q~ [47] 2 2 and is scaled by the spatial variance s (i.e. the scaled power spectrum is evaluated as ) [35]. Red color indicates higher values of the scaled power spectrum, . 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 565 sub-matrices). Bottom row: scale-dependent feedback model (data set 3). doi:10.1371/journal.pone.0092097.g005 transitions. To demonstrate the methods, we employed data 4. Probing statistical significance generated from simulations of spatially-explicit ecological models We make further comments on the statistical significance for with stochasticity that showed abrupt transitions from one state to indicators of spatial data, especially in the context of null models an alternative state. It is increasingly being recognized that spatial for spatial variance and spatial skewness. As an illustration of the dynamics pose challenges and provide opportunities for both basic coarse-graining method proposed in the ‘‘Methods’’ section, we science as well as management [58]. Research on spatial dynamics compared the trends in the generic leading indicators on coarse- in ecology is continually uncovering new patterns and mecha- grained matrix obtained by different dimensions of submatrix, nisms. Some of these processes are likely related to regime shifts. In s|s with s~1 (in this case the coarse-grained matrix is identical to this context, the main objective of this manuscript was to provide a the original matrix), 2 and 5 but starting from the same original methodological guide that can stimulate the application of spatial matrix of size 100|100 from our data set 2. As shown in the first indicators for ecological transitions on empirical data sets from real row of Fig. 8, the values of spatial variance along the degradation case studies. gradient differ substantially from the ones of the null model only in In recent years, much effort has been devoted to the search for the coarse-grained case (second and third column). The same ‘generic’ indicators based on the idea that there are some common result seems to be true, although the differences are not as behaviors across a range of complex systems as they approach a pronounced, for spatial skewness (second row). In contrast, the bifurcation point [59]. Taking into account the spatial organiza- coarse-graining method does not offer a good null model for tion of natural systems has revealed that many indicators may not spatial correlation (last row, Fig. 8). In summary, a random matrix behave in spatially-structured systems as they would in other that has the same dimensions and average as the original spatial systems. More specifically, recent theoretical studies have suggest- data can act as a null model for computing spatial correlation only. ed that trends of generic leading indicators (spatial variance, On the other hand, a comparison between indicators of coarse- spatial skewness, and spatial correlation between nearby sites) grained matrices of both original and random matrix data can act could be different in self organized patterned systems [12]. For as a null model for spatial variance and spatial skewness. such ecosystems, system-specific indicators may be more appro- priate. In particular, when spatial patterns are periodic or regular Discussion (Turing-like), the shape of the patterns may give an idea of the In this manuscript, we presented a systematic methodology for proximity to the threshold where the system may undergo a applying spatial early warning signals of abrupt ecological regime shift [22]; specifically, spots could warn of approaching PLOS ONE | www.plosone.org 8 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 colum (as in Fig. 1)n. First row: local facilitation model (original data not transformed). Bottom row: scale-dependent feedback model. Gray areas corresponds 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 desertification [16,22]. If the patterns are non-periodic or irregular Clearly, the usefulness of the indicators presented here depends (in particular, if the patch-size distribution is described by a heavy- greatly on the underlying mechanisms driving the change in the ecosystems studied. In all our analyses, we assumed that a single tailed distribution), the patch size distribution may contain information about the degradation level of the ecosystem but driver is monotonically changing, while the underlying environ- mental conditions are assumed relatively stable and the environ- more needs to be known about the underlying ecological mental stochasticity relatively weak. This may not always be the mechanisms to interpret the changes in the shape of this distribution [15,20,46,60]. In other words, the nature of spatial case. For example, spatial correlation can increase due to changes in the underlying spatial heterogeneity of the environment, or due organization of ecosystems seems to be a key factor in determining to alterations in local ‘mixing’ (or diffusion) levels in the landscape what type of indicators may be employed to detect an impending [5]. Effects of increasing spatial variance near a critical point could ecological transition in spatially-structured systems. Therefore, a be confounded by various intrinsic factors such as demographic first and essential step when starting to analyze a spatial data set is noise arising from changes in population sizes, state-dependent/ to get an idea of the type of spatial organization that one is dealing multiplicative noise. The nature of the dispersal processes between with. A good knowledge of the system and its underlying ecological patches in fragmented landscapes may affect expected trends in mechanisms (specifically those responsible for the spatial structure) the indicators. In addition, extrinsic factors such as environmental are required to know which indicator to use and how to interpret fluctuations that vary in space and time can complicate our the changes. Theoretical studies have started developing method- interpretations of early warning signals. In a similar way, system- ologies for inferring underlying mechanisms from a limited specific conditions may affect the behavior and thereby interpre- number of spatial snapshots [61,62]. Such knowledge will facilitate tation of the indicators. For example, in Mediterranean drylands, to assess the risks of ecological transitions while accounting for degradation is accompanied by a change in the patch-size potential false and failed alarms. distribution toward less large patches, whereas the opposite was PLOS ONE | www.plosone.org 9 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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: h-spectrum. Gray areas corresponds 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 observed in salt marshes because the underlying mechanisms have been extensively studied, may be a fertile ground for further driving the formation of the patterns differ in the two ecosystems research. [15,46,60]. Therefore, our interpretations of leading indicators are More importantly, despite a few recent studies, we still lack prone to both false positives and false negatives. empirical tests of spatial early warnings. Researchers have studied To prevent such errors, knowledge of the underlying heteroge- spatial warnings in laboratory populations of microbial organisms neity, repeated observations of the system [60], or a sufficiently such as Daphnia and yeast [25,31]. In these studies, individual well-described system so that it can be modelled are necessary to populations are maintained in locally well-mixed small beakers or know what to expect along a degradation trend. As this is not petri-dishes and they are ‘connected’ to other populations by usually the case, general models have been suggested to be fitted to ‘controlled dispersal’ where the researcher transfers a fraction of time series data[63,64] to simulate surrogate scenarios for the local population to its nearest neighbors. In the field, comparing trends obtained from the original data sets. In this researchers have employed space-for-time substitution; in this case the model provides an expectation of i) whether and when a approach, it is assumed that spatial patterns at locations with shift is likely, and ii) what trends should look like as the system is different values of stressors (e.g. grazing or rainfall) are equivalent approaching a shift. Such approaches, however, are yet to be to the dynamics of spatial patterns where the stressor is changing developed for spatially-explicit systems. Although we presented a with time. This is a widely used approach in ecology as an couple of ways to develop null models, these are rather simplistic alternative to long term ecological studies, for example to as they entirely neglect any underlying spatial structuring. investigate ecological succession or how ecosystems may respond Therefore, the design and selection of null models for spatial to climate change [67]. In the context of alternative stable states data, which is an area of active research in ecology, is another and tipping points, this method has been employed to establish the important avenue for further research [65]. There are other existence of (multiple) stable states in savanna ecosystems as a promising avenues of research to be pursued. Composite metrics function of rainfall [57,68], and to forecast how spatial self- that combine spatial patterns with their temporal dynamics could organization of semi-arid vegetation may respond to increasing potentially offer new and potentially more reliable indicators of stressors such as grazing [15]. However, the application of space- imminent transitions [66]. Ecological transitions and leading for-time method is not without limitations [69]. One needs to be indicators in the context of metapopulation dynamics, where cautious about the possibility of existence of various other sources factors that stabilize or make species more vulnerable to extinction of heterogeneity, both biotic and abiotic, along a spatial gradient of stressor. As we have argued before, the interpretation of the PLOS ONE | www.plosone.org 10 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 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 mimicks 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 262 sub-matrices. Middle: original data transformed using 565 sub-matrices. First row: spatial variance. Second row: spatial skewness. Third row: spatial correlation at lag one. In each panel, Kendall’s t 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 trends of early warnings are crucially dependent on the nature of have been used in models, experiments or retroactively). We hope local ecological interactions and the patterns they produce. that this work will stimulate further development, testing and application of spatial indicators in a broad range of ecosystems. 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 Supporting Information space-for-time substitutions. Appendix S1 Spatial indicators. In conclusion, both the theory and the application of the spatial (PDF) indicators lag behind the development of the temporal ones. Spatial patterns may however offer advantages for anticipating or Appendix S2 The three models used to generate the test detecting ecological transitions of the types studied here [11,32]. data sets. Unlike temporal indicators which require long, unbroken time (PDF) series of frequent observations, spatial indicators may be evaluated Appendix S3 Results compared to the null model based even if measurements are irregular or infrequent over time. While on white noise. spatial pattern measurements require intensive data collection at (PDF) each time point, in many cases this may be easier than high- frequency time series sampling [5,11,13]. In both space and time, Acknowledgments more empirical validation of the indicators currently proposed in the literature is needed. We do not yet have an example where This paper was based on the workshop ‘Practical Methods for Analysis of early warning signals were used to avert an upcoming shift (they Early Warnings for Regime Shifts’ organized by VD, SRC and MS in the PLOS ONE | www.plosone.org 11 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings Santa Fe Institute, in Santa Fe, New Mexico, USA, from 10 to 12 of Author Contributions October 2011. We thank Vincent Deblauwe for useful comments on Conceived and designed the experiments: SK VG VD SRC DAS VNL. previous versions of the manuscript and for sharing his r - and h - spectrum Performed the experiments: SK VG. Analyzed the data: SK VG. code. Contributed reagents/materials/analysis tools: SK VG VD. Wrote the paper: SK VG WAB SRC AME VNL DAS MS EHvN VD. References 1. Scheffer M, Carpenter SR, Foley JA, Folke C, Walker B (2001) Catastrophic 32. Dakos V, Carpenter SR, Ellison AM, Guttal V, Ives AR, et al. (2012) Early shifts in ecosystems. Nature 413: 591–596. warning signals of critical transitions: Methods for time series. PLOS ONE 7: 2. Scheffer M, Bascompte J, Brock W, Brovkin V, Carpenter S, et al. (2009) Early- e41010. warning signals for critical transitions. Nature 461: 53–59. 33. Legendre P, Legendre L (1998) Numerical ecology. Amsterdam: Elsevier 3. Wissel C (1984) A universal law of the characteristic return time near thresholds. Science. Oecologia 65: 101–107. 34. Guttal V (2008) Applications of nonequilibrium statistical physics to ecological 4. van Nes EH, Scheffer M (2007) Slow recovery from perturbations as a generic systems (phd thesis). The Ohio State University 3: 39–41. indicator of a nearby catastrophic regime shift. American Naturalist 169: 738– 35. Couteron P (2002) Quantifying change in patterned semi-arid vegetation by 747. fourier analysis of digitized aerial photographs. International Journal of Remote 5. Dakos V, van Nes E, Donangelo R, Fort H, Scheffer M (2010) Spatial Sensing 23: 3407–3425. correlation as leading indicator of catastrophic shifts. Theoretical Ecology 3: 36. Sethna JP (2007) Statistical Mechanics: Entropy, Order Parameters, and 163–174. Complexity. Oxford: Oxford University Press. 6. Held H, Kleinen T (2004) Detection of climate system bifurcations by 37. Dakos V, van Nes EH, Scheffer M (2013) Flickering as an early warning signal. degenerate fingerprinting. Geophysical Research Letters 31: L020972. Theoretical Ecology : DOI:10.1007/s12080-013-0186-4. 7. Carpenter SR, Brock WA (2006) Rising variance: a leading indicator of 38. Rietkerk M, van de Koppel J (2008) Regular pattern formation in real ecological transition. Ecology Letters 9: 308–315. ecosystems. Trends in Ecology & Evolution 23: 169–175. 8. Guttal V, Jayaprakash C (2008) Changing skewness: an early warning signal of 39. Manor A, Shnerb M (2008) Facilitation, competition, and vegetation patchiness: regime shifts in ecological systems. Ecology Letters 11: 450–460. From scale free distribution to patterns. Journal of Theoretical Biology 253: 9. Ke´fi S, Dakos V, Scheffer M, Van Nes E, Rietkerk M (2013) Early warning 838–842. signals also precede non-catastrophic transitions. Oikos 122: 641–648. 40. von Hardenberg J, Kletter A, Yizhaq H, Nathan J, Meron E (2010) Periodic 10. Obo´rny B, Mesze´na G, Szabo´ G (2005) Dynamics of populations on the verge of versus scale-free patterns in dryland vegetation. Proceedings of the Royal Society extinction. Oikos 109: 291–296. B: Biological Sciences 2077: 1771–1776. 11. Guttal V, Jayaprakash C (2009) Spatial variance and spatial skewness: leading 41. Scanlon T, Caylor K, Levin S, Rodriguez-Iturbe I (2007) Positive feedbacks indicators of regime shifts in spatial ecological systems. Theoretical Ecology 2: 3– promote power-law clustering of Kalahari vegetation. Nature 449: 209–212. 12. 42. Foti R, del Jesus M, Rinaldo A, Rodriguez-Iturbe I (2013) Signs of critical 12. Dakos V, Ke´fi S, Rietkerk M, van Nes E, Scheffer M (2011) Slowing down in transition in the everglades wetlands in response to climate and anthropogenic spatially patterned ecosystems at the brink of collapse. The American Naturalist changes. Proceedings of the National Academy of Sciences 110: 6296–6300. 177: E153–E166. 43. Deblauwe V, Couteron P, Lejeune O, Bogaert J, Barbier N (2011) 13. Carpenter S, Brock W (2010) Early warnings of regime shifts in spatial dynamics Environmental modulation of self-organized periodic vegetation patterns in using the discrete fourier transform. Ecosphere 1: art10. sudan. Ecography 34: 990–1001. 44. Newman M (2005) Power laws, pareto distributions and zipf’s law. Contempo- 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 rary physics 46: 323–351. ecosystems. Journal of Ecology 94: 537–547. 45. White E, Enquist B, Green JL (2008) On estimating the exponent of power-law 15. Ke´fi S, Rietkerk M, Alados CL, Pueyo Y, Papanastasis VP, et al. (2007) Spatial frequency distributions. Ecology 89: 905–912. vegetation patterns and imminent desertification in mediterranean arid 46. Ke´fi S, Rietkerk M, Roy M, Franc A, De Ruiter P, et al. (2011) Robust scaling in ecosystems. Nature 449: 213–217. ecosystems and the meltdown of patch size distributions before extinction. Ecology letters 14: 29–35. 16. Ke´fi S, Eppinga M, de Ruiter P, Rietkerk M (2010) Bistability and regular spatial patterns in arid ecosystems. Theoretical Ecology 3: 257–269. 47. Mugglestone MA, Renshaw E (1998) Detection of geological lineations on aerial 17. Maestre F, Escudero A (2010) Is the patch size distribution of vegetation a photographs using two-dimensional spectral analysis. Computers and geosci- suitable indicator of desertification processes? reply. Ecology 91: 3742–3745. ences 24: 771–784. 18. Ke´fi S, Alados C, Chaves R, Pueyo Y, Rietkerk M (2010) Is the patch size 48. Clauset A, Shalizi CR, Newman ME (2009) Power-law distributions in empirical distribution of vegetation a suitable indicator of desertification processes? data. SIAM review 51: 661–703. comment. Ecology 91: 3739–3742. 49. Seekell D, Carpenter S, Cline T, Pace M (2012) Conditional heteroskedasticity 19. Pueyo S (2011) Desertification and power laws. Landscape ecology 26: 305–309. forecasts regime shift in a whole-ecosystem experiment. Ecosystems 15: 741– 20. Lin Y, Han G, Zhao M, Chang S (2010) Spatial vegetation patterns as early 747. signs of desertification: a case study of a desert steppe in inner mongolia, china. 50. Shnerb NM, Sarah P, Lavee H, Solomon S (2003) Reactive glass and vegetation Landscape ecology 25: 1519–1527. patterns. Physical Review Letters 90: 038101. 21. Litzow M, Urban J, Laurel B (2008) Increased spatial variance accompanies 51. Guttal V, Jayaprakash C (2007) Impact of noise on bistable ecological systems. reorganization of two continental shelf ecosystems. Ecological applications 18: Ecological Mod-elling 201: 420–428. 1331–1337. 52. van Nes E, Scheffer M (2005) Implications of spatial heterogeneity for 22. Rietkerk M, Dekker SC, de Ruiter PC, van de Koppel J (2004) Self-organized catastrophic regime shifts in ecosystems. Ecology 86: 1797–1807. patchiness and catastrophic regime shifts in ecosystems. Science 305: 1926– 53. Keitt T, Lewis M, Holt R (2001) Allee effects, invasion pinning, and species 1929. borders. the American Naturalist 157: 203–216. 23. Hastings A, Wysham D (2010) Regime shifts in ecological systems can occur with 54. Rietkerk M, Boerlijst MC, van Langevelde F, HilleRisLambers R, van de no warning. Ecology Letters 13: 464–472. Koppel J, et al. (2002) Self-organization of vegetation in arid ecosystems. 24. Boettiger C, Hastings A (2013) Tipping points: From patterns to predictions. American Naturalist 160: 524–530. Nature 493: 157–158. 55. Livina V, Kwasniok F, Lenton T (2010) Potential analysis reveals changing 25. Drake J, Griffen B (2010) Early warning signals of extinction in deteriorating number of climate states during the last 60 kyr. Climate of the Past 6: 77–82. environments. Nature 467: 456–459. 56. Livina V, Kwasniok F, Lohmann G, Kantelhardt J, Lenton T (2011) Changing 26. Carpenter S, Cole J, Pace M, Batt R, Brock W, et al. (2011) Early warnings of climate states and stability: from pliocene to present. Climate dynamics 37: regime shifts: A whole-ecosystem experiment. Science 332: 1079–1082. 2437–2453. 27. Veerat AJ, Faassen EJ, Dakos V, van Nes EH, Lurling M, et al. (2012) Recovery 57. Hirota M, Holmgren M, Nes EV, Scheffer M (2011) Global resilience of tropical rates reflect distance to a tipping point in a living system. Nature 481: 357–359. forest and savanna to critical transitions. Science 334: 232–235. 28. Dai L, Vorselen D, Korolev K, Gore J (2012) Generic indicators for loss of 58. Walters C, Pauly D, Christensen V (1999) Ecospace: Prediction of mesoscale resilience before a tipping point leading to population collapse. Science 336: spatial patterns in trophic relationships of exploited ecosystems, with emphasis 1175–1177. on the impacts of marine protected areas. Ecosystems 2: 539–554. 29. Wang R, Dearing J, Langdon P, Zhang E, Yang X, et al. (2012) Flickering gives 59. Gilmore R (1981) Catastrophe Theory for Scientists and Engineers. New York: early warning signals of a critical transition to a eutrophic lake state. Nature 492: John Wiley and Sons. 419–422. 60. Weerman E, Van Belzen J, Rietkerk M, Temmerman S, Ke´fi S, et al. (2012) 30. Dakos V, Scheffer M, Van Nes E, Brovkin V, Petoukhov V, et al. (2008) Slowing Changes in diatom patch-size distribution and degradation in spatially self- down as an early warning signal for abrupt climate change. Proceedings of the organized intertidal mudflat ecosystem. Ecology 93: 608–618. National Academy of Sciences 105: 14308–14312. 61. Marcos-Nikolaus P, Martin-Gonza´lez J, Sole´ R (2002) Spatial forecasting: 31. Dai L, Korolev K, Gore J (2013) Slower recovery in space before collapse of Detecting determinism from single snapshots. International Journal of connected populations. Nature 496: 355–358. Bifurcation and Chaos 12: 369–376. PLOS ONE | www.plosone.org 12 March 2014 | Volume 9 | Issue 3 | e92097 Detecting Spatial Early Warnings 62. Mocenni C, Facchini A, Vicino A (2010) Identifying the dynamics of complex 66. Corrado R, Cherubini AM, Pennetta C (2013) Signals of critical transitions in ecosystems associated with fluctuations of spatial patterns. In: Noise and spatio-temporal systems by spatial recurrence properties. Proceedings of the Fluctuations (ICNF), 2013 22nd International Conference on. IEEE, pp. 1–4. National Academy of Sciences of the United States of America 107: 8097–8102. 67. Blois JL, Williams JW, Fitzpatrick MC, Jackson ST, Ferrier S (2013) Space can 63. Boettiger C, Hastings A (2012) Quantifying limits to detection of early warning substitute for time in predicting climate-change effects on biodiversity. for critical transitions. Journal of the Royal Society Interface 9: 1–29. Proceedings of the National Academy of Sciences 110: 9374–9379. 64. Lade S, Gross T (2012) Early warning signals for critical transitions: a 68. Staver AC, Archibald S, Levin SA (2011) The global extent and determinants of generalized modeling approach. Plos computational biology 8: e1002360. savanna and forest as alternative biome states. Science 334: 230–232. 65. Fortin M, Dale M (2005) Spatial analysis: a guide for ecologists. Cambridge, 69. Johnson EA, Miyanishi K (2008) Testing the assumptions of chronosequences in UK: Cambridge University Press. succession. Ecology Letters 11: 419–431. PLOS ONE | www.plosone.org 13 March 2014 | Volume 9 | Issue 3 | e92097

Journal

PLoS ONEPubmed Central

Published: Mar 21, 2014

References