Predicting discrete-time bifurcations with deep learning

Many natural and man-made systems are prone to critical transitions—abrupt and potentially devastating changes in dynamics. Deep learning classifiers can provide an early warning signal for critical transitions by learning generic features of bifurcations from large simulated training data sets. So far, classifiers have only been trained to predict continuous-time bifurcations, ignoring rich dynamics unique to discrete-time bifurcations. Here, we train a deep learning classifier to provide an early warning signal for the five local discrete-time bifurcations of codimension-one. We test the classifier on simulation data from discrete-time models used in physiology, economics and ecology, as well as experimental data of spontaneously beating chick-heart aggregates that undergo a period-doubling bifurcation. The classifier shows higher sensitivity and specificity than commonly used early warning signals under a wide range of noise intensities and rates of approach to the bifurcation. It also predicts the correct bifurcation in most cases, with particularly high accuracy for the period-doubling, Neimark-Sacker and fold bifurcations. Deep learning as a tool for bifurcation prediction is still in its nascence and has the potential to transform the way we monitor systems for critical transitions.


Introduction
Many systems in nature and society possess critical thresholds at which the system undergoes an abrupt and significant change in dynamics [1,2].In physiology, the heart can spontaneously transition from a healthy to a dangerous rhythm [3]; in economics, financial markets can form a 'bubble' and crash into a recession [4]; and in ecology, ecosystems can collapse as a result of their interplay with human behaviour [5,6].These events, characterised by a sudden switch to a different dynamical regime, are referred to as critical transitions.
Critical transitions can be better understood with bifurcation theory [7,8], a branch of mathematics that studies how dynamical systems can undergo sudden qualitative changes as a parameter crosses a threshold (a bifurcation).Many bifurcations are accompanied by critical slowing down-a diminishing of the local stability of the system-which results in systematic changes to properties of a noisy time series, such as its variance, autocorrelation and power spectrum [9][10][11].These properties can be approximated analytically in the presence of different bifurcations [10,[12][13][14], and a corresponding observation in data can be used as an early warning signal (EWS) for the bifurcation [11].Systematic changes in variance and lag-1 autocorrelation have been observed prior to transitions in climate [15][16][17], geological [18], ecological [19,20] and cardiac [21] systems, suggesting the presence of a bifurcation.However, these EWS have limited ability in predicting the type of bifurcation [14,22] and can fail in systems with nonsmooth potentials [23] or noise-induced transitions [24].
More recently, deep learning techniques have been employed to provide EWS for bifurcations [25][26][27].This involves training a neural network to classify a time series based on the type of bifurcation it is approaching, as well as appropriate controls [25][26][27].Unlike many applications of deep learning, this approach does not require abundant data from the study system, which, in the context of critical transitions, is often unavailable.(Unfortunately we do not have data from thousands or more ecosystems or climate systems that went through bifurcation.)Instead, the approach generates a massive library of simulation data from generic models that possess each type of bifurcation.The neural network then learns generic features associated with each type of bifurcation, that can be recognised in an unseen time series of the study system.This is enabled by the existence of universal properties of bifurcations that are manifested in time series as a dynamical system gets close to a bifurcation [7,9].In our previous work, we trained a deep learning classifier to provide an EWS for continuous-time bifurcations, and found it was effective at predicting transitions for real thermoacoustic, climate and geological transitions [25].
Bifurcations can be partitioned according to whether they occur in continuous or discrete-time dynamical systems [7,8].This distinction is important, since discrete-time dynamical systems (difference equations) can display very different behaviour to their continuous-time counterparts (differential equations).As an example, consider the logistic model for population growth.When set up in continuous time (appropriate for populations with overlapping generations e.g.humans), the population grows smoothly as the reproduction rate increases.Whereas, when set up in discrete time (appropriate for populations with non-overlapping generations, e.g.insects), the population undergoes a sequence of period-doubling bifurcations to chaos [28].This involves transitioning from equilibrium dynamics to an oscillation of period 2, to an oscillation of period 4, 8, 16 etc.until the dynamics become chaotic.It is therefore important to develop EWS suitable for both continuous and discrete-time bifurcations.While indicators like variance and lag-1 autocorrelation can provide EWS for discrete-time bifurcations, the ability of a deep learning classifier at this task has not been investigated.
As well as in ecology, discrete-time bifurcations arise naturally in physiology [3], epidemiology [29], and economics [30], where events can take place on a discrete timeline.To illustrate our approach, we will use model simulations from ecology, physiology and economics, as well as experimental data from spontaneously beating chick heart aggregates [21,31].Following administration of a drug, in some aggregates the time interval between two heart beats begins to alternate i.e. there is a period doubling bifurcation (Figure 1).Such transitions can also occur for the human heart in the form of T-wave alternans, which increases a patient's risk for sudden cardiac death [32].The period-doubling bifurcation is accompanied by critical slowing down, so systematic changes in variance and lag-1 autocorrelation are expected and have been shown to provide an EWS in this system [21].The chick heart aggregates serve as a good study system to test the performance of EWS since we have multiple recordings, not all of which underwent a transition, allowing us to test for false positives.
Among discrete-time bifurcations, there are many types, each with an associated change in dynamics [7].For this study, we focus on the five local bifurcations of codimension one (Box 1).In being 'local', these bifurcations are accompanied by critical slowing down, so systematic changes and variance and autocorrelation are expected.However, not all of these bifurcations result in a critical transition [22].They can instead involve a smooth transition to an intersecting steady state (transcritical) or to oscillations with gradually increasing amplitude (supercritical Neimark-Sacker).Predicting the type of bifurcation provides information on the nature of the dynamics following the bifurcation, something variance and autocorrelation alone do not provide.
Here, we train a deep learning classifier to provide a specific EWS for bifurcations of discrete-time dynamical systems.We train the classifier using simulation data of normal form equations appended with higher-order terms and noise.We then test the classifier on simulation runs of five discrete-time models used in cardiology, ecology and economics, and assess its performance relative to variance and lag-1 autocorrelation.We vary the noise amplitude and rate of forcing in model simulations to assess robustness of the EWS.Finally, we test the classifier on experimental data of spontaneously beating chick-heart aggregates that go through a period-doubling bifurcation.A reproducible run of all analyses may be performed on Code Ocean (https://codeocean.com/capsule/2209652/tree/v1)where the code is accompanied by the necessary software environment.

Performance of classifiers on withheld test data
As in previous work [25], we train two different types of classifiers and use their ensemble average to make predictions.Classifier 1 is trained to recognise bifurcation trajectories based on middle portions of the time series, whereas Classifier 2 is trained on end portions (see Methods).In this way, Classifier 1 provides an earlier signal of a bifurcation and Classifier 2 provides a more specific signal, as more information is revealed closer to the bifurcation.To quantify performance of the classifiers, we use the F1 score, which is a combined measure of sensitivity (how many of the true positives were predicted correctly) and specificity (how many of the positive predictions were actually true positives).On the withheld test data, Classifier 1 and 2 achieved an F1 score of 0.66 and 0.85, respectively.On the simpler, binary classification problem of predicting whether or not there will be any bifurcation, the classifiers achieved an F1 score of 0.79 and 0.97, respectively.Classifier 2 has a higher performance as it has the easier task of classifying data closer to the bifurcation where fluctuations are more pronounced.Performance on individual bifurcation classes is shown by confusion matrices (Sup.Fig. 1).The period-doubling, Neimark-Sacker and fold bifurcations are correctly classified with high sensitivity and specificity.On the other hand, the transcritical and pitchfork bifurcations are often mistaken for one another, likely due to having very similar normal forms (identical linear terms).Despite this, Classifier 2 can distinguish them at better than random, suggesting it is capable of recognising the different higher order terms in the data.From here onward, we report results using the ensemble prediction of the two classifiers, referred to collectively as the deep learning classifier.

Performance of EWS on theoretical models
We monitor variance, lag-1 autocorrelation and the deep learning classifier as progressively more of the time series is revealed.Variance and lag-1 autocorrelation are considered to provide an EWS if they display a strong trend, which we quantify using the Kendall tau statistic.For each of the five theoretical models (Fig. 2(a-e)), we observe an increasing trend in variance (Fig. 2f-j), and an increasing or decreasing trend in lag-1 autocorrelation (Fig. 2k-o).The direction of the trend in lag-1 autocorrelation prior to a bifurcation depends on the frequency of oscillations (θ) at the bifurcationequivalently the angle of the dominant eigenvalue from the positive x-axis (Box 1).For θ ∈ [0, π/2) lag-1 autocorrelation increases, whereas for θ ∈ (π/2, π] it decreases-insights that can be obtained from analytical expressions of the autocorrelation function [10,14].The period-doubling bifurcation is characterised by θ = π, and the Neimark-Sacker bifurcation shown here has θ ≈ π/4.The trends in variance and lag-1 autocorrelation therefore behave as expected and can be used as an EWS. The deep learning classifier assigns a probability to each of the six possible outcomes (null, perioddoubling, fold, Neimark-Sacker, transcritical and pitchfork).It is considered to provide an EWS when there is a heightening in the sum of the bifurcation probabilities (blue line, Fig. 2p-t).The

Box 1. Local discrete-time bifurcations of codimension one.
A discrete-time dynamical system has the form where ⃗ x t is a vector representing the state of the system and f is a function that maps to the next state.The system is stable about a steady state ⃗ x * if all eigenvalues (λ) of the Jacobian matrix (gradient of f at ⃗ x * ) have a magnitude less than one i.e. they lie within the unit circle on the complex plane (depicted in gray).A local bifurcation occurs when one or a pair of these eigenvalues (the dominant eigenvalue(s)) crosses the unit circle.The type of bifurcation depends on where the unit circle is crossed.A bifurcation is said to be of codimension one if it occurs upon changing only a single parameter of the system.There are five types of local, discrete-time bifurcations of codimension one: (a) period-doubling (b) Neimark-Sacker (c) fold (d) transcritical (e) pitchfork Bifurcation diagrams are shown in the panels below [7].The x-axis represents the bifurcation parameter and the y-axis a state variable.Solid lines show stable states or limit cycles, and dashed lines show unstable states.Arrows indicate direction of motion at different points in state space.The white points on the unit circle show the value of the dominant eigenvalue(s) at the bifurcation.The Neimark-Sacker bifurcation has a pair of complex conjugate eigenvalues, the angle (θ) being the frequency of oscillations at the bifurcation.type of bifurcation predicted is then taken as the highest individual bifurcation probability.For each simulation, the classifier is becomes more confident of an approaching bifurcation as time goes on, and its assigned bifurcation probability for the true bifurcation increases.The period-doubling, Neimark-Sacker and fold bifurcations are identified with high confidence well before the transition.The transcritical and pitchfork bifurcations are assigned roughly equal probability on their respective time series, suggesting they are difficult to tell apart-an observation consistent with the classifier's performance on its within-sample test data.
To obtain a measure of performance for the EWS, we need to test their predictions on both 'forced' time series (where a bifurcation is approached) and 'null' time series (where no bifurcation is approached).For the theoretical models, we generate null time series by keeping the bifurcation parameter fixed.We also test the robustness of the EWS to the rate of forcing and the noise amplitude of the simulations-two factors that have been shown to influence the performance of variance and lag-1 autocorrelation as an EWS [33,34].To this end, we simulate 100 forced and null time series at five different values of noise intensity and five different values of rate of forcing, resulting in a total of 5,000 time series for each theoretical model.Sample trajectories illustrating the different noise amplitudes and rates of forcing are shown in Sup.Fig. 3.We compute the probabilities assigned by the classifier and the Kendall tau values for variance and lag-1 autocorrelation at 80% of the way through the pretransition time series, and use these values as discrimination thresholds to construct ROC curves (Fig. 4a-e).Using the AUC score (area under the ROC curve) as a measure of performance, we find that the classifier outperforms variance and lag-1 autocorrelation for each theoretical model.When evaluated for each combination of noise amplitude and rate of forcing separately, the classifier has the highest AUC score in 100% of cases for the Neimark-Sacker, fold, and pitchfork models, 84% of cases for the period-doubling model, and 80% of cases for the transcritical model (Sup.Fig. 4).Similar to variance and lag-1 autocorrelation, the performance of the classifier is lower at higher rates of forcing.Noise amplitude affects performance differently depending on the model.In terms of predicting the correct bifurcation, the classifier typically performs better at slower rates of forcing (Sup.Fig. 5) and was able to classify the period-doubling and Neimark-Sacker bifurcations to high accuracy at all noise amplitudes and rates of forcing considered.

Performance of EWS on chick heart data
In the chick heart data, we mostly observe an increasing trend in variance and a decreasing trend in lag-1 autocorrelation prior to the period-doubling bifurcation, as previously reported [21].This is consistent with analytical approximations for variance and lag-1 autocorrelation prior to a perioddoubling bifurcation in a noisy dynamical system [10,14].The 46 records and their EWS are shown in Sup.Fig. 6-9, and a sample of five period-doubling records are shown in Fig. 3.The classifier correctly predicts a period-doubling bifurcation in 16 of the 23 period-doubling records.In other cases, it incorrectly predicts a Neimark-Sacker bifurcation (e.g.Fig. 3e).This seems to be linked to an early increase in lag-1 autocorrelation, perhaps caused by a non-monotonic approach to the period-doubling bifurcation.For predictions made at 60-100% of the way through the chick heart data, the classifier obtains the highest AUC score (Fig. 4f), a slight improvement on variance, with the advantage of also providing the bifurcation type.

Discussion
Many systems that evolve on a discrete timeline can undergo a sudden change in dynamics via a discrete-time bifurcation.We have found that a deep learning classifier is an effective tool for predicting discrete-time bifurcations in systems with a range of noise levels and rates of approach to the bifurcation.The classifier provides better sensitivity and specificity than variance and lag-1 autocorrelation-two commonly used EWS for bifurcations.Moreover, the classifier provides early indication of the type of bifurcation-an important piece of information given the qualitatively different dynamics associated with each bifurcation.A reliable early warning signal that specifies the type of bifurcation will help us prevent harmful bifurcations (e.g.dangerous heart rhythms [3]) and promote favourable transitions (e.g.ecosystem recovery [39]).
It may be possible to design a deep learning classifier that achieves a higher performance on our  [37].(e) Lorenz model going through a pitchfork bifurcation [38].(f) Chick heart aggregates going through a period-doubling bifurcation [31].Predictions are made 80% of the way through the pretransition data for the models, and 60-100% of the way for the experimental data.PD: period-doubling.NS: Neimark-Sacker.TC: transcritical.PF: pitchfork.
test data.First, there are many neural network architectures that could be investigated.For example, transformers, which are the current state-of-the-art for language models like GPT [40], may also be useful for time series classification [41].Second, the hyperparameters of the classifier could be systematically tuned to optimise performance.Third, there may be benefit to reframing bifurcation prediction as a hierarchical classification problem [42].One classifier could address the binary problem of flagging an approaching bifurcation, and a second classifier could address the multi-class problem of classifying the type of bifurcation given that a bifurcation is approaching.(The same way one might want to distinguish images of dogs and cats before attempting to classify dog breeds.)For a classifier to be effective, it must be trained on sufficiently diverse training data.As such, the method by which training data is obtained needs careful consideration.Our previous work on continuous-time bifurcations obtained training data from randomly generated dynamical systems with polynomial terms [25] and labelled the data using the bifurcation continuation software AUTO.This approach is appealing as it imposes relatively few restrictions on the models that are generated, and may include features associated with higher-order terms.Here, we opted for a more restricted approach that uses normal form models to generate the training data.This method has the advantage of being faster computationally, since the bifurcation of the model is known a priori.This study shows that even with this more restricted training data, a classifier can generalise to detecting bifurcations in more complex model and empirical systems.
We trained a classifier to provide EWS for a subset of bifurcations, namely local, codimensionone discrete-time bifurcations.While these bifurcations are present in many systems of interest, the real world presents many other classes of bifurcation in both continuous and discrete-time, including global bifurcations (e.g.homoclinic and heteroclinic), codimension-two bifurcations (e.g.cusp and Bogdanov-Takens), and bifurcations of attractors.For systems on attractors that explore a large portion of their phase space, empirical dynamical modelling [43], reservoir computing [44,45] and deep neural networks [46] can be used to make forecasts that may help predict critical transitions.In cases where spatial information is available, concepts from statistical physics may be useful [47], particularly in combination with deep learning [27].
One limitation of the present classifier is that it is only trained to predict discrete-time bifurcations.Therefore, one needs to know ahead of time whether continuous or discrete time is a better description for the system.In the case of the chick heart cells, we had prior knowledge that they are well described by a discrete-time dynamical system [48], and therefore appropriate for the classifier.An interesting avenue for future research is to build a classifier that works for both continuous and discrete-time bifurcations.This may be achieved by generating a training library from models with a range of discretised timesteps, from very large steps that generate discrete-time bifurcations, down to the limit of a discrete timestep of zero, where continuous-time bifurcations occur.With a large enough training set, one would not need to assume ahead of time whether continuous or discrete time is a better description for the system.
Our results demonstrate that combining dynamical system and deep learning methodologies can provide EWS for critical transitions that are both more reliable and more descriptive than non-hybrid approaches.This study has set a baseline for prediction performance across a variety of popular discrete-time models and an experimental data set.In providing a code capsule that reproduces this study, we hope to facilitate the development, testing and comparison of related methods.Building a universal predictor for critical transitions is not a job for a single research team [46], and will benefit from a variety of approaches and open source code.Depending on context, critical transitions can be devastating or highly desirable.Improved EWS would allow us to better prevent or promote such transitions.

Generation of training data for the deep learning classifier
Training data consists of simulation data from a library of 50,000 models.The models are generated at random from five different model frameworks, each possessing one of the bifurcations studied (perioddoubling, Neimark-Sacker, fold, transcritical, pitchfork).The models are composed of the normal form of the bifurcation [7], higher-order polynomial terms up to degree 10 with coefficients drawn from a normal distribution, and additive Gaussian white noise (ϵ t ) with amplitude (σ) drawn from a uniform distribution.In each case, the bifurcation occurs at µ = 0.
The model for the period-doubling bifurcation is where α i ∼ N (0, 1).The positive (negative) cubic term yields a supercritical (subcritical) bifurcation, and is chosen at random.The model for the Neimark-Sacker bifurcation is where α ij , β ij ∼ N (0, 1), R(θ) is the rotation matrix and θ ∼ U[0, π] is the angular frequency of oscillations at the bifurcation.The positive (negative) cubic term yields a subcritical (supercritical) bifurcation, and is chosen at random.The model for the fold bifurcation is where α i ∼ N (0, 1).The model for the transcritical bifurcation is where α i ∼ N (0, 1).Finally, the model for the pitchfork bifurcation is where α i ∼ N (0, 1).The positive (negative) cubic term yields a subcritical (supercritical) bifurcation, and is chosen at random.The library is composed of 10,000 models from each framework.For each model, we run a 'forced' simulation where the bifurcation parameter µ is increased linearly across the interval [µ 0 , 0], and a 'null' simulation where µ is fixed at µ 0 .The initial value for the bifurcation parameter µ 0 is drawn from a uniform distribution across all values that correspond to |λ| < 0.8, where λ is the eigenvalue of the Jacobian matrix in the model.This ensures that the training data contains simulations that start close to and far from a bifurcation.For the period-doubling, Neimark-Sacker, transcritical and pitchfork models, this means drawing µ 0 from U[−1.8, −0.2].For the fold bifurcation, this means drawing µ 0 from U[−0.9, −0.1].After a burn-in period of 100 iterations, we simulate each model for 600 iterations and keep the last 500 data points, or the 500 data points immediately preceding a transition if one occurs.We define a transition as a time when the deviation from equilibrium exceeds ten times the noise amplitude σ.We simulate one forced and one null simulation from each model, resulting in 50, 000 forced and 50, 000 null trajectories.To balance the number of entries for each class, we take 10, 000 null simulations at random, resulting in a total of 60, 000 entries in the training data set.Example trajectories for each class are shown in Sup.Fig. 2.

Architecture and training of the deep learning classifier
We use a neural network with a CNN-LSTM architecture and hyperparameters as in [25].This consists of a single convolutional layer with max pooling followed by two LSTM layers with dropout followed by a dense layer that maps to a vector of probabilities over the six possible classes.For training, we use Adam optimisation with a learning rate of 0.0005, a batch size of 1024, and sparse categorical cross entropy as the loss function.We use a training/validation/test split of 0.95/0.025/0.025.We found 200 epochs was sufficient to obtain optimal accuracy on the validation set.
To expose the classifier to time series of different lengths, we censor each time series in the training data.We train two classifiers independently using different censoring techniques.Classifier 1 is trained on time series censored at the beginning and the end, forcing it to learn from data in the middle of the time series.Classifier 2 is trained on time series only censored at the beginning, allowing it to learn from data right up to the bifurcation.The length for each censored time series L is drawn from U [50,500].Then, for Classifier 1, the start time of the censored time series t 0 ∼ U [0, 500 − L] and for Classifier 2, t 0 = 500 − L. The censored time series are then normalised by their mean absolute value and prepended with zeros to make them 500 points in length.We report results using the average prediction of the two classifiers.

Theoretical models
To test the deep learning classifier on out-of-sample data, we simulate a variety of nonlinear, discretetime models, each containing one of the studied bifurcations.To account for stochasticity, we include additive Gaussian white noise.We run forced simulations, where the bifurcation parameter is increased linearly up to the bifurcation point, and null simulations, where the bifurcation parameter remains fixed.To create a diverse set of test data, we vary the noise amplitude (σ) and rate of forcing (rate of change of the bifurcation parameter).We run 100 forced and null simulations at five different noise amplitudes and five different rates of forcing, resulting in 5000 simulations of each model.Values for the noise amplitude are on a logarithmic scale and values for the rate of forcing result in time series of length 100, 200, 300, 400, and 500.Sample simulations for each model at different noise amplitude and rate of forcing are shown in Sup.Fig. 3.The transition time for each forced simulation is taken as the moment when the bifurcation parameter crosses the bifurcation, or the moment when the state variable crosses a threshold, if specified.
(1) Fox model.To test detection of a period-doubling bifurcation, we use a model of cardiac alternans [35] with additive Gaussian white noise.This is given by where D n is the action potential duration of the nth beat, M n is a memory variable, I n is the rest duration following the action potential, T is the stimulation period, τ is the time constant of accumulation and dissipation of memory, α is the influence of memory on the action potential duration, and A, B, C and D are parameters governing the shape of the restitution curve.Following [35], we take A = 88, B = 122, C = 40, D = 28, τ = 180, α = 0.2, which give dynamics in good agreement with a complex ionic model.This yields a period-doubling bifurcation at approximately T = 200.Forced simulations are run with T decreasing linearly on the interval [300, 150] and null simulations are run with T = 300.Values for noise amplitude are 0.
(2) Westerhoff model.To test detection of a Neimark-Sacker bifurcation, we use a simple model of business cycles based on consumer sentiment [30] with additive Gaussian noise.This is given by where Y t is the national income at time step t, a is the level of autonomous expenditures of agents, b and c govern a curve that determines the fraction of income consumed by the agents, and d is the policy-maker's control parameter to offset income trends.Following [30], we take b = 0.45, c = 0.1, and d = 0, which yields a Neimark-Sacker bifurcation at a = 20 corresponding to the onset of business cycles.Forced simulations are run with a increasing linearly on the interval [10, 22.5] and null simulations are run with a = 10.Values for noise amplitude are 0.1 × {2 0 , 2 −1 , 2 −2 , 2 −3 , 2 −4 }.
(3) Ricker model.To test detection of a fold bifurcation, we use the Ricker model [36] with a sigmoidal harvesting term and additive Gaussian noise.This is given by where x t is the population size at time step t, r is the intrinsic growth rate, k is the carrying capacity, F is the harvesting rate, and h governs the steepness of the sigmoidal harvesting term.We take r = 0.75, k = 10, h = 0.75, which yields a fold bifurcation at F = 2.36.Forced simulations are run with F increasing linearly on the interval [0, 3.54] and null simulations are run with F = 0. We define a transitions as the time when x t drops below 0.45.Values for noise amplitude are 0.2×{2 0 , 2 −1 , 2 −2 , 2 −3 , 2 −4 }.
(4) Lotka-Volterra model.To test detection of a transcritical bifurcation, we use the discrete-time analogue of the Lotka-Volterra model, first studied by Maynard Smith [37].This system is especially relevant to arthropod predator-prey and host-parasitoid interactions.The rescaled equations [49] are r relates to the growth rate of the prey (x t ), and c relates to the foraging efficiency of the predator (y t ).We take r = 0.5, which yields a transcritical bifurcation at c = 1, the critical foraging efficiency beyond which the predator population can sustain themselves.Forced simulations are run with c increasing linearly on the interval [0.5, 1.25] and null simulations are run with c = 0.5.We look for early warning signals in the prey population.Values for noise amplitude are 0.01 × {2 0 , 2 −1 , 2 −2 , 2 −3 , 2 −4 }.
(5) Lorenz model.To test detection of a pitchfork bifurcation, we use the reduced discrete Lorenz system, which was first introduced as a demonstration of computational chaos [38].This is given by where state variables and parameters are derived from the full Lorenz equations [38].We take h = 0.5, which yields a pitchfork bifurcation at a = 0. Forced simulations are run with a increasing linearly over the interval [−1, 0.25] and null simulations are run with a = −1.We look for early warning signals in x t .Values for noise amplitude are 0.01

Experimental data
We use data from spontaneously beating aggregates of chick heart cells, which can undergo a transition to an alternating rhythm with the addition of the potassium-channel blocker E-4031 [31].To study the dynamics of these aggregates, we focus on the interbeat intervals-the time between consecutive beats.The transition to an alternating rhythm may be associated with a period-doubling bifurcation.Trends in variance and autocorrelation of the data have been shown to agree with theoretical expectations for a period-doubling bifurcation [21].We define the onset of the period-doubling bifurcation as the first time when the slope of a linear regression of the return map composed of a sliding window of interbeat intervals is below -0.95 for the next 10 beats.We have 23 time series that undergo a period-doubling bifurcation (Sup.Fig. 5, 6), and 18 time series that do not, from which we extract 23 segments at random with a random length between 100 and 500, to serve as null time series (Sup.Fig. 7, 8).For details on the chick heart experiments, see [21,31].
Computing and assessing the performance of EWS EWS are computed using the Python package ewstools [50].This involves first detrending the (pretransition) time series.For the model simulations, we use a Lowess filter with a span of 0.25 the length of the data.For the heart cell data, we use a Gaussian filter with a bandwidth of 20 beats.Variance and lag-1 autocorrelation are then computed over a rolling window of 0.5, which had higher performance than a rolling window of 0.25.The deep learning predictions at a given point in the time series are obtained by taking the preceding data from the time series, normalising it, prepending it with zeroes to make it 500 points in length, and feeding it into the classifier.
To compare performance of variance, lag-1 autocorrelation and the deep learning classifier, we use the AUC (area under curve) score of the ROC curve.The ROC curve plots the true positive rate vs. the false positive rate as a discrimination threshold is varied.We use the Kendall τ value as the discrimination threshold for variance and lag-1 autocorrelation, and the sum of the bifurcation probabilities for the discrimination threshold of the deep learning classifier.performed on Code Ocean at https://codeocean.com/capsule/2209652/tree/v1where the code is accompanied by the necessary software environment.

Figure 1 :
Figure 1: Period-doubling bifurcation in a spontaneously beating aggregate of embryonic chick heart cells following treatment with a potassium channel blocker (E-4031, 1.5µmol).(a) Interbeat intervals (IBI) for consecutive beats.A period-doubling bifurcation occurs at approximately beat 230.Arrows correspond to traces plotted in lower panels.(b-d) Normalised signal from optical imaging of the aggregate's motion.Traces are from a section well before (b), just before (c), and after (d) the perioddoubling bifurcation.

Figure 2 :
Figure 2: Trends in indicators prior to five different bifurcations in the theoretical models.(a-e) Trajectory (gray) and smoothing (black) of a simulation going through a period-doubling, Neimark-Sacker, fold, transcritical and pitchfork bifurcation, respectively.(f-j) Variance of residual dynamics after smoothing, computed over a rolling window (arrow) of size 0.5 times the length of the pretransition data.(k-o) Lag-1 autocorrelation.(p-t) Probabilities assigned by the deep learning (DL) classifier when given all preceding data.Orange line shows probability assigned to the true bifurcation.Gray lines show probabilities assigned to the other bifurcations.Blue line shows the sum of the five bifurcation probabilities.

Figure 3 :
Figure 3: Trends in indicators prior to a period-doubling bifurcation in five chick heart aggregates treated with a potassium channel blocker.(a-e) Inter-beat interval (gray) and smoothing.(f-j) Variance of residual dynamics after smoothing, computed over a rolling window (arrow) of size 0.5 times the length of the pre-transition data.(k-o) Lag-1 autocorrelation.(p-t) Probabilities assigned by the deep learning (DL) classifier to the period-doubling bifurcation (orange) and the other bifurcations (gray).Blue line shows the sum of the five bifurcation probabilities.

Figure 4 :
Figure 4: ROC curves for predictions of an upcoming transition in model and experimental data.ROC curves compare the performance of the deep learning classifier (DL, blue), variance (Var, red) and lag-1 autocorrelation (AC, green).For models (a-e), performance is assessed on 5,000 simulations with different noise amplitudes and rates of forcing.For experimental data (f), performance is assessed on 46 experimental runs.The area under the curve (AUC), abbreviated to A, is a measure of performance.Insets show the probabilities assigned by the classifier to each type of bifurcation (orange being the true bifurcation) among the trajectories approaching a transition.(a) Fox model going through a period-doubling bifurcation [35].(b) Westerhoff model going through a Neimark-Sacker bifurcation [30].(c) Ricker model going through a fold bifurcation[36].(d) Lotka-Volterra model going through a transcritical bifurcation[37].(e) Lorenz model going through a pitchfork bifurcation[38].(f) Chick heart aggregates going through a period-doubling bifurcation[31].Predictions are made 80% of the way through the pretransition data for the models, and 60-100% of the way for the experimental data.PD: period-doubling.NS: Neimark-Sacker.TC: transcritical.PF: pitchfork.

Figure S2 :
Figure S2: Sample simulations and corresponding bifurcation diagrams for each class in the training data.Panels show the trajectory of x n (blue), the corresponding bifurcation diagram (black) and the section used for training in each case (between the green vertical lines).The six possible classes are (a) null, (b) period-doubling, (c) Neimark-Sacker, (d) fold, (e) transcritical, and (f) pitchfork.These samples show supercritical versions of the period-doubling, Neimark-Sacker and pitchfork bifurcations, though subcritical examples are also present in the training data.The end time for the extracted data is taken as the time the bifurcation is crossed (t = 600) or the first time the deviation from equilibrium reaches ten times the noise amplitude, whichever comes first.

Figure S3 :
Figure S3: Sample simulations of each theoretical model approaching a bifurcation under a range of noise amplitudes (sigma) and rates of forcing (colour).The rates of forcing are chosen such that the bifurcation is crossed after 100 (green), 300 (red) and 500 (blue) units of time.The Fox model goes through a period-doubling bifurcation, the Westerhoff model a Neimark-Sacker bifurcation, the Ricker model a fold bifurcation, the Lotka-Volterra model a transcritical bifurcation, and the Lorenz model a pitchfork bifurcation.Model parameters are provided in Methods.

Figure S4 :
Figure S4: AUC score (area under the ROC curve) at different values of noise amplitude (sigma) and rate of forcing (RoF) for each indicator (columns) and theoretical model (rows).At each combination of RoF and sigma, we run 100 forced and null simulations, resulting in a total of 5000 simulations for each model.ROC curves are computed using predictions at 80% of the way through the pretransition time series.The bifurcations associated with each model are period-doubling (Fox), Neimark-Sacker (Westerhoff), fold (Ricker), transcritical (Lotka-Volterra), and pitchfork (Lorenz).

Figure S5 :
Figure S5: Proportion of predictions from the classifier that favoured the correct bifurcation for forced trajectories of different noise amplitude (sigma) and rate of forcing (RoF).For each combination of RoF and sigma, 100 forced trajectories are simulated, and predictions are made 80% of the way through the pretransition time series.The bifurcations associated with each model are period-doubling (Fox), Neimark-Sacker (Westerhoff), fold (Ricker), transcritical (Lotka-Volterra), and pitchfork (Lorenz).

Figure S6 :
Figure S6: Trends in indicators prior to a period doubling bifurcation in chick heart aggregates (a-l).See Fig. 3 caption for details.

Figure S7 :
Figure S7: Trends in indicators prior to a period doubling bifurcation in chick heart aggregates (m-w).See Fig. 3 caption for details.

Figure S8 :
Figure S8: Trends in indicators for chick heart aggregates that do not undergo a period-doubling bifurcation (a-l).See Fig. 3 caption for details.

Figure S9 :
Figure S9: Trends in indicators for chick heart aggregates that do not undergo a period-doubling bifurcation (m-w).See Fig. 3 caption for details.