Early Warning Signals for Critical Transitions: A Generalized Modeling Approach

Critical transitions are sudden, often irreversible, changes that can occur in a large variety of complex systems; signals that warn of critical transitions are therefore highly desirable. We propose a new method for early warning signals that integrates multiple sources of information and data about the system through the framework of a generalized model. We demonstrate our proposed approach through several examples, including a previously published fisheries model. We regard our method as complementary to existing early warning signals, taking an approach of intermediate complexity between model-free approaches and fully parameterized simulations. One potential advantage of our approach is that, under appropriate conditions, it may reduce the amount of time series data required for a robust early warning signal.


Introduction
Critical transitions are sudden, long-term changes in complex systems that occur when a threshold is crossed [1].Many systems are known to be at risk of such transitions, including systems in ecology [2], climate research [3], economics [4], sociology [5] and human physiology [6].Examples of critical transitions in ecology include shifts in food web composition in shallow lakes [7], degradation of coral reefs [8], degradation of managed rangelands [9], and desertification [10].
Warning signals for impending critical transitions are highly desirable, because it is often difficult to revert a system to the previous state once a critical transition has occurred [2,11].If an accurate mathematical model of the system is available then critical transitions can be predicted straight-forwardly, either by numerical simulation or by direct computation of the dynamical thresholds.For real world complex systems, however, sufficiently accurate models are in general not available, and predictions based on models of limited accuracy face substantial difficulties [12].Recent research has therefore focused on model-free approaches that extract warning signals from observed time series [13].Two of the most widely used approaches are increasing variance [14] and autocorrelation [15], both of which are caused by critical slowing down [16].Other approaches consider warning signals based on skewness [17], flickering [18] and spatial correlation [19].
In practice, robust calculation of early warning signals such as these is often limited by the need to acquire sufficient data.Especially in ecological applications, acquiring time series with the desired accuracy and sampling rate poses a significant challenge.One strategy for reducing the demand for time series is to utilize other information that may be available.This highlights the need for intermediate approaches, which can efficiently incorporate available information without requiring a fully specified mathematical model.
In the present Letter, we propose an approach for the prediction of critical transitions based on the framework of generalized modeling [20,21].The approach allows available information to be used to reduce the amount of time series data required or to increase the quality of the prediction.We demonstrate the applicability of the proposed approach by considering a simple one-population model, a previously studied fisheries model and a tri-trophic food chain.

Generalized modeling
Suppose that a system has been identified as being at risk of a critical transition.Even if very little specific information is available, the dynamics can generally still be captured by a so-called generalized model [20].Such a model captures the structure of the system, without restricting it to specific functional forms.
To formulate a generalized model we first identify important system variables (say, abundance or biomass of the populations in the system) and processes (for example, birth, death, or predation).As a first step, the generalized model can then be sketched in graphical form, such as in Fig. 1 below.This graphical representation is sometimes called a causal loop diagram [22].
To obtain a mathematical representation of the model we write a dynamical equation for each variable X i .In these equations we represent the processes by general functions.For instance we can model a single population X 1 subject to gains G and losses L by an ordinary differential equation or as a discrete-time map Note, that in contrast to conventional models, we do not attempt to describe the processes G and L by specific functional forms.Instead, we use unspecified functions G() and L() as formal placeholders for the (unknown) relationships realized in the real system.

Calculation of early warning signal
We assume that before the critical transition, the system can be considered in equilibrium.We emphasize that this does not require the system to be completely static or closed in a thermodynamic sense, but that, on the chosen macroscopic level of description, the system can be considered at rest.Additionally, the system is subject to a slowly changing external parameter that puts it at risk of undergoing a critical transition.The system is therefore at equilibrium only on a certain timescale.In the following we refer to this timescale as the fast timescale, while the dynamics of the whole system, including the slow change of the external parameter, takes place on the slow timescale.Using the definitions above critical transitions can be linked to instabilities (bifurcations) of the fast subsystem [23].For detecting these instabilities we construct the Jacobian matrix, a local linearization of the system around the steady state [24].A system of ordinary differential equations (ODEs) is dynamically stable if all eigenvalues of the Jacobian have negative real parts, whereas a discrete time map is stable if all eigenvalues have an absolute value less than one.Critical transitions are thus signified by a change in the external parameter causing at least one of the eigenvalues to cross the imaginary axis (ODE) or a unit circle around the origin (map).
To warn of impending critical transitions we monitor the eigenvalues of the Jacobian of the fast subsystem, which usually change slowly in time.A warning is raised if at least one of the eigenvalues shows a clear trend toward the stability boundary (for ODEs, zero real part; for maps, absolute value of one).The Jacobian itself can be computed directly from the generalized model, but will contain unknown terms reflecting our ignorance of the precise functional forms in the model.Previous publications [20] have shown that these unknowns can be treated as well-defined parameters with clear ecological interpretations.
In the present applications we estimate the unknowns in the Jacobian matrix from short segments of time series data.Thereby, a pseudo-continuous monitoring of the eigenvalues of the fast subsystem is possible.
We note that with given time series data estimating the generalized model parameters is simpler than estimating the entries of the Jacobian matrix directly, because the generalized model already incorporates structural information about the system.Further, many of the parameters in the generalized model may already be available in a given application, because they refer to well-studied properties of the species, such as natural life expectancy or metabolic rate.

