Early warnings of regime shifts in spatial dynamics using the discrete Fourier transform

Abstract

In ecosystems with spatial dispersion of nutrients, organic matter, or organisms, dispersal may dampen variability that could provide an early warning of regime shifts. The discrete Fourier transform (DFT) of spatial pattern rescales system dynamics into spatial frequencies where early warnings may be sharpened. We analyzed four spatial ecological models with the DFT. Different models represented single species in discrete and continuous time, and two different prey‐harvester systems. In all four systems, the DFT of transient data exhibited substantial increases prior to the critical transitions. The DFT adds to the arsenal of early warning indicators for spatially structured ecosystems. In addition it provides information about the spatial frequencies where destabilization first begins.

Introduction

Ecological regime shifts caused by critical transitions in nonlinear dynamics can be announced in advance by early warnings (Scheffer et al. 2009). As the transition point is approached, time series statistics exhibit systematic changes which include slowing down of the return rate after disturbance (van Nes and Scheffer 2007), increasing variance (Carpenter and Brock 2006), shift of variance spectra toward lower frequencies (Kleinen et al. 2003), and increasing skewness (Guttal and Jayaprakash 2008). All of these indicators are computed from time series observations, ideally with rather large sample sizes.

Spatial dynamics in ecology also undergo regime shifts (Sole et al. 1996, Peters et al. 2004, Rietkerk et al. 2004, van Nes and Scheffer 2005, Kefi et al. 2007). Indicators of spatial regime shifts are computed from measurements of spatial pattern, e.g., line transects through vegetation or two‐dimensional maps (Wardwell and Allen 2009). Several early warning indicators have been described for spatial ecological data. As grazing increases toward a critical point in arid grasslands, scale‐invariant power relationships become distorted and vanish with the transition to desert (Kefi et al. 2007). Other early warning indicators of spatial regime shifts include spatial variance (Oborny et al. 2005), skewness (Guttal and Jayaprakash 2009) and correlation (Dakos et al. 2010). Exploited fish populations exhibit reduced spatial heterogeneity and stronger response to climate fluctuations than unexploited species, and this difference may be related to changes in stability (Hsieh et al. 2008). It is important to emphasize that the spatial statistics are not computed by time series methods. They are computed for spatial patterns at a series of discrete times. Spectral analyses of time series (Kleinen et al. 2003) analyze patterns in the time domain whereas the analyses in this paper address the space domain.

In ecosystems subject to multiple regime shifts, early warnings can be muffled or magnified by interactions (Brock and Carpenter 2010). Muffling (magnifying) is a decrease (increase) of variance that diminishes (enlarges) the increase in variance that would be expected as a regime shift is approached. Muffling or magnifying can be caused by interactions or correlated responses to environmental shocks among components of the ecosystem. In spatial ecosystem dynamics, dispersal of populations, nutrients or organic matter is likely to cause muffling of variance because sharp changes in one patch could be compensated by flows to or from from neighboring patches. In systems subject to alternate states, flows from neighboring patches could either cause regime shifts or buffer local variation.

The discrete Fourier transform (DFT) of spatial pattern may be a sensitive indicator of impending regime shift in spatial ecological systems (Appendix). As a regime shift is approached, the DFT is likely to increase sharply at certain spatial frequencies whereas other indicators may be muffled by dispersal. The DFT is used in the context of optimal control problems (Bamieh et al. 2002), for example to study the stability of infinite horizon optimally controlled spatial systems in ecology and economics (Brock and Xepapadeas 2008, 2010). However these papers do not consider early warnings of critical transitions.

Here we investigate the use of the DFT as an early warning indicator for regime shifts in four spatial ecological systems. These include discrete‐time dynamics of a population subject to period‐doubling bifurcations and continuous‐time dynamics of a prey population subject to alternate states. We also consider two harvested systems, one where both prey and harvester disperse (as in a marine fishery) and one where harvesters disperse among patches but prey do not (as in a fishery for many lakes over a landscape). We show that early warning signals of impending bifurcations are more clearly evident in the DFT domain than in the spatial domain. As far as we know, this is a new contribution to the rapidly emerging literature on early warnings of regime shifts.

