Percolation-based precursors of transitions in extended systems

Abstract

Abrupt transitions are ubiquitous in the dynamics of complex systems. Finding precursors, i.e. early indicators of their arrival, is fundamental in many areas of science ranging from electrical engineering to climate. However, obtaining warnings of an approaching transition well in advance remains an elusive task. Here we show that a functional network, constructed from spatial correlations of the system's time series, experiences a percolation transition way before the actual system reaches a bifurcation point due to the collective phenomena leading to the global change. Concepts from percolation theory are then used to introduce early warning precursors that anticipate the system's tipping point. We illustrate the generality and versatility of our percolation-based framework with model systems experiencing different types of bifurcations and with Sea Surface Temperature time series associated to El Niño phenomenon.

Percolation-based precursors of transitions in extended systems

OPEN

Vctor Rodrguez-Mndez, Vctor M. Eguluz, Emilio Hernndez-Garca & Jos J. Ramasco

Abrupt transitions are ubiquitous in the dynamics of complex systems. Finding precursors, i.e. early indicators of their arrival, is fundamental in many areas of science ranging from electrical engineering to climate. However, obtaining warnings of an approaching transition well in advance remains an elusive task. Here we show that a functional network, constructed from spatial correlations of the systems time series, experiences a percolation transition way before the actual system reaches a bifurcation point due to the collective phenomena leading to the global change. Concepts from percolationtheory are then used to introduce early warning precursors that anticipate the systems tipping point. We illustrate the generality and versatility of our percolation-based framework with model systems El Nio phenomenon.

The occurrence of sharp transitions to dierent states or regimes during the evolution of complex systems is a phenomenon of major importance both from the fundamental point of view and for practical implementations of control and management. Examples of such abrupt changes can be found in ecology1, economy2, electrical engineering3, physiology4 or climate57. Detecting with sufficient anticipation the approach to a critical or tipping point thus becomes an important issue. Early-warning signals have been introduced and tested in recent works810,

including experimental verication in living and environmental systems9,11,12. These methods rely on the loss of resilience occurring generically when dynamical systems approach most (although not all) types of bifurcation points6,10,13. Recovery rates from perturbations become small, leading to critical slowing down of the dynamics, increased memory, long temporal autocorrelations, and to the growth of the temporal variance6,10,13. From a

dynamical viewpoint, these critical slowing down phenomena appear when the eigenvalue of the Jacobian matrix describing the rate of relaxation towards the attractor approaches zero close to bifurcation points. In many cases, particularly when dierent spatial parts of the system are coupled by diusion-like processes, the increase in temporal correlation is accompanied by the growth of spatial correlations1416. This, with the associated increase in spatial variance and response functions, is actually a standard method to characterize phase transitions in thermodynamic physical systems17.

The consideration of spatial correlations has led to a novel perspective for nding transition precursors through the use of correlation or functional networks1821. These networks are built by identifying spatial units as nodes in a graph, measuring the correlation among all pairs of them, and keeping the most signicant ones as link weights. Several network-based precursors have been proposed. Specically the values of the degree (number of connections per node), assortativity (degree-degree correlations), clustering (average density of triangles) and kurtosis rise when approaching a tipping point2224.

Here we show an important additional property of these functional networks: As internal correlations increase, networks evolve from a low to a high connectivity state and, before reaching its maximum link density at the bifurcation point, a percolation-like transition occurs in the network topology. Concepts from graph percolation theory25,26 can thus be imported to characterize this transition. Importantly, metrics can be dened that act as early warnings for this percolation transition, which itself is a precursor of the dynamic transition. The validity of these general ideas is tested by analyzing model systems displaying dierent types of bifurcations: steady and oscillatory, continuous and discontinuous. In all cases, we observe the occurrence of percolation transitions before the dynamical one, and we characterize it with quantities that can be used as early-warning signals. Our approach uses only time series from the elements of an extended system, without the need of specic knowledge about the underlying dynamics. Thus, it is suited to analyze observational data for which little or no modeling

Instituto de Fsica Interdisciplinar y Sistemas Complejos IFISC (CSIC-UIB), Campus Universitat de les Illes Balears,

A

SCIENTIFIC REPORTS

1

www.nature.com/scientificreports/

insight is available. This property of our framework is illustrated by applying it to temperature data from the Pacic Ocean associated to El Nio phenomenon.