Results
We applied the proposed approach to three case studies, focusing on a generic population with an Allee effect, a fisheries example, and a tri-trophic food chain.

Simulation with Allee effect
Allee effects, that is, positive relationships between per-capita growth rate and population size, are postulated in many populations and have been conclusively demonstrated in some [25].A population with an Allee effect can suddenly transition from a stable, non-zero population size to unconditional extinction [26].We constructed the Jacobian of such a system from a suitable generalized model (see Appendix S1 in Supporting Information).This Jacobian contains unknown parameters that depend on the specific system under consideration.An early warning signal is obtained by estimating these parameters from time series data of the population size and birth rate, and monitoring the eigenvalues of the Jacobian as external parameters change.
For testing the early warning signal, described above, we simulated a simple model based on Yeakel et al. [27].We included a small additive noise term σξ(t), with standard deviation σ and autocorrelation ξ(t)ξ(t ′ ) = δ(t − t ′ ).Throughout the simulation we slowly changed the mortality rate m according to m = 7.5 + 0.2t.A critical transition occurred, causing subsequent extinction of the population (Fig. 2).
The challenge we address is predicting the critical transition from a limited number (here, fifteen) of observations of population size and birth rate.We emphasize that we do not utilize any information on the functional forms of processes employed in the simulation, so that the prediction is based solely on the 15 observations and the assumed structural information (that is, one population subject to gains and losses).By estimating the parameters of the generalized model we determined the eigenvalues of the Jacobian as a function of time (Fig. 2b).A clear increase in the eigenvalue is detectable well before the critical transition, giving ample warning of the impending collapse.
Due to a phenomenon called bifurcation delay [23], the population size did not start to change rapidly until well after (t ≈ 13.5) the bifurcation point (t = 12.5).As previously observed by Biggs et al. [28], management action to reverse the change in bifurcation parameter may successfully avert the critical transition even after the fast subsystem's bifurcation has occurred, if still within the range of the bifurcation delay.In the case of Fig. 2(b), the eigenvalue trend is directed more towards the last possible time that successful management action can be taken than towards the time of the actual bifurcation.

Fishery simulation of Biggs et al.
Our second case study focuses on an example from fisheries.Increased harvesting of piscivores can induce a shift from the high-piscivore low-planktivore regime to a low-piscivore high-planktivore regime [29].Many fisheries are suspected to have undergone such transitions [30,31].Based on the causal loop diagram (Fig. 1), we formulated a discrete-time generalized model, describing the piscivore and planktivore populations at the end of each year, in the spirit of stock-assessment modeling (see Appendix S1).Thereby detailed modeling of the intra-annual dynamics was avoided.
To test the warning signal, we generated time series data with a detailed fishery model by Biggs et al. [28], which was developed from a series of whole-lake experiments [32].We note that this model differs significantly from our generalized model by a) accounting for the intra-annual dynamics and b) containing an additional state variable denoting the juvenile piscivore population.These discrepancies were intentionally included to reflect the limited information that would be available for the formulation of the generalized model in practice.In simulations the detailed model showed a transition to a low-piscivore high-planktivore regime as the harvesting rate was increased (Fig. 3a).
From this simulation, we recorded the simulated piscivore and planktivore density and piscivore catch at the end of each year.Because the simulated data was very noisy we estimated the Jacobian's eigenvalues after smoothing the recorded data (see Appendix S1).In addition to the time series data, the information on the natural adult piscivore mortality and reproduction rate and the planktivore influx from refugia were required (see Appendix S1).This type of information can be reasonably well estimated for most systems.We confirmed that our predictions (reported below) are not sensitive to the specific values used.Indeed, a simple approach for estimating these parameters is to recognize that the initial state, before the critical transition, is stable.In a number of test trials we confirmed that any reasonable combination of parameters used that corresponded to an initially stable steady state provided an early warning signal comparable to the results reported below.
An estimate of the Jacobian eigenvalues for the fisheries example is shown in Fig. 3b.As the system approaches the critical transition we observe that an eigenvalue approaches one, which signifies a critical transition for discrete time systems.This result is compared to the variance-based early warning signal computed by Biggs et al. [28], which uses a much more densely sampled time series including intra-annual dynamics.The comparison shows that the approach proposed here produces a signal of similar quality (although possibly too early), while requiring significantly less time series data.Further, comparison with a variance-based signal using the same amount of time-series data as the generalized model shows that the generalized model-based signal is a much clearer early warning signal in this case.In particular, the variance signal only rises during or after the transition.