Methods

We studied transient dynamics of four different spatial ecological models (Table 1). The models represent different types of spatial dynamics that have been used to represent prey populations subject to harvest. Models 1 and 2 represent single prey populations in discrete and continuous time, respectively. In Model 1, the prey undergo a period‐doubling bifurcation when the growth parameter a exceeds a critical value. In Model 2, the prey undergo a fold bifurcation, collapsing from a high‐biomass state to a low‐biomass state, when the parameter H (representing number of harvesters) exceeds a critical value. Models 3 and 4 build on Model 2 by adding dynamics of human harvesters. In Model 3 both prey and the human harvesters are mobile in space. Prey are capable of schooling behavior, and harvesters can adjust to movements of prey. In each spatial cell, prey can collapse from the high state to the low state if the number of harvesters exceeds a critical value. The collapse can be buffered by movement of prey from other cells, or intensified by aggregation of the harvesters in patches of high prey density. In Model 4, prey dynamics are local with no dispersion among patches, and only human harvesters are mobile. The dynamics of humans follow a model used in economics to represent dynamics of cities. In this model, humans tend to form aggregations that move slowly in space. As humans aggregate to high densities, prey collapse from the high state to the low state where the number of harvesters exceeds a critical value. Human harvesters can then move to nearby patches where prey are more abundant, perhaps leading to a spreading pattern of prey collapse. Further details of all models are provided below.

Models analyzed in this study. In all models x is the prey population, H is the harvester or predator population, t is time, and z is the spatial dimension.

For each model, we gradually changed a parameter over time until a critical transition point was reached, and studied the behavior of the system and its discrete Fourier transform as the critical transition was approached. Each model was analyzed on a one‐dimensional spatial coordinate system. The spatial coordinates were placed on a ring, or clockface, of N positions with arithmetic modulo N. For simulations reported here N = 128. Experiments with N = 32 and N = 512 yielded similar results. Dispersal terms in the spatial dimension z of the form D(∂²x/∂z²) were solved numerically following equation (A.2) in the Appendix (Levin 1976). All simulations were computed in R 〈http://www.r‐project.org/〉 by S.R.C. At each time step, discrete Fourier transforms (DFT) for the one‐dimensional spatial gradient were computed using the ‘fft' function and smoothed with a cubic spline.

Discrete time: Model 1

Model 1 represents discrete‐time dynamics of a prey population subject to density‐dependence and dispersal after reproduction. Parameters are a growth parameter a, a survival parameter α, dispersal parameter D, and shock standard deviation σ. Shocks are independent in space and time.

As the growth parameter a increases, the system eventually reaches a critical threshold where period‐doubling bifurcations begin, the point where the eigenvalue λ = 1 – ln(α) – a falls to −1 from above. We studied the behavior of the DFT before the critical threshold is reached. A study of a similar model did not find early warnings (Hastings and Wysham 2010). However, Hastings and Wysham (2010) focused on time series statistics for individual patches and not the overall spatial pattern.

We present two simulations with Model 1. In simulation 1, σ is small, following the value in Table 1. Simulation 2 adds relative large shocks from a uniform distribution as used by Hastings and Wysham (2010). Uniform shocks to a were drawn from a range ±0.001 of the mean, as in Hastings and Wysham (2010). The variance of a in Simulation 2 is more than 30 times greater than in Simulation 1. In addition, Simulation 2 added shocks to α at each time step, drawn from a uniform distribution with range [0.4717, 0.5282]. This range gives a variance 8 times greater than the single‐patch variance used by Hastings and Wysham (2010) who studied studied passive dispersal in blocks of 8 patches. Like Hastings and Wysham (2010) we drew α shocks from below (above) the mean when shocks to a were below (above) the mean. This strong correlation can cause the eigenvalue λ = 1 – ln(α) – a to switch randomly above or below −1 as a rises close to the deterministic threshold. Therefore Simulation 2 may challenge the ability of the DFT to provide an early warning.