Methods

Functional networks and precursors. The time evolution of extended dynamical systems is described by time-dependent spatial elds. Let (x, t) be one of such elds, and consider a suitable discretization of it {(xl, tk)}l,k, dening a time-series at discrete times tk, k=1, , R, from each spatial location xl, l=1, , N. The construction of functional networks implies to compute the Pearson correlation from the time series at every pair of locations:

= p t p t

( ) ( ) ,(1)

k a k b k

1 is the deviation of the eld from its temporal mean at each location. A network in which nodes are the spatial locations xl is dened by assigning links between pairs of nodes (xa, xb) for which the Pearson correlation in equation(1) is higher than a predened threshold : ab>.

To study percolation in these functional networks, as the control parameter p approaches transition points pd,

we measure several metrics: (a) S1, the relative size of the largest connected component, i.e. the fraction of nodes that are in the largest cluster. It abruptly changes from a value close to 0 to a value close to 1 at the percolation point. (b) The average size of the clusters excluding the largest one. This quantity is maximal at the percolation point25. It can be calculated as

where

where the sum runs over all cluster sizes s excluding the largest one, and cs is the fraction of nodes belonging to clusters of size s, cs=sns/N (ns is the number of clusters of size s present in the system26) thus cs also gives the probability that a randomly chosen node pertains to a cluster of size s.

The standard percolation indicators, S1 and s2 are not, however, the best precursors. We show below that better anticipation can be obtained by exploiting the interplay between the probabilities cs and the coming transition. In random graphs it is possible to perform analytic calculations on the behavior of cs. As new links are added at random, the percolation occurs when the mean degree k, which acts as the control parameter p, equals one.

Before this, the probabilities cs that a randomly chosen node belongs to a component of size s can be written as

= .

c e s k s

( ) / ! (3)

s

s k s 1

These probabilities have a maximum when the mean degree, which in this case is the parameter p, is

= =

k ps

s

ab

p t t t

x x

( ) ( , ) ( , )

l k l k R k l k

1 ,

2 2

=

sc N s n s

s

s

s

s

(2)

1 . The succession {ps} of the location of the maximum of cs converges to the percolation point p= 1 for increasing component size s, but for low s these maxima could be quite far from the percolation point, and in fact they anticipate it. We show below that the early-warning character of the maxima of cs holds also true in non-random functional networks obtained from systems undergoing very dierent dynamic transitions. This is due to the generic increase in the dynamical correlations of the elementary units of the system that the approach of a global tipping point brings. The sequence of peaks in cs oers thus a general and versatile tool to predict potential changes in the dynamics at a global scale.

Percolation and transitions in model systems. We analyze here three dierent extended dynamical systems displaying dierent types of bifurcations. The rst two examples experience steady bifurcations (a discontinuous saddle-node and a continuous pitchfork). The third case, the Lorenz96 system, experiences a variety of transitions being the rst one an oscillatory Hopf bifurcation between a steady state and traveling waves. Further transitions occur when changing the control parameter leading to low-coherence spatio-temporal chaos. In all cases we add random noise to the deterministic model. This represents the unavoidable stochastic uctuations to which real systems are always subjected and provides the necessary statistics to have well-dened spatial correlation functions. In all cases percolation occurs in the associated functional networks, providing robust early-warning signals of the approaching transition.

Results

A lake eutrophication model. As a rst example, we consider a lake eutrophication (LE) model which is a spatial version of a description of phosphorous recycling in a lake27. It suers a paradigmatic abrupt transition associated to a saddle-node (SN) bifurcation, namely a transition between two contrasted states for the phosphorus concentration (x, t) in the lake: low concentrations leading to clear water, and excess of phosphates leading to turbid water. The state of the system is given by the two-dimensional eld (x, t), representing the amount of phosphorus in the lake, evolving according to

= + + + .

s

= +

f ( ) /( 1)

x x x

( , ) ( , ) ( ) ( , ) ( , ) (4)

t

x t p b t rf t t

2

8 8 is a nonlinear response of the lake sediments to phosphorus, p is the nutrient input rate, taken here as the control parameter. is the strength of diusive spatial coupling, and represents an additive stochastic perturbation uncorrelated in space and time. We take b=r=1 and = 1.2. Space is discretized as a square lattice of N=7070= 4900 grid points separated by dx= 1. These will be the nodes of the functional

The function