Tri-trophic food chain
For our final example we consider a tri-trophic food chain.In ecological theory food chains play a role both as a prominent motif appearing in complex food webs and as coarse-grained models.Using generalized models, a general Jacobian for a continuous-time model of the tri-trophic food chain can be derived (see Appendix S1 and Gross et al. [33]).
For generating example time series data we used a conventional equation-based model in which a producer biomass, X 1 , predator biomass, X 2 , and top predator biomass, X 3 , follow In these equations we assumed Holling type-3 predator-prey interaction, logistic growth of the producer and linear mortality of the top predator and ξ i (t) are additive noise terms with If any biomass decreased to zero, we suppressed the noise term so that the corresponding population remained extinct.We simulated the Eqs.( 2) while increasing the mortality rate of the top predator.The resulting time series, for the chosen combination of parameters, show a slowly changing steady-state followed by a sudden transition to large oscillations, and a sudden collapse of all three populations (Fig. 4).
To provide an early warning signal for the transition we recorded time series of the three biomasses and the top-predator's death rate, and estimated the parameters of the generalized model from smoothed segments of these time series.Even for the smoothed data we find that one of the eigenvalues is very noisy and sometimes positive.We believe that the presence of this eigenvalue reflects the response of the prey to fluctuations on the higher trophic levels and therefore exclude this value from our analysis.As m is increased toward the critical transition at the onset of oscillations, two eigenvalues show a clear increase toward zero real part (Fig. 4).The two eigenvalues approach zero as a complex conjugate eigenvalue pair, which is indicative of the system undergoing a Hopf bifurcation [24].This bifurcation generally implies a transition from stationary to oscillatory dynamics.The large oscillations combined with stochastic fluctuations then lead rapidly to extinction.The early warning signal constructed here thus indicates the transition to an oscillatory state significantly before the transition occurs.

