Systematic analysis and optimization of early warning signals for critical transitions using distribution data

Summary Abrupt shifts between alternative regimes occur in complex systems, from cell regulation to brain functions to ecosystems. Several model-free early warning signals (EWS) have been proposed to detect impending transitions, but failure or poor performance in some systems have called for better investigation of their generic applicability. Notably, there are still ongoing debates whether such signals can be successfully extracted from data in particular from biological experiments. In this work, we systematically investigate properties and performance of dynamical EWS in different deteriorating conditions, and we propose an optimized combination to trigger warnings as early as possible, eventually verified on experimental data from microbiological populations. Our results explain discrepancies observed in the literature between warning signs extracted from simulated models and from real data, provide guidance for EWS selection based on desired systems and suggest an optimized composite indicator to alert for impending critical transitions using distribution data.


Classification as critical transitions
What are the best eary warning indicators near a critical point? The results are tested on a wide class of models and empirical data

INTRODUCTION
The dynamics of many complex systems is characterized by critical thresholds (tipping points) and abrupt shifts between alternative regimes. 1,2 Various examples have been observed in diverse research fields and include collapses of ecosystems, 3,4 sudden climate shifts 5,6 or financial crashes. 7,8 Abrupt regime shifts have particularly been theorized and observed in systems biology and medicine, [9][10][11] at the onset of certain disease states like atrial fibrillation 12 or epileptic seizures, 13 as well as in biological processes like regulation of gene networks 14,15 and cell fate decisions, 16,17 including epithelial-mesenchymal transitions. 18 Correctly detecting and alerting for these critical changes allows us to better understand complex developments and to anticipate dangerous outcomes. However, many such complex systems have not been fully characterized with mechanistic models, thus requiring simpler and more generic approaches to support data-driven estimates.
The critical transitions (CT) framework have been proposed to address tipping points using low-dimensional systems descriptions 19 and associated early warning signals (EWS), computed from statistical indicators extracted from data like increasing variance, autocorrelation or coefficient of variation. 20,21 These signs and derived indexes, 22-24 in principle generic for broad classes of systems, have been tested and applied on biological, epidemiological and medical data with alternate success. 25-28 Therefore, recent studies have recommended caution when attempting predictions based on EWS. [29][30][31] Since there is an increasing interest for EWS in systems biology and biomedicine, it is thus compelling to provide a unified framework for the analysis and interpretation of such indicators, to determine in which cases they can be safely applied and to understand their limitations, in particular when considering the type of data that are usually available from systems biology experiments. In addition, going beyond univariate indicators will improve their performance in detecting and alerting for impending critical transitions.
In this work, we provide a systematic analysis of the CT framework and its associated EWS, to define their range of applicability and understand why discrepancies have been observed between theoretical predictions and experimental data. 32,33 Systems biology is characterized by two main paradigms 34 : one investigating the single details of molecular combinations or regulatory networks, alike to ''microstates'' in statistical mechanics, 35 and another looking for general analytical models, built upon kinetic theories, to Box 1. Classification of critical transitions Consider a dynamical system whose state (or regime) is usefully characterized by a set of dynamic variables x˛ℝ m , whose relations to each other are modeled by a set of parameters p˛ℝ n : where F : R n+m /R n is a system of sufficiently smooth functions. If p is not explicitly dependent on time, the system is termed autonomous; if p = pðtÞ, the system is called non-autonomous. The distinction between autonomous and non-autonomous can be supported when considering naturally fixed parameters, 38 or when addressing timescale separation (''slow-fast system'') between biochemical processes, like mRNA transcription versus protein degradation times. 39 This results in sets of dynamical (for variables) and algebraic (for parameters, termed at quasi-steady state) equations. 40 Together, variables and parameters define and shape a state space (or ''landscape'') that, if Fðx; pÞ has elements of non-linearity, can be characterised by multiple attractors, 41 i.e., region of stability for system's states. If parameters are allowed to change (either non-autonomously, or at quasi-steady state), the state space is dynamic and attractors can change, as opposed to static landscapes like Waddington's.
The state space can be multidimensional. However, near bifurcation points, it can be aptly described using lowdimensional models associated with critical thresholds in the values of leading parameters (usually corresponding to the largest eigenvalues 42 ). Such models are termed ''normal forms'' of a dynamical system, simplified minimal-order forms that determine the system's behavior and retain universal properties of generic bifurcations (see Kuehn et al. 43 and STAR Methods). Normal forms can be inferred from bistability properties 14 or deduced from network models, if they are available for the considered systems. 44,45 In addition to bifurcation points, noise can characterize the system's dynamics. Noise is ubiquitous in biology 46,47 and can correspond to stochasticity in intrinsic biochemical processes or cell-cell variation. 48 Mathematically, noise variables can be modeled as fast degrees of freedom augmenting system ((1)), which is a dualistic representation to stochastic processes. 49 Noise can push the system out of original attractors onto new ones, therefore causing random switches between phenotypic states even in the absence of dynamical bifurcations.
We propose to use the relative timescales between dynamical variables, parameters and noise to develop a systematic classification of transitions between system states. This way, we synthesize and improve the contributions of Thompson et al., 50 Kuehn et al., 19 Ashwin et al., 51 Shi et al. 52 toward the establishment of a theory on critical transitions in real systems. To do so, let us extend and disentangle Equation 1 to explicit the dependencies on state variables xR m and system parameters p˛R n , on the introduced stochastic variables x˛R l and on the relative timescales modeled by time parameters t i ; i = fx; p; xg. This results in a multiscale slow-fast system 8 > > > > > > > < > > > > > > > : Using this representation, tipping systems can be classified into three main classes of critical transitions on the basis of relative timescales: bifurcation-induced (''b-tipping''), noise-induced (''n-tipping'') and rate-induced (''r-tipping''), following the nomenclature introduced by Ashwin et al. 51 If t x > t x , the system becomes ergodic and visits the full state-space uniformly without displaying transitions. 52 The b-tipping class thus encompasses all those transitions primarily driven by bifurcations, i.e., slow changes in control mechanisms modeled as quasi-steady approaches of leading parameters to their threshold values. They modify the attractor landscape, in the presence of low noise-to-signal ratios, and can be further sub-classified according to dimension m and co-dimension n. 50 In this work, we only consider low-dimensional ones, commonly found in cell dynamics studies. Examples include toggle-switch mechanisms for the lac-operon, 53 population collapses of microbiological colonies past threshold concentrations of stressors or nutrients, 54 or epithelial-mesenchymal determination. 55 Higher m and n yield more complex bifurcations associated with, e.g., neural network activity. 56 The n-tipping class groups various transitions driven by stochastic fluctuations on fixed landscapes, including large, impactful and unexpected events (sometimes called ''dragon kings'' 57 ). Example range from enzymes crossing activation chemical barriers via ''promoting vibrations'', 58  iScience Article understand complicated biochemical processes in simpler and general terms. 36 The latter allows us to construct classes of systems according to universal routes of dynamical development, regardless of the microscopic details. We leverage this paradigm to make sense of critical transitions and identify the most relevant classes pertaining to biological systems. Box 1 provides such classification, comparing critical transition to smooth transitions, illustrating a connection between mathematical modeling and empirical observations in the field. 37 In the rest of the article, we also provide guidance for EWS selection and optimization, depending on realistic noise properties and other notable features of classes of complex systems, developing new composite indicators.
Our work bridges mathematical insights and observations of real systems to classify various tipping mechanisms. There are ongoing debates whether regime shifts in biological systems, like cell-fate decision, are primarily driven by deterministic bifurcations 64,65 or by random fluctuations, 35,66 which prompted several authors to question the old ''Waddington landscape'' interpretation. 37 By systematically analyzing known regime shifts, we classify the mathematical models to address various types of critical transitions, subject to combinations of bifurcations and noise, 49 and to develop a method to extract systems' robustness proxies from data. For clarity, this topic is covered in Box 2.
For the subsequent analysis of EWS in settings that are typically observed in systems biology studies, we first employ a framework based on dynamical manifolds, underpinning universal routes to explosive transitions. 43 This characterizes the warning signals associated with ''noisy'' bifurcations, and studies their dependency on noise properties and other dynamical features like rapid approaches to threshold values. This way, we provide general results about EWS robustness and sensitivity to dynamical features, to guide applications on various systems, understand their limitations and promote future developments.
Then, we focus on a critical transition sub-class of high biological relevance, the stochastic saddle-node bifurcation. 36 For this tractable, yet realistic model of complex biological processes, we develop a composite EWS indicator to optimize the leading time of the alerts, i.e., how much in advance reliable signals are triggered, with respect to an impending transition. The new indicator is optimized over realistic noise types using the common genetic toggle switch model, 15 as representative of the considered CT class. This way, we overcome the limitations of other EWS from literature, which have mostly been developed over Gaussian noise whilst biological systems usually feature correlated and state-dependent noise. 48,76,80 Thanks to this extension, the indicator also provides additional insights about the systems under investigation, such as inference of noise type from data. The theoretical results are finally tested and verified on publicly available experimental data, demonstrating their potential for monitoring and interpreting diverse systems.