SCIENTIFIC REPORTS

2

www.nature.com/scientificreports/

Figure 1. Transition precursors for the LE model, equation (4). In (a) the steady homogeneous phosphorus concentration (orange) and the numerically obtained spatial average (black, further averaged over R=1000 temporal snapshots for each value of p) as a function of the control parameter p which is slowly increased from low to high values. In (b) , the spatial standard deviation (black) of used in refs 1315 as a transition precursor is displayed averaged over R= 1000 temporal snapshots. In orange, the clustering of the functional network built with threshold = 0.21 also used as a precursor in refs 22,23. In (c) the relative size of the giant component, S1 (orange), and the average size of the leover clusters (s2, black). (d) The probabilities c2 (black), c5 (orange), and c9 (purple) are shown. In all the panels, the vertical green line marks the position of the observed abrupt transition (Pd= 0.658). The grey area indicates the anticipation in parameter space gained over previous methods by using the peak of c2 as precursor of the percolation transition. All curves have been further averaged over 100 realizations of the random noise and initial conditions.

network. The Laplacian is discretized with the simplest nite dierences scheme and the deterministic terms in equation(4) are integrated with a 4th order Runge-Kutta method of time step dt= 0.05 aer which the term is implemented by adding an independent random number uniform in [a, a] to each lattice site (we use a=0.125). We approach from the left the SN bifurcation occurring at p = pSN = 0.660, above which the clear-water low-phosphorus state existing for p<pSN ceases to exist and the lake jumps to a eutrophicated high-phosphorus state. The stochastic perturbation makes the jump to occur at a value of p, pd, slightly below the SN value. On average (see Fig.1), we nd pd 0.658. This value is sufficiently close to pSN as to display the enhancement of correlations and slowing down which are at the basis of our method and of other early-warning methodologies. In this paper we show only results obtained when slowly increasing the control parameter p. When decreasing p from higher values hysteresis occurs and a dierent SN bifurcation is encountered at a lower .

p 0 389

d . The

sequence of precursors encountered when approaching this lower transition point is similar to the one shown here.

The increase of the spatial variance was used in refs 14,15 as transition precursor. This variance is dened for the discretized eld (xl, tk) in the asymptotic statistically steady state as:

N

=

1 1 is the average over the nodes. When the system approaches pd, the raise of spatial variance can be observed for the LE model (Fig.1b).

As explained above, we constructed functional networks by assigning links between locations among which correlations (as measured by Pearson correlation) are larger than a threshold for which we use = 0.21. We take averages over R= 1000 temporal snapshots. The spatial correlations, computed from equation(1), increase and lead to a growth of the link density in the vicinity of the critical point as indicated by the precursors proposed in refs 2224. In fact one of these network precursors, the clustering coefficient, is plotted in Fig.1b (orange) and

N t t

=

1 ( ( , ) ( )) ,

x

2

l

l

1

(5)

N

where

= =

t t

x

( ) ( , )

N l

l

SCIENTIFIC REPORTS

3

www.nature.com/scientificreports/

Figure 2. Role of the correlation threshold . c2 values, as given by the color bar, as a function of the control parameter p and the threshold used to build the functional network for the LE model. The white dots give the locations of the c2 maxima, while the black dots mark the maxima of s2 (percolation transition). The dynamical sudden jump occurs at pd= 0.658. The horizontal blue line identies the value of = 0.21 used in Fig.1.

has a peak at the dynamical transition. But we will show that there is also a percolation transition, and to capture it we studied how the size of the giant component, S1, and the average size of the leover clusters s2 change with p (see Fig.1c). Note that s2 has a peak at p = 0.648, which identies the occurrence of percolation in the network, way before the SN transition has happened, and thus it may be used as a signal that a dynamical transition is coming. The distance in the parameter space between that signal and the transition depends on which threshold is used to build the network, but there is an optimal value (see below).

The probabilities c2, c5 and c9 for the networks built from the LE correlations are shown in Fig.1d. Their maxima clearly anticipate the percolation transition signaled by S1 and s2, which gives itself an early warning of the

SN bifurcation. The peak in c2 appears at p2 = 0.635, anticipating pd more than twice as early as the percolation transition.

Building functional networks requires setting the correlation threshold above which two elements are considered as linked. Figure2 shows how the percolation transition and the values of c2 depend on both p and . If