In real‐world applications of early warning indicators, the true process generating the data is unknown and the DFT would be computed from observed spatial data. We simulate this situation for all four models. Additionally, for Model 1 we also consider the case in which the analyst knows the true model and can therefore linearize the DFT near steady‐state (Appendix). Even though this situation is not likely to occur in practice, the linearized calculation is useful for understanding how the DFT works as an early warning. Expanding the dispersal term, Model 1 is [Image Omitted. See PDF]The linearization of (1) is [Image Omitted. See PDF]As an early warning indicator, we present below the modulus of the growth term for[IMAGE OMMITED. SEE PDF.] which is[IMAGE OMMITED. SEE PDF.].

Continuous time: Models 2, 3 and 4

Models 2–4 represent continuous‐time dynamics of a prey population with alternative stable states depending on the number of harvesters in the system. Model 2, which is also the core of Models 3 and 4, is similar to a spatial model studied by van Nes et al. (2005). Models 3 and 4 generate more complex spatial dynamics through schooling of prey and pursuit by predators (Model 3) or social aggregation of harvesters (Model 4). To Models 2–4, small shocks were added as a Weiner process with standard deviation σ. Stochastic dynamics were solved by the Euler method using Ito calculus, as in previous studies of early warnings in continuous time (Carpenter and Brock 2006, Carpenter et al. 2008).

Only the prey population is subject to dispersal in Model 2. Harvester density H is the same in each patch at each time step. H is gradually increased to the critical threshold. Parameters common to Models 2–4 include coefficients for growth (r), capture (c), dispersal of prey (D_x), and carrying capacity of prey (K).

Both prey and harvesters are subject to dispersal in Model 3. Prey exhibit schooling behavior in this model (Tyutyunov et al. 2004), and harvesters move among patches in response to abundance of prey and presence of other harvesters. As in Model 2, the total number of harvesters across all patches is increased gradually to a critical transition. Model 3 could represent the dynamics of a commercial fishing fleet pursuing a fish species in an ocean or large lake, or a population of hunters pursuing a population of prey in an extensive terrestrial system. Velocity of prey movement v depends on a rate coefficient k₀, maximum crowding coefficient x_c, and dispersal rate D_v (Tyutyunov et al. 2004). Harvesters move in response to the perceived utility of each patch (J) which depends on local density of prey and harvesters through the weights p₁, p₂ and p₃. Interpretations of the other parameters follow Model 2.

In Model 4, harvesters disperse but prey cannot disperse. The total number of harvesters is increased gradually to a critical transition. Model 4 represents spatial dynamics of human fishers harvesting fish on a lake‐rich landscape where the spatial locations correspond to individual lakes (Carpenter and Brock 2004, Brock and Carpenter 2007). At each time step human harvesters choose lakes based on the perceived utility (J) of fishing in each lake (Krugman 1993, Fujita et al. 1999). J for a given lake is directly related to prey numbers and amount of information (I) available to the harvesters, and inversely related to the number of harvesters currently using the lake, through the parameters p₁, p₂, p₃ and p₄. Harvesters have more information about the lake they currently occupy and its neighbors, according to the spatial kernel k with spread parameter s_I that appears in the information term. Interpretations of the other parameters follow Model 2.

Results

Model 1: early warning of period‐doubling bifurcation

Simulation results for Model 1 are presented as surface plots (Figs. 1, 2) where the colors denote population or the DFT ordinate on the z‐axis, the x‐axis is spatial location or spatial frequency, and the y‐axis is the bifurcation parameter a. Colors mimic a topographic map, with blue for lowest values rising through green and yellow to tan for the highest values (note key to the right of each surface plot). Early warnings, if present, would appear as increases in the DFT ordinate (colors on the z‐axis) at certain spatial frequencies (x‐values) as the values of a (y‐axis) increase toward the critical threshold at a ≈ 2.693 (corresponding to an eigenvalue of −1).

Ricker population dynamics with dispersal (Model 1) and small shocks. A. Prey versus a and location. B. Discrete Fourier Transform (DFT) for prey versus a and frequency. C. Rate of increase computed from the Linearized DFT for prey versus a and frequency. D. Spatial variance across all sites versus a.