Robustness of EWS for noisy bifurcations
Within the class of critical transitions induced by bifurcations characterized by small fluctuations, discussed in Box 1 and 2, we study the EWS associated with impending tipping points, considering different noise types that are better representative of biological dynamics that additive Gaussian noise (see Box 2). Albeit any type of multiplicative noise can in principle be considered, we focus on the cases hðx;pÞ = x, hðx; pÞ = ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi f ðx; pÞ p and hðx; pÞ = f ðx; pÞ as the most commonly employed in the field, see Box 2 for explanation.
B-tipping and n-tipping directly link to the aforementioned debates in systems biology about deterministic or stochastic drivers of critical changes. R-tipping refers to critical ramping of control parameters, not coped by the system, which has been so far observed in climate 60 and engineering 61 systems. The heat-shock response of plants to ramping temperature conditions 62 may fall within this class, but further studies are required. The critical transition classes can be visualized on bifurcation diagrams or using quasi-potential landscapes, 63 which can be obtained as integrals of vector fields like Equation 1 or inferred from data. Figure 1 shows the classification between the transition classes, with illustrative examples of what can happen to systems within simplified attractors. Note that the hard-cut classification derives from the mathematical assumptions in Equation 3: gradients between the transition classes may exist and call for deep investigation. In particular, our work focuses on ''noisy bifurcations'', i.e., dynamics characterized by bifurcation points and the presence of low to moderate noise-to-signal ratio.