is very high, the network is never connected and there is no signal. Similarly, if is very low the network is always fully connected and there is no hint of the bifurcation. However, as shown in the gure, there is a range of values of where the percolation transition and its associated early warning signals appear. For a xed value of the peak of c2 (white dots) always occurs earlier than the percolation transition (black dots). As decreases, the peak of c2 occurs at earlier values of p. The curve of c2, however, widens and the resolution in the location of the peak gets poorer. The lowest (optimal) value of at which the peak can be distinguished marks thus the earliest warning signal that can be obtained for the bifurcation. Note that this does not imply that the method only works for a xed value of . In an empirical situation, one may need to explore this parameter, but there is a range of values of over the optimal that will provide valid early warning signals for the transition.

To prove the generality of the precursors, we analyze a dierent system, the time-dependent Ginzburg-Landau (GL) equation or model A28 that describes, for example, transitions in anisotropic ferromagnets. It is a paradigmatic model experiencing a continuous transition, namely a supercritical pitchfork bifurcation. We study the one-dimensional version for the magnetization (x, t):

= + + .

x t

( t p x t x t x t x t

3 2

As before, is diusive coupling and is an additive uncorrelated noise uniform in [a, a]. We discretize equation(6) into N = 5000 nodes, and take = 1.5, a = 0.01. The integration parameters dt and dx are as for the LE model. A continuous transition from zero magnetization occurs when increasing the control parameter p above pd=0 (Fig.3a).

We compute spatial correlations and build functional networks from R = 1000 snapshots using a threshold = 0.25. The corresponding percolation quantiers are plotted in Fig.3. Here correlations continuously build-up when increasing p towards the critical point at pd= 0 and, unlike the previous discontinuous transition, continuously decrease aer crossing it. Therefore, there are now two percolation transitions in the functional networks: one at each side of pd, as seen by the indicators S1 and s2 in Fig.3b. Figure3c depicts the probabilities c2, c5 and c9, with maxima giving a clear warning further away from pd. Figure3d displays the values ps corresponding to these maxima in cs. There are two successions of peaks converging to the percolation transitions occurring before and aer the bifurcation.

Figure4 shows what happens to the percolation transition when building the networks for dierent thresholds . The gure displays the maxima of s2 (which locate the percolation transitions) and the values of c2, in the

, ) ( , ) ( , ) ( , ) ( , ) (6)

SCIENTIFIC REPORTS

4

www.nature.com/scientificreports/

Figure 3. The Ginzburg-Landau model has a continuous transition at the critical value pd=0 marked by the green line. (a) The theoretical (yellow) and the numerically obtained (black dash line, further averagedover R= 1000 temporal snapshots) homogeneous value of the eld . (b) S1 (orange) and s2 (black). (c) The quantities c2 (black), c5 (orange), and c9 (purple). (d) Circles indicate the values of p, ps with s= 2, 5, 9, 16 and 21, for which cs attains its respective maximum. Functional networks were built using a threshold of =0.25.

Figure 4. The c2 values, as given by the color bar as a function of the control parameter p and the threshold used to build the functional network for the Ginzburg-Landau system. The continuous pitchfork bifurcation occurs at pd= 0. Black dots indicate the maxima of s2, which locate the percolation transition.

White dots locate the maxima of c2. We see how the percolation transition and its precursor c2 anticipate in dierent amounts (when increasing p from the low p state, or when decreasing p from the high p state) the dynamical bifurcation. The horizontal green line identies the value = 0.25 for which Fig.3 of the main text was constructed.

(, p) parameter space. We see that, when increasing p, the maximum in the precursor c2 anticipates the percolation transition, which itself anticipates the pitchfork bifurcation at p= 0. At the other side of the transition, when decreasing p from the high p state, the maximum in c2 also occurs before the percolation transition, which also anticipates the dynamical transition. As in the LE case the amount of anticipation is larger for lower , until the signal disappears.

Coupled chaotic oscillators display a large variety of dynamical regimes. Thus, due to the dierent bifurcations present in those models, they are an excellent test bed to prove the generality of the network-based percolation precursors. Here we consider the Lorenz96 model29. It was proposed by E. Lorenz as a simplied framework to investigate atmospheric predictability. It reads:

SCIENTIFIC REPORTS

5

www.nature.com/scientificreports/