Discussion
In this Letter we proposed an approach for anticipating critical transitions before they occur.In particular we showed that generalized modeling of the system can facilitate the incorporation of the structural information that is in general available.
We demonstrated the proposed approach in a series of three case studies.The first example showed that in simple systems even very few time points can be sufficient for clean prediction of the critical transition.The second example posed a hard challenge, where test data was generated by a model that differed considerably from the generalized model.Yet even in this case the generalized model significantly reduced the amount of data needed for predicting the transition.The third and final example demonstrated the ability of the proposed approach to distinguish between different types of critical transitions.
In all case studies we found that the proposed approach can robustly warn of critical transitions in the presence of noise.We believe that the performance of the approach under noisy conditions can be further improved by subsequent refinements.Such refinements could include combination with dynamic linear modeling [14], utilization of a parameter transformation (to 'scale' and 'elasticity' parameters) previously proposed for generalized models [20], or the use of optimized sampling procedures.
An advantage of the proposed approach is that it is generally becomes more reliable closer to a critical transition, where rates of change of state variables and other observables are generally larger, leading to better sampling.In contrast, methods based on variance or autocorrelation become statistically more difficult to estimate as the time series becomes less stationary.In this sense the proposed approach provides a tool complementary to established variance-based methods.
An important assumption in our present treatment was that the dynamics of the fast subsystem could, at least at some level of description, be considered as stationary.Let us emphasize that this is not a strong assumption because even systems that are primarily non-stationary, such as the fisheries example, can be modeled as stationary if a suitable modeling framework is chosen.Furthermore, ongoing efforts aim at extending the framework of generalized modeling to non-stationary dynamics, which may lead to a further relaxation of that assumption in the future [34].
In summary, we used generalized models to efficiently incorporate available information about a system without requiring detailed knowledge about that system.This led to the construction of early warning signals for critical transitions that required significantly less time series data than variance-based approaches.Thereby the proposed approach can contribute to the warning of critical transitions from a realistic sampling effort.(a) Results for the adult piscivore (blue line, left axis), juvenile piscivore (red line, left axis) and planktivore (green line, left axis) populations of the model of Biggs et al. [28], with the harvesting rate they specified of 1.5 + 0.005t per year.From these results, data were sampled from A and F , once per year, after maturation and mortality for that year had been computed (black markers/thick black line).The annual piscivore catch data used in the early warning signal calculation are also shown (thin black line, right axis).(b) Early warning signals: the eigenvalues (points, left axis) estimated by the generalized modeling approach; the intra-annual variance of the planktivore population (solid line, right axis) as calculated by the method of Biggs et al. [28]; and the planktivore variance using only year-end data (dot-dashed line, right axis).Only one of the eigenvalues is visible at this scale.The eigenvalues were always real.The year-end variance warning signal was calculated using a sliding window of the previous 15 year-end planktivore populations and a quadratic detrending within that window, which yielded generally stationary data for those windows that were before the critical transition.Early warning signal for a critical transition in a tri-trophic food chain.(a) Time series of X 1 (blue circles, left axis), X 2 (green crosses, left axis) and X 3 (red triangles, left axis) generated by Eq. ( 2).Only the data subsequently used in the early warning analysis were plotted.Some of the data during the oscillations between t = 40 and 45 were outside the scale of this graph, with X 3 exceeding 30.The estimates of top-predator mortality used to calculate the early warning signal are also shown (black dots, right axis).Parameters were A n = 4, B = 4, A p = 2, K p = 5, K n = 10, K 3 = 2, (σ 1 , σ 2 , σ 3 ) = (5, 0.02, 0.01) and m(t) = 0.4 + 0.02t.(b,c) Real and imaginary parts of eigenvalues estimated from the data in (a), using the procedure described in the text.Eigenvalues are denoted by dots and circles; a dot within a circle indicates two eigenvalues had the same real value.The markers and colors used in (b,c) have no correspondence to those used in (a).A third (purely real) eigenvalue was not plotted, for reasons described in the main text.In (b), the horizontal dashed line indicates the stability boundary at zero real part, while the vertical dashed line indicates the time that the (Hopf) bifurcation occurred in the fast subsystem of Eqs.(2).

Fishery simulation of Biggs et al.
We assumed the following knowledge: • The populations of adult piscivores A i and planktivores F i at the start of year i • The number of adult piscivores harvested Q i = Q(A i ) during year i • That adult piscivores die from other causes at a per-capita rate m • That adult piscivores predate on the planktivores • That there is a net movement of planktivores into (out of) the foraging arenas from (to) refugia [35] of R i in year i • That per-capita, and in the absence of interactions with other modeled populations, r piscivores would be born and mature into adult piscivores the following year • That however juvenile piscivores are predated on by planktivores [29] and are also controlled, for example by accidental predation, by the adult piscivores [32].
This knowledge is represented schematically in Fig. 1 of the main text.
We constructed a generalized model of this fishery by modeling only the year-end piscivore populations A i and planktivore populations F i with a discrete map.As a consequence we could avoid explicitly modeling the complications of intra-annual dynamics such as piscivore reproduction and maturation.We wrote: The first term of Eq. (5a) gives the number of surviving adult piscivores after harvesting (Q i ) and other causes of mortality (m).The second term rA i models reproduction without losses and the third C i = C(A i , F i ) gives the predatory or controlling effects of planktivores and adult piscivores on the number of new adult piscivores.Since C depends twice on the adult piscivore population, first through the number of juvenile piscivores born and again through the controlling effect of the adults, we assume C is quadratic in A i .
In Eq. (5b), R i models the net influx of planktivores into the foraging arenas.Since the planktivore population was initially small, we assumed there was negligible outflow of planktivores, and that the inflow depended on a refuge size that was constant.Therefore R i = R was constant.[In the other scenario considered by Biggs et al. [28], where habitat destruction causes the refuge size to decrease, R i would not be constant.]The remaining term D i = D(A i , F i ) models all other factors changing the planktivore population, most importantly predation by adult piscivores.Since the planktivore population was initially small, we assumed D to be linear in F i .
We generated data to test this early warning signal with the previously published fishery model of Biggs et al. [28].Their model is a hybrid discrete-continuous system and explicitly models the adult piscivore, juvenile piscivore and planktivore populations.The intra-annual dynamics are continuous, modeling: the harvest of adult piscivores; the predation and control of planktivores and juvenile piscivores, respectively, by adult piscivores; predation on juvenile piscivores by planktivores; and the movement of both planktivores and juvenile piscivores between the foraging arenas and refugia.At the end of each year discrete update rules were applied that controlled the mortality of adult piscivores, maturation of the surviving juvenile piscivores into adult piscivores, and birth of more juvenile piscivores.There was an additive noise term on the planktivore population dynamics.For further details of the model including its mathematical formulation see Biggs et al. [28].
From the output of this simulation we recorded annual adult piscivore density, planktivore density, and piscivore catch.These observations were sufficient to calculate the Jacobian of Eqs.(5), where we use the notation ∂H/∂x ≡ H ′(x) .The eigenvalues λ of the Jacobian, which may be complex numbers, are the solutions of det(J − λI) = 0 where I is the (2 × 2) identity matrix.
Our procedure was as follows.From the assumption of linearity in A, Q can be calculated immediately.The terms C i and D i can be found by solving Eqs.(5).From the assumptions that C is quadratic in A i and D linear in Finally, by solving linear Taylor expansions of the form the derivatives C ′(F ) i and D