OPEN ACCESS
iScience 26, 107156, July 21, 2023 3 iScience Article Analytic expressions for key summary statistics indicators can be obtained from Equation 6 using standard approaches for stochastic processes. 74,79 Their behavior as the control parameter changes provides EWS for approaching noisy bifurcations. 1 The lag-t autocorrelation function does not depend on h 2 ðx; pÞ but only on jv x f ðx s ; pÞj = k: Hence, the common indicator lag-1 autocorrelation (AC(1), with t = 1) only depends on the dampening rate. The power spectrum of the Fourier transforms and the variance, two common indicators, explicitly depend on h 2 ðx; pÞ: As discussed by Bury et al., 81 the power spectrum can provide an increasing signal over time-series data thanks to its quadratic scaling fk À 2 , but it requires high sampling frequency that is often not available. x and p may also correspond to network combinations of observable variables and parameters. 37,44 (A) Illustration of b-tipping and n-tipping using a bistable system with saddle-node bifurcations (unstable branch in red; saddle-node template shown in inset). Hysteresis can occur, i.e., asymmetric routes to tipping from one stable state or from the other (orange, from up to down with increasing p; black, from down to up with decreasing p). In b-tipping, the system (depicted as a red ball) approaches the bifurcation point. The associated landscape is molded by p and the basin of attraction becomes shallower (as visualized by the bars) until disappearing; there, the system tips. N-tipping: if subject to strong fluctuations (dashed arrow), the system ''fluctuates'' (dashed ball) and can overcome the barrier onto an alternative attractor, even if far from the bifurcation point.
(B) Illustration of r-tipping using the landscape representation: rapid ramping of the control parameter makes it as if the landscape shifts and the systems does not manage to move along, therefore tipping onto another attractor ''sliding'' underneath. Here, hysteresis may not be present, e.g., if the landscape is symmetric. See Ashwin et al. 51 for formal definitions.
(C) For contrast, an example of ''smooth'' transition without hysteresis is provided using a dynamical system close to a pitchfork bifurcation (inset) as template.
that is, the noise-to-signal ratio is not negligible but the slow-fast condition between variables and parameters still applies.
For this class, normal forms can be used to analytically study systems' robustness and derive EWS for impending tipping points. 19 Normal forms are general and low-dimensional models _ x = f ðx; pÞ that describe topologically equivalent systems within a bifurcation class, in the vicinity of critical points. 42 They allow to extract analytical and generic results for wide classes of systems, 43 at the price of neglecting homeostatic dynamics far from tipping points. As a result, they allow to focus on critical transition mechanisms across various systems, instead of studying the full evolution of a single system. Details about topological equivalence and construction of normal forms are in STAR Methods. Figure 2 shows an example of reduction to normal forms for two simple models.
Here, we consider those normal forms of primary biological interest. The saddle-node bifurcation, often associated with population collapses 1,26 or biological state transitions, 67 is defined by f ðx;pÞ = GpGx 2 . At p = 0, a stable (x s = ffiffiffi p p ) and unstable ðx u = À ffiffiffi p p Þ branch collide and vanish, resulting in a critical transition to an alternative branch (if it exists).
Transcritical bifurcations f ðx; pÞ = px À x 2 are characteristic, for instance, of epidemic outbreaks. 28 Here, the two equilibria x 1 = 0 and x 2 = p meet at p = 0 and exchange stability. Finally, the family of pitchfork bifurcations f ðx; pÞ = px + l x 3 describe branching processes from one to two states (or vice versa); l > 0 identifies subcritical bifurcations, associated with critical transitions, while l < 0 defines the supercritical case, with a continuous transition over mean values. This mechanism is identified in cell regulation processes. 37 Stochastically forced systems, associated with ''noisy b-tipping'', can be written in the Itô form 50 where dW is a Wiener process with variance s and f ðx; pÞ is a suitable normal form from those described above. The term hðx; pÞ allows us to represent different noise types, to reflect modern knowledge of stochastic processes occurring in biological systems. Additive Gaussian noise with hðx; pÞ = 1 is usually associated with extrinsic cell-cell variability. State-dependent (multiplicative) noise hðx; pÞsconst represents intrinsic stochasticity determined by, e.g., reaction rates, timescales or species concentrations of the underlying biochemical processes. 68 Combinations of additive and multiplicative noise, with various ratios depending on different systems, are more realistic 69,70 and fit experimental data better than Gaussian noise. 71 Another interesting case is colored (time-dependent) noise 32,72 ; however, combining additive and multiplicative noise processes is usually considered a valid modeling alternative to reproduce biological stochastic dynamics. 69 Hence, this paper focuses on the noise combinations, while the reader is referred to the previous publications for the specific case of colored noise. If the microscopic kinetics is known, the noise terms can be exactly derived from the Master equation using Gillespie formalism. 73 Alternatively, a diffusion approximation 74,75 derives noise terms proportional to system state (hðx; pÞ = x), or to the drift term of Equation 5, hðx; pÞff ðx; pÞ. Here, for multiplicative noise, we consider hðx; pÞ = ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi f ðx; pÞ p 68 and hðx; pÞ = f ðx; pÞ, to reflect modeling of biological regulatory circuits. 76 This way, mechanistic and stochastic normal-form bifurcation models are examined to study the effects of intrinsic and extrinsic noise on statistical patterns of variability and related EWS.
Following the procedure detailed in STAR Methods, Equation 5 is analyzed by solving the slow dynamics, linearizing around a trajectory inside the stable (attracting) manifold and changing the coordinates to highlight the residuals yðtÞ around the linearization. This procedure gives wherex s corresponds to the attracting part of the critical manifold (stable solutions). The linearised drift term corresponds to the leading eigenvalue of the deterministic normal form. Its magnitude is the asymptotic decay rate of a perturbation. It corresponds to the concept of engineering resilience, 77 which is akin to that of robustness. 78 A change of notation jv x f ðx s ðtÞ; tÞj = k makes explicit that Equation 6 corresponds to a (possibly non-autonomous) Ornstein-Uhlenbeck process, with critical k given by k 0 = 0. It is a well-studied problem in stochastic processes theory, with analytical solutions for its statistics in different regimes. 74,79 Equation 6 can be regarded as a first order autoregressive model. However, its derivation from normal forms allows more nuanced interpretation: rather than being hypothesized as a statistical model to capture simple relationships, it is general for all models that can be reduced to normal forms. where Pðy 0 Þ is the probability density function from the associated Fokker-Plank equation 79 and m is the expected average value. Skewness and kurtosis, sometimes suggested as indicators for EWS, 82 can be easily extracted from Equation 11 as third and fourth moments (n = 3 and 4). Entropy-based indicators are more challenging to derive in case of multiplicative noise, as their defining integrals may not be solvable. Their derivation in case of Gaussian noise is described in STAR Methods; for the other cases, their behavior is estimated below using computer simulations.
In all cases, the analytical results for each normal form can be obtained by substituting the corresponding dependency of the drift term to the control parameter: for the saddle-node, k = 2 ffiffiffi p p , for the transcritical k = p and for the pitchforks k = 2p. In Figure 3, the effect of multiplicative noise on the trends of common indicators is shown using hðx; pÞ = x and hðx; pÞ = x 2 . Figure 3 shows expected trends of common statistical indicators, for the three main normal forms and different noise types. The figure is derived by solving the equations above and substituting the corresponding noise functions h 2 ðx; pÞ. Although the scaling induced by differs, the qualitative trends are conserved across the bifurcations. This observation suggests genericity of EWS, but also difficulties to infer the existence of one or another bifurcation using statistical indicators alone (see also Kefi et al., 83 Boettiger et al. 84 ). Other methods (e.g., Angeli et al. 14 ) are recommended to complement the inference.
For Gaussian noise, EWS are associated with increasing trends of statistical indicators. 1,85 However, multiplicative noise may alter or completely disrupt them (as also noted by O'Regan et al. 68 ), resulting in no early warnings prior to tipping points. Equation 8 shows that even power spectrum trends can be subject to alterations from expected patterns, potentially resulting in spurious signals.
A preliminary investigation on ramping parameters 86 can also be conducted. In this case, t x x t p : the quasisteady-state (stationary) assumption is relaxed, but r-tipping may not yet occur. Let us consider linear ramping as k = k 0 À at, where k 0 is any initial condition, a is a small rate coefficient and the ramping stops at the critical value k = 0. Both coefficients are set to 1 to represent commensurable time scales. Only Gaussian noise is considered. This is a particular case of inhomogeneous processes 79 Figure 3 shows that trends of common indicators may be modified by commensurable time scales of parameters evolution. Hence, raising reliable alerts becomes more challenging.
Overall, this analysis demonstrates that theoretical EWS because of increasing trends of summary statistics are sensitive to the ''dynamical context'', i.e., noise properties and reciprocal time-scales. Hence, if the dynamical context is not carefully accounted for, spurious signals may be extracted from data, as observed in early findings from single systems. 28,87 If the context is known, the current results suggest which indicators to use to obtain robust early warnings. The autocorrelation is robust to whether noise is additive or multiplicative; the variance is more sensitive to multiplicative noise, but maintains its expected trends in case of ramping parameters. The coefficient of variation, calculated by using its definition Equation 10 together with Equation 13, and solving it using the Mathematica software, is also robust in case of commensurable time scales and copes well in case of certain types of multiplicative noise. Overall, what matters is the competition between changes in noise and changes in resilience: depending on which one is more rapid, the indicators and their associated EWS may perform as expected or fail to anticipate an impending critical transition.
Measurement processes or details of realistic models may further influence EWS. Measurement uncertainties, assumed as Poisson processes associated with measuring instruments or procedures and thus Here, other multiplicative noise forms may alter its behavior and shadow possible early warnings. Finally, skewness and kurtosis calculated from Equation 11 display increasing trends when Pðy 0 Þ is symmetric (STAR Methods). However, this may not be true in case of multiplicative noise, 15

Optimization of EWS
Having assessed in which cases the proposed EWS are expected to work for noisy b-tipping transitions, we now optimize their performance to provide significant and as-early-as-possible alerts, in a range of dynamical contexts and for the most common transitions observed in systems biology. To this end, we focus on multistable systems, 55 develop and solve an optimization problem using computer simulations to go beyond the first-order approximation from Equation 6 (see STAR Methods for details), and study a wide range of noise levels and types, to establish a composite indicator that is robust and performing across multiple systems.
Multistable systems are systems whose deterministic landscape features at least two attractors, 88 and usually undergo either saddle-node bifurcations or n-tipping. Bistability means local multistability across two attractors. Angeli et al. 14 provide necessary and sufficient conditions for bistability in a wide range of biological systems. Among them, a feedback model with three-points I/O characteristic curves suffices. A simple linear system with monotonic sigmoidal feedback can do the job, in a range of parameters (Figure 4). As a case study, the autocatalytic positive feedback loop derived from Michelis-Menten kynetics 15 _ x = f ðx; cÞ + hðtÞ = K + c x k 1+x k À x + hðtÞ : ( Equation 15) satisfies the bistability conditions, and can thus display transitions between attractors, if 0 < K < 1=ð3 ffiffiffi 3 p Þ for k = 2. 89 In Equation 15, x is the concentration of a transcriptional factor activator, activating its own transcriptions when bound to a responsive element; K is the basal expression rate, c is the maximum production rate, k is the Hill coefficient and hðtÞ accounts for the stochastic terms. Equation 15 comes from a two-variable genetic toggle switch, assuming slow-fast timescale separation between the two variables 90 and after a-dimensionalizing the chemical details to retain the dynamical scaffold. Notably, networks of Michelis-Menten regulators can be reduced to Equation 15 after dimension reduction techniques. 44 Equation 15 displays bistability for a range of values c (the exact range depends on K and k 91 ) and, in particular, a saddle-node bifurcation between two alternative steady states at a critical value c 0 of the parameter c, such that vf =vxj ðx;c0Þ = 0: ) where x 0 is the tipping value for the system state. Therefore, system (15) can be used as a paradigmatic example of biological systems, within the saddle-node b-tipping class, to perform optimization studies that go beyond the local and low-noise-to-signal-ratio approximation provided by normal forms. where a weights the additive or multiplicative noise component (a = 1 corresponds to additive Gaussian noise, a = 0 to solely multiplicative); like above, hðxÞ = x or hðxÞff ðxÞ = x k =ð1 + x k Þ 76 and dW is a Wiener process with variance s. Without loss of generality, 91 we set k = 2.
As EWS are associated with increases of statistical indicators, we need to establish a measure of statistically significant increase, to rule out false positives and false negatives because of random fluctuations in the indicators. To do so, we employ the p value analysis used in Proverbio et al. 91 (see STAR Methods for details). It allows us to measure at which value of the control parameter c, before c 0 , a significant signal is triggered, thus obtaining a ''lead-parameter'' c I sig ðs; aÞ depending on noise properties and the considered indicator I (see STAR Methods for details). c I sig ðs; aÞ is first computed for each indicator individually. Figure 5A shows the results in case of additive noise, while various functionals of multiplicative noise hðxÞ (with a = 0) are reported in Figure S3. Each indicator yields various c sig ; in Figure 5A, Var, AC(1) and H S maximize c sig over various noise levels, while other indicators like skewness and kurtosis perform poorly, as anticipated by the analytical results. CV and ID are also rather poor, likely because of fluctuations of mean values and anticipating n-tipping (cf. also Figures S2 and S3). For the case of multiplicative noise (Figure S3), H S keeps performing well while Var, as expected from the theoretical analysis, decreases its performance despite being still better than Skew and Kurt.
Complementing the analysis of the lead parameter requires understanding how many noise-induced tipping events occurred before it and assessing whether the increasing indicators alert for impending collapses or reflect transitions that have already happened. The analysis thus interprets warning indicators as ''anticipating'' or ''just-on-time detecting'' the tipping events. To do so, a counter C quantifies, for each parameter value c and for each noise level s, how many trajectories tip onto the alternative stable state. The results are in Figure 5B: as s increases, more n-tipping events occur before the bifurcation point. In  iScience Article particular for s > 0:42, several noise-induced transitions occur at cxc sig . Hence, as noise increases, the indicators capture ongoing critical transitions but are not able anymore to provide much earlier alerts. This likely explains the remarks from Dudney et al., 31 that EWS could not anticipate several transitions in realworld systems, in particular those characterized by high noise-to-signal ratios.
The previous results are also employed to define an optimization problem to maximize c I sig ðs; aÞ for varying a. To do so, we define a composite indicator as linear combination of indicators S = X k w k I k (Equation 18) and look for a set of weights w = fw k g that maximizes all c I sig ðs; aÞ as s increases (to guarantee robustness against noise levels), for the various a: where S are scores composed by sums of c S sig ðw; s l Þ over all s. In the set I , we include those indicators that are expected to be robust and performing, first and foremost in the additive noise case. Leveraging on the previous results, we therefore select Var, AC(1) and H S . As the problem is non-convex ( Figure 5C), we perform a grid search for all combinations of w k , with a stride 0.1 and such that P k w k = 1. See Fig. Figure 5C for the considered combinations to construct S. iScience Article Figure 5D reports the results of the optimization procedure. Combinations of Var and AC(1) make up for optimal indicators in case of additive noise, b w = ½0:9; 0:1; 0 for Var, AC(1) and H S , respectively, in case of a = 0. In this case, H S is log-proportional to Var (see Equation 30) and does not add much information. In turn, combining the indicators maximizes c S sig ðs; aÞ in case of mixed noise types. Finally, when multiplicative noise is prevalent in the system, using Shannon entropy is preferred ( b w = ½0; 0; 1 for a = 0). Note that, as the problem is non-convex, there may be more than one combination to create the optimal S. However, changes in weights w k are always within Dw k $ G10%w k and the trends are conserved (see dashed lines in Figure 5D). By following the computational procedure to perform a sensitivity analysis, we observe that such small Dw k yield changes of G4% on the scores S, on average over all a (DS˛½1:8; 6:5%), while offsetting w k by more than 50% (e.g., using full variance in case of multiplicative noise) worsens S (and consequently the optimal lead parameter) up to more than 20%.