Figure 5. (Panel a) shows the bifurcation diagram of the Lorenz96 model, equation(7), constructed from a 2-oscillator Poincar section (see text). The transition to traveling waves at p1= 8/9 and to spatiotemporal chaos (p2 4.1) are shown as vertical green lines. Functional networks were built using a threshold of =0.16. The percolation indicators S1 (orange) and s2 (black) are displayed in (panel b). (Panel c) shows c2 (black),c5 (orange), and c9 (purple). In (d), the circles indicate the value of p, ps with s= 2, 5, 9, 16 and 21, for whichthe cs curves attain their respective maxima. These indicators reveal a phase of percolated correlations in a parameter region which includes the interval [p1, p2], anked by two percolation transitions. The curves have been further averaged over 100 realizations of the random noise and initial condition.

d tdt t t t t t p

= + + .

+

1 2 1

k k k k k

k N

1 , arranged in a one dimensional ring around the globe (and thus having periodic boundary conditions). The structure of equation(7) contains some of the main elements of uid dynamics, namely dissipation, external forcing, and quadratic non-linearity through an advection-like term. In addition to the constant forcing p, which will be our control parameter, we include an additive stochastic perturbation k(t) uncorrelated in space and time. It is implemented here by adding independent random numbers uniform in [a, a] (a= 0.1) aer each time step (dt=1/64) of a fourth-order Runge-Kutta method which is used to integrate the rest of the terms. We focus in the behavior for N = 2500 elements or oscillators. In the absence of the random forcing, three dynamical regimes are easily identified30: (i) if p < p1 = 8/9, the system stabilizes in the homogeneous fixed point =

k(t) is meant to represent the values of some atmospheric variable at dierent locations k, =

t p k t

( ) , ,

k ;(ii) a Hopf bifurcation occurs at p = p1 so that for intermediate values of p, 8/9 < p < 4.1, the system is in a traveling-wave state; (iii) for p>p2 4.1, the system becomes spatiotemporally chaotic.

To display the bifurcations observed when integrating equation(7) for N= 2500 elements we have calculated a Poincares transversal section in the subspace of two contiguous oscillators: one of the oscillators, say 1, is monitored and when it crosses the value 1= 1 in the increasing direction the value of the contiguous oscillator, say = 2 is recorded and displayed. The bifurcation diagram showing the values of these sections is plotted in Fig.5a as a function of p. The three regimes described above for the deterministic system are readily identied here also.

Functional networks were constructed from R= 1000 snapshots by interpreting the locations k as nodes and assigning links between pairs of nodes when the Pearson correlation is larger than a threshold y=0.16. Figure5b shows the quantities S1, the size of the largest cluster, and s2, the mean cluster size excluding the largest one, for such network. Panels (c,d) display the behavior of s2, the probability that a randomly chosen node belongs to a cluster of size s. These gures clearly identify the presence of a percolated phase at intermediate values of p, started and ended by two percolation transitions. Figure6 shows the quantity c2 in the (, p) parameter plane for this model. We see that for increasing threshold the maxima of c2 approach the locations p1 and p2 at which traveling waves are born via a Hopf bifurcation and at which they destabilize into chaotic behavior, respectively. Thus, the percolating phase is a manifestation of the long-range coherence of the traveling wave state, whereas correlation length remains small in the homogeneous and in the chaotic regime. The quantity c2 (and indeed the other cs) clearly anticipates the rst bifurcation when increasing p. It also largely anticipates the occurrence of a chaos-order transition when decreasing p from large values.

k ( ) [ ( ) ( )] ( ) ( ) ( ) (7)

SCIENTIFIC REPORTS

6

www.nature.com/scientificreports/

Figure 6. The c2 values, as given by the color bar, as a function of the control parameter p and the threshold used to build the functional network for the Lorenz96 system. The dynamical transitions to traveling waves and to chaos are indicated by the vertical black lines. The horizontal green line identies the value =0.16 for which Fig.5 was constructed.