′(A) i
were estimated.This generalized modeling approach required estimates of the adult piscivore mortality m and reproduction rate r, and planktivore influx R from refugia.We assumed the fishery manager could estimate m = 0.4 and r = 1.1 from knowledge of the fish species and the fishery.Effective values m = 0.5 and r = 1 can be derived from the actual simulation's equations [28].The planktivore influx R would be harder to estimate in practice.In the absence of such knowledge, the ecosystem manager could make use of the fact that the fish populations were initially stable.This knowledge constrains R to the range 0.1 to 1.5, as values outside these ranges produce an eigenvalue with magnitude greater than one between times t = 0 to 10.We used R = 0.9.Indeed, this knowledge could also be used to constrain m and r to the ranges 0.1 to 0.7 and 0.95 to 1.5, respectively.We found that any combination of m, r and R that gave an initially stable system robustly predicted the later critical transition.
Since the ecosystem approached the critical transition more slowly than the model with Allee effect above, to remove noise we found it necessary to smooth estimates of the gradients.For each new observation at time t i , we fitted a second-order polynomial (since the method eventually results in differences of up to second order) to the record of each observed time series up to that time.Since it is the gradients at the current time that are most important, in the fit we applied a weight to each data point j of exp[(t j − t i )/τ ] with τ = 10.Gradients were then extracted directly from the coefficients of the polynomial fit.This simulated the calculation of early warning signals 'on the fly', as new observations become available.

Tri-trophic food chain
To obtain an early warning signal for this transition, we constructed a generalized model of the tri-trophic food chain where X 1 , X 2 and X 3 are the biomasses of the three species in the food chain, A(X 1 ) is the primary production rate of X 1 , G(X 1 , X 2 ) and H(X 2 , X 3 ) are the rates at which X 2 and X 3 consume X 1 and X 2 , respectively, and M (X 3 , µ) is the total mortality per unit time of the top-level predator X 3 , with the argument µ denoting it is subject to a changing external influence.

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Schematic of the fishery knowledge that was incorporated into the generalized model.

Figure 4 .
Figure 4. Early warning signal for a critical transition in a tri-trophic food chain.(a) Time series of X 1 (blue circles, left axis), X 2 (green crosses, left axis) and X 3 (red triangles, left axis) generated by Eq.(2).Only the data subsequently used in the early warning analysis were plotted.Some of the data during the oscillations between t = 40 and 45 were outside the scale of this graph, with X 3 exceeding 30.The estimates of top-predator mortality used to calculate the early warning signal are also shown (black dots, right axis).Parameters were A n = 4, B = 4, A p = 2, K p = 5, K n = 10, K 3 = 2, (σ 1 , σ 2 , σ 3 ) = (5, 0.02, 0.01) and m(t) = 0.4 + 0.02t.(b,c) Real and imaginary parts of eigenvalues estimated from the data in (a), using the procedure described in the text.Eigenvalues are denoted by dots and circles; a dot within a circle indicates two eigenvalues had the same real value.The markers and colors used in (b,c) have no correspondence to those used in (a).A third (purely real) eigenvalue was not plotted, for reasons described in the main text.In (b), the horizontal dashed line indicates the stability boundary at zero real part, while the vertical dashed line indicates the time that the (Hopf) bifurcation occurred in the fast subsystem of Eqs.(2).