Verification on experimental data
The theoretical predictions are verified and used to interpret experimental data from a previous publication. 26 The data are sampled from controlled experiments of budding yeast population collapse. Budding yeast cooperatively breaks down the sucrose necessary for its survival, thus inducing a density-dependent dynamics that realizes the Allee effect of bistable population dynamics (cf. Figure 2B). Repeated experiments empirically reproduced a saddle-node bifurcation by measuring population density (state variable) as a function of dilution factors (DF, control parameters) affecting the sucrose environment. Various EWS for population collapse can be estimated using distribution data. More details about data collection and analysis are in STAR Methods. Testing our theoretical results on a different system than Equation 15, yet still belonging to the saddle-node driven b-tipping class, would thus assess their generic applicability within this class. Figure 6 shows trends of each indicator individually, as function of the dilution factor (with critical value at 1600). The error bars are estimated from bootstrapping (STAR Methods). Figure 6 reproduces the results from Dai et al. 26 and includes the additional indicators considered in this paper. The mean is used to reconstruct the upper stable branch of a saddle-node bifurcation diagram (see Figure 1), reconstructed from data (the full diagram can be found in the original publication). However, it cannot be used as proper EWS as decreasing mean values could signify smooth changes rather than critical transitions, if the transition type and critical parameter are not known. Skewness and kurtosis fluctuate around 0 and 3, respectively, without providing EWS, as one expects in case of symmetric potentials (see Equation 36 and 0.18). AC(1) and the autocorrelation time (defined as À 1=log½ACð1Þ 26 ) first drop before increasing sharply just before the critical value. Comparing it with Figure 3, we speculate that there are commensurable time scales between the intake of sugar by yeast cells and their evolution in density. Further experiments are suggested to check for this intriguing hypothesis.
Even in this case, as expected, Var, Entropy (H S ), CV and ID display monotonous increasing trends close to the bifurcation point. The increases are thus assessed using the p value test (cf. STAR Methods) to check whether they are significant or associated with fluctuations. To trigger a significant early warning signal, we require a conservatory significant p value < 0:01. This way, we estimate the significant dilution factor DF sig for each indicator. For variance, DF sig = 1133, for the others DF. Comparing with the optimization results (from the previous section and Figure S3), we infer the presence of multiplicative noise in the system's dynamics. Note that entropy showcases the smallest p value at DF sig = 1000; it is also the most robust when changing the repetitions in the bootstrapping procedure (STAR Methods).
To test the hypothesis of association between EWS performance and noise type, we test combined indicators with H S and Var. According to the optimization above, the higher the variance content in the mixture, the lower the significance of the increasing trend. This is what is observed in Figure 7: having a balance between Var and H S yields DF sig = 1000, but with a higher p value than when reducing the ratio Var/ H S or when comparing with the case of entropy alone (from Figure 6).
Finally, we test combining CV and, since Figure 6 suggests that CV could perform well. Indeed, the new combined indicator yields DF sig = 750 (Figure 7, right), one dilution step before the others. This is not in contrast with the optimization analysis: CV is, in fact, expected to be as performing as H S if the noise levels are relatively high (see Figure S3). We recall that CV was not included in the optimization analysis to be ll OPEN ACCESS iScience 26, 107156, July 21, 2023 11 iScience Article generic and robust across noise types and levels. However, if high s in state-dependent noise is known, constructing a composite indicator using both CV and H S may improve the alerting performance.