Ricker population dynamics with dispersal (Model 1) and large shocks. A. Prey versus a and location. B. Discrete Fourier Transform (DFT) for prey versus a and frequency. C. Rate of increase computed from the Linearized DFT for prey versus a and frequency. D. Spatial variance across all sites versus a.

In Simulation 1 (small shock variance) with Model 1, the prey population increases by about 5% as the bifurcation parameter a rises gradually toward the critical transition at a ≈ 2.693 (Fig. 1A). At each value of a there is variability but no particular pattern is evident. Over the same range of a, ordinates of the DFT increase more than five‐fold, especially at low frequencies (Fig. 1B). The spatial variance computed across all patches increases more than five‐fold prior to the critical transition (Fig. 1D). The rising variance and rising DFT ordinates, especially at low frequencies, provide advance warning of the regime shift.

The modulus of the growth term obtained from the linearized DFT of Model 1 also increases especially at low frequencies (Fig. 1C). This growth term rises above 1 well before the critical transition at low frequencies, indicating that the system becomes unstable for these frequencies. This pattern is a powerful early warning that demonstrates the sensitivity of the linearized DFT to nearby regime shifts. Note, however, that it is necessary to know the true data generating process in order to compute the linearized DFT. The DFT in Fig. 1B and the spatial variance across all patches (Fig. 1D) are computed directly from the observed patterns. Although they are less powerful early warnings than the linearized DFT, they do not require prior knowledge of the true process that generates the data.

The outcome for Simulation 2 (large shock variance) with Model 1 is presented in Fig. 2. Population increases about 30%, and there is notable variability before a reaches the critical value of 2.693 (Fig. 2A). The DFT for direct observations increases more than five‐fold, especially at low frequencies (Fig. 2B). The growth term from the linearized DFT also increases with a, and at most frequencies it becomes unstable (exceeds 1; green tones) well before the critical threshold is reached. The spatial variance computed across all patches also rises more than five‐fold in advance of the critical transition. Therefore, early warnings are evident in the DFT and total variance for both simulations (low and high variance shocks) of Model 1.

Models 2–4: early warning of over‐harvest

For Models 2–4 (Figs. 3–5), we switch to a different graphical presentation because the shift between alternate states yields large changes in population and DFT ordinates that are difficult to portray in surface plots. Prey densities versus location, and DFT ordinates versus spatial frequency, are shown as colored lines, where a line is plotted for every 100 time steps of the simulation from step 1 to step 5000. At each point in time where a distribution is plotted, the color for the population distribution is identical to the color for its DFT. Colors ranging from blue to red to depict the gradient from “safe” values of H well below the critical threshold (blue) to “dangerous” values of H near the critical threshold (red) (see color key on each plot). The transition is easily seen as a jump from high to low prey values. An early warning would appear as a rise in DFT before the jump. For Models 2–4, the DFT is computed from observed patterns and does not require prior knowledge of the true data‐generating process.

Prey response in predator‐prey model with prey dispersal (Model 2) as H0 is gradually increased over time toward the critical value (from blue to red, see legend; color sequence is identical in each panel). A. Prey versus location. B. Discrete Fourier Transform of Prey versus frequency.

Harvester‐prey model with spatially‐fixed prey and mobile harvesters (Model 4) as H0 is gradually increased over time toward the critical value (from blue to red, see legend; color sequence is identical in each panel). A. Prey distribution versus location. B. Log (base 10) of the Discrete Fourier Transform (DFT) for prey versus frequency. C. Predator distribution versus location. D. Log (base 10) of the DFT for the predator versus location.

In Model 2, like Model 1, spatial dynamics are generated by prey dispersal (Fig. 3). However, in Model 2 the critical transition is a collapse of the prey population from a higher stable point to a lower stable point. There is a sharp decline in the prey population as the number of harvesters H rises to the critical transition at H_C ≈ 2.604 (Fig. 3A). The ordinates of the DFT rise substantially before the critical transition (Fig. 3B). The DFT increased by more than 100‐fold by 200 time steps before the transition. The rise in the DFT provides a clear early warning of the impending regime shift.