Percolation in sea temperature networks during El Nio events. To test the behavior of our precursors in observed real situations, we analyze sea surface temperature data from the region of the Pacic used to compute the NINO3.4 index31. El Nio-Southern Oscillation32,33 is the dominant variability mode in present-day climate, characterized by rather irregular (with average period of about 4 years) warm (El Nio) and cold (La Nia) episodes departing from the long-term mean temperature in the equatorial Pacic. These oscillations are related to the presence of a Hopf bifurcation in the coupled atmosphere-ocean system32,33. The bifurcation can be crossed or just approached, being then the oscillation excited by noise. In both cases there should be a build-up of correlations that would become visible in functional networks constructed from temperature time series. In this case there is no control parameter to x, but rather the equatorial Pacic evolves in time, coupled to the seasonal cycle, leading to changing spatial correlations. We will see that, despite this lack of control, and without using any information on the underlying dynamics, our approach is able to nd precursors of the relevant El Nio-La Nia events.

Sea Surface Temperatures were obtained from the ERA-interim reanalysis of the European Centre for Medium-Range Weather Forecasts34, with daily temporal resolution and a spatial resolution of x = 0.125, in the range of years 19792014 (Fig.7). Daily functional networks at day t were built from these time series computing the Pearson correlation with a time window of R = 200 days (100 days before and 100 days aer time t). The quantities plotted in Fig.7 are further averaged over 5 days.

In Fig.7, we have focused on three dierent periods: 19871989, during which a strong La Nia occurred, 19961998, featuring an El Nio-La Nia pair, and a recent El Nio in 2009. In contrast to the previous examples, this system is empirical and the contribution of the noise is more difficult to assess than in a model equation. Therefore, the systematic search for the optimal has not been performed. Nevertheless, the space of values of and how they aect c2 have been explored in Supplementary Figure S1. Interestingly, it seems that the range of values of necessary to observe peaks in c2 have moved toward lower in the early 2000s. We have xed two values of the threshold to produce Fig.7: = 0.99992 for the events of 19871989 and 19961998, and = 0.9986 for the event of 2009. In the two cases, this is where a nice compromise between signal-and-noise is found. In a practical situation, the selection of is not a post-hoc process: one can have a clear idea of the range of values to use from the previous events. Once is set at a xed value in this range, if c2 shows a peak, followed by a sequence of peaks of c3, c5, etc., the system is very likely going towards a new El Nio/Nia event. The panels a, b and c of Fig.7 display the variations of the ocean supercial temperature and also the moments at which an El Nio or La Nia event are officially declared are marked with a vertical orange line. The panels (df), on the other hand, depict the time evolution of the size of the largest connected component S1, which peaks before or on the arrival of the event. c2, represented in the lower panels (gi), also shows maxima way before the corresponding peak of S1. The anticipatory period since the c2 peaks to El Nio (La Nia) event is marked in gray in every plot. It corresponds to 240 days in 1988 (Fig.7g), 125 days in 1997 (Fig.7h), 175 days in 1998 (Fig.7h) and 115 days in 2009 (Fig.7i).

Conclusions

In summary, we have shown that consideration of the percolation transition in functional networks constructed from spatial correlations in extended systems provides powerful anticipatory tools for their dynamical regime shis. Precursors of the percolation transition itself, such as the probabilities cs for random nodes to belong to small clusters of size s, add extra anticipatory range. This is done by introducing a mesoscopic view of the system, instead of using global or local perspectives as the ones used in previous methods using the system dynamics slowing down, the spatial variability of the order parameter or the clustering and degree distribution of the functional network. Furthermore, the sequence of peaks of cs provides extra information on the distance still remaining to the (percolation) transition. We note that and s are methodological parameters, so that they can be explored even when far from the dynamical transition. Despite the fact that there exists an optimal for

SCIENTIFIC REPORTS

7

www.nature.com/scientificreports/

Figure 7. Application to El Nio phenomenon. The upper panel shows a map of the area over which themean Sea Surface Temperature T is monitored in the NINO3.4. Points denote locations used here as nodes in a functional network. The map was generated with Cartopy 0.11.036. Four events, two La Nia (cold) and two El Nio (warm), are shown in the time axis of (panels ai). Conventional starting dates of the events are marked by vertical orange lines. In (ac), the sea surface temperature T is shown as a function of time. A functional network is constructed from correlations at = 0.99992 for 19871989 and 19961998, and = 0.9986 for 2009. The size of the giant component S1 is shown in (df), showing percolating phases at a plateau, anked by two percolation transitions, which occurs before each of the events. Panels (gi) show c2 in the same time frame. Peaks in c2 ank both sides of the percolation plateaux, in a manner similar to the Ginzburg-Landau case shown in Fig.3. The time by which the peaks of c2 anticipate the conventional starting date of the event is marked in gray.