DISCUSSION
The paper provides a systematic classification of tipping mechanisms, highlights their underlying modeling assumptions, and bridges mathematical insights and observations of real biological systems to classify various tipping mechanisms, toward quantitative understanding and prediction of such relevant phenomena. The work shifts the focus from studying specific systems, that may undergo some transitions, to studying transitions, along with their classes and properties, which can accommodate the behavior of different systems. An interesting question for future studies will be to develop data-driven methods to classify each system within its corresponding class, much like those developed to distinguish stochastic or chaotic signals. 92 This will dramatically help the understanding of biological processes and guide the selection of EWS or other methods to anticipate critical transitions, as well as informing methods to reconstruct cell developmental trajectories like those proposed by Eugenio et al. 93 Moreover, we systematically investigate EWS associated with noisy bifurcation-induced transitions, key dynamical routes for the regulation and control of many natural processes. So far, EWS have been mostly studied in highly controlled computational settings, or checked on empirical data with alternate success.
Our results make sense of previous observations, help to define their range of applicability to reliably predict systems' behaviors, and allow to understand why spurious signals may be triggered in certain cases. We also assess whether and when EWS can be interpreted as anticipating or just-on-time detecting critical iScience Article transitions in the presence of noise. By carefully analyzing noise types and parameter dynamics, we also extend previous results to more realistic settings, to guide real-world applications.
Using both analytical and computational methods, we observe that the variance -a highly employed indicator for EWS -may be sensitive to state-dependent noise, while AC(1) can be skewed by ramping control parameters. Both are good indicators in case of quasi-steady-state dynamics and Gaussian noise, with the ability to provide information about augmented risk of tipping events. In the other cases, Shannon entropy is the most robust and performing indicator and is suggested for applications in case of uncertain settings.
If precise information about noise type and intensity are available, constructing composite indicators can improve the early-alerting performance, e.g., by combining CV and H S .
After preliminary suggestions by Drake et al., 20 composite indicators have been subsequently tested, mostly based on system-specific traits. 94 Contrarily, this paper introduces, demonstrates and optimizes a novel indicator based on summary statistics from normal forms, in principle extendable to numerous systems sharing similar dynamical features and noise properties. The optimization of composite indicators further points to the use of machine learning methods when abundant data are available, 95 but also opens important caveats for their application in real life: feature combinations may be optimized for certain settings (e.g., noise intensity or type) but may be hardly generalizable for others. Our results remark that training should be performed considering all possible combinations, or by first assessing which critical transition class is being considering. Otherwise, misleading signals may be triggered and wrong conclusions reached. On the other hand, our results can be used for feature selection of more interpretable machine learning algorithms that leverage the proposed composite indicators, insofar defined for a-priori assessment of systems that lack big data.
This work provides results and guidelines for the application of EWS from the critical transitions framework, but some points should still be covered by future studies. They include more refined analytical derivations of indicators in case of inhomogeneous processes as well as closed formulae for entropy in exotic settings. Further investigations on realistic systems, including non-autonomous transitions currently understudied in systems biology, are thus suggested as extensions of our work. Similarly to the case of univariate indicators, which may fail to distinguish critical and non-critical transitions, 83,84 multivariate indicators may be limited in capturing such differences. Other methods (e.g., Angeli et al. 14 ) may be used to complement the inference, but further studies are suggested toward the inference of transition types using EWS. Another limitation of the present study is the restriction to low dimensional systems. In principle, they are representative of any system after dimension reduction techniques are applied (see e.g. Gao et al. 44 and Laurence et al., 96 or Heino et al. 97 for Principal Component Analysis-based techniques), but it is necessary to assess if and how the latter induce performance drops. Extending the analysis to high dimensional systems, e.g., by testing multivariate indicators 98 or further refining EWS performance when multiple independent variables can be observed, is thus suggested to future studies. Moreover, our theoretical results have been verified on empirical data from literature, but we acknowledge the need of performing additional experiments to continuously validate our predictions. In particular, we suggest to design new experiments to test the quantitative predictions about lead parameters and to assess what happens in case of rapidly ramping parameters. Finally, although our results are in principle extendable to other fields, additional testing and validation is necessary if the data sampling differs and the ergodic equivalence between distribution iScience Article data and time series may be challenged, as it is the case, e.g., for ecology or climate. Future targeted studies, that are beyond the scope of the current article, will likely cover this point, e.g. on data from Dakos et al. 99 or Chen et al. 100 Our results can be readily tested and applied on real-world monitoring systems and can inform the development of new indicators to address specific problems like cancer onset, much like previous studies 22 did using less performing measurements. Our results can be further applied to a range of other diseases to detect their onset, from diabetes 101 and other complex diseases 22 to epidemics. 28 In addition, leveraging the sensitivity of indicators' trends to noise type and parameter dynamics can provide new methods to infer the latter from empirical data. For instance, comparing Figure 6 with Figure 3 supports hypothesis of commensurable time scales between intake of sucrose (affected by the dilution factor) and cells' growth in yeast experiments 26 ; such hypothesis, to be confirmed using controlled experiments and further computational studies, could advance our knowledge beyond the current slow-fast approximations. 40 Similarly, the prevalence of certain noise types can be inferred by comparing data and theory. Overall, we connect theory and data, such that knowledge about the dynamical settings allows optimizing EWS, and analysis of statistical indicators enables inference of dynamical properties.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Materials availability
This study did not generate new unique reagents.
Data and code availability d This paper analyzes existing, publicly available data. These accession numbers for the datasets are listed in the key resources table.
d All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table.
d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request

