Abstract Critical transitions – abrupt, often irreversible changes in system dynamics – arise across human and natural systems, often with catastrophic consequences. Realworld observations of such shifts remain scarce, preventing the development of reliable early warning systems. Conventional statistical and spectral indicators, such as increasing variance, tend to fail under realistic conditions of limited data and correlated noise, whereas existing deep learning classifiers do not extrapolate beyond their training data distribution. In this work, we introduce TipPFN, an in-context learning (ICL) framework that uses a prior-data fitted network to infer a system’s proximity to a critical transition. Trained on our novel synthetic data generator, which is based on canonical bifurcation scenarios coupled to diverse, randomized stochastic dynamics, TipPFN flexibly capitalizes on contexts of various sizes, complexity and dimensionalities. We demonstrate robust, state-of-the-art early detection of critical transitions in previously unseen tipping regimes, sim-toreal examples, and real-world observations in both ICL and zero-shot settings.

∗ Equal contribution

Preprint.

Introduction

Tipping points occur when variation in forcing or system parameters triggers an abrupt, potentially irreversible transition to a qualitatively different dynamical regime [1, 2]. Such transitions are central to climate [3], ecological [4], infrastructural [5], and collective biological and social systems [6, 7, 8], where local shifts can cascade through interacting components with farreaching consequences. Classical theory has largely focused on bifurcation (b-tipping) classes of tipping mechanism. Here, the loss of stability of an attracting state produces critical slowing down (CSD) [9], where recovery from perturbations becomes progressively weaker [10, 11]. This has motivated widely used early warning signals (EWS), including increasing lag-1 autocorrelation (AR1) [12], variance [13], and skewness [14]. However, these indicators rely on restrictive assumptions, including near-equilibrium linearization, stationarity, sufficiently long records, and simple bifurcation structure [15, 16]. They may therefore fail for rate-induced (r-tipping) and noise-induced (n-tipping), where transitions can occur without pronounced CSD [7, 17, 18].

1.0 0.9 0.8

AUROC

0.7 0.6 0.5 0.4 80

tcrit -10

Lead time ¢ TipPFN TabPFN 2.6

EWS ¿var Bury

Huang Zhuge

Figure 1: TipPFN detects critical transitions earlier and more accurately than classical early warning signals (EWS), state-ofthe-art deep learning and ICL-based methods across 14 semi-real, sim-to-real, and real-world systems spanning climate, engineering, and biology. Lines show acrosssystem mean and standard deviation.

From a machine learning (ML) perspective, characterizing critical transitions can be viewed as a form of out-of-distribution (OOD) task, where the system leaves one regime and enters into a dynamically distinct region [19, 20, 21]. This makes tipping prediction fundamentally different from standard regression tasks, since models must infer the onset and type of critical transitions from sparse, noisy data and partially observed trajectories, often without examples of the event itself. Recent ML approaches have shown promise by learning features of canonical tipping systems directly from data [22, 23, 7, 24, 25]. However, many of these methods are trained on particular classes, limiting their ability to generalize across observations and different critical mechanisms. A promising alternative has recently emerged in the form of prior-data fitted networks (PFNs), which are trained on large ensembles of synthetic datasets and learn to perform inference via incontext learning (ICL) [26]. Rather than fitting a model to a single dataset, PFNs are trained across millions of simulated tasks, enabling them to infer structure and make predictions directly from limited observations. Thus, PFNs have already transformed tabular ML through models such as TabPFN [27, 28, 29], which, despite being trained exclusively on synthetic data, achieve state-ofthe-art performance on real-world benchmarks. More broadly, PFNs have demonstrated remarkable generalisation across diverse tasks, including time series forecasting [30], causal inference [31, 32], and robotic controls [33]. Our key contributions are (see Figure 2): • Model: TipPFN, a transformer-based architecture for joint prediction of proximity to a critical transition from short, noisy time series. • Training and Benchmark: TipBox, our benchmarking suite and training data generator of stochastic dynamical systems across diverse tipping regimes and complexities. • Generalization: Robust performance across unseen tipping regimes, sim-to-real cases, and real-world systems in both ICL and zero-shot settings (see Fig. 1).

Background

Critical transition in dynamical systems Critical transition, specifically tipping points in dynamical systems, can arise through multiple distinct mechanisms, broadly classified as bifurcation-induced (b-tipping), noise-induced (n-tipping), and rate-induced tipping (r-tipping) [2]. Consider a general nonautonomous stochastic differential equation (SDE) of the form: 2

Context

tcrit

TipBox

Query Training Loop

RDTC

Training

tcrit

Dynamics

Time

? Time

Classes of tipping driver N-tipping

Unobserved data-generating mechanism

R-tipping

tcrit

Time

Observations

Inference

B-tipping

TipPFN

RDTC

Driver

LRDTC, Dynamics

Tipping driver nodes Dynamical state nodes

Time

Initial regime → transition regime

Random observed nodes

Figure 2: Overview. TipPFN is trained on synthetic data, based on embedding canonical b-tipping systems into high-dimensional, randomized stochastic dynamics. Primary training target is the relative distance to criticality (RDTC), a measure of how close a system is to a critical transition. Although the synthetic dynamics were only driven by b-tipping systems, TipPFN successfully generalizes to other classes of tipping systems, and a manifold of real-world systems.

dx = f (x, λ(t)) + σ ξ(t), (1) dt where x ∈ Rd denotes the system state, λ(t) is a time-dependent forcing parameter, σ > 0 controls the noise amplitude, and ξ(t) represents Gaussian white noise. Bifurcation-induced tipping In b-tipping, a transition occurs when a quasi-static equilibrium x∗ (λ) loses stability at a critical parameter value λcrit . Linearizing the dynamics around x∗ (λ) yields the Jacobian J(λ) = ∂f /∂x x∗ (λ) , whose eigenvalues {µi (λ)} govern local stability. A bifurcation occurs when the real part of the leading eigenvalue crosses zero, i.e., max Re(µi (λcrit )) = 0. i

(2)

Rate-induced tipping In r-tipping, a transition occurs when the rate of change of the forcing, given by λ̇, is sufficiently large such that the system fails to track the moving equilibrium x∗ (λ(t)), leading to a transition without any local loss of stability, i.e., the linearized system satisfy maxi Re(µi (λ)) < 0. Thus, these tipping mechanisms cannot, in general, be inferred solely from equilibrium stability analysis commonly used to check b-tipping (see Appendix A for more details). Noise-induced tipping In n-tipping, a transition occurs primarily through stochastic perturbations σ ξ(t) that induce escape from the basin of attraction of a stable equilibrium. This process can occur in conjunction with b-tipping and r-tipping in SDEs, a setup which form the basis for this work. Prior-data fitted networks PFNs are transformer-based neural networks trained on synthetically generated data drawn from a predefined prior, enabling approximate Bayesian inference via ICL [26]. In this framework, inference is amortized at the dataset level. The model learns a mapping from a dataset D to a posterior predictive distribution p(y|x, D) by training on large collections of simulated tasks (Di , xi , yi ). Rather than fitting a model to each new dataset, PFNs directly infer predictions conditioned on the observed data, performing Bayesian inference in a single forward pass. Related work In addition to the classical EWS based on statistics (AR1, variance, skewness), there is an emergence of dynamics-based approaches extracting additional spectral structure of the systems, including tracking of dominant eigenvalues [34, 35]. Recently, fully data-driven ML approaches, including reservoir computing [25, 23], convolutional, and recurrent-based neural networks [22, 24], have been proposed to detect tipping points directly from time series data. These models can outperform both classical and spectral EWS when trained on simulated data, but are tailored to specific tipping regimes (e.g., b-tipping [22, 24] or r-tipping [18]). They also show limited generalization capabilities beyond the training distribution, which we further showcase throughout. 3

TipPFN: Predicting tipping behaviors with PFNs

TipPFN is designed to infer approaching critical transitions directly from observed trajectories. To this end, we train a model on synthetic dynamical systems in which a control parameter λ is gradually varied to produce trajectories that progress from stable regimes to towards a critical transition. By embedding controlled, stereotypical tipping systems within larger, randomized systems of nonlinear SDEs, the synthetic data provide a scalable source of diverse transition patterns and allow the model to learn transferable structure rather than system-specific rules. This is essential for our setting, where the goal is generalization to previously unseen systems. Next, we define the prediction target, the relative distance to criticality. We then describe synthetic data generation, and finally the TipPFN training and inference setup; also see Fig. 3 and Appendix B. Relative distance to criticality (RDTC) Direct supervision on critical events is often ill-posed, since the observed transition time can depend on noise, finite-time effects, and threshold choices. For the bifurcation-induced tipping systems used in training, however, proximity to criticality is well defined through the forcing parameter λ. We therefore introduce the relative distance to criticality (RDTC) Λ as a continuous supervision target. Let λ̃ ∈ [0, 1] denote the normalized control parameter, where λ̃ = 0 corresponds to a stable regime and λ̃crit = 1 to the forcing value at which the underlying deterministic system undergoes a bifurcation. We define Λ = 1 − λ̃. (3) For Λ = 1, the system is maximally far from the bifurcation within the normalized forcing range; for 0 < Λ < 1, it remains in the sub-critical regime, reaching criticality at Λcrit = 0. RDTC serves as the primary training target for synthetic bifurcation systems. Furthermore, it is intended to act as a transferable measure of proximity to criticality across otherwise distinct systems. Prior structure TipPFN’s prior p(ψ) over generative processes ψ factorises over a canonical tipping system M , a random potentially cyclic interaction graph G and an auxiliary nonlinear SDE system Z. G specifies the interaction structure of Z, and how it is driven by the state variables zM (t) = (xM (t), yM (t)). This driving timeseries is generated by sampling from a flexible family of low-dimensional nonlinear dynamical systems designed to exhibit canonical tipping behavior [22] (see Appendix B.1 for details). We define the prior over the generative process ψ = (M, G, Z) as: p(ψ) = p(M ) p(G) p(Z | M, G),

(4)

where M ∼ Uniform({fold, Hopf, transcritical}) determines the bifurcation class inducing qualitatively different signatures in the pre-tipping dynamics. G is generated by drawing from a distribution

tcrit

context observed window forecast

¤¤

¤¤ features signal

moving windows

time

nc:1-med=0.69

truth nowcast forecast

signal

W = 128

fc:1-med=1.05

¿var =0.293 + 15 features

100

200

300

400

500

sequence index

Figure 3: Tipping Prediction with TipPFN. (a) Example multi-variate time series from a critical episode and underlying RDTC Λ. Green bars mark candidate observation windows ending ∆ time steps before the critical time tcrit . (b) Example query episode at ∆ = 30 and context composed from a critical (red) and non-critical (blue) episode. TipPFN and TabPFN are conditioned on the shaded context and observed window, including the signal and all available features, and then predict relative distance to criticality, Λ∗ , for the nowcast and forecast region. These predictions form the basis of the tipping-risk scores used for early-warning assessment. Baseline scores (e.g. EWS τvar ) are computed only from the observed signal window. 4

over directed, potentially cyclic graphs with power-law distributions over incoming and outgoing node degrees [36] to yield the interaction structure within the auxiliary dynamic nodes, followed by random attachment of two driver nodes representing M ’s state variables. Based on this dynamical causal structure, nonlinear SDE interaction terms are drawn from a distribution over multi-layer perceptrons. Additional linear stabilization terms and a global timescale parameter are added to the auxiliary nodes and also randomized. The hyperparameter distribution of the auxiliary system is chosen to keep it in a stable regime, preventing additional, spontaneous bifurcations due to intrinsic auxiliary system dynamics, which would confound the signal provided by the driving system. Prior-data generation As illustrated in Fig. S1a, to fit the prior, we specify a sampling scheme over context-query pairs D of the form p(D) = Eψ∼p(ψ),S∼p(S|ψ) p(D | S) , (5) which first samples a generative process ψ = (M, G, Z) ∼ p(ψ) and then a synthetic episode ensemble S = {Ek }K k=1 ∼ p(S | ψ) which fixes the canonical tipping system M , the interaction k graph G and SDE system Z. Each episode Ek = (z1:T , uk1:T , λ̃k1:T ) is based on an independent k noise realization and contains the trajectories of the driving system z1:T and the auxiliary nodes uk1:T k created by the independent linear forcing schedule λ̃1:T . Forcing schedules are constrained to produce trajectories spanning multiple dynamical regimes: critical (tipping), approaching towards or receding from criticality, constant, and equilibrium (non-critical); see Appendix B.7. From each sample query S = {Ek }K is held out as the query and a context set C ⊂ E \ E query including k=1 , one episode E |C| ∈ {0, .., 3} context trajectories is sampled uniformly from the remaining K − 1 episodes. A subset of feature dimensions of the driving system z and the auxiliary system u is randomly selected to yield the observed features o with dim(o) ∈ {1 : 16}, rendering some the resulting prediction problems only partially observed. By resampling the time series with different resolutions, we realize variable-context training tasks in which query and context episodes share the same underlying system class and coupling structure, but differ in forcing and stochastic realization. This yields a query-context pair D = ((Λ1:T , o1:T ) ∪ C). As a result, the model learns to use context when informative without depending on its presence. In total, we trained on ∼1M context-query pairs D from ∼120k ψ samples with K = 6 episodes and M ∈ R2 , Z ∈ R14 . Input representation Each training task is encoded as a tabular multi-episode sequence with episode and time identifiers, observed variables, target variables, and prediction masks. To avoid target leakage, variables are normalized using only observed values. Observed and target variables are centered and transformed with asinh scaling, which reduces the effect of large amplitude differences across heterogeneous signals. RDTC is transformed separately as Λ∗ = tanh (5Λ), the nonlinearity providing finer resolution near criticality. In the query episode, RDTC is always fully masked, whereas the remaining target variables are only partially masked. Architecture TipPFN is a transformer-based architecture with a similar structure to TabPFN [27, 28]; the main adaptations for tipping prediction lie in the task construction, masking scheme, and training objective described above. Transformer attention is structured, such that all tokens from the context set C can attend to all other tokens of C. Tokens from the query trajectory can attend to all tokens from C, themselves, and all tokens from the query trajectory with smaller time features, i.e., tokens representing preceding time steps (causal attention), as illustrated in Fig. S1b. This ensures that no future information leaks into the estimate, as required for early warning signals. More details on the architecture of PFNs is provided in Appendix B. Training TipPFN is trained as a TabPFN-style transformer on tabularized multi-episode trajectories. The model predicts quantile distributions for masked targets and is optimized with a weighted pinball loss. Along the temporal dimension, attention is causal, so each prediction can depend only on past and present observations. Context episodes are fully visible. The parameters θ of a transformer model qθ are optimised using the following training loss: " T # X L = ED=((Λ1:T ,o1:T )∪C)∼p(D) − log qθ Λt | o1:min(t,tnc ) , C .

(6)

t=1

Here, Λt denotes RDTC at time step t in the query episode, and o1:min(t,tnc ) denotes query observations up to step t or a cut-off tnc < t, up to which observational data is available, to train nowcasting 5

and forecasting capabilities. The RDTC prediction loss is augmented with a similar loss, scaled by a factor of 0.2 relative to the RDTC loss, which requires to predict a randomly selected subset of the observational features from the remaining set of at least 1 feature. This loss was added to encourage the model to learn a more general, task-agnostic representation of the nonlinear stochastic dynamics of the synthetic training data. Posterior Predictive Distribution for RDTC By minimising training loss given in Equation 6, the transformer model qθopt approximates the true Bayesian posterior predictive distribution (PPD) [26]: qθopt Λt | o1:min(t,tnc ) , C ≈ p(Λt | o1:min(t,tnc ) , C) Z = p(Λt | o1:min(t,tnc ) , C, ψ) p(ψ | o1:min(t,tnc ) , C) dψ. Ψ

Thus, TipPFN approximates this intractable integral in a single transformer forward pass by amortising inference over a large ensemble of synthetic prediction problems sampled from p(ψ) during training, yielding a predictive distribution approximating Bayes-optimal uncertainty estimates. When no context episodes are available, the prediction relies entirely on the query trajectory and the prior p(ψ). Increasing the available context sharpens the posterior p(ψ | o1:t , C) over generative processes compatible with the observed data, depending on the information content of the added context samples. Inference TipPFN supports two inference modes depending on the operational setting: For nowcasting, it estimates the current RDTC at time t, conditioned on all observations up to that point: qθopt (Λt | o1:t , C) ,

(7)

where Λt denotes the RDTC at the current time step t. This mode is appropriate when the primary goal is to monitor the system’s instantaneous distance to a critical point in real time, issuing an updated estimate as each new observation arrives. In the forecasting setting, TipPFN predicts the future trajectory of RDTC, conditioned on observations up to the current time step: qθopt (Λt:T | o1:t , C) , (8) where Λt:T denotes the sequence of RDTC values from the current time step t to a future horizon T . This mode is appropriate when the goal is to anticipate if the system will approach a tipping point, enabling proactive intervention before the critical transition occurs.

Experiments

