Early warning signs for saddle-escape transitions in complex networks

Many real world systems are at risk of undergoing critical transitions, leading to sudden qualitative and sometimes irreversible regime shifts. The development of early warning signals is recognized as a major challenge. Recent progress builds on a mathematical framework in which a real-world system is described by a low-dimensional equation system with a small number of key variables, where the critical transition often corresponds to a bifurcation. Here we show that in high-dimensional systems, containing many variables, we frequently encounter an additional non-bifurcative saddle-type mechanism leading to critical transitions. This generic class of transitions has been missed in the search for early-warnings up to now. In fact, the saddle-type mechanism also applies to low-dimensional systems with saddle-dynamics. Near a saddle a system moves slowly and the state may be perceived as stable over substantial time periods. We develop an early warning sign for the saddle-type transition. We illustrate our results in two network models and epidemiological data. This work thus establishes a connection from critical transitions to networks and an early warning sign for a new type of critical transition. In complex models and big data we anticipate that saddle-transitions will be encountered frequently in the future.

The low-dimensional systems that are commonly investigated in the context of critical transitions (or tipping points) [1,2] typically settle to stable, but not necessarily stationary, states [3].The defining features of such states is that the system returns exponentially to the state after sufficiently small perturbations [4].When environmental parameters change, a critical transition may occur when thresholds are crossed, where the system becomes unstable to certain perturbations.Low dimensional systems then typically depart quickly from the original state before approaching another, possibly distant, state.The thresholds at which such transitions occur are called bifurcation points [5] and the corresponding transitions are local bifurcations [6,3].Warning signs for such bifurcation-induced transitions detect the loss of exponential restoring dynamics either through its impact on the statistics of noise-induced fluctuations [1,7,8] or by direct measurement of recovery rates [9,10].Thus, bifurcation-induced critical transitions are well understood, and the corresponding warning signs have been analyzed mathematically [11] and tested in experiments [12,13].
For understanding the alternative mechanism of critical transitions that is the focus of the present paper, consider that many real-world complex systems do not reside in stable states.An intuitive example is provided by the outbreak of an epidemic invasion in a population occuring under stable enviromental conditions.Although the precise mechanism for spreading after invasion can be very complex [14], many examples show that introducing a certain new or previously extinct pathogen into an unprepared population can have drastic consequences.This illustrates that the original state of the system, in which the pathogen is absent, was unstable with respect to a specific perturbation, corresponding to the introduction of the pathogen.Note that even after its introduction, the pathogen may be extremely rare for an extended amount of time, residing in small subpopulations or animal vectors, such that on a macroscopic level a disease-free state is still observed for significant time.We may only see a very small rise in the number of infected individuals for a long time under fixed environmental conditions but eventually a drastic jump to an endemic regime occurs.
The example above differs fundamentally from bifurcation-induced transitions, because the qualitative change is not induced by a change of environmental parameters, but rather by a specific 'rare' perturbation.A system is susceptible to such perturbation-induced transitions if it resides in a saddle point, a state that is stable with respect to some perturbations, but unstable with respect to others.In the following we refer to critical transitions caused by the departure from saddle points far from bifurcations as saddle-escape transitions; see also (Fig. 1) for a mathematical normal form example for passage near a saddle; we remark that this example is generic in the sense that mathematical theory guarantees that other systems with nondegenerate saddles show the same dynamics up to coordinate changes and by using a suitable notion of equivalence for the dynamics.
In simple low-dimensional, modular or symmetric (i.e.effective few-variable) systems that are typically studied in the context of early-warnings signs, saddle-escape transitions may occur [15,16] but are relatively rare in practice, particularly in comparison to saddles in high-dimensional (many-variable) systems.To understand this difference, first, consider the abundance of available saddle states.Whether a steady state is an attractor, repeller, or saddle is determined by the eigenvalues of the systems Jacobian matrix, which provides a linearization of the system around the state in question [4,5].A state is a repeller when all of these eigenvalues have positive real parts, an attractor if all eigenvalues have negative real parts, and a saddle when there are eigenvalues with positive real part and also eigenvalues with negative real-part.In a complex heterogeneous system the eigenvalues can, to first approximation, be considered as random variables [17] that have negative real parts with a certain probability q.Since the total number of eigenvalues increases with the number of variables N , the proportion of attractors decreases as q N , the proportion of repellers decreases as (1 − q) N , whereas the proportion of saddles increases rapidly with increasing system size (see Supplementary Information, Section 1).We can therefore expect to find an abundance of saddle states in generic complex high-dimensional heterogeneous systems.
Given the existence of saddles we may also want to to ask whether systems can even approach a single saddle.This is possible, e.g., if a trajectory exists that connects the saddle state to itself or several trajectories connect between different saddles.Such homoclinic or heteroclinic orbits [4,5] exist in lowdimensional systems only if certain conditions are met exactly, e.g.parameters are tuned exactly right.However, if we allow parameters to change dynamically, the dimension of the system increases naturally and long-time homoclinic and heteroclinic dynamics can appear robustly (see Supplementary Information, Section 2, for an epidemic model example).A simple example is the fold-homoclinic (or square-wave) bursting mechanism observed in many neurons [18], where additional macroscopic slow gating variables are frequently introduced to capture the complex processes in the neuron in addition to the usual voltage dynamics.This leads to a robust repeated passage near a homoclinic structure.Furthermore, if there are certain symmetries in the system, then heteroclinic behavior may occur generically as well [19,20].
Saddle-escape can only be considered a critical transition if the system resides in the saddle for a sufficiently long time such that the saddle is perceived as the natural state.The residence time of a system near a saddle can be long as the system moves slowly near steady states (see Supplementary Information, Section 1, for a standard calculation of the residence time).When the system is subject to a small perturbation, its response can be decomposed into different fundamental modes.Some of these modes decay quickly, restoring the system to the steady state.However, at least one mode exists that once excited initially grows exponentially, and thus leads to an escape from the saddle.Here we consider the state where the initial perturbations from those growing modes are very small or entirely absent, since otherwise the system would depart very quickly from the saddle.
In full mathematical generality we would usually expect that a typical perturbation of the system excites all modes to some degree, and thus triggers departure from the saddle point.Based on this reasoning we would expect that already the first perturbation of a system residing in a saddle launches it into exponential departure.However, in real high-dimensional systems the situation is not so simple.For example, hidden symmetries or constraints may exist that effectively decouple certain variables from perturbation of other variables.The most important constraint is the absolute nature of the zero line.Consider again the example of a rare pathogen introduced into a population.Here it is immediately apparent that a perturbation cannot excite the exponential growth of the infected population, unless it already involves the introduction of at least some infected individuals.Furthermore, the perturbation direction to the saddle may initially be very weak in comparison to the strong stable directions, which also leads to long residence times near saddle points (see Supplementary Information, Section 2, for a simple compartmental epidemic model illustrating this effect).Finally, in many models stochastic effects cannot be ignored close to the steady state (see Supplementary Information, Section 2).For instance it is well known that small populations that are feasible in the deterministic limit may still go extinct due to fluctuations in populations size, a phenomenon known as demographic extinction [21].These mechanisms, i.e., micro-level stochastic extinction, slowness of departure and decoupling of perturbation modes, can effectively stabilize the saddle state for a long time until a particular perturbation, or series of perturbations, launches the system on an escaping trajectory.Considered together, the well known-arguments presented above suggest that saddle-escape may play a role as a critical transition, particularly in high-dimensional systems.
Let us now ask whether warning signs for saddle transitions can exist.Because saddle-escape is triggered by the occurrence of a certain perturbation, which is inherently unpredictable in the context of the model, we cannot hope to detect the saddle escape far before the exponential departure from the saddle starts.However, consider that in contrast to bifurcation-induced transitions, saddle escape will generally occur at states that are far from bifurcations.Such states are said to be hyperbolic, and perturbations to these states grow or decline exponentially depending upon the perturbation direction.The departure form the saddle is initially slow and dynamics around the saddle are often perceived as meta-stable.Phenomenologically, the transition therefore appears very similar to a bifurcation-induced critical transition, showing first a slow drift, followed by a sharp spike.Arguably, in many applications it should thus be possible to restore the system to the saddle during the initial phase where the intrinsic dynamics are still slow.Our aim is thus to detect the saddle-escape after the critical perturbation has occurred, but before the system has moved so far from the saddle that the dynamics has accelerated too much.In fact, if we know, a priori, the system is approaching a saddle, not just a fully stable state, then the warning sign described below may help to detect a potential critical transition already at the very beginning of the metastable phase; otherwise, we can apply it during the metastable phase.
In our mathematical treatment we consider a scenario where the system starts from some initial point, then approaches the saddle and stays for a significant time in the vicinity of the saddle point before departing again.For detecting the onset of the departure we exploit the exponential form of the departure trajectory.Between two time points t 1 and t 2 the logarithmic distance along a trajectory close to the saddle point x * will eventually be dominated by a scaling of the form where λ u > 0 is the real part of the largest eigenvalue of the Jacobian and k 0 is a constant (see Supplementary Information, Section 2), i.e., the logarithmic distance increases linearly in forward time.Although saddles generically have positive eigenvalues, the influence of λ u , which can be estimated by Eq. 1, on the overall system dynamics is small as long as the system is approaching the saddle sufficiently close to its stable directions.In this case, the dynamics does not yet involve a significant component in the direction of unstable eigenvectors or the system would not approach the saddle at all.Let us illustrate this again by the epidemics example.While the pathogen is absent the variable that captures the density of the infected population remains fixed to zero, since there are no dynamics in the infected population variable we only measure a stable eigenvalue with real part λ s < 0 via an analogous logarithmic scaling relation as Eq. ( 1) (see Supplementary Information, Section 2), where the logarithmic distance decays linearly in forward time.Only after the pathogen is introduced and the infected population   2).The red/blue linear interpolants yield two approximations for the stable eigenvalue λ s ≈ −1.10, −0.99 and the black lines for the important unstable eigenvalue λ u ≈ 0.51, 0.57; the true values are (λ s , λ u ) = (−1, 0.5).The black squares in (b) can be obtained from x 1,2 ∼ e λut .Note that the choice of B is a choice of sliding window length (or lead time) for prediction as in the case for bifurcation-induced tipping.
starts growing, we start to pick up the positive eigenvalue λ u associated to the exponential growth.
Hence, we may use logarithmic distances between points as a measure to determine, which eigenvalue λ is currently dominating as shown in (Fig. 1 (b)-(c)).If we find negative λ this signals that we are in a regime dominated by a stable direction and approach the saddle, whereas a positive λ signals a departure.Therefore, the emergence of positive λ beyond a certain threshold can provide a warning signal for saddleescape transitions.If we know, a priori, that we are approaching a saddle, not just a fully stable state, then the warning sign already detects the saddle when the logarithmic distance reduction decays linearly.If it could be a stable or saddle state, we can only apply the logarithmic distance reduction warning sign during the metastable phase near the saddle.
To illustrate the saddle-escape mechanism for critical transitions and to test the proposed early-warning signal we consider two recent adaptive-network models in which critical transitions as well as saddle-type behavior play key roles in the dynamics and an epidemiological data set.The first system is a model from evolutionary game theory, in which the evolution of cooperative behavior in a network of interacting, self-interested agents is studied [22].
The model describes a network of agents, connected by social contacts.Each agent pursues one of two possible strategies, which we call cooperate and defect.The agents engage in pairwise interactions with their neighbors, which are modeled as a snowdrift game [23].In this game the highest social payoff in produced by mutual cooperation.However, for the individual agent, defecting yields a higher payoff when the agent is interacting with a cooperator.In time both the agents' strategies and the network of interactions change as agents switch to the strategy that performs optimally in the population, and also rewire their connections to other agents following the more successful strategy (see Supplementary Information, Section 3).
Previous work [22] has shown that in the limit of infinite population size the system robustly approaches the state of full cooperation, where the probability that a randomly drawn agent follows the cooperative strategy is one.In large, but finite, networks a state of almost full cooperation is reached that is disrupted by large outbreaks of defection (see Supplementary Information, Section 3).(red) and the density of defectors x 2 (blue); note that we slowly increase the parameter p in time at a constant rate, i.e., p can be viewed as a time variable.In (a1)-(a2) the minima and maxima of a moving average are shown whereas (a3) shows the actual time series.The vertical dashed curve (thin black) in (a1) indicates the theoretically-predicted transition to oscillations; see (Supplementary Information, Section 3).In (b1) the variances for x 1,2 are calculated using a moving window technique up to the gray vertical lines; note that the scaling of the V 1,2 -axis is 10 −5 .Observe that x 1 does not show a clear scaling law while the scaling of x 2 can be used for predicting the transition from steady state to oscillations using classical variance-based warning signs.The predicted transition point from extrapolating the increasing variance scaling law [11] is marked as vertical dashed line (black) in (a2); note that there is a delay in the Hopf bifurcation point so the predicted critical transition to matches, from a practical viewpoint, the data better than the second-order moment closure theory [22].In (b2) the period T of the oscillation is measured and 1/T is linearly interpolated to approximate the period blow-up [5] point (yellow).This is used to predict the transition point (green) from a periodic to a saddle-type/homoclinic regime; note that this period blow-up is not the saddle-mechanism we focus on in this paper but another new warning sign we just note as an interesting related result.The predicted transition is marked by the dashed vertical line (green) in (a3).
Then we also show the logarithmic distance reduction measured from (a3) in (b3).The ellipses in (b31) indicate the regime where the decay-scaling for the saddle-approach breaks down.Note that the ellipses are there to guide the eye.If one would want to give an explicit warning sign, a threshold for λ u has to be specified, which is not done in this qualitative example.A detailed quantitative analysis of thresholds is carried out for a data set below using ROC analysis.Here we just want to point out the existence of saddles and the qualitative change in the distance reduction near the saddle, i.e. the parts (a3) and (b3) illustrate the main ideas for saddle-escape warning signs.
We now explore whether the logarithmic scaling discussed above is capable of providing early-warning of these outbreaks.We build on full agent-based simulation but our analysis focuses on a pair of observables.
In particular, we study the density of cooperators, x 1 , i.e., the proportion of agents whose current strategy is cooperation, and the density of links (per agent) that exist between cooperators and defectors, x 2 .Note that per-capita densities of populations are also frequently the only natural variables available in data, so they are a natural choice to detect critical transitions and warning signs.
First, let us establish that the outbreaks of defection are indeed triggered by a saddle-escape transition.Although bifurcation-induced transitions are not the main focus of this paper, we provide a brief description of the known bifurcation structure here as a parameter is varied and consider the dynamics as a function of the rewiring p, which measures the relative time scale of structural changes of the network to internal changes in the agents strategy.This rewiring rate was previously identified as a key parameter of the system (see also Supplementary Information, Section 3).
Increasing p at a constant rate from a small initial value, we observe two main dynamical changes.First, a stationary solution turns into stable oscillations, which increase in amplitude and period.Here we observe the classical increase in variance before a Hopf bifurcation [1,11] (Fig. 2 (b1)), which is a good predictor for the transition from random fluctuations to small deterministic oscillations.Second, critical transitions occur at larger rewiring rates, where the cooperator density x 1 (t) drops sharply to lower values before rising again slowly.These transitions are associated to the presence of a homoclinic loop in the system [22], which is attached to the fully cooperative state x * 1 = x * 2 = 1.In (Fig. 2 (b2)), we show that the period T of the oscillations may grow rapidly in finite time, indicating a global bifurcation that gives rise to the homoclinic loop (see Ref. [5] for further background on homoclinic bifurcations).Trajectories close to the homoclinic loop remain near the saddle point for a long time before making a fast excursion.Locally, near x * 1 = x * 2 = 1, this is precisely the situation of saddle-escapes.The logarithmic distances d 1,2 are shown in (Fig. 2 (b3)).As a warning sign to predict the rapid drops in x 1 in the range of p ∈ [0.93, 0.96], we assume that we know an eventual saddle-instability will happen so we just have to look for a change of linear scaling induced by stable directions for d 1 during the phases when x 1 is gradually increasing.These changes can be seen for p ≈ 0.944 and p ≈ 0.957 as predicted by the theory see (Fig. 2 (b3)).Note that the complex saddle-escape dynamics occurs in the regime of relatively large re-wiring rate, when the dynamical process of node update and re-wiring act on similar time scales, while the bifurcation-induced warning signs worked for low re-wiring in a quasi-stationary scenario.
As a second example, we consider a susceptible-infectious-susceptible (SIS) epidemiological model on an adaptive network [24,25].The network consists of susceptible (S) and infectious (I) agents.An infection spreads along S-I-links with probability p per unit time, infectious agents recover with probability r, and susceptible agents try to avoid infectious ones by rewiring S-I-links to S-S-links with a probability that is proportional to the total number of infectious agents (see Supplementary Information, Section 4).
In simulations of this system the number of infectious agents shows distinguished peaks in time, which can be interpreted as epidemic outbreaks (Fig. 3).In (Fig. 3 (b)) the logarithmic distances between consecutive points are shown for the density of I-nodes and S-I-links.Despite the strong fluctuations away from the peaks, both warning signals show three phases after a peak: (1) strong stabilization, (2) plateauor noise-type behavior and (3) a trend towards instability before the next spike (see Fig. 3 (a2)).In fact, similar phases can also be observed for the first model in (Fig. 2 (b3)).
Furthermore, we note that the logarithmic distance increases much earlier for the S-I-links than for the I-nodes, so that monitoring the links between infectious and susceptible agents provides an earlier warning signal for epidemic outbreaks (Fig. 3).This is in accordance with the intuitive idea that knowledge about the contact dynamics among infectious and susceptible agents should allow to predict epidemic outbreaks more easily.
The two models suggest that there are regions in parameter space where saddle-escapes play a key role.However, it is also known that epidemic network models can exhibit bifurcation-induced critical transitions with classical warning signs critical transitions [11,26].However, in those cases one usually assumes that a parameter, e.g. the infection rate, is very slowly varying and this causes the critical transition.This (3) directly motivates the question whether there are data sets available where epidemics are driven by saddle escapes.
Here we focus on repeated measles outbreaks documented biweekly between 1944 and 1966 in 60 cities in the United Kingdom [27,28].The time series of the proportion of infected individuals shows long repeated periods of low disease prevalence interspersed with large, but very short, outbreaks.This strongly indicates that a saddle-type mechanism may be at work (Fig. 4 (a)).We use a receiver-operating-characteristic (ROC) curve [29,30] to quantify the performance of our warning sign for saddle escapes.We briefly recall, how ROC curves are calculated.First, one defines a scalar precursory variable X to be computed from observations before the transitions and considers a threshold δ such that if X > δ an alarm is given.Then the two ratios of correct alarms to the total number of actual tipping events and false predictions to the total number of non-tipping events are calculated for various thresholds δ.This provides a quantitative indicator for the ability of the precursor variable to detect tipping points (see Supplementary Information, Section 6, for more background on ROC curves and their interpretation).
For the ROC curves in our case, we compute a least squares fit of the parameter λ u =: X from equation (1) within a time window of k data points as a precursory variable.We give an alarm for an imminent outbreak at the end of the time window when λ u > δ for some threshold δ.This prediction is considered correct if the disease prevalence exceeds 0.1 within the next 5 data points, which corresponds to an outbreak of at least ten percent of the maximum outbreak coming up within the next 2.5 months.In (Fig. 4  threshold values δ.Both curves lie above the diagonal, indicating that our method is better than a purely random prediction.Moreover, the predictions using positive threshold values δ are, on average, better than the ones using δ < 0, as the former produce a ROC curve farther away from the diagonal.This is a natural result, because δ > 0 corresponds to the detection of an actual instability while 0 > δ ≫ −1 only corresponds to the detection of a weakly stable direction.Clearly, very short prediction times (small k) are problematic as they increase the error in λ u leading to false alarms.Long prediction times (large k) lead to very few correct predictions because the exponentially stable approach to the saddle strongly dominates on long time scales.Although the ROC results show that we can potentially improve the prediction of epidemic outbreaks, the situation is far from ideal.The ROC curve is still quite far from the top-left corner in (Fig. 4 (b)), which would be perfect prediction.Hence, there is substantial work to be carried out to try to increase the practical performance of the warning sign.
Note that, although we demonstrated many factors which strongly indicate a saddle-escape mechanism in measles data, it is unlikely that a test exists which can guarantee detailed knowledge of the underlying dynamical mechanisms based on just a relatively short uni-variate time series of the epidemic.However, there is evidence from several epidemic models, which do display saddle states [31,32].It then remains as an open problem to link our warning sign analysis of measles data to particular models of measles [33,34], which we leave as a challenge for future work.Furthermore, we remark that although the indicator we have used seems reasonably efficient, its robustness will strongly depend on the noise level in the system.
As a general comment, we note that the exact mechanism by which a critical transition occurs may well lie in the eye of the beholder.A given critical transition observed in nature may appear as a bifurcationinduced transition in one model and as a saddle-escape in a different model describing the same phenomenon with different variables.
In summary, we have identified a new type of critical transitions that could be relevant for real-world complex many-variable saddle-type systems.We proposed a simple, intuitive early warning signal for these transitions, and demonstrated the application of this warning signal in two adaptive network models and real-world data.We believe that this warning signal will be useful for detecting saddle-escape transitions in network models and data sets, for instance for ecological and socio-economic systems.Because the underlying mechanism considered here differs from the ones typically studied in critical transitions, the proposed early warning sign is complementary to existing approaches.
It is crucial to point out that significant additional research is necessary to make warning signs, for saddle-escape as well as the classical bifurcation-induced tipping scenarios, more applicable for many realworld applications.The first main issue is to connect results to a more detailed statistical analysis, which is currently work in progress by several groups [30].A second challenge is to link qualitative modelling results better to quantitative warning signs.Indeed, both necessarily depend on each other for saddleand bifurcation-mechanisms.In both cases the predictions are tremendously improved if one knows a priori that the current dynamical evolution of the system may be towards a destabilizing tipping scenario and one just does not know where the system tips, on what time horizon it happens, and which precise mechanism occurs.A third challenge is to keep new potential applications in perspective and there has been significant recent progress in this direction in realistic systems for bifurcation-induced transitions, for example in ecology [10].Various models in neuroscience [18] and in ecology [16] suggest that saddle points can occur, which are bound to be preceeded by the generic logarithmic distance reduction we use as a warning sign.In an ecological context a natural future test case could be bio-invasions [35,36].
It is extremely important to note that the argument here is based upon certain mathematical assumptions to make it rigorous.However, it is strongly expected that if we weaken the assumptions in various ways, we still find saddle points frequently in high-dimensional systems.In fact, let us point out that the idea to characterize instability in large-scale systems using random matrix theory is well-known [38] but is still a topic of very recent interest [39].However, previously one only had the semi-circular law available that required the symmetry of A n or one had to rely on structured matrices, for example certain types of food webs [17].These assumptions are usually too strong for complex dynamical networks, which can be highly heterogeneous and yield unstructured, non-symmetric ODEs.This makes the recent progress on proving the full circular law conjecture important in our context.As discussed above, one expects that the results can be generalized even further to include even larger classes of complex systems.Furthermore, note carefully that even if n is small, one may still have saddles, which may be relevant for the dynamics, i.e., the assumptions we make are sufficient to prove saddle existence with high probability but the assumptions may not be necessary to find saddles at all.

Residence Times
Another main point of our argument is that the systems can spend a much longer time near saddle points than away from them.This leads to metastable behavior near saddle points.This can be illustrated with the simplest two-dimensional case given by the ODEs where λ s < 0 < λ u .Note carefully that it is justified to reduce the dimension of the system after the large dimensionality of the system has led to saddle points; the mathematically rigorous reduction just follows from center manifold theory [4].The system (3) decouples with solution Suppose we start with some x 1 (0) = κ > 0 and want to reach a small neighborhood of the origin with x 1 (T ) = δ ≪ κ.This takes a time T = λ −1 s ln(δ/κ).Viewing T as a function of κ shows that the time increases logarithmically.Therefore, a trajectory spends a much longer time near the equilibrium in comparison to the approach towards the equilibrium.Similarly, we can require a trajectory to start in a small neighborhood of the saddle with x 2 (0) = δ > 0 and end at x 2 (T ) = κ > 0. Then T = λ −1 u ln(κ/δ) and the same arguments apply to show that the initial time spend near the equilibrium is much longer than the escape time.Although we have only worked with a linear system (3), similar conclusions apply for the nonlinear case as long the passage near the hyperbolic saddle occurs sufficiently close to its stable and unstable manifolds [4].

Basics
Locally near a hyperbolic saddle point, which we can assume without loss of generality to be at x = (0, 0, . . ., 0) =: 0, we can work with the linearization so that for some matrix A ∈ R n×n with eigenvalues λ i ∈ C for i ∈ {1, 2, . . ., n} and associated eigenvectors v i .We are going to assume that the eigenvalues are distinct which is generic within the space of matrices.Standard linear algebra gives a coordinate transformation P : R n → R n , x = P y, so that where the k matrices B 1 , . . ., B k are the usual Jordan blocks and P maps the standard basis vectors to the basis {v i }.Since the eigenvalues are distinct we have B j ∈ R or B j ∈ R 2×2 .It is straightforward to observe that the escape near saddles is governed by the weakest stable and the strongest unstable directions.
More precisely, we will only consider at most four eigenvalues λ s , λ s , λ u , λ u where overbar denotes complex conjugation so that It is extremely important to highlight again the logic in the previous derivations: First, we start with a large-dimensional system, where it can be shown that saddle points are frequent.For each hyperbolic saddle point, there are many eigenvalues.However, for the dynamical approach or departure of the saddle point, the dynamics is locally governed by the weakest stable and strongest unstable directions, i.e., stronger stable directions damp out very quickly, while weak unstable directions are generically dominated by the strongest unstable mode.Hence, we may develop a local theory for high-dimensional hyperbolic saddles by focusing on the leading directions in the stable and unstable manifolds, which are generically low-dimensional.

The Planar Saddle
We start with the case n = 2 and λ s,u ∈ R. Setting x = P y gives y ′ = P −1 AP y = By with solution y(t) = y(0)e tB or x(t) = P y(t) so that for vectors v 1,2 that can be calculated explicitly.Since we want to approach the saddle point and stay near it for some significant amount of time we must have that |y 2 (0)| = 0 is small so that a solution starts close to the stable manifold W s (0) = {ρv 1 : ρ ∈ R}.Let • denote the usual Euclidean norm.Observe that the term y 1 (0)e λst v 1 → 0 as t → ∞ and, if y 2 (0) = 0, y 2 (0)e λut v 2 → ∞ as t → ∞.Hence, initially x(t) decreases exponentially until a unique minimum and then x(t) increases exponentially.For times t j > 0 such that y 2 (0)e λut j is small it follows for some where ) is a constant that will not be of relevance here.Observe that (6) allows us to estimate λ s from data.Then we can consider it a warning sign when the logarithm of the distance between points starts to deviate from the linear fit (6).In the regime where |y 1 (0)|e λst v 1 is small a similar procedure allows us to estimate the strongest unstable eigenvalues since then we find where k 2 = ln(|y 2 (0)| v 2 ).From the knowledge of λ u we can predict how rapidly x(t) is expected to grow.In practice, we can estimate the eigenvalues λ s,u from a uni-variate coordinate time series x i (t) by looking at a fixed time point T and a set of K previous times Computing d i (T ) for different times T gives that in different regimes (stable/unstable) we have d i (T ) ∼ λ u,s T + K 2 for some constant K 2 .