Both prey and the harvesters are mobile in space in Model 3 (Fig. 4). As in Model 2, the critical transition is a collapse of the prey population from a higher stable point to a lower stable point as the total number of harvesters in the system, H₀, rises above a critical value near 2.604 (Fig. 4A). As the critical transition approaches, the spatial pattern of harvesters becomes more patchy (Fig. 4C). The DFT ordinates for both prey (Fig. 4B) and harvesters (Fig. 4D) increase substantially prior to the critical transition, especially at lower frequencies. Spatial transects of either prey or harvesters provide a clear early warning through the DFT, even though changes in density of the harvesters (Fig. 4C) are rather subtle.

Predator‐prey model with prey dispersal and mobile predator (Model 3) as H0 is gradually increased over time toward the critical value (from blue to red, see legend; color sequence is identical in each panel). A. Prey distribution versus location. B. Discrete Fourier Transform (DFT) for prey versus frequency. C. Predator distribution location. D. DFT for predator versus location.

In Model 4, only harvesters are mobile in space (Fig. 5). As in Models 2 and 3, the critical transition is a collapse of the prey population from a higher stable point to a lower stable point as the total number of harvesters in the system, H₀, rises above a critical value (Fig. 5A). The critical transition is evident from the sharp drop in prey density. However, the critical value of H₀ cannot easily be computed for Model 4, in part because of the self‐organization of the harvesters (Fig. 5C). Over time the harvesters tend to aggregate, as has been shown for other spatial economic models that formulate human movement in similar ways (Krugman 1993, Fujita et al. 1999, Carpenter and Brock 2004). DFT ordinates for both prey and harvesters rise over time (Figs. 5B, 5D) and the change was so extreme as to require base‐10 log transform of the y‐axes for visualization. The increase in the DFT is steepest at low frequencies, though all frequencies show increases especially near the critical transition.

Discussion

The spatial DFT of transient data is a sensitive leading indicator of impending critical transition. We analyzed transients from diverse ecological models, in discrete and continuous time, for single species and for prey‐harvester systems. In all cases the DFT exhibited clear early warnings in advance of regime shifts. Even if we do not know the dynamical system and, hence, do not know how to perform a linearization, early warning signals tend to be earlier and sharper in the DFT domain compared to the spatial domain. The DFT complements other early warning indicators for spatial systems such as the spatial variance, skewness, autocorrelation, or power relationships (Oborny et al. 2005, Kefi et al. 2007, Guttal and Jayaprakash 2009, Dakos et al. 2010). In addition, the DFT provides information about the spatial frequencies (ω values) where the regime shift first occurs.

In earlier analyses of a spatial Ricker system, Hastings and Wysham (2010) did not discern early warnings of regime shift. Their analyses focused on the time series from individual patches. We used the information in the overall spatial pattern and found that the aggregate spatial variance and the DFT increased steeply in advance of the regime shift. The early warning occurred even when we added large correlated shocks to a and α as used by Hastings and Wysham (2010).

Spatial pattern may offer advantages for anticipating or detecting ecological regime shifts of the types studied here (Guttal and Jayaprakash 2009, Dakos et al. 2010). Unlike temporal indicators which require long unbroken time series of frequent observations, spatial indicators may be useful even if measurements are irregular or infrequent over time. In our analyses of Models 2–4, for example, early warnings were evident from samples separated by more than 200 time steps. While spatial pattern measurements require intensive data collection at each time point, in many cases this may be easier than high‐frequency time series sampling.

Nonetheless, dispersion in spatially‐structured systems could dampen fluctuations and thereby weaken early warnings of regime shifts. The Fourier transform addresses the dynamics in terms of spatial frequencies, which are decoupled (Appendix). We speculate that this decoupling accounts for the high sensitivity of the DFT observed in our computer experiments.