Topological equivalence and normal forms
Bifurcations model drastic changes in the qualitative behaviour of dynamical systems, such as shifts in equilibria and regimes. 42,43 Before delving into bifurcations and their representation as normal forms, recall the concept of topological equivalence.
Local topological equivalence between two dynamical systems fT ; R n ; 4 t g and fT ; R; j t g is established if there exist a homeomorphism h : R n /R n that maps orbits of the first system to orbits of the second one, and the direction of time is preserved. Local topologically equivalence near an equilibrium b u is, in turn, established between a dynamical system fT ; R n ; 4 t g and a dynamical system fT ; R; j t g near an equilibrium b y if there exist a homeomorphism h : R n /R n that is defined in a small neighborhood U˛R n of b u, satisfies b y = hðb uÞ, and maps orbits of the fT ; R n ; 4 t g˛U onto orbits of fT ; R; j t g˛V = hðUÞ3R n while preserving the direction of time.
A bifurcation consists in the appearance of a topologically non-equivalent phase portrait under variation of parameters. The difference between the dimension of the parameter space and the dimension of the corresponding bifurcation boundary is called ''codimension''.
To determine a system's behaviour near bifurcations, minimal-order forms, called ''normal forms'', can be employed. In fact, the normal form of the bifurcation is locally topologically equivalent near an equilibrium to all systems exhibiting that certain type of bifurcation. 102  The sign ''À'' in ''Àk'' is included so that Equation 28 is interpreted as the associated Langevin equation to a Ornstein-Uhlenbeck process. 79 The term multiplying the deterministic drift can thus be interpreted as À vV =vx where V ðxÞ is the potential governing the drift of the particle subjected to random noise. In our case, thanks to the choices made,