which the anticipatory power is largest, there is typically a wide range of values of for which the method works. Furthermore, the sequence of peaks of cs can give a good hint on the distance to the approaching transition.

The tools presented here work in a variety of transition types, the condition being the increase of spatial correlations when approaching the transition point. This happens generally at least for systems close to bifurcations characterized by critical slowing down and with spatial locations coupled by diusion. Most of the local bifurcations types satisfy the critical slowing down criterion (see discussion in refs 6,10). Although diusive coupling is sufficient to provide increasing spatial correlations when combined with critical slowing down, it is by no means necessary, as the example of the Lorenz96 model (for which spatial coupling is of advection type rather than diusive) shows.

Through this paper we have focused in spatially embedded complex systems, in which network nodes are associated to spatial locations. Since the only information needed to apply our framework is a set of time series coming from dierent network nodes, we expect our approach to be also useful in more general network systems experiencing regime transitions35, beyond the spatial ones. In addition we have shown that the percolation-based precursors can be used even in cases, such as El Nio events, where very little or no information on the underlying dynamics is available. This is, therefore, a fresh perspective on a known phenomenon with the bonus of oering a method that can become instrumental in the monitoring and management of complex systems.

SCIENTIFIC REPORTS

8

www.nature.com/scientificreports/

References

1. Scheer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. Catastrophic shis in ecosystems. Nature 413, 591596 (2001).2. Yan, W., Woodard, R. & Sornette, D. Diagnosis and prediction of tipping points in nancial markets: Crashes and rebounds. Physics Procedia 3, 16411657 (2010).

3. Dobson, I., Carreras, B. A., Lynch, V. E. & Newman, D. E. Complex systems analysis of series of blackouts: Cascading failure, critical points, and self- organization. Chaos 17, 026103 (2007).

4. van de Leemput, I. A. et al. Critical slowing down as early warning for the onset and termination of depression. Procs. Natl. Acad. Sci. USA 111, 8792 (2014).

5. Dakos, V., Scheer, M., van Nes, E. H., Brovkin, V., Petoukhov, V. & Held, H. Slowing down as an early warning signal for abrupt climate change. Procs. Natl. Acad. Sci. USA 105, 1430814312 (2008).

6. Thompson, J. M. T. & Sieber, J. Predicting Climate Tipping As a Noisy Bifurcation: a Review. International Journal of Bifurcation and Chaos 21, 399423 (2011).

7. Lenton, T. M., Livina, V. N., Dakos, V., van Nes, E. H. & Scheer, M. Early warning of climate tipping points from critical slowing down: comparing methods to improve robustness. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 11851204 (2012).

8. Scheer, M. et al. Anticipating Critical Transitions. Science 338, 344348 (2012).9. Wang, R. et al. Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature 492, 419422 (2012).10. Dakos, V., Carpenter, S. R., Van Nes, E. H. & Scheer, M. Resilience indicators: prospects and limitations for early warnings of regime shis. Philosophical Transactions of the Royal Society of London B: Biological Sciences 370, 20130263 (2015).

11. Veraart, A. J. et al. Recovery rates reect distance to a tipping point in a living system. Nature 481, 357359 (2012).12. Quail, T., Shrier, A. & Glass, L. Predicting the onset of period-doubling bifurcations in noisy cardiac systems. Procs. Natl. Acad. Sci. USA 112, 93589363 (2015).

13. Scheer, M. et al. Early-warning signals for critical transitions. Nature 461, 5359 (2009).14. Guttal, V. & Jayaprakash, C. Spatial variance and spatial skewness: leading indicators of regime shis in spatial ecological systems, Theoretical Ecology 2, 312 (2009).

15. Dakos, V., van Nes, E. H., Donangelo, R., Fort, H. & Scheer, M. Spatial correlation as leading indicator of catastrophic shis. Theoretical Ecology 3, 163174 (2010).

16. Dakos, V., K, S., Rietkerk, M., van Nes, E. H. & Scheer, M. Slowing down in spatially patterned ecosystems at the brink of collapse. The American Naturalist 177, E153E166 (2011).

17. Chaikin, P. M. & Lubensly, T. C. Principles of condensed matter physics (Cambridge University Press, Cambridge, UK, 1995).18. Tsonis, A. A. & Roebber, P. J. The architecture of the climate network, Physica A 333, 497504 (2004).19. Eguluz, V. M., Chialvo, D. R., Cecchi, G. A., Baliki, M. & Apkarian, A. V. Scale-free brain functional networks. Phys. Rev. Lett. 94, 018102 (2005).