Researchers have noted several reasons that early warnings might not be detected prior to a regime shift (Kleinen et al. 2003, Carpenter and Brock 2006, Carpenter et al. 2008, Biggs et al. 2009, Contamin and Ellison 2009, Scheffer et al. 2009, Brock and Carpenter 2010, Hastings and Wysham 2010): (1) large process error—the magnitude of shocks to the dynamics may overshadow the early warnings or trigger a shock before the warning can be detected; (2) large observation error—imprecise observations may make it impossible to discern the early warnings; (3) rapid transition through the critical transition—if the system is forced across the critical threshold too fast, there may not be time for the early warning to be detected; (4) muffling—interactions of multiple thresholds and correlated process errors may weaken the signal of the early warning; (5) higher‐order terms, neglected in linearizations around the critical transition, may not damp off as expected and may interfere with early warnings; (6) certain kinds of critical transitions, such as global bifurcations, may not generate the usual early warning signals.

In the settings analyzed in this paper, limitations 1, 4 and 5 have little effect on the DFT. Like all early warning indicators proposed to date, the DFT is subject to limitations 2, 3 and 6. Simulations suggest that the DFT may be less sensitive to limitation 1 than some other indicators (compare Figs. 1 and 2), although very large process error could obliterate early warnings. Mathematical arguments (see Appendix, “Discrete Fourier Transform of Linearized Systems”) suggest that muffling (limitation 4) has little effect on the DFT, and this is a major advantage in the analysis of spatial data. The effect of higher‐order terms (limitation 5) on the DFT may be minor as shown in the argument below.

Where the underlying dynamics are known or good approximate models exist, the DFT of the linearized system provides especially powerful early warnings. In our simulations, the growth term of the linearized DFT rose above 1, indicating explosive change, well before the leading eigenvalue left the unit circle. In many cases good approximating models will not be known. Yet even when linearization is not possible the DFT provides early warnings.

To understand why the linearization is not essential, consider what happens if we add the higher‐order terms [Image Omitted. See PDF]to the linear approximation (equation (A.3) of the Appendix). Then the equation is not an approximation. It is exact. Here o_i denotes a term that goes to zero faster than the distance of x(t, z) from[IMAGE OMMITED. SEE PDF.] goes to zero. Hence the difference between the DFT of the observed pattern {x(t, z)} and the DFT of the centered raw data,[IMAGE OMMITED. SEE PDF.] is the term [Image Omitted. See PDF]Because [Image Omitted. See PDF]the difference in the DFT's will not obscure any early warning signals. The difference between the DFT of the raw data and the linearization's data is the DFT of the o_i which is small. While this heuristic argument is certainly not a proof of why the DFT domain is a better place to look for early warnings than the z‐domain, it is suggestive.

In applications to field data, it may be possible to improve the sensitivity of the DFT method by fitting models, such as polynomials, to the observed data, and then computing the DFT for the resulting polynomial forecasts. Following the same reasoning as above, if we expand the original dynamics around[IMAGE OMMITED. SEE PDF.] and keep the first few higher order terms beyond the linear terms, then a polynomial fit to the raw data might do a fair job of approximating this Taylor series. The higher‐order terms may also lead to early‐warning signals from higher moments such as skewness or kurtosis (Guttal and Jayaprakash 2008, 2009).

Several important issues require for further research. (1) We analyzed pattern‐forming or reaction‐diffusion systems (Murray 1983) using a discrete ring of interacting sites. The critical transitions in our models were generated by period‐doubling or fold bifurcations driven by slow changes in certain parameters. It remains to be seen how the DFT will perform with other types of spatial critical transitions. (2) Seasonal change and other cyclic phenomena affect ecosystem dynamics, and could interfere with spatial as well as temporal early warnings. Further study, in models and in the field, is needed to understand how early warnings are affected by seasonality. (3) Fourier transforms can be defined for very general spaces (Bamieh et al. 2002), whereas we focused on one spatial dimension as might be observed in a line transect through an ecosystem. Further work is needed to see how Fourier transforms perform in detecting early warnings in two or three dimensional data. (4) Also, research is needed on the impact of different boundary conditions on the detection of early warnings. We assumed a periodic structure (the ring) in our pattern formation systems because we wanted to focus on the ability of the DFT to detect early warnings that were not forced by boundary conditions. However, many important pattern formation systems in biology and ecology have different boundary conditions, some of which play a role in producing patterns (Murray 1983). (5) Finally, ecologists use a wide range of models and methods to study spatial pattern in ecosystems (Fortin and Dale 2005). Other statistics for spatial pattern (e.g., wavelets) may be as effective, or more effective, than the DFT for detecting critical transitions.