REAGENT or RESOURCE SOURCE IDENTIFIER
that is, a quadratically shaped adjoining potential typical of an overdamped oscillator under noise, of which k represents the depth. The working hypothesis is that boundary of the ideal potential V can grasp the boundary of the attracting basin of the original model after sufficiently long time. Equation 28 is analytically tractable to understand the main qualitative features of more complicated critical transitions. However, it requires ad hoc extensions when studying system-specific quantitative details like observability boundaries and lead times. Gardiner 79 also extends Equation 28 to inhomogeneous processes with ramping parameters, used in Equation 13. Panel (c) shows the bifurcation diagram, over an unfolded supercritical pitchfork bifurcation, of _ x = q + pðx À 1Þ À ðx À 1Þ 3 , which corresponds to the bifurcation normal form, shifted (to better visualize the diagram) and modified by a small perturbing term q = 0:01 unfolding the bifurcation 50 into a smooth branch. In brief, an unfolding of a dynamical system under static equivalence is one that exhibits all possible bifurcations of the equilibrium (rest) points, up to topological equivalence of the set of equilibria. 42 In other terms, it investigates what happens when small terms are added to the original bifurcation, mimicking extra parameters, small offsets or ''impurities''.

Reproduce Figure 1
The illustrative attractors in panel (a) and (b) are two-well potentials associated, e.g., to the cusp bifurcation (aka ''organising center'' 50,93 ), a generic bifurcation described by _ x = a + bx À x 3 , where the combination of a and b determine bistability and the route to a saddle-node bifurcation.