In this section, we evaluate the performance and predictive skill of TipPFN for tipping behavior across different mechanisms. We benchmark against a comprehensive set of competitive baselines on datasets ranging from synthetic systems to real-world observations. Evaluation is based on classification performance (critical vs. non critical) using receiver operating characteristic (ROC) curves [37, 38] and Area Under the ROC curve (AUROC) scores [39, 40], as well as lead-time analysis; more details in Appendix C. Validation datasets We generated and collected cross-domain validation datasets covering 12 model families (canonical and semi-real) and 9 observational systems (sim-to-real & real-world); see Appendix C.1 for the full list and descriptions. Canonical datasets are held-out realizations from the synthetic prior. Semi-real datasets are reduced-order models of real-world dynamics that are more complex than the canonical systems and include unseen tipping behavior, e.g., r-ipping. Sim-to-real datasets contain both simulated and observed trajectories, allowing us to test whether matched simulations can provide context for empirical transition prediction. Real-world observation datasets consist of empirical time series from systems without fully specified generative dynamics, including neurological seizure recordings, ocean circulation (AMOC), power-grid blackout, and cyanobacteria population dynamics; depending on the available data, they support either context-based transfer or strict zero-shot RDTC prediction. 6

Baselines We compare TipPFN against classical EWS, including AR1, variance, and skewness in rolling windows, summarize their temporal trends with Kendall-τ [41], and obtain ROC curves by thresholding these trend scores, following standard practice [9]. ML baselines include Bury [22], Huang [18], and Zhuge [24], operating on univariate time series. TabPFN2.6 [28] provides a state-ofthe-art ICL baseline; more details are provided in Appendix C.2. Querying procedure & scoring All methods are evaluated on matched moving query windows ending ∆ time steps before the critical event. Uni-variate baselines were evaluated on the driving time series, zM (t). PFN-based models additionally receive context episodes, whereas classical and ML baselines use only the query window; see Fig. 3. The original window’s time is re-indexed to start at 0 such that no positional information is leaked. At test time, no method has access to future observations, hidden system parameters, or ground-truth RDTC, except for Zhuge, which also observes the forcing parameter and therefore receives privileged query-episode information. As the ML baselines require a window size W ≥ 500 steps while TipPFN uses W ≤ 128, we construct length-adapted baseline inputs by resampling, back-filling, or forward-filling the query window. For each model, we pick the best-performing variant when reporting results (Appendix C.5.1). TipPFN and TabPFN receive the same context episodes and query window. For each query trajectory, they predict the (transformed) RDTC Λ∗ over the observed window and future time points; see Fig. 3. TabPFN is cast as tabular regression by treating trajectory time points as rows and context episodes as in-context training examples. Predicted RDTC distributions are then converted into tipping-risk scores, e.g. 1−median(Λ) or CDF(Λ < 0.05) at the final nowcast or forecast time point, which rank queries by predicted risk. For TabPFN and TipPFN, we report the final forecast-time 1 − median(Λ) score for all methods. As illustrated in Fig. 3b, this score is the parameter-free estimate of forecast remaining RDTC at the forecast horizon, and therefore aligns directly with the early-warning task. Note that the forecast is always 192 − W steps into the future, so shorter windows make a longerhorizon forecast, thus also increasing uncertainty. We report score sensitivity analyses in Table S14 and Appendix C.5.1. 4.1

Generalization across synthetic transition systems

We start with the canonical bifurcation setting, evaluating on held-out realizations from the fold, Hopf, and transcritical families used to define the training prior. TipPFN already performs strongly

Table 1: Average pre-tip AUROC by system and dataset. Balanced AUROC for detecting critical episodes, averaged over positive lead times (∆ > 0). Each column uses one fixed evaluation setting across all listed rows that yields highest AUROC across all these systems (see Table S13 for details). Dataset

TipPFN 0c

TipPFN 1c

TipPFN 2c

TabPFN 1c

TabPFN 2c

EWS τvar

Bury [22] Ptip

Huang [18] Ptip

Zhuge [24] ∆± tip

Canonical B-Fold B-Hopf B-Trans.

.876 .861 .864

.929 .935 .925

.955 .952 .950

.574 .640 .581

.786 .820 .807

.633 .635 .543

.697 .842 .722

.462 .496 .453

.522 .665 .458

Semi-real B-Harv. B-RM TC B-RM Hopf B-SEIRx B-AMOC R-Bautin R-SN R-Compost R-AMOC

.846 .815 .769 .582 .853 .308 .644 .767 .257

.958 .951 .829 .822 .898 .543 .802 .623 .782

.974 .963 .904 .905 .940 .700 .892 .665 .913

.685 .494 .499 .495 .616 .404 .335 .472 .521

.810 .741 .689 .685 .825 .489 .569 .493 .654

.914 .604 .557 .498 .925 .451 .440 .481 .270

.967 .784 .638 .468 .948 .308 .354 .727 .102

.873 .604 .484 .507 .509 .538 .273 .678 .712

.734 .570 .422 .381 .679 .180 .080 .506 .112

.361 .512 .663

.816 .519 .833

.486 .530 .502

.560 .474 .641

.515 .565 .540

.674 .476 .483

.385 .522 .404

– – –

.769

.859

.522

.670

.572

.613

.527

.443

Real-world & Sim-to-Real SWEC-iEEG .674 TAC .521 DaphniaExt .527 All datasets

.678

ctx=0

a 0.55

0.60

0.74

ctx=1 0.77

0.75

0.73

0.89

ctx=2 0.90

0.81

0.80

0.94

b 0.95

real context

synthetic context 20

0.60

0.75

0.78

0.66

0.80

0.83

0.68

0.70

0.87

0.90

0.54

0.56

0.66

0.68

0.51

0.65

0.68

0.55

0.56

0.77

0.82

128

0.51

0.54

0.58

0.56

0.75

0.80

0.51

0.55

0.58

0.56

0.76

0.81

0.54

0.58

0.60

0.57

0.78

0.83

0.5

ctx

TPR

AUROC

TabPFN

window size

TipPFN

128

0.5

FPR

0.5

FPR

features

Figure 4: Context matters. (a) Multi-parameter sweep over number of context episodes, observed feature channels, and window size W for TipPFN and TabPFN. Color denotes AUROC averaged over all ∆ > 0 and all datasets in Table 1; stars mark the best configurations. Note that TabPFN cannot be used in the zero-context setting. (b) TipPFN ROC of Daphnia dataset for ∆ = 2, W = 16 and varying number of real or simulated context episodes.

without additional context trajectories, with robust AUROC across all three systems and clear gains over classical early-warning baselines (Table 1, Fig. S14). Importantly, the binary task is stricter than a simple critical-versus-noncritical split: negative examples include not only equilibrium trajectories, but also trajectories that remain subcritical while approaching or receding from criticality. Strong performance therefore indicates that TipPFN learns reusable dynamical signatures rather than memorizing individual trajectories. Beyond canonical bifurcations We next test generalization. On semi-real reduced-order systems, TipPFN typically matches or exceeds the evaluated baselines, showing transfer from simple synthetic bifurcations to more realistic unseen dynamics (Table 1; Fig. S14). A stronger test is r-tipping, which lies outside the training mechanisms altogether. Despite this mismatch, TipPFN remains competitive on Bautin, compost-bomb, saddle-node, and AMOC r-tipping systems. These results use one fixed configuration per method, selected for strong aggregate performance across all systems. Under a more permissive best-per-system comparison, several baselines improve on individual tasks, but their gains are often system-specific (Table S10). TipPFN remains bestperforming on many datasets, and where it is outperformed, it typically is close to the best baseline. 4.2

Leveraging context

A central motivation for TipPFN is that context provides system-specific information at inference time. This is particularly valuable for tipping prediction, where observations are often short, noisy, and only weakly informative on their own. Rather than reducing trajectories to one-dimensional early-warning statistics, TipPFN conditions on full multivariate context episodes and can therefore infer which feature patterns distinguish tipping from non-tipping behavior. Figure 4a shows that TipPFN benefits most from combining additional context, higher-dimensional observations, and sufficiently long windows. For both TipPFN and TabPFN, performance improves when context episodes are added even under a fixed input budget, suggesting that the gain comes not only from observing the query longer, but from comparing it against a broader range of system behavior. Increasing the number of observed features yields a further boost, highlighting the value of multivariate inputs. 4.3

Real-world transfer

We next evaluate TipPFN on sim-to-real and real-world datasets, where true RDTC labels are typically unavailable. We therefore construct surrogate RDTC targets, Λ = Λ(t), from experimentally observed or independently estimated transition markers, and test whether TipPFN can generalize using limited empirical context or, where available, matched simulations. 8

a AMOC

b Blackout

c Microcosm

0.50

signal

150 20

0.25

0.00

100

¤¤

1.0

0.5

2018-Feb

count

0.0 40

1.0

0.5

0.0 40

0/100 crossed

0.0 40

100/100 crossed median lead time ¢ = 145 s

2008

2012

2016

2020

2024

date

400

200

100/100 crossed median lead time ¢ = 6.79 days

Lead time ¢ [s] observed

bifurcation

¤ = 0:05 crossing

nowcast q50 §80% CI

Lead time ¢ [days] median crossing

first crossing

Figure 5: Zero-shot TipPFN RDTC nowcasts on three uni-variate real-world time series, without per-system fine-tuning. For each system, the top row shows the observed empirical signal, the middle row shows an ensemble of TipPFN predictions for Λ∗ (median and 10–90% predictive bands over 100 stochastic retained time points), and the bottom row shows the distribution of the first predicted crossings Λ∗thrs = tanh(5Λthrs ) ≈ 0.245 across the ensemble. If available, the red dashed line marks bifurcation time from spectral-based dominant eigenvalue analysis [35]. On the SWEC-iEEG dataset, we test TipPFN performance on neurological seizure prediction from multichannel iEEG bandpower trajectories. Since seizure dynamics vary substantially across patients, the benchmark probes whether context episodes enable patient-specific in-context adaptation. TipPFN benefits from this structure and performs strongest overall and in per-patient comparison (see Fig. S16 and Appendix C.4.1). On DaphniaExt, we analyse TipPFN on extremely sparse (W = 16) real observations and matched simulations [42]. As shown in Fig. 4b, simulated and observed context reach comparable performance but respond differently to context size. Adding observed episodes keeps improving performance and simulated context provides even better results, compensating for sparse observations. Furthermore, while TabPFN and TipPFN both perform well on predicting the extinction event, only TipPFN discerns the underlying bifurcation [43]; see Figs. S9, S21 and Appendix C.4.2. On TAC, we evaluate TipPFN on a noisy experimental system with a stochastic subcritical Hopf transition [44], using the known critical point to define surrogate RDTC targets. While the Hopfspecific Bury baseline performs strongly, TipPFN improves substantially with real and simulated context, especially when the supplied episodes are themselves critical (see Fig. S22, Appendix C.4.3). Zero-shot prediction on real-world systems We also evaluate TipPFN in a strict zero-shot setting, without context trajectories. The relevant output is the predicted RDTC trajectory, a time-resolved nowcast of distance to criticality from a single observed time series. This allows critical-point timing to be estimated even when no matched empirical or simulated context is available. Across six real-world observational time series (see Fig. 5, Appendix S23), TipPFN’s RDTC predictions are informative: In the AMOC example, we detect an increase in RDTC variance coinciding with an abnormal reversal in AMOC slow-down [45]. In a power-grid time series preceding a blackout, TipPFN signals hidden instability 2 minutes before the failure becomes apparent. In the microcosm series, it estimates a cyanobacteria population collapse approximately a week from the actual event. Because visible regime shifts may lag behind hidden loss of stability, zero-shot RDTC prediction can provide advance warning and create a window for mitigation or response.

Conclusion

We introduce TipPFN, a transformer-based prior-data fitted network designed to predict the relative distance to criticality (RDTC) directly from observed trajectories. TipPFN is trained exclusively on a novel synthetic data-generation framework combining canonical tipping systems with randomized high-dimensional nonlinear stochastic dynamics, enabling universal prediction across diverse transition regimes. We evaluate TipPFN against state-of-the-art machine learning and dynamical systems baselines on more than 20 datasets spanning synthetic, sim-to-real, and real-world settings across 9

multiple domains, noise levels, and tipping mechanisms. TipPFN consistently outperforms existing approaches while demonstrating strong generalization across previously unseen systems. Unlike existing tipping-predicting methods, TipPFN can leverage additional context trajectories, auxiliary features, and variable observation windows through ICL, while remaining effective even in the absence of context information. The framework naturally integrates real observations, simulated trajectories, or hybrid combinations thereof, enabling prediction in settings where observational data are sparse or expensive [46, 47]. By inferring hidden proximity to a transition rather than relying on explicit early warning indicators, TipPFN provides a universal, system-agnostic estimate of criticality without requiring system-specific retraining or manual selection of diagnostic metrics. Finally, TipPFN consistently and accurately provides early prediction of tipping characteristics, enabling real-world deployment where early intervention is critical. Limitations In real-world online settings, where observations are sequentially ingested as they become available, future work will explore combining ICL with near-real-time model adaptation, for example through low-rank updates [48] or test-time training [49, 50]. This would allow learned priors to evolve as new observations arrive, rather than remaining fixed after pretraining. Future work will also investigate the spatial dependence of critical transitions, which is central to many high-impact systems, including climate [3, 51] and ecological systems [14].

Acknowledgment JN, HF, PG acknowledge funding, computing, and storage resources from the NSF Science and Technology Center (STC) Learning the Earth with Artificial Intelligence and Physics (LEAP) (Award #2019625) and by NASA under award No 80NSSC25K0062.

References [1] Sheng Fang, Ziyan Wang, Jürgen Kurths, and Jingfang Fan. Tipping points and cascading transitions: Methods, principles, and evidences. arXiv preprint arXiv:2511.01168, 2025. [2] Peter Ashwin, Sebastian Wieczorek, Renato Vitolo, and Peter Cox. Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1962):1166–1184, 2012. [3] Anastasia Romanou, David Rind, Jeff Jonas, Ron Miller, Maxwell Kelley, Gary Russell, Clara Orbe, Larissa Nazarenko, Rebecca Latto, and Gavin A Schmidt. Stochastic bifurcation of the north atlantic circulation under a midrange future climate scenario with the nasa-giss modele. Journal of Climate, 36(18):6141–6161, 2023. [4] Annelies J Veraart, Elisabeth J Faassen, Vasilis Dakos, Egbert H Van Nes, Miquel Lürling, and Marten Scheffer. Recovery rates reflect distance to a tipping point in a living system. Nature, 481(7381):357–359, 2012. [5] Western Electricity Coordinating Council. Western systems coordinating council disturbance report, 1996. [6] Luis Gómez-Nava, Robert T. Lange, Pascal P. Klamser, Juliane Lukas, Lenin Arias-Rodriguez, David Bierbach, Jens Krause, Henning Sprekeler, and Pawel Romanczuk. Fish shoals resemble a stochastic excitable system driven by environmental perturbations. Nature Physics, 19(5):663–669, February 2023. [7] Paul DL Ritchie, Hassan Alkhayuon, Peter M Cox, and Sebastian Wieczorek. Rate-induced tipping in natural and human systems. Earth System Dynamics, 14(3):669–683, 2023. [8] Yunus Sevinchan, Petro Sarkanych, Abi Tenenbaum, Yurij Holovatch, and Pawel Romanczuk. Collective decision-making with heterogeneous biases: Role of network topology and susceptibility. Physical Review Research, 7(1), March 2025. [9] Vasilis Dakos, Egbert H Van Nes, Paolo d’Odorico, and Marten Scheffer. Robustness of variance and autocorrelation as indicators of critical slowing down. Ecology, 93(2):264–271, 2012. [10] C Wissel. A universal law of the characteristic return time near thresholds. Oecologia, 65(1):101– 107, 1984. [11] Egbert H Van Nes and Marten Scheffer. Slow recovery from perturbations as a generic indicator of a nearby catastrophic shift. The American Naturalist, 169(6):738–747, 2007. [12] Vasilis Dakos, Marten Scheffer, Egbert H Van Nes, Victor Brovkin, Vladimir Petoukhov, and Hermann Held. Slowing down as an early warning signal for abrupt climate change. Proceedings of the National Academy of Sciences, 105(38):14308–14312, 2008. [13] Stephen R Carpenter and William A Brock. Rising variance: a leading indicator of ecological transition. Ecology letters, 9(3):311–318, 2006. [14] Vishwesha Guttal and Ciriyam Jayaprakash. Changing skewness: an early warning signal of regime shifts in ecosystems. Ecology letters, 11(5):450–460, 2008. [15] Vasilis Dakos, Stephen R Carpenter, Egbert H van Nes, and Marten Scheffer. Resilience indicators: prospects and limitations for early warnings of regime shifts. Philosophical Transactions of the Royal Society B: Biological Sciences, 370(1659), 2015. [16] Sonia Kéfi, Vasilis Dakos, Marten Scheffer, Egbert H Van Nes, and Max Rietkerk. Early warning signals also precede non-catastrophic transitions. Oikos, 122(5):641–648, 2013. [17] Induja Pavithran and RI Sujith. Effect of rate of change of parameter on early warning signals for critical transitions. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(1), 2021. [18] Yu Huang, Sebastian Bathiany, Peter Ashwin, and Niklas Boers. Deep learning for predicting rate-induced tipping. Nature Machine Intelligence, 6(12):1556–1565, 2024. [19] Jiashuo Liu, Zheyan Shen, Yue He, Xingxuan Zhang, Renzhe Xu, Han Yu, and Peng Cui. Towards out-of-distribution generalization: A survey, 2023. [20] Kaiyang Zhou, Ziwei Liu, Yu Qiao, Tao Xiang, and Chen Change Loy. Domain generalization: A survey. IEEE transactions on pattern analysis and machine intelligence, 45(4):4396–4415, 2022. 11