While our results show that the DFT is a sensitive early warning indicator of some types of critical transitions, much work remains to be done. The problem of early warnings for critical transitions in spatially‐structured ecosystems seems likely to attract the attention of researchers for a long time into the future.

Acknowledgments

Both authors contributed equally to this paper. We thank M.‐J. Fortin, M. L. Pace, D. A. Seekell and two anonymous referees for helpful comments. The National Science Foundation and the Vilas Trust funded our work.

Appendix

Discrete Fourier Transform of Linearized Systems

This Appendix develops some technical details of the DFT that are useful for understanding why it can be used as an early warning indicator. The DFT of spatial dynamics linearized near steady state should be very sensitive to approaching critical transitions. We provide a brief description of this linearized DFT in three parts: (1) a derivation of the linearized DFT; (2) the identity of the eigenvalues of the dynamical system in spatial coordinates to the eigenvalues in the spatial frequency coordinates of the DFT; and (3) a general discrete‐time case to discuss the sensitivity of the DFT to impending regime shifts.

The DFT as an early warning signal

This subsection develops the linearized DFT for a general model of an ecosystem variable in space. Many spatial dynamical systems in ecology are modeled using partial differential equations (Levin 1976). To simplify the exposition, we will use one spatial dimension to represent a transect through an ecosystem. We will solve the models on a ring to avoid edge effects. Consider the model [Image Omitted. See PDF]where z is a one dimensional spatial coordinate denoting site z in the space of N sites arranged in a continuous circle with circumference N, x(t, z) is an n‐dimensional vector, D is a n by n diagonal matrix, and the operators dx(t, z)/dt and ∂²x(t, z)/∂z² are the continuous derivative operator and the second order partial derivative operator or their discrete counterparts. We work with the continuous version of dx(t, z)/dt and the discrete version of the second order partial derivative operator (Levin 1976) [Image Omitted. See PDF]Let w̄(θ) be a spatially homogeneous steady state for parameter vector θtheta;. Expand (A.1) at w̄(θ) to obtain the linearization [Image Omitted. See PDF]We put A(θ) ≡ δf(w̄(θ), θtheta;)/δx. In expressions (A.3)–(A.5) we use x(t, z) to denote the difference, x(t, z) − w̄(θ), in order to simplify notation. Let [Image Omitted. See PDF]denote the DFT of the ith element of the n‐vector, x(t, z), i = 1, 2, … , n. Then we have [Image Omitted. See PDF]Now let θtheta; slowly (relative to the speed of the dynamics in (A.2) and (A.3)) approach a critical value θtheta;_c. At θtheta;_c the real part of the leading eigenvalue (or pair of complex eigenvalues) hits zero from below for continuous time systems and leaves the unit circle in the complex plane for discrete time systems. It is important to note that the eigenvalues of the system are the same in the z‐space and in the Fourier‐transformed ω‐space (see next section).

Intuition suggests that the spatial dispersal in (A.3) might muffle a potential early warning signal of the approach to local instability at any site z by spreading the effects of the instability over many sites. In contrast, the DFT's product expression might sharpen or even magnify a potential early warning signal, especially at an ω where the [–1 + cos(2πω/N)] term in (A.5) magnifies the real part of the leading eigenvalue(s) of system (A.5). This effect could be especially sharp in discrete time systems, e.g., the discrete time version of (A.5), where overshoot bifurcations can occur when the leading eigenvalue (s) pass through −1 (see “A Discrete Time Example” below).

Eigenvalues are identical in z‐space and ω‐space