20. Timme, M. & Casadiego, J. Revealing networks from dynamics: an introduction. Journal of Physics A: Mathematical and Theoretical 47, 343001 (2014).

21. Donges, J. F., Zou, Y., Marwan, N. & Kurths, J. Complex networks in climate dynamics. The European Physical Journal Special Topics 174, 157179 (2009).

22. Van Der Mheen, M. et al. Interaction network based early warning indicators for the Atlantic MOC collapse. Geophysical Research Letters 40, 27142719 (2013).

23. Tirabassi, G. et al. Interaction network based early-warning indicators of vegetation transitions. Ecological Complexity 19, 148157 (2014).

24. Feng, Q. Y., Viebahn, J. P. & Dijkstra, H. A. Deep ocean early warning signals of an Atlantic MOC collapse. Geophysical Research Letters 41, 60086014 (2014).

25. Stauer, D. & Aharony, A. Introduction to percolation theory 2nd ed (Taylor & Francis Inc, Philadelphia, USA, 1994).26. Newman, M. E. J. Networks: An introduction (Oxford University Press, Oxford, UK, 2010).27. Carpenter, S. R., Ludwig, D. & Brock, W. A. Management of eutrophication for lakes subject to potentially irreversible change. Ecological Applications 9, 751771 (1999).

28. Hohenberg, P. C. & Halperin, B. I. Theory of dynamic critical phenomena. Rev. Mod. Phys. 49, 435479 (1977).29. Lorenz, E. N. Predictability: a problem partly solved. in Proceedings of the Seminar on Predictability Vol. I, ECMWF Seminar, edited by Palmer, T. & Hagedorn, R. (ECMWF, 1996) pp. 4058.

30. Karimi, A. & Paul, M. R. Extensive chaos in the Lorenz-96 model. Chaos 20, 043105 (2010).31. Ludescher, J. et al. Improved El Nio forecasting by cooperativity detection. Procs. Natl. Acad. Sci. USA 110, 1174211745 (2013).32. Sarachik, E. S. & Cane, M. A. The El Nio-Southern Oscillation Phenomenon (Cambridge University Press, Cambridge, UK, 2010).33. Dijkstra, H. A. The ENSO phenomenon: theory and mechanisms. Advances in Geosciences 6, 315 (2006).34. ECMWF. European Centre for Medium-Range Weather Forecasts, 2009: ERA-Interim Project. Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory, Boulder, CO. Accessed on September 15, 2015, [available online at http//apps.ecmwf.int/datasets/data/interim-full-daily] (2009).

35. Dorogovtsev, S. N., Goltsev, A. V. & Mendes, J. F. F. Critical phenomena in complex networks. Rev. Mod. Phys. 80, 12751335 (2008).36. Met Office. Cartopy a cartographic python library with a matplotlib interface, version 0.11.0, (Exeter, Devon). Installed in May 2015, [available at http://scitools.org.uk/cartopy/], (20102016).

Acknowledgements

V.R.-M. was supported by the European Commission Marie-Curie ITN program (FP7-320 PEOPLE-2011-ITN) through the LINC project (Grant no. 289447). We also acknowledge support from FEDER and Spanish Ministry of Economy and Competitiveness (MINECO) through projects INTENSE@COSYP (FIS2012-30634) and NOMAQ (FIS2014-60343-P), and from the European Commission through project LASAGNE (FP7-ICT-318132). J.J.R. acknowledges funding from the Ramn y Cajal program of MINECO.

Author Contributions

V.R.-M., V.M.E., E.H.-G. and J.J.R. (all the authors) conceived the method, V.R.-M. conducted the numerical simulations and data treatment, V.R.-M., V.M.E., E.H.-G. and J.J.R. (all the authors) analyzed the results, contributed in the writing and revised the manuscript.

Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing nancial interests: The authors declare no competing nancial interests.

SCIENTIFIC REPORTS

9

www.nature.com/scientificreports/

How to cite this article: Rodrguez-Mndez, V. et al. Percolation-based precursors of transitions in extended systems. Sci. Rep. 6, 29552; doi: 10.1038/srep29552 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

SCIENTIFIC REPORTS

10