[21] Xin Wu, Fei Teng, Xingwang Li, Ji Zhang, Tianrui Li, and Qiang Duan. Out-of-distribution generalization in time series: A survey, 2025. [22] Thomas M Bury, RI Sujith, Induja Pavithran, Marten Scheffer, Timothy M Lenton, Madhur Anand, and Chris T Bauch. Deep learning for early warning signals of tipping points. Proceedings of the National Academy of Sciences, 118(39):e2106140118, 2021. [23] Xin Li, Qunxi Zhu, Chengli Zhao, Bolin Zhao, Xue Zhang, Xiaojun Duan, and Wei Lin. Ultraearly prediction of tipping points: Integrating dynamical measures with reservoir computing. arXiv preprint arXiv:2603.14944, 2026. [24] Chengzuo Zhuge, Jiawei Li, and Wei Chen. Deep learning for predicting the occurrence of tipping points. Royal Society Open Science, 12(7), 2025. [25] Shirin Panahi, Ling-Wei Kong, Mohammadamin Moradi, Zheng-Meng Zhai, Bryan Glaz, Mulugeta Haile, and Ying-Cheng Lai. Machine learning prediction of tipping in complex dynamical systems. Physical Review Research, 6(4):043194, 2024. [26] Samuel Müller, Noah Hollmann, Sebastian Pineda Arango, Josif Grabocka, and Frank Hutter. Transformers can do bayesian inference. In International Conference on Learning Representations, 2022. [27] Noah Hollmann, Samuel Müller, Katharina Eggensperger, and Frank Hutter. Tabpfn: A transformer that solves small tabular classification problems in a second, 2023. [28] Noah Hollmann, Samuel Müller, Lennart Purucker, Arjun Krishnakumar, Max Körfer, Shi Bin Hoo, Robin Tibor Schirrmeister, and Frank Hutter. Accurate predictions on small data with a tabular foundation model. Nature, 637(8045):319–326, 2025. [29] Jingang Qu, David Holzmüller, Gaël Varoquaux, and Marine Le Morvan. Tabiclv2: A better, faster, scalable, and open tabular foundation model, 2026. [30] Shi Bin Hoo, Samuel Müller, David Salinas, and Frank Hutter. From tables to time: Extending tabpfn-v2 to time series forecasting, 2026. [31] Jake Robertson, Arik Reuter, Siyuan Guo, Noah Hollmann, Frank Hutter, and Bernhard Schölkopf. Do-pfn: In-context learning for causal effect estimation, 2025. [32] Vahid Balazadeh, Hamidreza Kamkari, Valentin Thomas, Benson Li, Junwei Ma, Jesse C. Cresswell, and Rahul G. Krishnan. Causalpfn: Amortized causal effect estimation via in-context learning, 2025. [33] David Schiff, Ofir Lindenbaum, and Yonathan Efroni. Gradient free deep reinforcement learning with tabpfn, 2025. [34] Yuta Miyauchi, Masahiro Ikeda, and Yoshinobu Kawahara. Generalized stochastic resilience for early warning signals based on koopman operator. Nonlinear Dynamics, 114(4):246, 2026. [35] Florian Grziwotz, Chun-Wei Chang, Vasilis Dakos, Egbert H van Nes, Markus Schwarzländer, Oliver Kamps, Martin Heßler, Isao T Tokuda, Arndt Telschow, and Chih-hao Hsieh. Anticipating the occurrence and type of critical transitions. Science Advances, 9(1):eabq4558, 2023. [36] Étienne Birmelé. A scale-free graph model based on bipartite graphs. Discrete Applied Mathematics, 157(10):2267–2284, 2009. [37] James A Hanley and Barbara J McNeil. The meaning and use of the area under a receiver operating characteristic (roc) curve. Radiology, 143(1):29–36, 1982. [38] Tom Fawcett. An introduction to roc analysis. Pattern recognition letters, 27(8):861–874, 2006. [39] WWTG Peterson, T Birdsall, and We Fox. The theory of signal detectability. Transactions of the IRE professional group on information theory, 4(4):171–212, 1954. [40] Andrew P Bradley. The use of the area under the roc curve in the evaluation of machine learning algorithms. Pattern recognition, 30(7):1145–1159, 1997. [41] Maurice G Kendall. A new measure of rank correlation. Biometrika, 30(1-2):81–93, 1938. [42] W. E. Ricker. Stock and recruitment. Journal of the Fisheries Research Board of Canada, 11(5):559–623, May 1954. [43] John M Drake and Blaine D Griffen. Early warning signals of extinction in deteriorating environments. Nature, 467(7314):456–459, 2010. 12

[44] Giacomo Bonciolini, Dominik Ebi, Edouard Boujo, and Nicolas Noiray. Experiments and modelling of rate-dependent transition delay in a stochastic subcritical bifurcation. Royal Society open science, 5(3), 2018. [45] Sang-Ki Lee, Dongmin Kim, Fabian A Gomez, Hosmay Lopez, Denis L Volkov, Shenfu Dong, Rick Lumpkin, and Stephen Yeager. A pause in the weakening of the atlantic meridional overturning circulation since the early 2010s. Nature Communications, 15(1):10642, 2024. [46] Juan Nathaniel, Jiangong Liu, and Pierre Gentine. Metaflux: Meta-learning global carbon fluxes from sparse spatiotemporal observations. Scientific Data, 10(1):440, 2023. [47] Siyeon Kim, Juan Nathaniel, Zhewen Hou, Tian Zheng, and Pierre Gentine. Spatiotemporal upscaling of sparse air-sea pco2 data via physics-informed transfer learning. Scientific data, 11(1):1098, 2024. [48] Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Liang Wang, Weizhu Chen, et al. Lora: Low-rank adaptation of large language models. Iclr, 1(2):3, 2022. [49] Yu Sun, Xiaolong Wang, Zhuang Liu, John Miller, Alexei Efros, and Moritz Hardt. Test-time training with self-supervision for generalization under distribution shifts. In International conference on machine learning, pages 9229–9248. PMLR, 2020. [50] Arnuv Tandon, Karan Dalal, Xinhao Li, Daniel Koceja, Marcel Rød, Sam Buchanan, Xiaolong Wang, Jure Leskovec, Sanmi Koyejo, Tatsunori Hashimoto, et al. End-to-end test-time training for long context. arXiv preprint arXiv:2512.23675, 2025. [51] Juan Nathaniel, Yongquan Qu, Tung Nguyen, Sungduk Yu, Julius Busecke, Aditya Grover, and Pierre Gentine. Chaosbench: A multi-channel, physics-based benchmark for subseasonal-toseasonal climate prediction. Advances in Neural Information Processing Systems, 37:43715– 43729, 2024. [52] Eusebius J Doedel, Alan R Champneys, Fabio Dercole, Thomas F Fairgrieve, Yuri A Kuznetsov, B Oldeman, RC Paffenroth, B Sandstede, XJ Wang, and CH Zhang. Auto-07p: Continuation and bifurcation software for ordinary differential equations. 2007. [53] Benjamin Herdeanu, Juan Nathaniel, Carla Roesch, Jatan Buch, Gregor Ramien, Johannes Haux, and Pierre Gentine. CausalDynamics: A large-scale benchmark for structural discovery of dynamical causal models. In Advances in Neural Information Processing Systems 38 (NeurIPS 2025), Track on Datasets and Benchmarks, 2025. arXiv preprint arxiv:2505.16620. [54] Juan Nathaniel, Carla Roesch, Jatan Buch, Derek DeSantis, Adam Rupe, Kara D Lamb, and Pierre Gentine. Deep koopman operators for causal discovery. Communications Physics, 8(1):513, 2025. [55] Robert M May. Thresholds and breakpoints in ecosystems with a multiplicity of stable states. Nature, 269(5628):471–477, 1977. [56] Hassan Alkhayuon, Peter Ashwin, Laura C Jackson, Courtney Quinn, and Richard A Wood. Basin bifurcations, oscillatory instability and rate-induced thresholds for atlantic meridional overturning circulation in a global oceanic box model. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 475(2225), 2019. [57] CM Luke and PM Cox. Soil carbon and climate change: from the jenkinson effect to the compost-bomb instability. European journal of soil science, 62(1):5–12, 2011. [58] Sebastian Wieczorek, Peter Ashwin, Catherine M Luke, and Peter M Cox. Excitability in ramped systems: the compost-bomb instability. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467(2129):1243–1269, 2011. [59] Alessio Burrello, Lukas Cavigelli, Kaspar Schindler, Luca Benini, and Abbas Rahimi. Laelaps: An energy-efficient seizure detection algorithm from long-term human iEEG recordings without false alarms. In 2019 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 752–757, 2019. [60] Viktor K. Jirsa, William C. Stacey, Pascale P. Quilichini, Anton I. Ivanov, and Christophe Bernard. On the nature of seizure dynamics. Brain, 137(8):2210–2230, 2014. [61] Patrick Mergell, Hanspeter Herzel, Thomas Wittenberg, Monika Tigges, and Ulrich Eysholdt. Phonation onset: vocal fold modeling and high-speed glottography. The Journal of the Acoustical Society of America, 104(1):464–470, 1998. 13

[62] Preston R Murray and Scott L Thomson. Vibratory responses of synthetic, self-oscillating vocal fold models. The Journal of the Acoustical Society of America, 132(5):3428–3438, 2012. [63] Ryo Shimamura and Isao T Tokuda. Effect of level difference between left and right vocal folds on phonation: Physical experiment and theoretical study. Journal of the Acoustical Society of America, 140(4_Supplement):3393–3394, 2016. [64] Stephan Wagner, Janina Steinbeck, Philippe Fuchs, Sophie Lichtenauer, Marlene Elsässer, Jos HM Schippers, Thomas Nietzel, Cristina Ruberti, Olivier Van Aken, Andreas J Meyer, et al. Multiparametric real-time sensing of cytosolic physiology links hypoxia responses to mitochondrial electron transport. New Phytologist, 224(4):1668–1684, 2019. [65] Ben I Moat, David Smeed, Darren Rayner, William E Johns, Ryan H Smith, Denis L Volkov, Shane Elipot, Tillys Petit, Jules B Kajtar, Molly O Baringer, et al. Atlantic meridional overturning circulation observed by the rapid-mocha-wbts (rapid-meridional overturning circulation and heatflux array-western boundary time series) array at 26n from 2004 to 2023 (v2023. 1). 2024. [66] Michael L Rosenzweig and Robert H MacArthur. Graphical representation and stability conditions of predator-prey interactions. The American Naturalist, 97(895):209–223, 1963. [67] Tamer Oraby, Vivek Thampi, and Chris T Bauch. The influence of social norms on the dynamics of vaccinating behaviour for paediatric infectious diseases. Proceedings of the Royal Society B: Biological Sciences, 281(1780), 2014. [68] A Demetri Pananos, Thomas M Bury, Clara Wang, Justin Schonfeld, Sharada P Mohanty, Brendan Nyhan, Marcel Salathé, and Chris T Bauch. Critical dynamics in population vaccinating behavior. Proceedings of the National Academy of Sciences, 114(52):13762–13767, 2017.

Appendix contents A Additional background

B Prior-Data Fitted Networks

B.1 Driver variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.2 Training Dataset Forcing Schedule . . . . . . . . . . . . . . . . . . . . . . . . . .

B.3 Augmenting SDE System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.4 Task Construction for TipPFN . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.5 TipPFN Model Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.6 Training Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.7 Linear Forcing schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.8 Compute Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C Experiment details

C.1 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.2 Baselines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.3 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.4 Further Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.5 AUROC Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Additional background

Critical transitions. As discussed in the main text, tipping phenomena can be broadly classified into three distinct types depending on the mechanisms driving state transitions, namely bifurcationinduced tipping (b-tipping), noise-induced tipping (n-tipping), and rate-induced tipping (r-tipping) [2]. Refer to Fig. 2 for illustration on the different mechanisms underlying each one of them. Bifurcation-induced tipping and critical slowing down. In b-tipping, critical transitions arise from local bifurcations, where an equilibrium x∗ (λ) loses stability as parameters vary. Linearizing the dynamics around x∗ (λ) yields dϵ = J(λ) ϵ, (S1) dt where ϵ = x − x∗ (λ) and the Jacobian matrix is Jij (λ) =

∂fi (x, λ) . ∂xj x=x∗ (λ)

(S2)

A bifurcation occurs when the linear stability of the equilibrium changes, typically because one or more eigenvalues µi (λ) of J(λ) cross the stability boundary in the complex plane, i.e., max Re(µi (λc )) = 0. i

(S3)

The manner in which eigenvalues approach this boundary constrains the bifurcation type: a real eigenvalue crossing zero is associated with steady-state bifurcations (e.g., fold, transcritical), whereas a complex-conjugate pair crossing the imaginary axis is associated with a Hopf bifurcation. In discrete-time systems, the corresponding stability boundary is the unit circle. As the system approaches such a bifurcation, perturbations decay increasingly slowly, a phenomenon known CSD [1, 12]. In particular, as λ → λcrit , the leading eigenvalue satisfies Re(µi (λ)) → 0, implying a vanishing local contraction rate and increasingly slow recovery from perturbations. 15

Rate-induced tipping. In contrast to b-tipping, r-tipping arises from the inability of the system to track a moving attracting state under time-dependent forcing [2]. Consider again the nonautonomous system with control parameter λ(t). For each fixed λ, assume the associated frozen system admits a locally stable equilibrium x∗ (λ) satisfying f (x∗ (λ), λ) = 0,

max Re(µi (λ)) < 0. i

(S4)

Defining the tracking error

ϵ(t) := x(t) − x∗ (λ(t)), and linearizing around x (λ(t)) yields dϵ d = J(λ(t)) ϵ + x∗ (λ(t)) +higher-order and stochastic terms. dt dt | {z }

(S5)

∗

(S6)

branch drift

d ∗ The additional term dt x (λ(t)) represents the motion of the equilibrium induced by the changing control parameter. R-tipping occurs when this control-induced drift is sufficiently large relative to the local contraction governed by J(λ(t)), so that the trajectory can no longer track the moving equilibrium, even though the instantaneous dynamics remain locally stable. Unlike in b-tipping, no eigenvalue needs to cross the stability boundary, and therefore classical indicators based on CSD may fail to provide reliable early warning signals.

Noise-induced tipping. In n-tipping, transitions occur through stochastic perturbations that induce escape from the basin of attraction of a stable equilibrium, even when Re(µi (λ)) < 0. This mechanism can act independently or interact with both b-tipping and r-tipping in stochastic systems. Early warning signals. CSD manifests in observable time series statistics [9, 13]. For a discrete trajectory {xt }Tt=1 , the lag-1 autocorrelation (AR1) is defined as PT (xt − x̄)(xt−1 − x̄) AR1 = t=2PT , (S7) 2 t=1 (xt − x̄) PT where x̄ := T1 t=1 xt denotes the empirical mean over the time window. To connect AR1 and variance to local stability, consider the scalar linearized stochastic dynamics xt+1 − x̄ = a(xt − x̄) + σξt ,

(S8)

where {ξt } is an i.i.d. zero-mean stochastic process and |a| < 1 ensures stability. Under stationarity and independence of the innovations σξt , AR1 ≈ a. If the system arises from a continuous-time linearization with leading eigenvalue µ < 0 and sampling interval ∆t, then a = eµ∆t , so AR1 → 1 as µ → 0− . Similarly, the stationary variance satisfies VAR(x) =

σ2 , 1 − a2

(S9)

which increases as a → 1. A detailed derivation of the relationship between AR1, variance, and the underlying spectral properties of the linearized dynamics is provided in e.g., [35, 34]. These signatures, rising AR1 and variance, form the basis of classical EWS for b-tipping. However, these indicators rely on local linearization near equilibrium and are therefore intrinsically tied to eigenvalue dynamics. In n-tipping, transitions occur due to stochastic escape from a basin of attraction even when Re(µi ) < 0, while in r-tipping, rapid parameter changes prevent the system from tracking its quasi-static equilibrium without any eigenvalue crossing. In both cases, critical slowing down may be weak or absent. While bifurcation theory provides a principled framework for characterizing tipping behaviour, its practical application is limited. Identifying tipping points requires access to the governing equations, parameters, and equilibrium states needed to compute the Jacobian and its eigenvalues. In real-world systems, these quantities are typically unknown or only partially observed, motivating data-driven approaches that infer tipping behaviour directly from observed trajectories without requiring explicit knowledge of the underlying dynamical system. 16

B a

Prior-Data Fitted Networks Inference

Prior-fitting

query

𝑞𝜃 (̇Λreal 𝑡 |𝑜1:3

Generate synthetic samples prior ψ ∼ p(ψ) from prior over generative processes Real dataset 𝑜t with 𝒞 = ∅

Sample episodes from generative process 𝒮 ∼ p(𝒮|ψ) based on individual noise realizations and forcing schedules Sample context-query pair 𝒟 from 𝒟 ∼ p(𝒮|ψ)

Full visibility

Real or “sim-toreal” contextquery pairs

con.,𝑗

query

𝑜2

Train TipPFN q𝜃 on contextquery pairs {𝒟1, … , 𝒟n}

Context 𝒞

, 𝒞)

Query