Entropy in case of Gaussian noise
Within a symmetric potential, elicited by a (locally) quadratic normal form, consider a Gaussian distributed variable y $ N ðm; VarÞ. Its entropy is:

Measurement noise
Consider a measurement process with uncertainties s 2 m , independent from system variance (Equation 9). The resulting expected error, obtained from summing the two standard deviation in quadrature, 105  . In principle, we can explicitly consider multiplicative noise like in the main text. However, the goal in this case is to compute if notable discrepancies exist between ideal measurements (no uncertainty) and realistic measurements (with some uncertainty, that can be filtered to correspond to additive noise). Hence, only the case of additive white process noise is currently considered. This results in: , we can immediately see that measurement uncertainties s m induce small scaling but do not alter the functional. Only relatively high measurement uncertainty levels change the absolute values of expected lag-1 autocorrelation, but maintain the increasing patterns close to critical points.

Skewness and kurtosis
For certain simulated systems, the third statistical moment (skewness) has been suggested to provide useful early warnings. 82 However, experimental results 26 were not able to confirm the expectations, estimating flat and fluctuating trends before a tipping point.
For a stochastic process with quasi-steady state parameter and small noise limit, its statistical moments are Consequently, the normal forms considered above, under small noise or in case of symmetric potential, are expected to display a flat skewness.
On the other hand, the integral (.16) may be non-zero, and even dependend on the drift parameter k, if ms 0 or if PðyÞ is asymmetrical. In the first case, solving Equation 35 yields (provided that Re½k > 0): ) In this case, as k/0, the skewness is expected to increase, potentially providing an early warning On the other hand, an asymmetric potential can be obtained in case of multiplicative noise. 15,79 Depending on the specific form, it may be possible to observe increasing trends associated to EWS, but they may be system-specific and not generalisable. In this sense, there is no ambiguity between the results of Guttal et al. 82 and Dai et al. 26 : they were studying systems with different properties, using an indicator that is not particularly performing and generalisable.