Critical transitions are identified by certain changes in the eigenvalues of the dynamical system. If the DFT is to be useful as an early warning, we must understand how the eigenvalues for the linear system (equation A.3) in z‐space are related to the eigenvalues of the Fourier transformed system in ω‐space. Here we show that these eigenvalues are identical. The discussion below adds to the findings of Bamieh et al. (2002, Table 1 and surrounding material) who present Fourier pairings between z‐spaces and ω‐spaces for linear dynamic settings.

Rewrite (A.3) in vector form, [Image Omitted. See PDF]where the notation x_i(t,.) denotes the vector with components {x_i(t, z), z = 0, 1, … , N – 1} and the N by N matrix J_i is the matrix that takes the N by 1 vector x_i(t,.) into the vector y_i(t,.) with components y_i(t, z) = D_ii(x_i(t, z+1) – 2x_i(t,z)+x_i(t, z – 1)), z = 0, 1, … , N – 1. With this notation we may write (A.4) in the vector form [Image Omitted. See PDF]where the N by N invertible matrix F is defined by (A.4). In order to define eigenvalues and eigenvectors, we view the system (A.6) as a system that takes nN dimensional vectors into nN dimensional vectors. An eigenvector and corresponding eigenvalue of the differential equation system (A.1) is defined as follows [Image Omitted. See PDF]where the nN vector of initial conditions is an eigenvector of system (A.8) corresponding to eigenvalue λ.

Because [IMAGE OMMITED. SEE PDF.], we may rewrite the Nn by Nn system of differential equations (A.8) thus, [Image Omitted. See PDF]

For each i = 1, 2, … , n, (A.9) is the same as (A.5). Thus the N by N matrix FJ_iF⁻¹ must be the same as the N by N diagonal matrix with (ω,ω) element 2D_ii[−1 + cos(2πω/N)], ω = 0, 1, … , N – 1. Therefore if λ is an eigenvalue of the original system in z‐space it is an eigenvalue of the Fourier transformed system in ω‐space. The reverse argument is similar. Hence the eigenvalues of both systems are the same.

A discrete‐time example

How does the DFT augment aggregated indicators of spatial pattern, such as the overall spatial variance or autocorrelation? This section shows that the DFT is especially sensitive to the frequency where the critical transition first occurs. Consider the discrete time version of equation (A.3), where x_t(z) is scalar and the N‐vector[IMAGE OMMITED. SEE PDF.] and I denotes the N by N identity matrix and e is a vector of random shocks: [Image Omitted. See PDF]Here C is the circulant matrix whose first row is (c₀, c₁, c₂ ,… , c_N-2, c_N-1) = (−2, 1, 0, … , 0, 1) and each further row is a cyclic shift of the row above it. A N by N circulant matrix C is a matrix where each row is a cyclic shift of the row above it. Formulae for the eigenvalues, eigenvectors, etc., of circulant matrices are given by Gray (Gray 2006) . We also use the convolution theorem and the shift theorem for DFTs (Chu 2008). The convolution theorem states that the DFT of a convolution sum is the product of the DFT's. The shift theorem states that the DFT applied to {x(t, z – k), z = 0, 1, … , N – 1} for any integer k = 0, 1, … , N – 1 is [IMAGE OMMITED. SEE PDF.][IMAGE OMMITED. SEE PDF.]

The Fourier transform of (A.3b) is [Image Omitted. See PDF]or, equivalently, applying the shift theorem, we have, [Image Omitted. See PDF]In equation (A.3b) the patches are coupled by the diffusion term but in equation (A.3d) the frequencies are uncoupled. The decoupling occurs because, applying the shift theorem to the special form of the circulant matrix in (A.3b) above, we have, [Image Omitted. See PDF]To put it another way, Fourier transformation take the nondiagonal system [Image Omitted. See PDF]into the diagonal system (A.3d). Since the N eigenvalues of equation (A.3d) are just the coefficients of [IMAGE OMMITED. SEE PDF.], our search for an early warning signal is directed to [Image Omitted. See PDF]At this frequency, the linearized DFT should be maximally sensitive to an impending regime shift.