real q𝜃 (Λreal 𝑜1:𝑡 , 𝒞 sim/real (nowcasting) 𝑡 real sim/real q𝜃 (Λreal (forecasting) 𝑡:𝑇 𝑜1:𝑡 , 𝒞

query

𝑜1

Λ𝑚

query

, 𝑡2

query

, 𝑡1

con.,j

, 𝑜𝑚

con.,j

, 𝑡𝑚

Λcon.,𝑘 , 𝑜𝑛con.,𝑘 , 𝑡𝑛con.,𝑘 𝑛 query

𝑜3

query

Λcon.,𝑙 , 𝑜𝑝con.,𝑙 , 𝑡𝑝con.,𝑙 𝑝

, 𝑡3

Causal attention within query episode

Full attention within context

Figure S1: (a) Prior-data fitting and inference and (b) attention structure of TipPFN. Schematics based on [26].

B.1

Driver variables

Following [22], the driver M is generated from randomly constructed two-dimensional dynamical systems of the form

ẋ =

10 X

ai pi (x, y),

ẏ =

i=1

10 X

bi pi (x, y),

(S10)

i=1

where (x, y) ∈ R2 and {pi (x, y)}10 i=1 denotes the set of all monomials up to third order: p(x, y) = (1, x, y, x2 , xy, y 2 , x3 , x2 y, xy 2 , y 3 ).

(S11)

The coefficients ai , bi are independently sampled from N (0, 1), after which a random subset (50%) is set to zero to induce sparsity. To encourage bounded trajectories, coefficients associated with cubic terms are constrained to be negative. We then perform integration for 104 timesteps and 10−2 discretization to select trajectories that converge to an equilibrium, defined as the final 10 points with difference of less than 10−8 . We then use AUTO-07P program [52] to identify bifurcation types and their corresponding critical values along the identified equilibrium branch as each nonzero parameter is varied within the interval [−5, 5]. Table S1 lists the relevant driver-generation hyperparameters. Table S1: Hyperparameters of the synthetic polynomial driver. Quantity

Value or distribution

Polynomial degree

All monomials up to degree 3; 10 coefficients per equation, 20 total. aj , bj ∼ N (0, 1) before sparsification. Exactly 50% of the 20 coefficients are set to zero uniformly at random. Cubic coefficients x3 , x2 y, xy 2 , y 3 in both equations are replaced by their negative absolute values. z0 ∼ N (0, 4I2 ). Forward Euler with ∆t = 0.01 for up to 100 time units. Norm of the range of the final 10 states below 10−8 and maximum trajectory amplitude below 103 . All Jacobian eigenvalues at the equilibrium have negative real part. rM = | maxi Re µi | at the stable equilibrium. Up to 100 random polynomial models per TipBox attempt. Up to 50 TipBox attempts per raw sample; AUTO output suppressed.

Coefficient prior Coefficient sparsity Boundedness bias Initial condition for equilibrium screening Equilibrium screening integration Convergence criterion Stability criterion Recovery rate Model retry budget Bifurcation retry budget

For a selected bifurcation point, let p0 be the initial coefficient value and pcrit the AUTO-07P bifurcation value. The schedule λ̃k (t) is applied as λM,k (t) = p0 + λ̃k (t) pcrit − p0 . (S12) The RDTC target is Λk (t) = 1 − λ̃k (t). (S13) Thus Λ = 0 marks the bifurcation point, Λ > 0 is subcritical, and Λ < 0 means the schedule has passed the detected critical value. Driver SDE simulation Given M and a forcing schedule, TipBox simulates the stochastic driver with Euler–Maruyama at integration step ∆t = 0.01 and output sampling interval 1.0. A 100-timeunit burn-in is run at the initial parameter before the recorded trajectory starts. The additive driver noise scale is √ σtilde = 0.01, ξ ∼ Triangular(0.75, 1.0, 1.25), (S14) σM = 2rM σtilde ξ, √ and Brownian increments are scaled by σM ∆t. B.2

Training Dataset Forcing Schedule

For the training dataset, forcing schedules are sampled independently per episode, while the bifurcation system M is shared across the K = 6 episodes of a raw sample. Table S2 summarizes the full schedule prior. All schedules start at λ̃ = 0. Examples of bezier-type schedules are shown in Fig. S2.

Figure S2: Bezier-type forcing schedules used in generation of the training dataset. Note that this partially non-linear and beyond-criticality forcing is unlike that used in validation, see Section B.7. During training, it is causally permissible to also show non-linear forcings, through which the model may learn to recognise a wider range of dynamical system responses to changes in forcing. B.3

Augmenting SDE System

The two-dimensional driver is embedded into a higher-dimensional stochastic system so that TipPFN does not only see canonical normal-form trajectories. The auxiliary variables uk (t) = (u1,k (t), . . . , u14,k (t)) solve an Ito SDE with diagonal diffusion, dui,k (t) = [(1 − ηi )fi (ri,k (t)) − γi ui,k (t)] dt + σi dWi,k (t),

(S15)

where ri,k (t) concatenates the parent auxiliary variables and external-input channels selected by G for variable i. The functions fi are sampled MLPs, one per dynamic variable. The same graph, MLPs, and SDE parameters are shared across the episodes of a raw sample, while initial states, Brownian noise, TipBox noise, and forcing schedules vary by episode. Directed graph sampler The graph G, which is structuring the interactions within the variables of the auxiliary system [53], is generated by a scale-free graph model based on bi-partite graphs [36]: For each raw sample, the code first samples a directed bipartite graph between variable nodes V and interaction nodes H, with |V | = 14 and |H| = 8. For every variable i, the desired outgoing and incoming interaction degrees Di+ and Di− are sampled independently from the truncated power law d−2 Pr(D = d) = P8 , −2 ℓ=1 ℓ 18

d ∈ {1, . . . , 8}.

Table S2: Forcing schedule prior. Final levels below 1 create subcritical approaches, level 1 reaches the bifurcation, and levels above 1 pass the bifurcation. Schedule

Probability

Parameters

Definition and effect

bezier_and_hold

0.4

start=0 smooth_hold=true h ∼ U (0.6, 0.8) L = 1 w.p. 1/2; else L ∼ U (0.6, 1.2) nc ∈ {3, 4, 5, 6} with ratios 2 : 3 : 2 : 1 control range [0.1L, 0.9L]

Bezier ramp from 0 to final level L by fraction h, followed by a hold at L. Interior control-point values are clipped to [0, L]; duplicating the final control point makes the join to the hold phase approximately smooth.

ramp_and_hold

0.4

h ∼ U (0.6, 0.8) L = 1 w.p. 1/2; else L ∼ U (0.6, 1.2)

Linear ramp from 0 to L until fraction h, followed by a constant hold at L.

constant

0.2

λ̃(t) = 0

Stationary no-forcing baseline at the initial stable equilibrium, with Λ(t) = 1.

Conditional on these sampled degrees, each candidate edge i → h is included with probability Di+ /8, and each candidate edge h → i with probability Di− /8, independently across interaction nodes. The bipartite graph is then projected to a directed variable graph: an edge i → j is added whenever there exists an interaction node h such that i → h and h → j. Multiple such paths collapse to one directed edge, and self-loops are retained. After this projection, the canonical drivers M are added as external input channels. Each external-input node connects to each auxiliary variable independently with probability pinp = 0.8. MLP hyperparameters are given in table S3 Table S3: Parameter hyperparameters for the auxiliary SDE. These are sampled once for the shared SDE system inside a raw sample. Quantity in Eq. S15

Value or distribution

MLP depth Hidden width Activation Weight and bias initialization scale

2 linear layers 32 hidden units tanh For each MLP, sw = 1.0 + HalfCauchy(1.0) resampled until sw ≤ 1.5; weights and biases use N (0, s2w ). For each MLP, sp = 0.5 − HalfCauchy(0.5) resampled until sp ≥ 0.1; hence sp ∈ [0.1, 0.5]. Hidden-layer weights only; retained weights are scaled by (1 − sp )−1/2 . With 2 layers, there is no interior hidden-to-hidden layer, so this setting is effectively inactive for variable MLPs. 0.0. γi = 20 |N (0, 1)|. ηi = 0.5 |N (0, 1)|; the code multiplies the MLP flow by 1 − ηi . σi = 0.05 |N (0, 1)| with independent diagonal noise.

Hidden-layer sparsity Sparsity application

Preactivation noise Self-loop coefficient Friction coefficient Diffusion coefficient

B.4

Task Construction for TipPFN

Each raw sample is transformed into a masked tabular time-series task. Row selection by resampling A subset of time steps is selected based on stratified jitter, to yield query time-series with lengths between 92 and 256 time steps, constrained to a total context size of 320 time steps or less, to reserve 192 of the 512 step budget for the query trajectory and potential forecasting time steps, with indices 1 : T ∗ referring to the resampled time steps. Using Λ1:T ∗ = 1− λ̃1:T ∗ , this yields the query-context pair D = ((Λ1:T ∗ , o1:T ∗ )∪C) with prior distribution p(D) (see Fig. 3), hyperparameters are given in Table S4. 19

Table S4: Episode partitioning and temporal subsampling. Quantity

Value or distribution

Maximum sequence length Query episodes Context episodes Episode selection Context length budget Context per-episode bounds Query length budget Temporal subsampling Overflow handling Identifier normalization

512 rows. Always 1 episode. Nctx ∈ {0, 1, 2, 3} with probabilities (0.2, 0.3, 0.3, 0.2). Deterministic-random selection; query and context episodes do not overlap. Nctx = 0 : [0, 0], 1 : [92, 256], 2 : [184, 320], 3 : [276, 320] total rows. Minimum 92 rows, maximum 256 rows. Target 192 rows, with per-episode bounds [1, 256]. Stratified jitter with endpoints retained. Shrink context; use available episodes if insufficient; error on empty query. Episode ID and time are min-max normalized to [−1, 1] over valid rows.

Column selection The dynamic variables are partitioned into target columns, which the model may be asked to predict, and observation columns. RDTC is always the first target column. Additional targets are sampled from the auxiliary time series u, and the drivers x and y. Feature columns are sampled from the remaining var. The number of target columns is sampled from {1, . . . , 16} with probabilities proportional to

ρact =

(1.00, 0.95, 0.90, 0.86, 0.82, 0.78, 0.74, 0.70, 0.67, 0.64, 0.61, 0.58, 0.55, 0.52, 0.49, 0.46).

The number of feature columns is sampled from {1, . . . , 20} with probabilities proportional to

ρfeat =

(1.00, 0.97, 0.94, 0.91, 0.88, 0.85, 0.82, 0.79, 0.76, 0.73, 0.70, 0.67, 0.64, 0.61, 0.58, 0.55, 0.52, 0.49, 0.46, 0.43).

Masking RDTC is always masked for every query row, making Λ(t) the primary hidden target. The additional query task is sampled with a 50/50 mixture. In the none task, no additional query mask is applied; these views train nowcasting of current distance to criticality from currently observed signals. In the forecast task, a cutoff fraction c ∼ U (0.2, 0.4) is sampled and further selected action columns are masked after that cutoff. Normalization and missing values Normalization is applied after partitioning and masking, using only information available under inference-time constraints. Raw RDTC is transformed to the model target Λ∗ = tanh(5Λ). (S16) All other columns are median-centered and transformed with an adaptive inverse-hyperbolic-sine map, x − median(x) x 7→ asinh , s = max(0.01, 0.5 IQR(x)). (S17) s The median and IQR are computed from the observed rows only, excluding masked prediction targets and excluding RDTC from the statistics. Training targets The model predicts 99 quantiles for every masked action value and is optimized with pinball loss. For distributional decoding and validation metrics, quantile outputs are sorted to enforce non-crossing quantiles, and exponential tails are used for extrapolation outside the represented quantile range. RDTC targets receive unit loss weight, while auxiliary signal reconstruction targets receive weight 0.2. This makes distance-to-criticality prediction the primary objective while using reconstruction of observed dynamics as an auxiliary regularizer. 20

B.5

TipPFN Model Architecture

The trained model is a transformer architecture. A task is represented as a padded table with row index r = 1, . . . , R and column index split into action columns and feature columns. Each row has two identifier values, episode id and time, denoted ιr = (er , tr ). Let yr,a denote the value in action slot a ∈ {1, . . . , A} and xr,f the value in feature slot f ∈ {1, . . . , F }. The first action slot is reserved for RDTC, yr,1 = Λ∗r after the tanh transform in Equation S16. A Boolean mask mr,a indicates which action values are prediction targets. The model uses separate scalar encoders for action and feature values, Ey : R → Rdmodel ,

Ex : R → Rdmodel ,

and fixed random positional buffers for rows, action columns, and identifier/feature columns. For active, unmasked slots the initial tokens are therefore row h(0) + pact r,a = Ey (yr,a ) + pr a ,

(0)

hr,A+g = Ex (vr,g ) + prow + pfeat r g ,

where a = 1, . . . , A, g = 1, . . . , F + 2, and vr = (er , tr , xr,1 , . . . , xr,F ). Masked action tokens are replaced by a zero vector before entering the transformer. Table S5: Resolved TipPFN architecture hyperparameters. Quantity

Value

Maximum rows R Identifier columns Maximum feature columns F Maximum action columns A Total token columns per row Hidden width dmodel Transformer blocks Attention heads MLP width Dropout Output dimension Quantile levels Row attention mode Causal time tolerance Learned parameters Fixed positional buffers

512. 2: episode id and time. 20 payload features, plus the 2 identifier columns in the feature stream. 16. 16 + 2 + 20 = 38. 1280. 8 row-column blocks. 8 heads; head dimension 160. Ratio 4; hidden width 5120. 0.1 in attention and MLP residual branches during training. 99 per action token. αj = j/100, j = 1, . . . , 99. Causal. 10−6 . 209,959,779, excluding fixed positional buffers. 704,000 scalar entries: 512 × 1280 row, 22 × 1280 id/feature-column, and 16 × 1280 action-column buffers.

Each row-column block applies column attention, row attention, and an MLP with pre-normalization and residual connections: H ← H + Dropout Attncol (LN(H)) , H ← H + Dropout Attnrow (LN(H)) , H ← H + Dropout MLP(LN(H)) . Column attention is applied independently within each row across the active action, identifier, and feature columns. Row attention is applied independently for each column across the temporal/episode rows. Attention uses fused scaled dot-product attention with bias-free query, key, value, and output projections. The MLP is a two-layer feed-forward network with GELU activation using the tanh approximation. Causal row mask Let ci indicate that row i belongs to a context episode. In causal mode, the set of row keys visible to query row i is {j : cj = 1} ∪ {i}, ci = 1, A(i) = (S18) {j : cj = 1} ∪ {j : ej = ei , tj < ti − 10−6 } ∪ {i}, ci = 0. Padded keys are removed from A(i). Thus context rows can exchange information bidirectionally within the context set, while query rows can attend to all context rows and only to their own past and present query rows. 21

B.6

Training Configuration

Training uses the synthetic task distribution described above. For an output vector q̂r,a,1:99 and target yr,a , the pinball loss at quantile level α is α(y − q̂), y ≥ q̂, ℓα (q̂, y) = (α − 1)(y − q̂), y < q̂. The per-entry loss is the mean of ℓαj over the 99 quantile levels. Let Mrdtc be the masked RDTC entries, and let Msig be all other masked target entries. The training loss is L = 1.0 mean(r,a)∈Mrdtc ℓ(q̂r,a , yr,a ) + 0.2 mean(r,a)∈Msig ℓ(q̂r,a , yr,a ), with empty groups skipped. Training hyperparameters are given in table S6. Table S6: Resolved optimization and trainer hyperparameters. Quantity

Value

Per-device batch size Global batch size Data-loader workers Precision Optimizer Learning rate Weight decay

40 transformed views. 80 transformed views across 2 DDP workers. 1 worker, persistent_workers=false, pin_memory=true. bf16-mixed. Muon for 2D trainable tensors, AdamW for remaining tensors. 0.01. 0.01 for decay groups; no decay for biases, 1D norms, and parameters matching no-decay keywords. Momentum 0.95; Nesterov enabled; Newton-Schulz steps 5; ϵ = 10−7 . β1 = 0.9, β2 = 0.999, ϵ = 10−7 . Linear warmup followed by cosine decay. 1% of 200,000 steps, i.e. 2,000 optimizer steps. 0.1, giving final LR 0.001 from the initial LR 0.01. 200,000. Global norm clipping at 1.0. Lightning seed_everything=true.

Muon parameters AdamW parameters Scheduler Warmup Minimum LR scale Maximum optimizer steps Gradient clipping Random seeding

B.7

Linear Forcing schedules

In the following we provide more details on the forcing schedule to generate the various tipping scenarios used for validation datasets, see AppendixC.1. Note that these differ from the non-linear forcing schedules used in the training dataset. For each episode in the synthetic and semi-real b-type systems, the a forcing schedule was randomly selected. To gain a balanced distribution of positive (critical) and negative (non-critical) episodes, an episode was grouped as critical with 50% probability. If not critical, one of the four non-critical groups described in Table S7 was uniform randomly selected. The goal of these forcing schedules was to provide a realistic evaluation task: instead of just distinguishing between critical and equilibrium trajectories, the evaluation task includes distinguishing non-critically forced b-type systems (approaching criticality, but not reaching it) from actual critical forcings. We employ exclusively linear forcing schedules for all experiments; changing the forcing at an unobserved future time would be causally intractable for a prediction model [54]. As figures S3 and S4 show, these normalized forcing schedules Λ(t) = 1 − λ̃(t) also varied in the initial starting point and the time of reaching criticality. To avoid strong perturbation, the tipping system was warmed-up by slowly and smoothly varying the control parameter from equilibrium to the parameter that is the start of the forcing schedule (not shown). R-tipping The notion of a normalized forcing schedule does not easily translate to r-type tipping systems, which require to surpass not a certain forcing value but a critical change in forcing λ̇. In particular, there is no straight-forward equivalent of the approaching or receding groups. To still 22

v3 synth | RDTC-slope histograms

a critical n = 4; 060 (49.6%)

equilibrium n = 1; 042 (12.7%)

flat n = 1; 034 (12.6%)

receding n = 1; 007 (12.3%)

approaching n = 1; 049 (12.8%)

¡1

1.0

0.5

0.0

b count

200

100

¡1

¢¤=episode

¤ = 0:05

¢¤=episode

80% interval

¢¤=episode

mean linear trend

Figure S3: Forcing schedules Λ(t) and overall group sizes (number of episodes) and composition for the canonical evaluation datasets (b_fold, b_hopf, b_transcritical). Refer to Table S11 for the actually used episode counts out of this larger dataset. v3 systems-bif | RDTC-slope histograms

a critical n = 4; 115 (50.2%)

equilibrium n = 999 (12.2%)

flat n = 1; 025 (12.5%)

receding n = 990 (12.1%)

approaching n = 1; 063 (13.0%)

¡1

1.0

0.5

0.0

b count

200

100

¡1

¢¤=episode

¤ = 0:05

¢¤=episode

80% interval

¢¤=episode

mean linear trend

Figure S4: Forcing schedules Λ(t) and overall group sizes (number of episodes) and composition for the semi-real bifurcation systems (see Table S8). Refer to Table S11 for the actually used episode counts out of this larger dataset. create some variability in schedules, we split the non-critical group into equilbibrium and flat groups, where the flat group has a higher but still non-critical maximum forcing rate.

Table S7: Trajectory classes defined by the evolution of the normalized control parameter λ̃(t) and their corresponding Λ(t) behavior. Trajectory type Definition of λ̃(t) Λ(t) behavior λ̃0 + αt, t < tc , Critical λ̃(t) = α>0 Λ(t) → 0. Critical value (tipping) 1, t ≥ tc , λ̃ = 1 is reached at tc . Approaching λ̃(t) = λ̃0 + αt, λ̃(t) < 1 − ϵ ∀t, ϵ > 0 Λ(t) → ϵ. The tra(non-critical) jectory remains bounded away from criticality. Receding λ̃(t) = λ̃0 − αt, λ̃0 < 1 ∀t Λ(t) → 1. The trajectory (non-critical) moves away from the critical point. Flat λ̃(t) = λ̃0 + δ(t), |δ(t)| ≪ 1 Λ(t) ≈ const. (non-critical) Equilibrium λ̃(t) = 0 ∀t Λ(t) ≈ 1. System remains far from criticality. B.8

Compute Resources

All experiments were conducted on rented GPU servers (IONOS cloud) and a single dedicated CPU server (Hetzner, Germany). We report hardware specifications, per-stage compute, and total project compute below. • Dataset generation. The synthetic training dataset was generated on a single CPU server (AMD EPYC 9454P, 48 cores / 96 threads, 125 GB RAM) over approximately 8 days at ∼ 80% CPU utilisation (∼ 7,400 CPU core-hours). • Model training (reported). The checkpoint used in all experiments was trained on 4× NVIDIA H200 NVL GPUs (141 GB HBM3e each, 1 TB system RAM) for approximately 47 hours (∼ 189 H200 GPU-hours) at near-full utilisation (∼ 100% SM occupancy, ∼ 493 W average draw per device). • Evaluation (reported). All evaluation runs used NVIDIA H200 NVL GPUs (1–4 per run). The paper-bound evaluation runs (7 datasets, including ablations and step-size variants) consumed approximately ∼ 150 GPU-hours in total. Individual runs ranged from ∼ 1 h (lean / fast configurations) to ∼ 20 h (full SWEC and TAC sweeps on 4× H200), with TabPFN and Bury baselines requiring considerable CPU and GPU time. • Total project compute. The full research project, including preliminary and failed experiments not reported in the paper, consumed substantially more compute than the final experiments. Over the course of development, the training project logged ∼ 1,400 GPUhours across 44 runs (spanning smaller NVIDIA RTX 6000 Ada, RTX A6000, RTX Pro class GPUs). For evaluation, we logged ∼ 750 GPU-hours across 347 runs (all on H200). Including dataset generation, we estimate the total project compute at ∼ 2,150 GPU-hours plus ∼ 7,400 CPU core-hours.

Experiment details

C.1

Datasets

Summary Table S8 summarizes the datasets included in TipBox as well as additional benchmark datasets for evaluation on real-world examples. Sample tipping and non-tipping trajectories, along with their specific instantiation of the bifurcating parameter or forcing rate schedules are shown in Figures S5 and S6. Table S8: Datasets used for evaluating tipping point prediction across b-tipping (B), r-tipping (R), and n-tipping (N), spanning canonical, semi-real, sim-to-real, and real-world systems. Dataset

System

Description

Type

References

Canonical

b_fold

[22]

Canonical Canonical

b_hopf b_transcritical

Fold (saddle-node) bifurcation Hopf bifurcation Transcritical bifurcation

B B

[22] [22]

Semi-real

b_harvesting

[55]

Semi-real

b_rosenzweig_macarthur_tc

[22]

Semi-real

b_rosenzweig_macarthur_hopf

[22]

Semi-real

b_seirx_tc

[22]

Semi-real Semi-real Semi-real

b_amoc r_bautin r_compost_bomb

B R R

[7, 56] [7] [57, 58]

Semi-real

r_saddle_node

[2]

Semi-real

r_amoc

[7, 56]

Sim-to-real

TAC

[44]

Sim-to-real

DaphniaExt

Thermoacoustic instability experiments Ecological extinction events

[43]

Real

SWEC-iEEG

[59, 60]

Real

microcosm

B/N

[4]

Real

voice

[61, 62, 63]

Real

cellular_atp

B/N

[64]

Real

greenhouse_earth

[12]

Real

blackout_frequency

[5]

Real

AMOC

Epileptic seizure onset (iEEG) Cyanobacteria under light stress Phonation onset under increasing flow (Hopf) Cellular energy status under hypoxia Calcium carbonate in the end of greenhouse Earth Bus voltage frequency before power grid failure Observed variability of the meridional overturning circulation (MOC)

R/B

[65]

May’s harvesting model Consumer-resource (transcritical) Consumer-resource (Hopf) SEIRx diseasevaccination AMOC 2-box model Bautin normal form Compost bomb instability Saddle-node normal form AMOC 2-box model

b_fold

The fold (saddle-node) bifurcation is described by the normal form

dx = µ − x2 , (S19) dt where two equilibria collide and annihilate at µ = 0, leading to an abrupt transition. This bifurcation is associated with catastrophic shifts, where the system loses stability and jumps to a distant attractor. b_hopf

The Hopf bifurcation involves a two-dimensional system of the form

dx = µx − y − x(x2 + y 2 ), dt (S20) dy = x + µy − y(x2 + y 2 ), dt where a pair of complex conjugate eigenvalues crosses the imaginary axis at µ = 0, leading to the emergence of a stable or unstable limit cycle. This transition results in oscillatory dynamics rather than a shift between steady states. b_transcritical

The transcritical bifurcation is given by

dx = µx − x2 , (S21) dt in which two equilibria exist for all parameter values and exchange stability at µ = 0. In contrast to the fold case, the transition is continuous and does not involve the disappearance of equilibria, but rather a reorganization of stability.

(a) May’s harvesting

(b) RosenzweigMacArthur

(d) Bautin

Figure S5: Representative B-tipping (a, b) and R-tipping (c, d) examples. In each panel, the upper plot shows examples of tipping (red) and non-tipping (blue) trajectories, and the lower plot shows the corresponding parameter schedule. Note that the validation datasets use a more diverse set of forcing schedules, see Figures S3 and S4, which vary the initial value and slope (for b-type systems only) as well as the timing of the critical transition (for all systems, including r-type). b_harvesting We consider May’s single-species harvesting model [55] describing the dynamics of a biomass variable x(t) ∈ R+ . The stochastic dynamics are given by dx x x2 = rx 1 − − h(t) 2 + σx ξ(t), (S22) dt k s + x2 where r is the intrinsic growth rate, k is the carrying capacity, s controls the saturation scale of harvesting, h(t) is the harvesting rate, and ξ(t) denotes Gaussian white noise with amplitude σx . The logistic growth term promotes population persistence, while the nonlinear harvesting term induces a destabilizing feedback. As the harvesting rate h increases, the system undergoes a fold (saddle-node) bifurcation at a critical value h⋆ , beyond which the stable positive-biomass equilibrium disappears, leading to an abrupt collapse of the population. 26

We use parameters r = 1, k = 1, s = 0.1, and σx = 0.01, for which the deterministic system exhibits a fold bifurcation at h⋆ ≈ 0.26. Tipping trajectories are generated by ramping the harvesting rate from h = 0.15 to h = 0.27, crossing the bifurcation threshold, while non-tipping trajectories follow the same initial ramp but are capped at h = 0.25, remaining below the critical point. Tipping is identified as the collapse of the biomass toward the low-density state. b_rosenzweig_macarthur_tc and b_rosenzweig_macarthur_hopf We consider the Rosenzweig– MacArthur consumer–resource model [66, 22] with state x(t) = (x(t), y(t))⊤ ∈ R2+ , where x(t) denotes the resource population and y(t) the consumer population. The stochastic dynamics are given by dx a(t) xy x − = rx 1 − + σx ξx (t), (S23) dt k 1 + a(t)hx dy e a(t) xy = − my + σy ξy (t), (S24) dt 1 + a(t)hx where r is the intrinsic growth rate of the resource, k its carrying capacity, a(t) the consumer attack rate, e the conversion efficiency, h the handling time, and m the consumer mortality rate. The noise terms ξx (t) and ξy (t) are independent Gaussian white noise processes with amplitudes σx and σy . As the attack rate a increases, the system exhibits two bifurcations that structure its dynamical regimes. At a ≈ 5.60, a transcritical bifurcation marks the consumer invasion threshold: for a < 5.60, the system converges to a consumer-free equilibrium (x∗ , 0), while for a > 5.60, a stable coexistence equilibrium (x∗ , y ∗ ) emerges. As a increases further, the coexistence equilibrium remains stable until a Hopf bifurcation at a⋆ ≈ 15.69, where a pair of complex conjugate eigenvalues crosses the imaginary axis, leading to the onset of predator–prey cycles. For a > a⋆ , the system exhibits stable limit cycle behaviour. Initial conditions are sampled near the equilibrium corresponding to the initial parameter value. In the transcritical regime, trajectories are initialized near the resource-only equilibrium, with x(0) close to carrying capacity and a small but positive consumer density y(0) > 0, since y = 0 is always an invariant equilibrium and would otherwise prevent consumer invasion even after crossing the transcritical threshold. In the Hopf regime, initial conditions are sampled near the stable coexistence equilibrium at the ramp start, ensuring trajectories begin on the equilibrium branch that later destabilizes. We use parameters r = 4, k = 1.7, e = 0.5, h = 0.15, and m = 2. For transcritical tipping, we generate trajectories by ramping a from below to above a ≈ 5.60, enabling the transition from a consumer-free to a coexistence state. For Hopf tipping, trajectories are generated by ramping a from 12 to 17, thereby crossing the Hopf threshold, while non-tipping trajectories are capped below the respective bifurcation points. Tipping is identified either as successful consumer invasion (transcritical) or as the transition from equilibrium dynamics to persistent oscillations (Hopf). b_seirx_tc We consider a stochastic SEIRx model [67, 68, 22] capturing the coupled dynamics of infectious disease transmission and vaccination behaviour. The state is given by x(t) = (S(t), E(t), I(t), x(t))⊤ , where S, E, and I denote the susceptible, exposed, and infectious populations, respectively, and x(t) represents the fraction of individuals with provaccine sentiment. The recovered population is given by R(t) = N − S − E − I. The dynamics are dS SI = µN (1 − x) − β − µS + σS ξS (t), dt N dE SI =β − (ϵ + µ)E + σE ξE (t), dt N dI = ϵE − (γ + µ)I + σI ξI (t), dt dx = κx(1 − x) −ω(t) + I + δ(2x − 1) + σx ξx (t), dt

(S25) (S26) (S27) (S28)

where µ is the birth/death rate, β the transmission rate, ϵ the exposed-to-infectious rate, γ the recovery rate, κ the social learning rate, δ the strength of injunctive social norms, and ω(t) the perceived risk 27

of vaccination relative to infection. The noise terms ξS , ξE , ξI , ξx are independent Gaussian white noise processes. As the perceived vaccination risk ω increases, the system undergoes a transcritical bifurcation at ω = δ, corresponding to a loss of stability of the disease-free, high-vaccination equilibrium. Below the threshold (ω < δ), the system remains in a regime with high vaccination uptake (x ≈ 1) and low infection prevalence. Above the threshold (ω > δ), vaccination sentiment declines, leading to the emergence of endemic infection. This transition from disease-free to endemic dynamics constitutes bifurcation-induced tipping. We use parameters N = 100,000, µ = 0.02/52, β = 10.5, ϵ = 0.7, γ = 0.7, κ = 0.007, and δ = 50, corresponding to a typical pediatric infectious disease. Noise amplitudes are σS = σE = σI = 5 and σx = 5 × 10−4 . Tipping trajectories are generated by ramping ω from 0 to 100, thereby crossing the bifurcation threshold, while non-tipping trajectories remain below ω = δ. Tipping is identified as a collapse in vaccination sentiment (x(t) decreasing) accompanied by a resurgence of infection (I(t) increasing). r_saddle_node We consider a stochastic saddle-node normal form with scalar state x(t) ∈ R and time-dependent forcing λ(t). The dynamics are given by √ dx = (x + λ(t))2 − 1 + 2σx ξ(t), dt

(S29)

where σx controls the noise magnitude and ξ(t) denotes Gaussian white noise. The external forcing evolves according to λmax λmax ϵt λ(t) = tanh +1 , (S30) 2 2 where λmax sets the forcing amplitude and ϵ determines the rate of change. For each fixed λ, the frozen system exhibits the characteristic saddle-node structure with a stable and unstable equilibrium that collide at a critical parameter value. For sufficiently slow forcing (ϵ small), trajectories track the attracting equilibrium branch. However, when the forcing rate exceeds a critical threshold, the system fails to follow this quasi-static equilibrium and undergoes a rapid transition, even though the instantaneous system remains locally stable. This behaviour constitutes rate-induced tipping. We use λmax = 3 and compare a tipping regime with ϵ = 1.25 to a non-tipping regime with ϵ = 0.625, with noise level σx = 0.008. For analysis, we consider the scalar observable x(t) and define tipping as the event x(t) ≥ 0. r_bautin We consider a stochastic Bautin system with complex state z(t) = x(t) + iy(t) and time-dependent forcing Λ(t). The dynamics follow a shifted Bautin normal form dz = (a + iω) z − Λ(t) − b |z − Λ(t)|2 z − Λ(t) + |z − Λ(t)|4 z − Λ(t) + σz ξz (t), (S31) dt where a controls linear growth, ω is the angular frequency, b is the cubic coefficient, σz sets the noise amplitude, and ξz (t) is complex-valued white noise. The external forcing evolves according to Λmax Λmax rt Λ(t) = tanh +1 , (S32) 2 2 where Λmax determines the forcing magnitude and r controls the rate of change. The forcing Λ(t) translates the underlying Bautin system through state space. For sufficiently large r, the system fails to track the quasi-static attractor and undergoes rate-induced tipping, even though the instantaneous dynamics remain locally stable. We use parameters a = 0.1, ω = 3, b = 1, σz = 0.2, and Λmax = 8. Tipping runs are generated with r = 0.10, while non-tipping runs use r = 0.05. For analysis, we consider the scalar observable p ρ(t) = x(t)2 + y(t)2 , and define tipping as the event ρ(t) ≥ 10. 28

(S33)

r_compost_bomb We consider a stochastic compost-bomb model describing thermally driven instability in a coupled temperature–carbon system [57, 58]. The dynamics are given by dTs = −λ(Ts − Ta (t)) + ACs r0 e(α(Ts −Tref )) + D3 ξ(t), dt dCs = Π − Cs r0 e(α(Ts −Tref )) , dt

(S34) (S35)

where Ts (t) denotes soil temperature and Cs (t) the soil carbon content. The external forcing is given by a linearly increasing atmospheric temperature Ta (t) = vt,

(S36)

where v controls the rate of forcing. The nonlinear term r(Ts ) = r0 e(α(Ts −Tref )) represents temperature-dependent microbial respiration, giving rise to a positive feedback loop in which increasing temperature accelerates respiration, releasing heat and further increasing temperature. For sufficiently small forcing rates v, the system tracks a quasi-static equilibrium corresponding to stable temperature–carbon states. However, when v exceeds a critical threshold, the system fails to track this equilibrium branch and undergoes rapid thermal runaway, despite remaining locally stable at each instantaneous forcing level. This behaviour constitutes rate-induced tipping. We use parameters µ = 2.5 × 106 , λ = 5.049 × 106 , A = 3.9 × 107 , Π = 1.055, r0 = 0.01, and α = log(2.5)/10, with noise magnitude D3 = 50. Tipping and non-tipping regimes are generated by varying the forcing rate v. Tipping is identified when Ts > 30, indicating entry into the runaway regime.

(a) AMOC box-model (B-tipping)

(b) AMOC box-model (R-tipping)

Figure S6: Examples of idealized AMOC box-model trajectories. (a) B-tipping and (b) R-tipping. Each panel shows left: phase-space evolution in (ST , SN ), center: AMOC strength Q, and right: the corresponding hosing or control schedule. Note that for B-tipping, we use a more diverse set of forcing schedules than shown here, see Figures S3 and S4, which vary the initial value and slope as well as the critical time. b_amoc and r_amoc We consider a reduced Atlantic Meridional Overturning Circulation (AMOC) box model derived from the five-box formulation in [7, 56]. By treating the Southern Ocean and bottom water salinities as slow variables and enforcing salt conservation, the system reduces to a 29

two-dimensional stochastic dynamical system with state x(t) = (SN (t), ST (t))⊤ ∈ R2 , where SN and ST denote North Atlantic and tropical Atlantic salinity, respectively. The AMOC strength is diagnosed as λ [α(TS − T0 ) + β(SN − SS )] Q= , (S37) 1 + λαµ and serves as the observable for tipping detection. The salinity dynamics are piecewise-defined depending on the sign of Q. For Q ≥ 0, dSN = Q(ST − SN ) + KN (ST − SN ) − FN (H(t))S0 + σN ξN (t), dt dST VT = Q γSS + (1 − γ)SIP − ST + KS (SS − ST ) dt + KN (SN − ST ) − FT (H(t))S0 + σT ξT (t),

(S38) (S39) (S40)

while for Q < 0, dSN = |Q|(SB − SN ) + KN (ST − SN ) − FN (H(t))S0 + σN ξN (t), (S41) dt dST = |Q|(SN − ST ) + KS (SS − ST ) + KN (SN − ST ) − FT (H(t))S0 + σT ξT (t). (S42) VT dt Here, FN and FT denote freshwater fluxes induced by a hosing forcing H(t), and ξN (t), ξT (t) are independent Gaussian white noise processes with amplitudes σN = σT = 1.0. VN

For quasi-static forcing, the system exhibits b-tipping (fold bifurcation) at a critical freshwater forcing H ⋆ ≈ 0.4236, corresponding to the collapse of the overturning circulation [56]. Tipping trajectories are generated by linearly ramping H from 0.0 to 0.4736, crossing the bifurcation threshold, while non-tipping trajectories follow the same ramp but are capped at H = 0.3736, remaining on the stable branch. To isolate r-tipping, we use a transient forcing protocol H0 + ∆H sech r(t − tcrit ) , t < tcrit , H(r, t) = H0 + ∆H, t ≥ tcrit ,

(S43)

such that tipping and non-tipping trajectories differ only through the forcing rate r. For sufficiently large r, the system fails to track the quasi-static equilibrium and undergoes collapse even though the instantaneous system remains locally stable. We use r = 0.017 for tipping and r = 0.005 for non-tipping, with tcrit = 500. Tipping is identified as a rapid decline in AMOC strength Q(t), indicating a transition from the strong to the weak circulation state. SWEC-iEEG We analyze continuous intracranial electroencephalography (iEEG) recordings from the SWEC-ETHZ database [59], which comprises 116 expert-annotated seizures across 2,656 hours of multi-channel iEEG from 18 patients with pharmacoresistant epilepsy. Signals were recorded intracranially via strip, grid, and depth electrodes, band-pass filtered between 0.5 and 150 Hz, and digitized at 512 Hz. Seizure onset and end times were determined by a board-certified EEG epileptologist. Seizure onset can be interpreted as a bifurcation-induced tipping event: in the Epileptor framework [60], the transition from interictal to ictal dynamics corresponds to a saddlenode bifurcation driven by a slow permittivity variable that modulates neural excitability. The system remains near a stable resting-state equilibrium during interictal periods, and seizure onset occurs when this equilibrium is annihilated through the bifurcation, triggering a rapid transition to high-amplitude ictal oscillations. For evaluation, we extract fixed-length segments preceding each annotated seizure onset as critical trajectories, paired with multiple interictal segments, at least 1h away from any seizure, as non-critical controls. We do not feed raw iEEG into our framework, but a multi-resolution band-power representation that summarises spectral content on a common 1 s grid, see Fig. S7. For each patient we down-sample to C=16 evenly-spaced electrodes and compute log-power in K=16 contiguous frequency bands whose edges are geometrically spaced from 0.5 to 45 Hz (geomspace(0.5, 45, 17)). Each band’s power is estimated with Welch’s method on a causal-trailing window of max(2 s, 6/flo ) seconds, where flo 30

ID01 episode=1 seizure RDTC

bandpower features, vertically offset (robust centered/scaled)

1.0 mrbp_33p97_45p00hz mrbp_25p64_33p97hz mrbp_19p35_25p64hz

0.8

mrbp_14p61_19p35hz mrbp_11p03_14p61hz mrbp_8p32_11p03hz

0.6 rdtc

mrbp_6p28_8p32hz mrbp_4p74_6p28hz mrbp_3p58_4p74hz

0.4

mrbp_2p70_3p58hz mrbp_2p04_2p70hz mrbp_1p54_2p04hz mrbp_1p16_1p54hz

0.2

mrbp_0p88_1p16hz mrbp_0p66_0p88hz mrbp_0p50_0p66hz

0.0 0

100 time_s

150

200

ID01 episode=10 baseline RDTC

100 time_s

150

200

bandpower features, vertically offset (robust centered/scaled)

mrbp_33p97_45p00hz mrbp_25p64_33p97hz mrbp_19p35_25p64hz mrbp_14p61_19p35hz mrbp_11p03_14p61hz mrbp_8p32_11p03hz mrbp_6p28_8p32hz mrbp_4p74_6p28hz mrbp_3p58_4p74hz mrbp_2p70_3p58hz mrbp_2p04_2p70hz mrbp_1p54_2p04hz mrbp_1p16_1p54hz mrbp_0p88_1p16hz mrbp_0p66_0p88hz mrbp_0p50_0p66hz

1.04

1.02

rdtc

1.00

0.98

0.96

100

150 time_s

200

250

300

100

150 time_s

200

250

300

Figure S7: Representative SWEC-iEEG traces around seizure onset and baseline. Right-hand side plot shows vertically stacked sub-rows of the K=16 multi-resolution log-power features, one band per row, with seizure onset at 160 s. The expert-annotated seizure time is used for the surrogate RDTC definition (top row). is the band’s lower edge (so low-frequency bands integrate longer windows for adequate spectral resolution); we cap nperseg at 1024 samples and zero-pad the FFT to 2048 to fix the frequency-grid spacing across bands. Band power is then obtained by trapezoidal integration of the PSD over the band edges, log10 -transformed, and aggregated across the 16 channels via median. The resulting K=16-dimensional time series, sampled at 1 Hz, is the input to all methods on SWEC-iEEG. This dataset provides a challenging real-world benchmark with high dimensionality (36–100 electrodes per patient or representative bandpowers), substantial inter-patient variability, and clinically defined transition times. TAC We consider a stochastic thermoacoustic system exhibiting a subcritical Hopf bifurcation, following the experimental and modelling study of Bonciolini et al. [44]. The system describes the evolution of acoustic pressure amplitude A(t) in a combustion chamber, where thermoacoustic feedback can lead to large-amplitude oscillations. Near the bifurcation, the dynamics can be approximated by a stochastic normal form dA = µ(t)A − αA3 + βA5 + σξ(t), dt

(S44)

where A(t) denotes the oscillation amplitude, µ(t) is a time-dependent control parameter, α, β > 0 determine the nonlinear saturation, and ξ(t) is Gaussian white noise. The subcritical nature of the bifurcation implies bistability between a low-amplitude (stable) state and a high-amplitude limit cycle. The control parameter µ(t) is ramped in time, corresponding to changes in operating conditions of the combustor. In the quasi-static case, the system undergoes a bifurcation-induced transition when the stable equilibrium disappears. However, for finite ramping rates, the system exhibits a rate-dependent delay: the transition to large-amplitude oscillations occurs beyond the quasi-static bifurcation point, 31

Critical vs non-critical Tac examples critical: episode 1 ramp_a_cycle_00 (critical, target=hopf) amplitude SN Hopf tip

non-critical: episode 9 stationary_1_window_006 (non-critical, baseline) 1.0

amplitude

1.04

0.6

60 0.4 40

100 1.00 0.98

50 0.2

0.96 0.0

0 2

6 time (s)

amplitude SN Hopf tip

RDtC

non-critical: episode 10 stationary_1_window_023 (non-critical, baseline) 200

1.0

amplitude

RDtC

1.04 0.8 0.6

150 0.4

RDtC

200

150

100

amplitude

250

time (s)

critical: episode 2 ramp_a_cycle_01 (critical, target=hopf) 300

1.00

100

0.98

0.2

1.02 RDtC

0.96 0.0

0 1

0 0

time (s)

critical: episode 3 ramp_a_cycle_02 (critical, target=hopf) amplitude SN Hopf tip

non-critical: episode 11 stationary_1_window_028 (non-critical, baseline) 1.0

amplitude

RDtC

1.04

150 0.8

1.02 0.6

150 0.4

100

RDtC

200

RDtC

amplitude

250

4 time (s)

100 1.00 0.98

0.2

0.0

RDtC

amplitude

RDtC

1.02

amplitude

RDtC

150 0.8 RDtC

amplitude

100

RDtC

amplitude

120

0.96 0

0 0

time (s)

Figure S8: Examples of critical and non-critical TAC datasets trajectories. Dashed lines indicate the varios transitions, of which we use the Hopf bifurcation as surrogate RDTC target.

with the delay increasing as the rate of change increases [44]. This phenomenon leads to dynamic hysteresis when the parameter is ramped forward and backward. For the TipPFN benchmark, the query trajectories are real experimental episodes extracted from the Bonciolini pressure recordings, rather than simulated test trajectories. We use the Mic1 pressure channel sampled at 2 kHz, extract a 100–200 Hz band-pass Hilbert amplitude envelope, and construct a balanced real-query pool of ramp-up tipping episodes and low-state stationary baseline episodes; see Fig. S8 for example trajectories. The additional three microphone channels are processed analogously and provided as features. We then provide real and simulated context data for up to four context episodes. RDTC is defined with respect to the fitted Hopf crossing; the actual forcing µ(t) is not provided to the models. This setting combines bifurcation structure, finite-rate forcing, and stochastic transition variability, making it a stringent sim-to-real test of whether TipPFN can use limited context to forecast criticality in noisy laboratory trajectories. DaphniaExt We use experimental observations of population collapse from controlled microcosm experiments reported in [43]. The system consists of replicate populations of Daphnia magna subjected to gradually deteriorating environmental conditions through sustained reductions in food availability. The resulting population dynamics can be interpreted through a minimal stochastic growth model with time-dependent parameters, dN = r(t) N + σξ(t), dt

(S45)

where N (t) denotes population size, r(t) is a time-varying growth rate reflecting environmental deterioration, and ξ(t) is Gaussian white noise [42]. 32

Experimental forcing schedule

food supply (uL)

400 300 200 100 0

100

150

deteriorating

200 day constant control

250

decline starts

300

350

400

bifurcation window

Figure S9: Experimental forcing schedule and estimated bifurcation range [43] for the DaphniaExt dataset. Observed prefix with 15 forced Ricker simulations (I10) 140

real prefix Ricker paths q10-q90

120

median simulation start day 154

population

100 80 60 40 20 0 0

100

150

200

250

300

350

day

Figure S10: Initial empirical DaphniaExt observations, extended with simulated stochastic growth model trajectories which are used for synthetic context trajectories. As environmental conditions worsen, r(t) decreases and crosses zero, corresponding to a transcritical bifurcation in which the stable positive population state loses stability and extinction becomes inevitable. The observed dynamics exhibit a transition from stationary fluctuations around a stable equilibrium to a persistent decline toward extinction. Tipping is identified as the onset of irreversible population decline.

Observed prefix with 15 unforced Ricker simulations (I10) 140

real prefix unforced Ricker paths q10-q90

120

median simulation start day 154 constant K approx 35.7

population

100 80 60 40 20 0 0

100

150

200

250

300

350

day

Figure S11: Like above, but without forcing, corresponding to a non-critical trajectory.

estimated_bif

real baseline (n=25)

real extinction (n=30)

population

RDtC 1.00

300 0.75 200 0.50 100

0.25

0.00

300

1.00 0.75

200

0.50 100 0.25 0.00

synthetic baseline (n=100) synthetic extinction (n=100) real control_extinct (n=5)

1.00 200

0.75 0.50

100

0.25 0

0.00 1.00

300 0.75 200 0.50 100

0.25

0.00

300

1.00 0.75

200

0.50 100

0.25 0.00

0 0

100

200 day

300

400

100

200 day

300

400

Figure S12: Real and synthetic DaphniaExt episodes with corresponding heuristic RDTC ramps, anchored at a random point within the estimated bifurcation range, see Fig. S9.

observed_ext

real baseline (n=25)

real extinction (n=30)

population

RDtC 1.00

300 0.75 200 0.50 100

0.25

0.00

300

1.00 0.75

200

0.50 100 0.25 0.00

synthetic baseline (n=100) synthetic extinction (n=100) real control_extinct (n=5)

1.00 200

0.75 0.50

100

0.25 0

0.00 1.00

300 0.75 200 0.50 100

0.25

0.00

300

1.00 0.75

200

0.50 100

0.25 0.00

0 0

100

200 day

300

400

100

200 day

300

Figure S13: Like above, but with the Λ = 0 point set to the actual extinction.

400

microcosm We analyze a cyanobacteria population in chemostats subjected to dilution events under gradually increasing light levels [4]. Population density is inferred from the light attenuation coefficient, computed from continuous measurements of outgoing light intensity. The full dataset comprises 7,784 observations over 28.86 days at a sampling interval of 5 minutes. The time series is divided into six segments separated by dilution events, which are external perturbations and not part of the intrinsic population dynamics. To ensure consistency, we restrict analysis to recovery phases between dilution events. Specifically, for each segment, we use the final 250 time points (approximately one day) prior to the next dilution event, yielding a total of 1,500 data points across all segments. voice We analyze experimental time series of phonation onset, corresponding to the emergence of vocal fold oscillations, which can be interpreted as a Hopf bifurcation [61, 62]. The data are obtained from a physical replica of human vocal folds (EPI model) with a multilayer body–cover structure that reproduces realistic tissue mechanics. Oscillations are induced by gradually increasing airflow, while subglottal pressure is recorded using a pressure transducer. The available time series (from [35] consists of 7,501 data points over 0.375 s with a sampling interval of 5 × 10−5 s. Tipping corresponds to the transition from a stable, non-oscillatory state to sustained oscillations as the flow rate crosses a critical threshold. We use the measured pressure signal as the observable for tipping detection. cellular_atp We analyze experimental time series of cytosolic ATP (i.e., the readily available energy pool) dynamics in living plant cells, measured via a genetically encoded FRET sensor sensitive to MgATP2− [64]. The data capture the response of leaf tissue under gradually increasing hypoxia, induced by sealing samples in a dark environment where respiration depletes available oxygen. ATP levels remain relatively stable during oxygen decline and then undergo a sudden collapse, indicating a critical transition in cellular energy state. The dataset consists of 271 observations over 3 minutes with a sampling interval of 0.011 minutes. We use the fluorescence ratio signal as the observable, with tipping identified as the abrupt drop in ATP concentration. greenhouse_earth We analyze a paleoclimate time series of calcium carbonate (CaCO3 ) associated with the end of greenhouse Earth, marking the transition from an ice-free state to the formation of polar ice caps [12]. The data exhibit increasing autocorrelation prior to the climate shift, consistent with critical slowing down. The dataset consists of 462 observations spanning 5.9 million years with a sampling interval of 0.013 million years. Tipping is identified as the transition to a glaciated climate state. blackout_frequency We analyze bus voltage frequency data preceding the Western Interconnect blackout of August 1996, measured within the Bonneville Power Administration network [5]. The time series captures system dynamics leading up to grid separation and exhibits critical fluctuations prior to the blackout, consistent with EWS of instability. The dataset consists of 11,301 observations over 565 seconds with a sampling interval of 0.05 s. Tipping is identified as the transition to systemwide failure, i.e., blackout. AMOC The RAPID-MOCHA-WBTS dataset [65] provides continuous measurements of the AMOC at 26.5°N, derived from current velocity, temperature, salinity, and pressure observations collected by a trans-basin mooring array across the Atlantic and cable measurements across the Florida Straits. These variables are used to compute oceanic volume and heat transports in both depth and density space, yielding a 12-hourly, 10-day low-pass filtered transport time series spanning April 2nd 2004 to March 27th 2024. C.2

Baselines

In addition to the classical indicators such as AR1 and variance, we provide additional details on the comprehensive set of DL-based baselines included in this work. TabPFN (v2.6) We include TabPFN [27, 28] as a state-of-the-art ICL baseline. TabPFN is a Prior-Data Fitted Network (PFN), a transformer model trained on large amounts of synthetically generated tabular data to approximate Bayesian inference via in-context learning. In our setting, we adapt TabPFN to time series by representing each trajectory segment as a tabular sample with 36

engineered features (e.g., recent observations or summary statistics), and evaluate its ability to predict RDTC. We evaluate TabPFN’s performance for an increasing number of context episodes (1, 2, 3). Zero-context datasets are skipped for TabPFN because the row-wise regressor requires observed target rows. The context information is transformed into a tabular regression problem suitable for TabPFN in the following way: Let the processed single-sample batch have row index r = 1, . . . , R, identifier columns ιr = (er , tr ), selected feature vector xr , action vector yr = (yr,1 , . . . , yr,A ), action mask mr = (mr,1 , . . . , mr,A ), context indicator cr , and valid-row indicator vr . The first action slot is transformed RDTC, yr,1 = Λ∗r . The row-wise TabPFN feature vector is XrTab = [ιr , xr , yr,2:A , mr,2:A , cr , vr ] ,

YrTab = yr,1 .

(S46)

Q = {r : vr = 1, cr = 0, mr,1 = 1}.

(S47)

The supervised fit set and prediction set are T = {r : vr = 1, mr,1 = 0, YrTab finite},

Non-finite entries in X Tab are passed as missing values. TabPFN is fit separately for each query episode on T and predicts only at Q. The model is evaluated using default parameters. It returns logits Pj ℓq,j over histogram bins and bin borders b0 , . . . , bB . With pq,j = softmax(ℓq )j , Cq,j = k=1 pq,k , and Cq,0 = 0, the inverse-CDF quantile is Q̂q (α) = bj−1 +

α − Cq,j−1 (bj − bj−1 ), max(pq,j , 10−12 )

j = min{k : Cq,k ≥ α}.

(S48)

Bury We include the deep learning approach of [22], which uses a CNN-LSTM architecture trained on simulated time series from canonical dynamical systems (Hopf, transcritial and fold). The model predicts the probability of an upcoming tipping point directly from observed trajectories in a single forward pass, without relying on handcrafted indicators. As it is trained on bifurcation normal forms, the method is primarily designed for b-tipping and has previously not been generalize to r- or n-induced transitions. Zhuge We include the deep learning approach of [24], which builds on [22] using a CNN-LSTM architecture trained on simulated time series from canonical bifurcation normal forms. The model learns a direct mapping from observed trajectories (and control parameters) to the predicted tipping point, and can handle irregularly sampled data. As it is trained on normal-form dynamics, the method (like [22]) is primarily designed for B-tipping and relies on shared features of local bifurcation structure. Huang We include the deep learning approach of [18], which develops a CNN-based classifier for predicting R-tipping under time-varying forcing and stochastic perturbations. Unlike Bury and Zhuge, this method is explicitly designed for r-tipping and predicts the probability that an observed trajectory will tip, rather than relying on CSD indicators. The original study trains and evaluates separate models on three prototypical R-tipping systems: a saddle-node system, a Bautin system, and a compost-bomb model. It further considers extrapolation to out-of-sample forcing rates, showing that models trained at one forcing rate can generalize to unseen rates in some systems, while performance is system-dependent. In our experiments, we use this method in a stronger cross-system setting: the model is trained only on the saddle-node r-tipping system and then evaluated on other R-tipping systems. This differs from the original setup, where separate system-specific training was used. We adopt this protocol because the publicly released code and data primarily provide the saddle-node training workflow/checkpoints, limiting direct replication of the full multi-system training procedure for our benchmark. C.3

Metrics

We describe the metrics used in this work. For tipping point detection, we consider a binary classification setting where each sample is labeled as either tipping (y = 1) or non-tipping (y = 0). Let ŝ : D → R denote a scoring function that assigns a confidence score to each observed trajectory, where higher values indicate a higher likelihood of an impending tipping event. A prediction ŷ = 1 is made whenever ŝ > τ for a threshold τ . 37

For each threshold τ , we define the true positive rate (TPR) and false positive rate (FPR) as TPR(τ ) =

|{i : yi = 1, ŝi > τ }| , |{i : yi = 1}|

FPR(τ ) =

|{i : yi = 0, ŝi > τ }| . |{i : yi = 0}|

(S49)

The receiver operating characteristic (ROC) curve [38, 37] is obtained by plotting TPR(τ ) against FPR(τ ) as the threshold τ varies over all possible values. The ROC curve characterizes the trade-off between correctly detecting tipping events and incorrectly classifying non-tipping trajectories. The Area Under the ROC Curve (AUROC) [39, 40] is defined as Z 1 AUROC = TPR FPR−1 (u) du,

(S50)

which corresponds to the probability that a randomly chosen tipping trajectory is assigned a higher score than a randomly chosen non-tipping trajectory. An AUROC of 0.5 indicates random performance, while 1.0 corresponds to perfect discrimination.

C.3.1

Scores

Table S9 gives an overview of computed scores for each model and baseline. Table S9: Score definitions across all methods. TipPFN and TabPFN expose six heads per scope: two distributional summaries of the predicted RDtC quantile distribution (1-mean, 1-median) and four CDF probabilities at thresholds τ ∈ {0.05, 0.1, 0.2, 0.3}. The model emits RDtC in tanhbounded coordinates; raw thresholds τ in physical RDtC units are mapped by τ̃ = tanh(τ /0.2). Both scopes use the last query position in their respective scope (forecast horizon aggregation if a horizon row is present). Baseline score names follow <method>:<short_mode>:<scope>:<head> where <scope> is wN for a fixed query window of length N or whist for the whole available history; <short_mode> ∈ {pad, backfill, resample} controls how short signals are extended to the model’s input length. EWS does not consume short_mode; Zhuge additionally carries a <control_source> ∈ {param, rdtc} slot. Higher score ⇒ closer to tipping for every entry except where noted. Score name

Definition

TipPFN, TabPFN — nowcast scope (last nowcast query position) <m>:nc:1-median 1 − median(Qnc ) <m>:nc:1-mean

1 − E[Qnc ]

<m>:nc:P(rdtc<τ ) Pr(Qnc < τ̃ ) TipPFN, TabPFN — forecast scope (forecast query position(s)) <m>:fc:1-median 1 − median(Qfc ) <m>:fc:1-mean <m>:fc:P(rdtc<τ )

1 − E[Qfc ] Pr(Qfc < τ̃ )

EWS[12, 9] — rolling-window indicator with Kendall-τ trend test dews:wN :var_tau τ Varw (xt ) dews:wN :ar1_tau dews:wN :acf_tau dews:wN :skw_tau dews:wN :lambd_tau

τ AR(1)w (xt) τ ACFw (xt ) τ Skeww (x t) τ λ̂w (xt )

Bury [22] — CNN classifier over {fold, hopf, transcritical, null} bury:<sm>:wN :prob_tip Pr(tip) = 1 − Pr(null) bury:<sm>:wN :prob_fold bury:<sm>:wN :prob_hopf bury:<sm>:wN :prob_transcritical

Pr(fold) Pr(hopf) Pr(transcritical)

Huang [18] — CNN classifier (tip vs. no-tip) huang:<sm>:wN :prob_tip

Pr(tip)

Notes Qnc : predicted RDtC quantile distribution at the last nowcast query Trapezoidal integration of the quantile function over p ∈ [0, 1] τ̃ = tanh(τ /0.2), τ ∈ {0.05, 0.1, 0.2, 0.3} Last forecast position; multi-horizon runs aggregate per-horizon scores As above; same scope rules as fc:1-median τ̃ , τ as in the nowcast row Kendall-τ of rolling variance over window of length N Kendall-τ of rolling lag-1 autocorrelation Kendall-τ of rolling autocorrelation function (lag 1, window) Kendall-τ of rolling skewness Kendall-τ of rolling local return rate λ̂ Sum of fold + hopf + transcritical class probabilities Per-class softmax output Per-class softmax output Per-class softmax output Single-output sigmoid probability of an upcoming tipping event

Zhuge [24] — neural regressor for the tipping bifurcation parameter zhuge:<cs>:<sm>:wN :signed_tip_margin ∆± tip = −p̂rem /|p∆ |

Normalized signed margin to the bifurcation parameter; identical for <cs> ∈ {param, rdtc}

C.4

Further Results

In this section, we provide additional information to the conducted experiments and their results. AUROC over lead time by system/data row Column-wise best: one fixed table setting per method column, shared across all rows

1.0

B-Fold

B-Hopf

B-Trans.

B-Harv.

B-RM TC

B-RM Hopf

B-SEIRx

B-AMOC

R-Bautin

R-SN

R-Compost

R-AMOC

AUROC

0.8 0.6 0.4 0.2 0.0

1.0

AUROC

0.8 0.6 0.4 0.2 0.0

1.0

AUROC

0.8 0.6 0.4 0.2 0.0

SWEC-iEEG

1.0

TAC

AUROC

Lead time

DaphExt

TipPFN / 0c TipPFN / 1c TipPFN / 2c TabPFN / 1c TabPFN / 2c EWS / ¿var

0.8 0.6 0.4

Bury / ptip Huang / ptip

0.2 0.0

Zhuge / ¢ § tip

Lead time

Figure S14: AUROC over lead time for all systems, selecting each model’s score variant that performed best across all these systems (same as in Table 1 for ∆ > 0.

AUROC over lead time by system/data row Row-wise best: table setting re-selected separately within each system/data row B-Fold

B-Hopf

B-Trans.

B-Harv.

B-RM TC

B-RM Hopf

B-SEIRx

B-AMOC

R-Bautin

R-SN

R-Compost

R-AMOC

SWEC-iEEG

TAC

DaphExt

1.0

AUROC

0.8 0.6 0.4 0.2 0.0 1.0

AUROC

0.8 0.6 0.4 0.2 0.0 1.0

AUROC

0.8 0.6 0.4 0.2 0.0

1.0

AUROC

Lead time TipPFN / 0c TipPFN / 1c TipPFN / 2c TabPFN / 1c TabPFN / 2c EWS / best Bury / best Huang / best Zhuge / best

0.8 0.6 0.4 0.2 0.0

Lead time

Figure S15: AUROC over lead time for all systems, selecting each model’s score variant that performed best on that particular system (corresponding to Table S10) for ∆ > 0.

Table S10: Row-wise best pre-tip AUROC by system and dataset. Balanced AUROC for criticalvs-all tipping detection, averaged over positive lead times (∆ > 0). Each cell re-selects the best score/configuration within that system/data row and method column. The small line below each AUROC is the picked configuration: w<N> for query-window length, res = resample, bf = backfill, pad = pad; for TipPFN/TabPFN it also lists the selected score head (e.g. fc:1-median or fc:P(rdtc<0.05)). Feature budget is fixed at 16 within the TipPFN/TabPFN columns; baseline score heads match the Table S13. For the TAC dataset, higher performance can be achieved by selective context composition, see S22. System/data

TipPFN 0c

TipPFN 1c

TipPFN 2c

TabPFN 1c

TabPFN 2c

EWS best

Bury [22] best

Huang [18] best

Zhuge [24] best

Canonical B-Fold B-Hopf B-Trans.

.935

.929

.955

.584

.786

.655

.697

.577

.522

w96 / fc:P(rdtc<0.2)

w128 / fc:1-median

w64 / fc:P(rdtc<0.3)

w64 / fc:1-median

w128

res / w128

pad / w128

res / w128

.893

.935

.952

.640

.820

.755

.842

.497

.665

w128 / fc:P(rdtc<0.1)

w128 / fc:1-median

w128 / fc:P(rdtc<0.05)

w64 / fc:1-median

w128

res / w128

res / w96

res / w128

.874

.925

.950

.601

.824

.650

.722

.490

.458

w128 / fc:P(rdtc<0.3)

w128 / fc:1-median

w128 / fc:P(rdtc<0.3)

w128 / fc:1-median

w128

res / w128

bf / w64

res / w128

Semi-real B-Harv. B-RM TC B-RM Hopf B-SEIRx B-AMOC

.970

.971

.980

.691

.847

.991

.967

.873

.734

w128 / fc:1-median

w128 / fc:1-mean

w96 / fc:1-median

w96 / fc:P(rdtc<0.3)

w128

res / w128

res / w64

res / w128

.901

.952

.963

.569

.741

.733

.784

.604

.570

w96 / fc:P(rdtc<0.2)

w128 / fc:P(rdtc<0.05)

w128 / fc:P(rdtc<0.1)

w128 / fc:P(rdtc<0.3)

w64 / fc:1-median

w128

res / w128

res / w64

res / w128

.770

.845

.912

.529

.689

.627

.638

.484

.422

w128 / fc:P(rdtc<0.3)

w128 / fc:P(rdtc<0.1)

w128 / fc:P(rdtc<0.05)

w128 / fc:P(rdtc<0.3)

w64 / fc:1-median

w128

res / w128

res / w64

res / w128

.723

.831

.907

.583

.699

.536

.500

.507

.426

w128 / fc:P(rdtc<0.3)

w128 / fc:1-mean

w64 / fc:P(rdtc<0.05)

w96

pad / w128

res / w64

bf / w96

.906

.969

.972

.632

.829

.935

.948

.509

.679

w96 / fc:P(rdtc<0.3)

w128 / fc:1-median

w128 / fc:P(rdtc<0.05)

w64 / fc:P(rdtc<0.3)

w128

res / w128

res / w64

res / w128

.580

.643

.776

.468

.501

.451

.541

.538

.283

w64 / fc:P(rdtc<0.05)

res / w64

bf / w96

R-Bautin

w64 / fc:1-median

w96 / fc:P(rdtc<0.3)

w64

R-SN

.754

.864

.931

.501

.592

.440

.520

.959

.242

w64 / fc:1-median

w96 / fc:P(rdtc<0.3)

w128 / fc:P(rdtc<0.3)

w64 / fc:P(rdtc<0.3)

w64

bf / w96

pad / w128

bf / w96

.927

.763

.793

.500

.505

.481

.727

.737

.506

w128 / fc:P(rdtc<0.05)

w128 / fc:1-median

w64 / fc:P(rdtc<0.3)

w64

res / w128

.584

.798

.924

.616

.699

.270

.519

.982

.350

w64 / fc:P(rdtc<0.3)

w128 / fc:P(rdtc<0.3)

w128 / fc:P(rdtc<0.05)

w64 / fc:P(rdtc<0.3)

w64

bf / w96

bf / w128

bf / w96

R-Compost R-AMOC

Real-world & Sim-to-Real .679 SWEC-iEEG w128 / fc:P(rdtc<0.2) .536 TAC w96 / fc:P(rdtc<0.05) .556 DaphniaExt w16 / fc:1-median All datasets

.773

w64 / fc:P(rdtc<0.1)

.624

.925

.486

.687

.539

.696

.393

w128 / fc:P(rdtc<0.3)

w64 / fc:1-median

w128 / fc:P(rdtc<0.3)

w96

res / w64

res / w128

.526

.554

.533

.476

.568

.479

.522

w128 / fc:P(rdtc<0.3)

w64 / fc:P(rdtc<0.1)

w128 / fc:P(rdtc<0.1)

w128

res / w64

– –

.736

.833

.514

.672

.580

.651

.404

w16 / fc:P(rdtc<0.05)

w16 / fc:1-median

w16 / fc:P(rdtc<0.1)

w8 / fc:1-mean

w12

pad / w12

res / w8

–

.821

.889

.563

.691

.614

.682

.605

.488

Table S11: Source-trajectory support for the AUROC evaluation. For each system or empirical benchmark row, we report the number of unique source samples and query episodes represented in the result tables, split by real versus synthetic source origin and by critical versus non-critical labels. Counts are deduplicated over lead-time offsets and model-configuration sweeps. Canonical and semi-real systems are synthetic simulation benchmarks by construction. nneg split

samples

episodes

real npos

real nneg

synth. npos

synth. nneg

B-Fold

164

B-Hopf

184

B-Trans.

164

B-Harv.

B-RM TC

124

B-RM Hopf

132

B-SEIRx

B-AMOC R-Bautin R-SN R-Compost R-AMOC

27 34 35 31 28

108 136 140 124 112

0 0 0 0 0

54 66 71 63 56

54 70 69 61 56

9 equilibrium, 13 flat 7 receding, 11 approaching 15 equilibrium, 18 flat 16 receding, 13 approaching 17 equilibrium, 16 flat 15 receding, 19 approaching 11 equilibrium, 7 flat 5 receding, 10 approaching 14 equilibrium, 12 flat 13 receding, 15 approaching 38 equilibrium, 32 flat 41 equilibrium, 28 flat 30 equilibrium, 31 flat 25 equilibrium, 31 flat

Others SWEC-iEEG TAC DaphniaExt

18 2 2

810 320 110

116 80 30

694 80 25

0 80 30

0 80 25

694 baseline 80 real baseline, 80 synth. baseline 25 real baseline, 25 synth. baseline

406

2,776

226

799

877

874

799 real baseline, 105 synth. baseline, 265 synth. equilibrium 241 synth. flat, 133 synth. receding, 130 synth. approaching

System/data Canonical

19 equilibrium, 15 flat 31 receding, 19 approaching 18 equilibrium, 21 flat 29 receding, 23 approaching 28 equilibrium, 17 flat 17 receding, 20 approaching

Semi-real

Summary Total

C.4.1

SWEC-iEEG

Seizure onset prediction is a challenging and clinically relevant test case for transition forecasting, as it requires anticipating a rapid shift from interictal to ictal dynamics. In SWEC-iEEG, the model operates on multichannel iEEG bandpower trajectories. We interpret seizure onset as the critical event and evaluate prediction on preictal segments relative to interictal control segments from the same patient, drawn from periods at least one hour away from any seizure episode. TipPFN performs especially well with two context episodes (Table 1), outperforming TabPFN and remaining stronger than the classical baselines even after selecting the best bandpower time series for each baseline method. The corresponding lead-time analysis is shown in Figures S15 and S17. With neurological EEG patterns being highly patient-specific, one interesting use case for TipPFN is to evaluate individual seizure risk by including only same-patient observations into the context. Fig. S16 shows resulting patient-specific AUROC scores, with TabPFN reaching highest median score. However, the TipPFN prediction seems not to work well for some patients, specifically those with fewer recorded seizures. While the performance of Bury is lower overall, we find a narrow between-patient variance.

1.0

SWEC-iEEG per-patient AUROC (left: all n=18; right: ¸ 4 seizures, n=12)

AUROC over ¢ > 0

0.8

0.6

0.4

0.2

0.0

all patients (n=18) ¸ 4 seizures (n=12)

TipPFN 1c

TipPFN 2c

TabPFN 1c

TabPFN 2c

DEWS best

Bury best

Huang best

Figure S16: SWEC per-patient AUROC distribution. For each method column we show two boxes: the full n = 18 patient pool (lighter fill, left) and the n = 12 subset of patients with ≥ 4 reported seizures (darker fill, right). Each dot is one patient’s mean AUROC over ∆ > 0 at the setting that maximizes the patient-macro AUROC for that column. For the uni-variate baseline models, this includes a per-signal-band selection: the bands mrbp_0p50_0p66hz (delta), mrbp_3p58_4p74hz (theta), and mrbp_19p35_25p64hz (beta) are scored separately and the best band is used, matching the choice in Table 1) and Fig. S14.

Figure S17: SWEC-iEEG dataset AUROC over lead time for different number of features shown. For zero features, TipPFN and TabPFN scores are degenerate. C.4.2

Daphnia Extinction

The controlled Daphnia magna extinction experiment of Drake and Griffen [43] is a valuable real-world transfer benchmark because it combines empirical observations with experimentally controlled deterioration, a well-characterized transcritical transition, and allows to match a Rickermap simulation [42] to generate simulated context. To construct context, we consider two surrogate RDTC choices: either setting Λ = 0 at the observed extinction event, or sampling it from the bifurcation-time interval estimated in [43]. Both observed and simulated contexts can yield comparable predictive performance, suggesting that matched simulations can act as a practical surrogate when repeated real trajectories are limited. Simulated context tends to peak already with fewer episodes, see Fig. S18 and S19. The two Λ target constructions also reveal a qualitative difference between PFN-style models: TabPFN scores en par with TipPFN when it comes to distinguishing extinction; however, the arguably more difficult detection is that of the underlying bifurcation, which TipPFN discerns more clearly than other baselines, see Fig. S21.

Daphnia ROC: Observed extinction anchor, TipPFN forecast (context hue) ¢=4

¢=2

¢ = ¡2

¢=1

0.5

synthetic context TPR

0 1

context count

real context TPR

5 0.5

0 0

0.5 FPR

Figure S18: DaphniaExt dataset ROC curves for W = 16, different ∆ and increasing context episodes, using real observations as context (top row) or purely simulated context (bottom row). Surrogate Λ was condition on Λ = 0 at the extinction time of the respective population. Daphnia ROC: Estimated bifurcation anchor, TipPFN forecast (context hue) ¢=4

¢=2

¢ = ¡2

¢=1

0.5

synthetic context TPR

0 1

context count

real context TPR

5 0.5

0 0

0.5 FPR

Figure S19: Like above, but surrogate Λ was conditioned on Λ = 0 at a random time within the bifurcation-time interval [271, 316] days estimated in [43], aiming to estimate the unterlying transition.

Daphnia: AUROC by context (w=16) daphnia_extinction@observed_ext

context=0

1.0

context=1

context=2

context=5

context=10

context=20

context_source_kind=real AUROC

0.8 0.6 0.4 0.2 0.0

context_source_kind=synthetic AUROC

1.0

TipPFN nc:1-median TipPFN fc:1-median dews:w16:ar1_tau dews:w16:var_tau bury:resample:w16:prob_transcritical bury:resample:w16:prob_tip huang:resample:w16:prob_tip TabPFN nc:1-mean TabPFN fc:1-median

0.8 0.6 0.4 0.2 0.0

Lead time

Figure S20: DaphniaExt dataset AUROC curves for W = 16, different context sizes, using real observations as context (top row) or purely simulated context (bottom row). Surrogate Λ was conditioned on Λ = 0 at the extinction time of the respective population. Daphnia: AUROC by context (w=16) daphnia_extinction@estimated_bif

context=0

1.0

context=1

context=2

context=5

context=10

context=20

context_source_kind=real AUROC

0.8 0.6 0.4 0.2 0.0

context_source_kind=synthetic AUROC

1.0

TipPFN nc:1-median TipPFN fc:1-median dews:w16:ar1_tau dews:w16:var_tau bury:resample:w16:prob_transcritical bury:resample:w16:prob_tip huang:resample:w16:prob_tip TabPFN nc:1-mean TabPFN fc:1-median

0.8 0.6 0.4 0.2 0.0

Lead time

Figure S21: Like above, but with surrogate Λ on the estimated bifurcation, the systems hidden critical transition.

C.4.3

TAC

The thermo-acoustic combustor dataset provides a particularly noisy real-world benchmark built around a stochastic subcritical Hopf transition [44]. We formulate the task as prediction of the Hopf bifurcation and assign surrogate Λ labels by anchoring a linear ramp at the known critical point. In the basic setting, the Hopf-specific Bury model performs very well, consistent with its close match to the target mechanism. TipPFN performance is higher for an increased context size of 4; it additionally increases when composing the context with exclusively critical context episodes, presumably because they are the more informative and distinguishing signal. Fig. S22 shows the effect on AUROC when filling the context with only non-critical, only critical, or random episodes, causing TipPFN to outperform Bury. TAC context-composition regimes features=3, ctx=[2,3,4], composition=random; baselines pooled globally all non-critical

mixed

all critical

1.0

0.8

0.6

0.4

0.2

AUROC

1.0

0.0

40 20 Lead time

TipPFN P(rdtc < 0:3) (¹ > 0 =0.60) TabPFN P(rdtc < 0:3) (¹ > 0 =0.50) EWS ¿var (¹ > 0 =0.57) Bury ptip (¹ > 0 =0.48) Huang ptip (¹ > 0 =0.51)

0.0

40 20 Lead time

TipPFN P(rdtc < 0:3) (¹ > 0 =0.61) TabPFN P(rdtc < 0:3) (¹ > 0 =0.55) EWS ¿var (¹ > 0 =0.57) Bury ptip (¹ > 0 =0.48) Huang ptip (¹ > 0 =0.51)

0.0

40 20 Lead time

TipPFN P(rdtc < 0:3) (¹ > 0 =0.65) TabPFN P(rdtc < 0:3) (¹ > 0 =0.50) EWS ¿var (¹ > 0 =0.57) Bury ptip (¹ > 0 =0.48) Huang ptip (¹ > 0 =0.51)

Figure S22: TAC dataset AUROC scores over lead time for contexts of different composition.

C.4.4

Zero-shot datasets

In addition to Fig. 5, we study three further zero-shot time series, depicted in Fig. S23. While we find early warning (∆ > 0) for the mitochondria and voice dataset, the RDTC prediction for the paleoclimate (greenhouse_earth) dataset only crossed the threshold after the critical transition. a Paleoclimate

b Mitochondria

c Voice 1.0

¤¤

signal

0.85 50

0.80 0.75

0.5

0 1.0

1.0

0.5

count

0.0 50 0

0.0 100/100 crossed median lead time ¢ = -0.4 Myr

5.0

2.5

50 0

0.0

Lead time ¢ [Myr]

0.0 100/100 crossed median lead time ¢ = 0.244 s

100/100 crossed median lead time ¢ = 0.0231 s

0.2

Lead time ¢ [s] observed

bifurcation

nowcast q50 §80% CI

¤ = 0:05 crossing

0.1

0.0

Lead time ¢ [s] median crossing

first crossing

Figure S23: Zero-shot TipPFN RDTC nowcasts on three additional uni-variate real-world time series (see Fig. 5). Note the crossing threshold Λ∗ = tanh(5Λ) ≈ 0.245.

C.5

AUROC Uncertainty

Table S12: Per-system standard errors for the macro AUROC over ∆ > 0 for the values reported P (c) (c) in Tables 1 and S10. Total samples is c (npos + nneg ) across all contributing (∆, system) points P (c) (c) pooled nmin is c min(npos , nneg ), the effective balanced-sample count entering equation (S51). SEHM is reported at A = 0.8; values at A = 0.5 are ∼ 25% larger. SEBS is the propagated balancedsubsample uncertainty from equation (S52); for TAC it is effectively zero because npos =nneg at every contributing point and no subsampling occurs. System / dataset total samples pooled nmin SEHM at A=0.8 SEBS Canonical B-Fold B-Hopf B-Trans.

911 1,019 910

407 465 418

0.020 0.019 0.020

0.0009 0.0008 0.0005

Semi-real B-Harv. B-RM TC B-RM Hopf B-SEIRx B-AMOC R-Bautin R-SN R-Compost R-AMOC

434 672 737 370 597 738 768 682 608

194 300 335 164 273 318 346 308 272

0.029 0.023 0.022 0.031 0.024 0.022 0.022 0.023 0.024

0.0008 0.0021 0.0013 0.0022 0.0019 0.0042 0.0017 0.0023 0.0029

Others SWEC-iEEG TAC DaphniaExt

14,580 21,600 4,400

2,088 10,800 2,000

0.009 0.004 0.009

0.0006 ∼0 0.0006

Summary All systems

49,124

18,688

0.003

0.0005

For each AUROC value reported, we compute a balanced AUROC by drawing K = 10 classbalanced subsamples without replacement from the available critical and non-critical query windows, evaluating sklearn.metrics.roc_auc_score on each, and reporting the mean. Each subsample contains nmin = min(npos , nneg ) examples per class. We attach a sample-size standard error using the Hanley–McNeil conditional-Gaussian approximation, s A (1 − A) SEHM = , (S51) npooled min P (c) (c) where npooled min = c min(npos , nneg ) sums per-point minority-class counts across all lead times ∆ > 0 contributing to the reported macro AUROC. Equation (S51) gives a slightly conservative bound because it ignores correlation between adjacent lead times. AUROC differences below 2 SEHM should be regarded as not statistically resolved. We additionally report the propagated standard error of the macro AUROC from the balanced(c) subsampling procedure. For each contributing (∆, system) point we have an empirical std σBS across the K=10 subsamples; pooling these across the Kmacro points contributing to a system’s macro AUROC gives v u Kmacro X (c) 2 1 u t1 SEBS = σ . (S52) Kmacro K c=1 BS This quantity captures sensitivity of the AUROC estimator to which nmin examples are chosen for a given imbalance ratio at each point. It does not capture sampling variation in the underlying episode population, and is therefore consistently 5–30× smaller than SEHM in our data. We report it as a 50

separate diagnostic rather than combining it with SEHM , since the two estimate related but distinct variances and the second-order contribution would not move the displayed precision. Per-system pooled nmin , SEHM , SEBS , and the underlying total sample count are shown in Table S12. Selection bias from best-of-K picking. For columns whose value is the best across a family of configurations (e.g. any column that selects across window size, feature-budget, or score-head combinations within a method), the reported AUROC is the maximum over K correlated noisy estimates and is therefore positively biased relative to the true expected AUROC at the chosen setting. Under the null hypothesis that all candidate settings have equal expected AUROC, the expected maximum exceeds the per-cell mean by approximately σ · Φ−1 (1 − 1/K), where σ is the per-cell sampling std and Φ−1 is the inverse normal CDF. For our typical K ∈ [10, 150] candidate cells per (method, system) and σ ≈ SEHM , the per-cell inflation is 0.02–0.06 in the iid-null limit; correlation across candidates pulls this back to a typical 0.01–0.03. An unbiased estimate would require holding out a separate validation split for cell selection, which the current evaluation pipeline does not expose. Differences between methods at the per-system level should therefore be interpreted as upper bounds on the true gap rather than unbiased point estimates. C.5.1

Score Analysis

The CNN baselines (Bury, Huang, Zhuge) require fixed-length input windows; signals shorter than the network’s expected length are extended via one of three short_mode schemes (pad, backfill, resample). Rather than pinning a single short-mode globally, each baseline column in the headline AUROC table independently picks its column-best short_mode, and the row-wise variant (Table S10) re-picks per system; the picked token is reported in every baseline cell. Across all baseline cells of the row-wise table, resample dominates (∼71%), with backfill (∼19%) and pad (∼10%) chosen on a minority of systems where the signal-to-padding scaling differs. Because the picked token is already disclosed at cell granularity and no global pin is asserted, we do not include a dedicated short-mode sensitivity SI table; the row-wise picks are self-documenting. A handful of baseline scores are anti-correlated with criticality on specific datasets, producing AUROCs significantly below chance (< 0.5). The clearest case is Zhuge’s signed_tip_margin on SWEC-iEEG, where the bifurcation-parameter regression sign convention is inverted relative to the labelling: critical episodes receive lower scores than non-critical ones, collapsing AUROC to ≈ 0. Huang’s prob_tip shows a milder version of the same effect on high-frequency SWEC bandpower bands (14–25 Hz, AUROC ≈ 0.18–0.20), and saturates at 1.0 on the lowest band (0.50– 0.66 Hz), where AUROC degenerates to 0.5. While excluding degenerate scores, we report all baseline AUROCs as-is without auto-flipping anti-correlated scores – flipping at evaluation time would amount to a per-cell sign-tuning. The per-band signal selection used on SWEC partially masks the saturation/inversion artefacts, since the best-AUROC band is picked per cell and pathological bands are dropped, but for methods whose inversion is uniform across bands (Zhuge on SWEC) no band selection can rescue the score.

Table S13: Selected configuration for each column of the headline AUROC table (Table 1). For TipPFN and TabPFN the configuration spans the score head, context size, query-window length, and feature budget; for the baseline methods only the parameters consumed by that method are listed (e.g. EWS reads window and indicator; Bury, Huang and Zhuge additionally read short_mode). Each baseline column independently picks its best query-window length (column-best over {64, 96, 128}) and, for SWEC where the score is computed across multiple bandpower bands, its best signal band. Daphnia uses an incomparable window range {8, 12, 16}, so every column reports the per-cell best-of{8, 12, 16} for daphnia (matching the row-wise convention). The row-wise best variant (Table S10) re-selects window, short_mode, and (on SWEC) signal band separately per system. Column

Score name

Configuration

TipPFN / 0c TipPFN / 1c TipPFN / 2c TabPFN / 1c TabPFN / 2c EWS / best Bury / best Huang / best Zhuge / best

tippfn:fc:1-median tippfn:fc:1-median tippfn:fc:1-median tabpfn:fc:1-median tabpfn:fc:1-median dews:w64:var_tau bury:resample:w128:prob_tip huang:resample:w64:prob_tip zhuge:rdtc:resample:w128:signed_tip_margin

head=fc:1-median, ctx=0, w=128, features=16 head=fc:1-median, ctx=1, w=128, features=16 head=fc:1-median, ctx=2, w=128, features=16 head=fc:1-median, ctx=1, w=64, features=16 head=fc:1-median, ctx=2, w=64, features=16 w64, indicator=var_tau short_mode=resample, w128, head=prob_tip short_mode=resample, w64, head=prob_tip control=rdtc, short_mode=resample, w128, head=signed_tip_margin

Table S14: Head sensitivity of the All-datasets summary row. Mean AUROC over ∆ > 0 across all 15 system/data rows, broken down by TipPFN/TabPFN context column (column groups), score scope+head (rows), and query-window length (sub-columns) at features=16. The headline col-wise table (Table 1) pins TipPFN/TabPFN to fc:1-median for every p column. The best cell per column 2 2 group is shown in bold with an asterisk (·∗ ); other cells within σbest + σcell of the best (crosssystem SE) are also bold (no asterisk), indicating statistically indistinguishable choices. TipPFN / 0c

TipPFN / 1c

TipPFN / 2c

TabPFN / 1c

TabPFN / 2c

Score head

w64

w96

w128

w64

w96

w128

w64

w96

w128

w64

w96

w128

w64

w96

w128

fc:1-median fc:1-mean fc:P(rdtc<0.05) fc:P(rdtc<0.1) fc:P(rdtc<0.2) fc:P(rdtc<0.3) nc:1-median nc:1-mean nc:P(rdtc<0.05) nc:P(rdtc<0.1) nc:P(rdtc<0.2) nc:P(rdtc<0.3)

.684 .698 .705∗ .705 .703 .694 .637 .652 .500 .500 .509 .618

.688 .692 .692 .696 .698 .701 .702 .701 .673 .682 .698 .700

.688 .700 .689 .690 .691 .690 .685 .697 .677 .684 .683 .678

.704 .727 .733 .734 .735 .739 .671 .693 .527 .551 .598 .663

.759 .777 .775 .775 .772 .774 .777 .793 .580 .609 .693 .773

.776 .791 .793 .792 .790 .796 .808 .817∗ .636 .713 .784 .815

.769 .790 .786 .787 .788 .802 .786 .811 .567 .613 .700 .773

.833 .844 .836 .837 .838 .842 .849 .865 .659 .714 .795 .855

.861 .868 .864 .865 .864 .865 .876 .885∗ .717 .782 .858 .884

.524 .519 .497 .513 .532 .548 .561∗ .545 .516 .519 .533 .541

.514 .511 .493 .506 .526 .541 .541 .524 .483 .485 .516 .537

.512 .507 .485 .500 .527 .543 .525 .509 .436 .455 .521 .550

.672 .652 .606 .620 .638 .664 .743∗ .725 .669 .670 .690 .717

.667 .648 .604 .617 .635 .660 .719 .700 .623 .633 .669 .699

.665 .647 .598 .614 .640 .663 .702 .682 .564 .601 .664 .694

Record · ID 178904 · SHA-256 0d2e2fb29adbb5ad
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.