Complex Eigenvalues
For the case n = 3 we will again assume that (5) holds and that the complex conjugate eigenvalue pair has negative real part i.e. we consider λ s = a s + ib s , λ s = a s − ib s , λ u with associated real eigenvectors v 1,2,3 .With x(t) = P y(t) the general solution is Writing the solution for x(t) in the basis of the eigenvectors v i gives As before, we are going to distinguish two regimes in the time domain, starting with the assumption that y 3 (0)e λut v 3 is small which yields exponentially decaying oscillations in time series for each coordinate x i .Let T be the time between successive maxima or minima then b s = 2π/T .Furthermore, if t 2 > t 1 > 0 as previously and for a positive constant k 3 .The last equation can then be used to estimate a s as shown in Section 2.2.The case of two complex conjugate eigenvalue pairs is similar and will not be discussed here.

Noisy Saddles
An important question is to consider the influence of noise as natural systems, and in particular the measurement of natural systems, are often well-described by an underlying deterministic system with additional random fluctuations.Consider the standard one-dimensional Ornstein-Uhlenbeck (OU) x = x(t) stochastic process generated by the stochastic differential equation (SDE where W = W (t) is a standard 1-dimensional Brownian motion and the equation is interpreted in the Itô-sense.It is well-known [40] that the solution to (9) and the resulting variance V (t) = Var(x(t)) can be calculated For a < 0 it follows that V (t) → −σ 2 /a as t → ∞ and for a > 0 one gets Generalizing (9) to the simplest possible saddle point yields dx = Ax dt + σ dW (11) where W = W (t) now denotes a standard 2-dimensional Brownian motion and A has two eigenvalues λ s < 0 < λ u .Then the same conclusion as before apply since the entries of the covariance matrix C(t) = Cov(x(t)) are generically linear combinations of two decoupled OU-processes, one stable with asymptotically constant variance for a = λ s and one with diverging variance for a = λ u .Hence one could also attempt to use the scaling (10) to get an estimate for λ u by considering a moving window analysis of the logarithm for the variance.This could be of particular interest in case the more straightforward logarithmic distance reduction method does not work.Most likely this will be the case only for very particular intermediate noise strengths.For small noise the estimator based on distances works quite well as shown in the main manuscript.However, if the noise is too large one never reaches a neighborhood of the saddle point with high-probability so that predictions become impossible anyway i.e. the events acquire a purely noise-induced character.

An Example from Epidemics
Saddle points also appear frequently in many applications.Here we briefly illustrate the dynamics near saddles for a compartmental epidemic model [31,41].The basic model is given by where S, I, R are the number of susceptible, infected and recovered individuals in the population (with T := S + I + R), d and b are per capita death and birth rates, γ is the per capita recovery rate, and ρ is the per capita loss of immunity constant.φ is the main bifurcation parameter of the model and controls the interaction strength for the nonlinear incidence function I 2 S. Assuming d = b to keep the total population constant, introducing the new variables s := S/T, i := I/T, r = R/T, and using the constraint 1 = s + i + r, one arrives at a two-dimensional ODE system [31] i The system (13) has a number of different dynamical regimes depending upon the parameter values.Figure 5 shows a simulation of the dynamics for parameter values b = 1, ρ = 0.132051, φ = 81.88,γ = 1.4.In the (i, r)-phase space plot in Figure 5(a) a saddle point (i, r) = (i, r) has been marked as a dot.The trajectory approaches a neighbourhood of the saddle several times, where it only evolves slowly as discussed above.These long periods near the saddle are then interspersed with several short periods consisting of large epidemic outbreaks.The dynamics in this parameter regime is transient and will settle after a very long time to a stable sink equilibrium.However, before this occurs, the saddle dynamics plays the main role.It is important to note that the observed effect can also occur for the case when the saddle point lies precisely on the zero line.Indeed, if we consider a coordinate change ĩ = i − i * , then the saddle point lies on the zero line {i = 0} and we observe the same saddle escape phenomenon as above.

Cooperation Games on Networks
In this section we give a more detailed technical description of the snowdrift game network.The evolutionary game is defined between agents (nodes) that interact on an adaptive network via the links between them [22].The number of nodes N and number of undirected links K is fixed.However, as described below, the adjacency matrix A = (a ij ) may change in time and interacts with the dynamics.This makes the system and adaptive, or co-evolutionary, network.An agent i interacts with an agent j at a given time step if there is a link between i and j.The interaction takes place via a game between the two nodes.In this game, an agent can have two possible strategies σ i , cooperation C or defection D, which are the two dynamical states of the nodes.The payoff agent i receives from agent j via an interaction is modeled via the snowdrift game [23] with interaction matrix where c represents the cost of cooperation and b the benefit.M 11 represents cooperation of both agents, M 22 defection of both agents and the off-diagonal entries correspond to the mixed cases where one agent tries to cooperate but the other agent defects.The total payoff π i for node i is The network is made adaptive by a probabilistic rule.After the game has been played, choose a link at random.With probability p re-wire this link and with probability 1 − p one of the linked agents adopts the other agent's strategy.The two events of re-wiring and adaptation have to specified in more detail.Define the performance φ(σ) of a strategy σ ∈ {C, D} as where n σ is the fraction of agents using strategy σ.If the strategy adoption event takes place agent j adopts the strategy of agent i with probability and i adopts j's strategy with probability f β (j, i) = 1 − f β (i, j).For a re-wiring event of a link between i and j, delete it and select a random node k.Then the link between k and i is generated with probability f α (i, j) and between k and j with probability f α (j, i).
The main dynamical variables we are interested in are the fraction of cooperators and defectors n C and n D , as well as the link densities l CC , l CD and l DD .Since the number of nodes and links is constant it suffices to restrict attention to a single node density and two link densities.Consider the parameter set α = 30, β = 0.1, b = 1, c = 0.8, N = 50000, K = 500000 where the re-wiring rate p is the primary bifurcation parameter that is varied between p = 0 and p = 1.There are three main dynamical regimes, for small p the densities n CD remain almost constant and there are just finite-size effect stochastic fluctuations.Note that for small p, the network topology is very close to being static.Increasing p gives rise to a supercritical Hopf bifurcation to oscillations in a system, where the population and link densities between different types of agents are taken into account [22], i.e. the linearization at the near-homogeneous steady state has a pair of complex conjugate eigenvalues, which cross the imaginary axis at nonzero speed upon variation of p.This leads to small-scale deterministic oscillations, which grow in amplitude upon increasing p further.For large p the periodic dynamics approaches a nearhomoclinic orbit with saddle-type escape dynamics and long periods of high cooperation values with n C near 1.More precisely, the period of the oscillations increases and long times are spend near a saddle steady state and eventually the periodic orbit limits onto a homoclinic orbit, which means trajectories come extremely close to the saddle steady state.Note carefully that high values of p mean that the network topology can change very quickly, i.e. quickly in comparison to the changes of the dynamical states of the agents.This leads to a highly heterogeneous complex system, which can drastically change its entire topology and dynamical state; for a more detailed description of the dynamics we refer to [22].

Epidemics on Networks
In this section we provide a more detailed overview of the epidemiological network model discussed in the main text.The total number of nodes N and links L is assumed to be fixed.Nodes can be in either in a susceptible (S) or infected (I) state.Links between nodes represent potential transmission routes of the disease.Loops and double-links are not allowed.At each time step an infected node recovers with a probability (or recovery rate) r into a susceptible node.For every SI-link the disease spreads with probability p so that the S node becomes an I node upon infection.This basic dynamics just represents the standard susceptible-infected-susceptible model [42].In addition, susceptibles can try to avoid contact with infected and this is modeled via a probability w of re-wiring an SI-link.In this case, the susceptible node S cuts its link to a node I and establishes a new link to another susceptible node.
As before, we are mainly interested in the node and link densities and their trajectories.It is natural to view the parameters (r, p, w) as bifurcation parameters.We briefly describe the dynamics that have been found in [24].It is observed that cluster formation and degree correlation depend on the re-wiring r, e.g. higher re-wiring rates lead to higher degree correlation.Furthermore, the main bifurcation point upon varying w is the epidemic threshold corresponding to a transcritical bifurcation.Hopf bifurcations, saddle-node bifurcations, oscillations and hysteresis can occur.
The re-wiring mechanism can be refined by introducing awareness of susceptibles to the disease.Let ρ = i/N ∈ [0, 1] denote the infected fraction of the population.Now define w = w 0 ρ where w 0 is a fixed constant.In addition to the bifurcation phenomena observed previously, a homoclinic bifurcation is found upon varying p.The large-amplitude oscillations that occur near the homoclinic bifurcation are of interest for our study of saddle escapes in the main text.
The parameter values used for full network simulation are p = 0.0058, r = 0.002 and w 0 = 0.6 with an initial susceptible density s t=0 = 0.98.The total number of nodes was fixed to N = 10 5 and the total number of links to L = 10 6 .These yield oscillations and their analysis using saddle-type escape dynamics gives the results shown in (Fig. 3).

Measles Data
The data set we use for testing our methods are measles epidemic data recorded in 60 UK cities from 1944 to 1966; the data has been downloaded from the website [43,44].A detailed description and analysis of this data set is given in the two papers [28,27].For each city the total number of cases has been reported biweekly.Depending on the size of the city, measles is expected to occur in endemic cycles (large cities) or in recurrent epidemics with local extinction.It is important to observe that either of those two phenomena depends on the fact that the virus is somewhere in the total host population.Therefore, a low infected density is generically expected to lie very close to the stable manifold of a saddle point when considered in a sufficiently large phase space.
An interesting aspect of the modeling of Grenfell et al. [28,27] is that, although their model is stochastic, they recognize that during the initial phase and throughout the epidemic, the underlying dynamical system seems to be behave deterministically.This is precisely the same behavior one can observe from a model such as the epidemiological adaptive-network model described in the last section.

ROC curves
Here, we briefly review the main idea of ROC (receiver operating characteristic), also called the ROC curve, and give the formal definitions for the general case.Denote the points in a given time series by I j := I(t j ) ∈ R for j = 1, 2, . . .and let Y m be a binary random variable with Next, one defines a precursory variable X m := pre(I k 1 ,k 2 ), where pre : R k 1 −k 2 +1 → R is a mapping, which computes out of the observations a scalar-valued precursor.Obviously the variable X m can depend upon the choice of k 1 , k 2 and on further parameters (this problem is currently being studied by the first author and several colleagues for the case of B-tipping).We give an alarm when X m > δ for the event some δ ∈ R. The precursory variable enables us to calculate the rate of correct predictions as well as the rate of false positives r c = #correct predictions #events/outbreaks and r f = #false positives #non-events .
Both rates will obviously depend upon k 1 , k 2 , δ and the dependence will be upon δ only if the window length is fixed.One may also write r c and r f by using aposterior probability density functions as [29] r c = {Xm>δ} P(X m |Y m = 1), The ROC-curve is a plot of the rates in the (r f , r c )-plane for different values of the threshold δ; see (Fig. 4) for an example.There are several important standard observations about ROC curves.Perfect prediction occurs if no false positive and all true events are detected.This implies that the ROC curve should consist just of the point (r f , r c ) = (0, 1).The point (0, 0) in the lower-left corner of an ROC curve represents a value of δ that is so high that no alarm is given at any point while the point (1, 1) at the upper right corner represents when alarms are given at every time step.The diagonal connecting (0, 0) and (1, 1) is precisely, where true positive rate equals the false positive rate, which is equivalent to making random guesses.A precursor with performance better than random guesses corresponds to a point in the upper triangle with r c > r f .

Figure 1 :
Figure 1: Dynamics near a planar saddle with small noise.(a) Phase space (x 1 , x 2 ) with a trajectory (black) passing near the saddle point.The stationary state (dark green dot) and a circle (gray) of radius r = 0.5 indicate a neighborhood of the stationary state outside of which the trajectory is shown as a dashed curve.(b) Time series for x 1 (red) and x 2 (blue).The gray vertical lines indicate entry and exit to the ball B = {x ∈ R 2 : x 2 1 + x 2 2 < r 2 }.The black squares are predicted values from the warning signals obtained inside B. (c) Plot of the logarithmic distance reduction d(T ) as crosses;(see Supplementary Information, Section 2).The red/blue linear interpolants yield two approximations for the stable eigenvalue λ s ≈ −1.10, −0.99 and the black lines for the important unstable eigenvalue λ u ≈ 0.51, 0.57; the true values are (λ s , λ u ) = (−1, 0.5).The black squares in (b) can be obtained from x 1,2 ∼ e λut .Note that the choice of B is a choice of sliding window length (or lead time) for prediction as in the case for bifurcation-induced tipping.

Figure 2 :
Figure 2: Critical transitions for an evolutionary game.(a) Time series for the density of cooperators x 1(red) and the density of defectors x 2 (blue); note that we slowly increase the parameter p in time at a constant rate, i.e., p can be viewed as a time variable.In (a1)-(a2) the minima and maxima of a moving average are shown whereas (a3) shows the actual time series.The vertical dashed curve (thin black) in (a1) indicates the theoretically-predicted transition to oscillations; see (Supplementary Information, Section 3).In (b1) the variances for x 1,2 are calculated using a moving window technique up to the gray vertical lines; note that the scaling of the V 1,2 -axis is 10 −5 .Observe that x 1 does not show a clear scaling law while the scaling of x 2 can be used for predicting the transition from steady state to oscillations using classical variance-based warning signs.The predicted transition point from extrapolating the increasing variance scaling law[11] is marked as vertical dashed line (black) in (a2); note that there is a delay in the Hopf bifurcation point so the predicted critical transition to matches, from a practical viewpoint, the data better than the second-order moment closure theory[22].In (b2) the period T of the oscillation is measured and 1/T is linearly interpolated to approximate the period blow-up[5] point (yellow).This is used to predict the transition point (green) from a periodic to a saddle-type/homoclinic regime; note that this period blow-up is not the saddle-mechanism we focus on in this paper but another new warning sign we just note as an interesting related result.The predicted transition is marked by the dashed vertical line (green) in (a3).Then we also show the logarithmic distance reduction measured from (a3) in (b3).The ellipses in (b31) indicate the regime where the decay-scaling for the saddle-approach breaks down.Note that the ellipses are there to guide the eye.If one would want to give an explicit warning sign, a threshold for λ u has to be specified, which is not done in this qualitative example.A detailed quantitative analysis of thresholds is carried out for a data set below using ROC analysis.Here we just want to point out the existence of saddles and the qualitative change in the distance reduction near the saddle, i.e. the parts (a3) and (b3) illustrate the main ideas for saddle-escape warning signs.

Figure 3 :
Figure 3: Epidemic outbreaks and prediction in an adaptive SIS model.(a) Normalized time series for the infected density I (red) and the susceptible-infected link density SI (blue).(b) Logarithmic distances for I and SI are shown as well (see Supplementary Information, Section 2).The linear interpolations (black) indicate the expected linear upward trend before a saddle-escape; the slopes of the four black lines (from left to right) are approximately 7.364, 12.516, 5.466 and 3.461 respectively.The three ellipses in (a2) highlight the three typical regimes between spikes discussed in the text and are there to guide the eye as in (Fig.2).Parameter values for this figure are p = 0.0058, r = 0.002 and w 0 = 0.6.

Figure 4 :
Figure 4: Measles epidemics in the UK between 1944 and 1966.(a) Typical time series of a city (here: Birmingham) for the infected population I; the series has been normalized by the maximum outbreak.(b)ROC curves for a five (dots) and ten (stars) data point prediction averaged over all cities.The diagonal is shown as well.Blue corresponds to a precursor volume with δ > 0 and red to a precursor volume with δ < 0. A prediction time window of 5 months is indicated by 'stars' while a time window of 2.5 months is indicated by 'dots'.

Figure 5 (
b) shows a time series of the infected population density i corresponding to the trajectory from Figure5(a).

Figure 5 :
Figure 5: Numerical simulation for the compartmental epidemic model (13).(a) Phase space plot with saddle point (green dot) and a trajectory segment (black curve) showing multiple epidemic outbreaks with passages near a saddle.(b) Time series for the infected population density i corresponding to the trajectory from (a).The i-value of the saddle point is indicated by a dashed green line.The long passages near the saddle between larger epidemic ourbreaks are clearly visible on this time scale.

Y m := 1 a
tipping/event occured at time m 0 no tipping/event occured at time m, i.e.Y m just records whether a tipping point occured at time m or not.Y m is taken as a random variable since we do not a priori when events occur.Consider some subset of previous observationsI k 1 ,k 2 := (I m−k 1 , I m−k 1 −1 , . . ., I m−k 2 ) ∈ R k 1 −k 2 +1 with 0 < k 2 < k 1 ,where k 1 − k 2 + 1 is also referred to as the sliding window length or observational window length.