Let's do the time-warp-attend: Learning topological invariants of dynamical systems
Authors: Noa Moriel
Matt Ricci
Mor Nitzan
Venue: ICLR 2024 poster (2024)
Keywords: dynamical systems
bifurcations
topological invariance
Hopf bifurcation
physics-informed machine learning
augmentation
single-cell RNA-sequencing
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Dynamical systems across the sciences, from electrical circuits to ecological networks, undergo qualitative and often catastrophic changes in behavior, called bifurcations, when their underlying parameters cross a threshold. Existing methods predict oncoming catastrophes in individual systems but are primarily time-series-based and struggle both to categorize qualitative dynamical regimes across diverse systems and to generalize to real data. To address this challenge, we propose a data-driven, physically-informed deep-learning framework for classifying dynamical regimes and characterizing bifurcation boundaries based on the extraction of topologically invariant features. We focus on the paradigmatic case of the supercritical Hopf bifurcation, which is used to model periodic dynamics across a wide range of applications. Our convolutional attention method is trained with data augmentations that encourage the learning of topological invariants which can be used to detect bifurcation boundaries in unseen systems and to design models of biological systems like oscillatory gene regulatory networks. We further demonstrate our method's use in analyzing real data by recovering distinct proliferation and differentiation dynamics along pancreatic endocrinogenesis trajectory in gene expression space based on single-cell data. Our method provides valuable insights into the qualitative, long-term behavior of a wide range of dynamical systems, and can detect bifurcations or catastrophic transitions in large-scale physical and biological systems.

Published as a conference paper at ICLR 2024
LET'S S DO THE TIME- TIME-WARP-ATTEND: LEARNING TOPOLOGICAL INVARIANTS OF DYNAMICAL SYSTEMS
Noa Moriel The Hebrew University noa www.meereimll.....mmm........aaa www.iaaa. ac.
Matthew Ricci
The Hebrew University
matthew ricci@mail hujji ac
Mor Nitzan
The Hebrew University
mor. nitzaanmaail..nnmmmmmmaaaaaaaiim huji.aa.. ac.1l
ABSTRACT
Dynamical systems across the sciences, from electrical circuits to ecological ner- works, undergo qualitative and often catastrophic changes in behavior, called bifurcations, when their underlying parameters cross a threshold. Existing methods predict oncoming catastrophes in individual systems but primarily time-series- based and struggle both to categorize qualitative dynamical regimes across diverse systems and generalize real data. To address this challenge, we propose a deep-learning framework for classifying dynami-
data-driven, thysically-informe cal regimes and characterizing bifurcation boundaries based on the extraction of topologically invariant features. We focus on the paradigmatic case of the super- critical Hopf bifurcation, which used model periodic dynamics across a wide range of applications. Our convolutional attention method is trained with data augmentations that encourage the learning of topological invariants which can be used to detect bifurcation boundaries in unseen systems and to design models of biological systems like oscillatory gene regulatory networks. We further demon- strate our method's use in analyzing real data by recovering distinct proliferation and differentiation dynamics along pancreatic endocrinogenesis trajectory in gene expression space based on single-cell data. Our method provides valuable insights into the qualitative, long-term behavior of a wide range of dynamical systems, and can detect bifurcations or catastrophic transitions large-scale physical and biological systems.
dynamical systems, bifurcations, topological invariance, Hopf bifurcation, physics-informed machine learning, augmentation, single-cell RNA-sequencing.
Keywords:
INTRODUCTION
Gradually changing the parameters of dynamical system will typically result in innocuous, quantita- tive alterations to the system's behavior. However, exceptional cases, parameter perturbations will result in qualitative, potentially catastrophic change in the system's 's dynamical behavior called bifurcation. Bifurcations may lead to extreme consequences, as when a beam buckles and breaks under increasing weight (Strogatz, 1994), when disastrous oscillations presage the collapse of poorly designed bridges, like the infamous Tacoma Narrows Bridge (Billah & Scanlan, 1991; Deng et al., 2016), or when airplane exceeds speed threshold, leading to dangerous self-excited oscillations in its tail (Weisshaar, 2012; Dimitriadis, 2017). In biological systems, bifurcations marking the transition away from oscillatory dynamics on cellular or organ-level scale, such those regulating
"quual contribution. Corresponding author.

---

Published as a conference paper at ICLR 2024
blood pumping and circadian rhythms, can have dire health consequences (Xiong & Garfinkel, 2023; Qu et al., 2014).
Understanding bifurcations which transition between dynamical regimes, be they in physics, engi- neering, or biology, is a subject of intense research in the study of nonlinear systems (Kuznetsov, 1998; Karniadakis al., 2021; Ghadami & Epureanu, 2022). One of the central challenges across all applications is that distinct dynamical regimes can be qualitatively identical (e.g. the emergence of oscillations) while being quantitatively very different (e.g. the amplitude or frequency of those oscillations); smoothly warping the state space has no qualitative effect on the dynamics. Put another way, a detector for different dynamical regimes and the transitions between them must be selective for topological structure of underlying dynamics but invariant to the geometric nuisances of those dynamics. Predicting the onset of such bifurcations often requires the explieit use of governing equa- tions, which must be hardwired or inferred. Hence, the analysis of real data, where such equations are often unknown or difficult to define, is particularly challenging (Kuznetsov, 1998). Applications to real data either use hardwired (Wissel, 1984; Wiesenfeld, 1985; Scheffer et al., 2009) or learned (Bury et al., 2023) time series features, but no approach, to our knowledge, explicitly addresses the selectivityiivaaiance invariance trade-off off inherent to the topological nature of the problem. Additionally, where vector field data is readily available and time-series tracing is not, e.g., in capturing gene expression of cells in high-throughput single-cell experiments, it remains unclear how extract topological information either from such noisy and sparse vector fields directly or from their integration.
We address this challenge by learning universal dictionary of dynamical features which is explicitly encouraged to be invariant to geometric nuisances but selective for the topological features which define a particular type of bifurcation (Fig. 1). This is achieved primarily by the use of topological data augmentation: each sample of our training data is a velocity field which comes bundled with numerous deformed but identically labeled versions. The features extracted from these augmented samples are then passed to a supervised classifier which must contend with the substantial geometric nuisances induced by the augmentation scheme. We demonstrate our method for the ubiquitous supercritical Hopf bifurcation whereby a stable fixed point (or "equilibrium") loses stability and transitions into a regime exhibiting a stable limit cycle. We train on synthetic, augmented data generated from a simple oscillator prototype. Despite being only trained on data from one simple equation, our architecture can generalize to radically different systems, such as the nonlinear Selkov glycolysis oscillator and the Belousov-Zhabotie Zhabotinsky model of chemical oscillations.
Our approach is geared towards realistic, high--dmensionall. noisy, and complex data. For instance, we apply our technique to the repressilator, a gene regulatory network which can generate oscillatory gene expression, and infer the probability for oscillatory behavior under conditions of varying gene transcription and degradation rates. Finally, we apply our method to the analysis of complex, large-scale biological data. By extracting features from single-cell gene expression velocity data during the differentiation of pancreatic cells, our method can distinguish between proliferation and differentiation dynamics with high fidelity. Our results indicate that our approach can be applied to broad range of dynamical systems, providing valuable insights into their dynamic characteristics and bifurcation behavior !
The main contributions of this paper are:
1. A novel approach to topologically invariant feature learning in dynamical systems. The approach uses warped vector field data as augmentations to generalize across systems exhibiting similar dynamics in terms of their topological representation in phase space, allowing for a more selective and invariant understanding of bifurcations in diverse real- world systems.
2. Invariant classification of complex and nonlinear synthetic data, including models of diverse chemical and biological systems, achieved by leveraging knowledge of a single prototypical system.
3. The identification of bifurcation boundaries using classifier robustness.
4. The recovery of distinct proliferation and differentiation dynamics along the pancreatic dif- ferentiation trajectory from single-cell gene expression data, providing practical application of our method to the analysis of biological systems.
Cooee available at: http://ggthhb.coommmiraaaambmmmt
2

---

Published as a conference paper at ICLR 2024
X = (xx;0) vector field, y(f) attractor label
augmentations
(1) angle transformation
(2) self-atiention layers
9() predicted
attractor
periodi
point
Xy space
Ny space
xy space
Ky space
CNN
Prodict (4) for unsecn system $
Infer bifurcation boundary
Predicting proliferating and differentiating regimes
space
parameter space
gene expression space
Figure 1: Overview: Our framework classifies long-term dynamical behavior of point vs periodic attractor across systems and regimes by leveraging pelloniccallyinvriinn augmentations of single prototypical system. To encourage generalization,
magnitudes and angular representations instead, and (2) deploy attention layers which focus on learning essential dynamical cues. Using this framework. we are able to classify and recover the bifurcation boundaries of complex, unseen systems. We apply this method to the problem of recovering cell-cycle scores which distinguish between proliferating and differentiating cell populations. All vector fields plotted as streamlines.
abstract system-specifie
RELATED WORK
Classical methods for detecting limit cycles. Solving non-linear differential equations is in general difficult, and it has been recognized since the time of Poincaré that reasoning about dynamical systems' characteristics can be made easier by qualitatively delineating solutions according their long-term behaviors. For autonomous (i.e., time-invariant) systems in two dimensions, if the dynamics do diverge, these behaviors are either fixed points limit cycles, and transitions between these dynamical regimes as a result of changing parameters are called bifurcations. More explicitly, bifurcation these cases corresponds to change in the topological structure of a phase portrait as parameter of the dynamical system is varied.
Full characterization of location and number of limit cycles in two-dimensional. polynomial systems is challenging open problem in dynamical systems theory (Hilbert, 2019; Smale, 2000). Early theoretical results guarantecing the existence or absence of a limit cycle in two-dimensional and Bendixson-Dula
autonomous systems the theorems of Further guarantees are possible, for example, when the system can be reduced to the form a so-called Liénard equation. None these results, however, is constructive or easy use in practice. Center manifold theory (Carr, 2012), for its part, analyzes dynamics on representative subspaces where systems reduce normal forms having known bifurcation behaviors, though this method requires knowledge of the underlying governing equations. The same limitation applies to Conley- Morse theory, which represents dynamics as a graph whose nodes correspond to long term-term& behaviors (Conley, 1978). Importantly, none of these approaches is easily adapted to real data, where equations are general unknown and system measurements can sparse, incomplete and noisy.
respectively.
Data-driven methods for detecting limit cycles. Data-driven methods fall into non-topological and topological classes. In the former class, several approaches use compressed sensing of trajectory data to reduce noise, followed by identification of system parameters and associated bifurcation diagrams (Wang et al., 2011; 2016), though this approach requires inferring governing equations for each new sample. Other work sought to reduce vector field measurements et al., 2007b;a) for which a database over parameter regimes can be constructed via combinatorial dynamics exploration and graph equivalence comparisons (Arai et al., 2009). Yet other approaches
Conley-Morse graph (Chen
3

---

Published as a conference paper at ICLR 2024
used estimates of Lyapunov exponents to delineate between periodic and stationary systems, although, as we show below, this approach can be sensitive to noise (Stoker, 1950; Zhang et al., 2009).
On the topologiaa front, there is a growing body of work using persistent homology to detect periodicity in time series data (Ravishanker & Chen, 2021). For example, Perea & Harer (2015) extracted persistence barcodes from delay embeddings of time series data to measure periodicity, although this process is computationally expensive and there is no standard approach to classifying barcodes. Persistent homology approaches often assume that the underlying bifurcation parameter is known (Tymochko et al., 2020), which greatly simplifies the problem at the cost of realism. To our knowledge, there is no data-driver method which is efficient (not requiring expensive barcode calculations for each sample), general (applicable to novel systems without retraining) and robust (not sensitive to noise).
Generalization across dynamical systems. Comparatively few approaches exist for the simultaneous analysis of diverse dynamical systems and their bifurcation behavior. Work by Brown & Sethna (2003); Transtrum & Qiu (2014) and Quinn et al. (2021) mapped families of dynamical systems to a manifold whose geometry encodes different qualitative behaviors, but these methods have no straightforward way of identifying bifurcations between these behaviors. More recent work from Ricci et al. (2022) uses a data-driven approach called phase2vec to learn manifolds of dynamical systems, leading to dynamical embeddings which encode qualitative properties of the underlying systems, though it is not tailored for the detection of distinct dynamical regimes and their associated bifurcations per se.
Diffeomorphic augmentations for generalization. Diffeomorphic transformations which preserve topology in the context of dynamical systems have been used principally in computer vision. In order to train a classifier of a limited image set of biological entities such as diatoms, Vallez et al. (2022) used diffeomorphic augmentations to expand its train data, however, since the labels of these data were not necessarily preserved by these transformations, a manual post-filtering was required. In a second set of works concerning dynamical systems specifically, Khadivar al. (2020) and Bevanda et al. (2022) fit diffeomorphic transformation into closed or linearized forms, and work by Bramburger et al. (2021) used an autoencoder approach to reduce a complex dynamical system to topologically conjugate (i.e., dynamically equivalent) form of discrete-time mapping. Here, we use diffeomorphisms to curate a set of consistently labeled dynamical augmentations.
METHODOLOGY
Topological equivalence in dynamical systems. Let à = {z; 0) for z E R² be parameterized dynamical system with parameter vector 8 O. A solution to this equation, F( t), with an initial condition, <x(1 = 0) = zo. is called trajectory through wo. A standard approach to understanding dynamical systems (Kuznetsov, 1998) is to characterize them by the topology of all of their trajectories, construed as functions of z and t or just For example, do trajectories extend indefinitely, stop at point (a fixed point) or loop back on themselves (a limit cycle)? More precisely, we say that f has a fixed point x. and 0 = 8r if f(="; 0²) = 0. We say that f has limit cycle at 0" if there exists trajectory of f which is closed and isolated (neighboring trajectories are not closed). Importantly, these properties are preserved when the state space is subjected
topology- prreggeviigeeeingggg map, making them topological invariants of the dynamical systems. However, these properties may change as 0 changes; i.e. when the system undergoes bifurcation. For example, in system that undergoes supercritical Hopf bifurcation, there exists a parameter 0 such that below a certain threshold of 0 the system exhibits a stable fixed point, and above it, the fixed point loses its stability and stable limit cycle emerges. See Appendix 7.1 for definitions.
Dynamical equivalence and topological augmentations. what follows, we will consider autonomous systems with either a single stable fixed point (i.e., equilibrium) or stable limit cycle with an unstable fixed point in its interior. We term these two regimes as dynamical classes, denoted v(f), and we denote the fixed point and limit cycle classes considered here the "pre-Hopf" and "post- classes, respectively. As discussed above, these classes are invariant to topology- preserving deformations of the state space. Here, will use diffeomorphisms, h E H which induce an equivalence relation on dynamical systems in a standard and well-studied way: we have v(s) = v(ʃ h) for all h. E H so that given H defines equivalence relation which partitions the
4

---

Published as a conference paper at ICLR 2024
space of dynamical systems into its dynamical classes. We refer to parameterized systems, f.g such that = f a h as being dynamically equivalent
In this light, classification becomes a problem of learning representations of dynamical systems which are invariant to all h E И but selective for the underlying dynamical class. To that end, we learn features on numerous, transformed versions of prototypical system, simple oscillator f defined for a point at radius r and angle & with parameters 0 = (a,w):
f r(a r²)
Ø w,
which has a sparse functional form and undergoes supercritical Hopf bifurcation for any w * 0 at 0- = (0, w) We validate these learned features by measuring how well they encode the dynamical class and bifurcation boundary of testing data, which, though superficially very different from the training data, having different functional forms and parameter regimes, are nevertheless in the same equivalence classes induced by H.
Our training data is constructed by first acquiring planar vector fields. {f(z, 0) E R2} evaluated at all points, x, on fixed, 64 x 64 lattice. Parameters, 0, are sampled uniformly a E [-0.5, 0.5), w 11,, and each system is assigned a label according to whether it is in the pre- or post-Hopf regime. Then, the dataset is augmented with diffeomorphisms generated using monotonic rational-quadratic splines. In particular, we use the spline construction from Durkan et al. (2019), where instead of training the neural spline flows proposed in that work, we use its random initialization which generates a diverse set of diffeomorphic augmentations. For details on data augmentation, see Appendix 7.2.1.
Convolutional attention training. We construe our full, augmented data set as collection of labeled arrays in R2 22x64x x4 64 which we use to train a convolutional network. Our architecture is adapted from the discriminator of self-attention generative adversarial network (SAGAN) of Zhang al. (2019), coupled with a multi-layered perception (MLP). For standardization of the data, we reduce each vector field to only its angular components, ignoring magnitudes which carry little information about eyclic behavior. The feature extraction module consists of four convolutional blocks with self-attention mechanism into the last two convolutional layers. Our reasoning for incorporating self self-attention is that only few dynamical keypoints in the vector field (i.e., the fixed point or the limit cycle) carry information about class label, and the rest of the array can be (nearly) ignored. Together, the data warping and attentional components of our pipeline comprise our eponymous method. Features from this network are passed to a multilayer perceptron classifier which is trained with a supervised, cross-entropy loss. Because of our augmentation scheme, this classifier must contend with data which, though radically different in appearance, belong to the same dynamical class, encouraging the classifier to learn the true dynamical invariances underlying the Hopf bifurcation. Training and architecture details are described in Appendix 7.2.2.
4 RESULTS
4.1 TRANSFERRING KNOWLEDGE OF DYNAMICAL CLASSES FROM SIMPLE PROTOTYPES TO COMPLEX TEST DATA
We evaluate our model and its baseline competitors on a diverse, synthetic dataset of noisy (zero- mean gaussian), real-world systems drawn from across the sciences. Test datasets include the simple harmonic oscillator ("SO"), its warped variant ("Augmented SO"), a system undergoing supercritical Hopf bifurcation as described in (Strogatz, "SuperrSuprriiiiccal Hopf"), examples based on I. Liénard equations ("Liénard Polynomial" and "Liénard Sigmoid")) and systems used to model oscillatory phenomena in nature and engincering applications, namely the van der Pol system, used for modeling electrical circuits in early radios ("Van der Pol"), the chemical oscillations ("BZ reaction"). and the Selkov oscillator, used in modeling glycolitic cycles (" ("Selkov"), see Appendix 7.3 and Strogatz (1994).
model, model of
Our baselines comprise: (1) "Critical Points", a heuristic based on critical points identification as derived from Helman & Hesselink (1989; 1991) algorithm, (2) "Lyapunov" which distinguishes pre-
"e choose the neutral term "dynamical equivalence" to distinguish our notion from the similar notions of topological conjugacy and topological equivalence (Kuznetsov, 1998). See Appendix 7.1.2 for details.
5

---

Published as a conference paper at ICLR 2024
and post-Hopf classes according to the value of the Lyapunov exponent of a time series integrated from the vector field, (3) "Phase2vee", vector field embeddings from Ricci et al. (2022), a state- of-the-art seep-leerning method for dynamical systems embedding, and (4) "Autoencoder", latent features of a convolutional vector field autoencoder, see Appendix 7.4. As with our method, for all baselines except Critical Points, classifications are produced by a linear classifier trained on features from augmented simple harmonic oscillator data (Augmented SO) in order to discriminate between point and cyclic attractors. The accuracy average and standard deviation for our model and the trained baselines are computed over 50 training runs.
We find that our framework generalizes across a wide range of dynamical systems and classifies distinct dynamical regimes in noisy (zero-mean gaussian) conditions better on average than baseline methods (see Table 1, Table A4, and examples in Fig. A8). Critical Points achieves nearly perfect classification on noiseless data but is easily disrupted by noise (e.g., the perfect classification of Supercritical Hopf drops to 56% with Gaussian noise
0.1). Likewise, Lyapunov suffers a substantial decrease in accuracy for classical oscillators (e.g.. 89% accuracy of Supercritical Hopf' reduces to 70% accuracy with Gaussian noise o = 0.1). Both with and without noise, Phase2vec extends mainly to SO and Supercritical Hopf. Autoencoder obtains competitive results but exhibits little generalization to Liénard Polynomial system, and essentially no generalization to Selkov system. We also test resilience to noise in the Augmented SO data (see Table 2). Despite the substantial distortion of the vector fields caused by noise, the accuracies of our model and that of autoencoder representation only gradually decrease, unlike the sharp decrease in accuracy for the other baselines.
Table 1: Test accuracy with added Gaussian noise of a = 0.1.
SO
Aug. SO
Supercritcaa Hupt
Liénard Puly
Läčmurd Signoid
Van der Pol
BZ Rescttion
Sellow
Average
Our Model
0.98:00.1
0.93-20.01
0.99-0.01
0.86=00113
0.92::0.07
0.83:20.15
0.82::.11
0.65:00.03
0.87
Critical Poinss
0.56
0.70
0.49
0.56
0.51
0.60
Lyspunov
0.60
0.71
0.67
0.87
0.85
0.57
0.72
Phase2vee
0.73=0008
0.71
0,88=003
0,49=20.06
0.48==0.0
0.49:00044
0.50000
0.49:0000
0.6
Autoencoder
0.95 20000
0.8720.01
0.7700011
0.63:20.16
0.882:0.99
0.811:0.16
0.82::0133
0.49::0.02
0.8
Table 2: Evaluation of Augmented SO with added Gaussian noise of scale a.
0.10
0,20
0.30
0:40
0.50
0.50
0.70
0.30
0.90
1.00
Model
0.93.10.11
0,93..001
0.9310.01
0,911.002
0.86..0.0
0.78.10.06
0.7.10.06
0.54.D004
06610.01
0,57.0003
0,51__0022
Criticl Points
0,93
0.53
0.507
0.51
0.51
0.51
0.50
Lyspanov
0.71
0,52
0.50
0.50
0.50
0.50
0.50
0.770003
0.7-0.05
0.59+0011
0.51+001
0.5+000
0,5+000
0.5+0.0
0.5+00
Autoencoder
0,36.1011
0.87.0001
0.4410022
0,810.04
0.7610055
0.72.10.05
0.6910.05
0,68
6.65.0.06
0.6410055
0.62.0005
In Appendix 7.4, we based on flattened vector of coefficients of degree three polynomial fit to each vector field with least squares, common method for bifurcation analysis, although this resulted in a poor average accuracy of 53%.
we report the results of another bascline, "Parameters" which trains a linear classifier
Furthermore, we find that all modules in our pipeline (topological augmentation, angular represen- tation of data, self-attention layers) contribute, on average, to accuracy in complex systems, see Appendix 7.5. Concretely, the average classification accuracy of the framework, trained on data without augmentations. using vector representation, lacking attention layers, or without any of these, 78%, 82%, 77%, 68%. respectively, and notably, except for the ablation of attention layers, other ablated models do not generalize to the complex BZ reaction and Selkov systems. We also found that topological augmentations are particularly necessary when generalizing to systems with fixed point far from the center of the grid. For example, "centered" systems such as Liénard Sigmoid and Liénard Polynomial, are well-characterized without augmentations, reaching 0.96 and 0.94 accuracy respectively, whereas systems that are not centered, e.g. BZ reaction whose fixed point is close to the edge of the grid, require augmentations for generalization.
We also evaluate the robustness of our method as a function of (1) the training dataset size, where we find good generalization even relatively small data sets (0.84 accuracy after training on ~ 2k vector fields), (2) vector field resolution, where performance is stable when doubling the vector resolution (0.86 accuracy), and (3) architectural and optimization hyperparameters (learning rate, training epochs, etc.), obtaining almost consistently an overall accuracy of between 0.85 and 0.88. See Appendix 7.5 and Table A4 for details.

---

Published as a conference paper at ICLR 2024
4.2 RECOVERING BIFURCATION BOUNDARIES
Classifier confidence is expected to decrease for systems near the bifurcation boundary since limit cycles in this regime shrink below the vector field resolution, making them indistinguishable from fixed points. This effect can be used to pinpoint the bifurcation boundary in parameter space. For example, when we plot each system's average prediction as function of its bifurcation parameter (for Augmented SO, the parameter values are chosen as the original a. w). we find an inverse relation between the parametric distance from the bifurcation boundary and classifier confidence (see Fig. 2, Fig, All and Fig. A12 for bifurcation diagrams of alternative baselines). An exception is the upper boundary of the Selkov system where instances of high b value are misclassified as exhibiting cyclic dynamics, and hence, the upper bifurcation boundary in the (Fig. 2, right plot). We explain the systematic misclassification as result of damping which makes trajectories move in slow, tight spirals towards the fixed point attractor in a manner numerically indistinguishable from limit cycles (Fig, A7). These overall suggest that our method can be geared towards recovering the bifurcation boundaries of unseen systems.
parameter plane is undetected
Augmenced SO 68
Supercribcal
Lienard Paly
Lionard Siamaid
BZ onaction
0.
DF
-1.
-6.8
0.0
-0.4
03
0.
-04
055
04
-0.8
0.0
2.
-0.8
0.
0.
-4.8
0.5
D
0.04
Figure 2: Inference of biffrrattion boundaries. Each system is plotted at coordinates of its simulated parameter values and colored according to its average cycle prediction across 50 training re-runs. The true parametric bifurcation boundary is plotted in black. Colors scaled between perfect cycle classification (red) and perfect point classification (blue) with white at evenly split classification.
4.3 PREDICTING OSCILLATORY REGULATION OF THE REPRESSILATOR SYSTEM
In order to demonstrate the viability of our method for higher-dimensional systems under realistic noise conditions, we apply our topological approach to computational model in well-studied biological setting, the so-called repressilator model (Elowitz & Leibler, 2000). The repressilator system is synthetic, genetic regulatory network that was introduced by Elowitz and Liebler into bacteria to induce oscillatory regulation of gene expression (Fig. 3). The system is modeled by the mRNA and protein counts of three genes (TetR, Lacl and À phage cl) that inhibit each other in cyclic manner (see Appendix 7.8 for governing equations). Theoretically and experimentally. the existence and characteristics of oscillations have been shown to depend on numerous factors like the rates of transcription, translation, and degradation, as well as intrinsic and extrinsic noise (Verdugo, 2018;
et al., 2016). Here we test our method in classifying a regime of the repressilator which has been shown to undergo a Hopf bifurcation, though, unlike the previous complex simulated systems, it spans six dimensions and is constructed from noisy cellular series.
To generate 2D vector field data, we consider trajectories of cellular gene expression measured by protein levels of TetR and LacI (Fig. 3, middle). Different instances of simulations share the initial condition but vary by the rate of transcription, rates, B. We sample 100 cell states along the temporal trajectory and add Gaussian noise of a 0.5 to the gene expression states (to protein and mRNA levels). For each set of measurements, interpolate an input 64x64 vector field, based on their velocities across Letl protein phase space, examples at Fig. A14 and Appendix 7.8 for further details. We predict for each vector field its dynamical class, achieving 87% classification accuracy and proximate bifurcation boundary to the yittttil-l---evveedd boundary (Fig. 3, right (Verdugo, 2018)). In Fig. A13 we show prototypes of trajectories across the a - 3 parameter space, where we see that where we fail to classify a point attractor, trajectories of cells loop closely in the TetR-Lacl protein plane and generate data samples that highly resemble trajectories with a periodic attractor rather than with a point attractor.
and by the ratio of protein and mRNA degradation
4.4 DISTINGUISHING DYNAMICAL REGIMES IN SINGLE-CELL DATA OF THE PANCREAS
Last, we leverage our method to characterize distinct dynamical regimes corresponding to different biological processes in real data. Specifically, we focus distinguishing proliferation vs. differ-
7

---

Published as a conference paper at ICLR 2024
Simulated oscillating trajectory example
Predicted oscillations
plac
Protsin Por Coll
Lacl
Aee
02
15
Figure 3: Oscillations in a repressitator xystem. Example of oscillations in the repressilator system viewed as a time-series (left) or projected onto TetR-Lac protein plain (middle). (right) Predicted oscillatory behavior by our model oulput upon varying transcription rate without repression (ai) and the ratio of protein and mRNA degradation rates (5)) The thcoretical boundary is marked in black. See details in Appendix 7.8.
entiation dynamics in single-cell RNA-sequencing data. A population of proliferating, precursor cells expresses genes associated with DNA duplication in a cyclic fashion (Tirosh al., 2016; Schwabe et al., 2020; Karin et al., 2023). Cell differentiation, on the other hand, and similarly to the Waddington landscape abstraction (Waddington, 2014), is often characterized by cellular flow in gene expression space towards localized, differentiated state (Sáez et al., 2022). We show that our model can detect the qualitative change in the flow of cells in gene expression space as they mature and differentiate out of their precursor cell state, moving from dynamical regime characterized by periodic behavior to that characterized by fixed points.
We focus the process of pancreas endocrinogenesis where Ductal and Ngn3+ precursor cells rapidly proliferate to give rise to functional cell types Alpha, Beta, Delta and Epsilon (Bastidas-Ponce al., 2019); In that study, single-cell RNA-sequencing dataset was collected where the expression of over 2,000 genes are measured per cell, for total of 3,600 cells (Fig. 4, left). While the RNA expression profile, providing a gene expression vector per cell, provides for each cell its location within gene expression space, RNA velocity (Bergen et al., 2021; Gorin et al., 2022) can be used to compute the effective velocity of each cell (based on differences in spliced and unspliced RNA counts). We use se Velo (Bergen et al., 2020) to infer cellular velocities in the high-dimensionall gene expression space and to project them to a 2D UMAP representation. Simultaneously, we also compute for each cell 'cell-cycle score', which captures the extent of expression of cell-cycle genes (in phase and in grouped G2 and M phases, Fig. 4, middle; see Appendix 7.9, Bergen al. (2020); Tirosh et al. (2016)). For each set of cells, we interpolate an input 64x64 vector field, based on their RNA velocities, over the 2D UMAP representation.
To generate multiple frames from the trajectory, we randomly split cells of each cell type into groups of at least 50 cells resulting in 69 (sparse, noisy) vector fields with corresponding cell-cycle score (Fig. 4, middle, Appendix 7.9). We then used our model (still trained solely on simulated augmented simple oscillator data) to characterize the dynamic regimes of the cellular behavior. Our model can successfully identify proliferative, cyclic behavior within the single-cell data, as reflected by 94% accuracy in predicting cyclic behavior matching high cell cycle scores (Fig. 4, right), substantially exceeding the prediction accuracy reached by alternative baselines (Critical Points: 54%, Lyapunov: 49%, Phase2vec: 30%, Autoencoder: 30%). As opposed to the cell-cycle score, our model does not use any prior knowledge regarding cell-cycle related genes, and is only based on the structural features of the cellular vector fields, and therefore can also be employed when such prior knowledge unavailable.
LIMITATIONS
Our framework learns to be invariant to geometric (as opposed to topological) properties of systems through topological data augmentation and convolutional attention, achieving competitive results in bifurcation detection across real-world systems. Nevertheless, we note here several important limitations of our method in its current form. The method is currently confined to the important but circumscribed setting of two-dimensional dynamics. Scaling up to much higher dimensions in these contexts presents technical challenges like the exponential growth of input size with dimension, as well as theoretical challenges such as the lack of structural stability in higher dimensions and the introduction of chaotic regimes. Since inherently low-dimensional dynamics often underlie high-
8

---

Published as a conference paper at ICLR 2024
Velocity gene expression space
Tree cell-cycle
Predicted cell-cycle score
Figure 4: Prolifration-to-iiffeeettttttete transition in pancreas development. (left) Gene expression velocity of cells undergoing pancreas endocrinogenesis in UMAP space. Colors distinguish between different cell types. (middle) Cell-cycle score hased the fraction of cells with high S-phase gene expression (with S-score 0.4). (right) Predicted cell-cycle score computed as the average classification prediction across 50 training re-runs. See details in Appendix 7.9.
dimensional measurements, especially in biology (Xiong & Garfinkel, 2023; Khona Ficte, 2022), high dimensions can be mitigated by hardwired, or optimized dimensionality reduction methods as we deploy for the -x-dimensional repressilator and for analyzing over 2,000 gene expression (dimensional) dynamics of the pancreas single-cell dataset. Additionally, we used a supervised, multi-class setup so that the classifier can predict that both limit cycle and fixed point co-exist. Empirically, we find that our method learns that the point and periodic attractors are exclusive, in fact, the learned logits consistently negate each other. For emphasis, we tested our trained model on samples representing subcritical Hopf bifurcation having regime with point and periodic attractors whose logits were of intermediate value but still negations, see Appendix 7.7. In the future, we aim to extend combinatorial annotations of invariant sets, as elaborated in the discussion.
DISCUSSION
The problem consider is the data-driven classification of diverse set of both real and synthetic dynamical systems into topological equivalence classes. This broad scope fundamentally distin- guishes our work from most existing methods in bifurcation analysis which deal exclusively with the prediction of oncoming catastrophes from time series data of single, isolated systems. The fact that our method can be simultaneously applied to numerous systems of diverse functional and parametric forms demonstrates that our model has learned an abstract notion of these equivalence classes: it can applied
to novel scenarios without the need for retraining. The integration of physics-infond informee data augmentation within the pipeline is crucial for attaining this high-level performance for complex systems that are substantially different from the prototype system underly- ing the training data. Our model captures visual uncertainty in systems near bifurcation boundaries, allowing for the extraction of an equation-free representation of these transitions. Furthermore, our model demonstrates robustness across samples disrupted by noise, withstands data produced from complex, high- digh-dimensional time-series that is prominent for real-world applications such as the design of synthetic gene regulation, and successfully detects in real single-cell data the qualitative shift in cellular flow topology, corresponding to shifts of biological behaviors, from precursor to differentiated cell states, without relying on gene semantics.
Our work suggests numerous next steps including further translating our prior knowledge of topo- logical properties, like the importance of critical points, into either explicit features or trainable, physics-informed network components (Karniadakis et al., 2021). Another goal is to broaden the scope to encompass a wider range of invariance types, combinatorial comparisons of invariant sets (e.g., via approximation of Conley-Mores graphs), and extensions to stronger notions of equivalence, such as topological conjugacy.
accommodate more versatile invariances and bifurcations, future work could extend our framework to the unsupervised case. For example, given that our data augmen- tations produce samples which are the same dynamical equivalence class, a contrastive loss could be used to shape the feature space so that dynamical classes are represented implicitly during training. Or, while challenging, following Skufca & Bollt (2008), a purely unsupervised approach could be used explicitly construct warpings between systems in order to characterize equivalence classes. Our work demonstrates the potential for scaling up to applications such as identifying dynamical regimes and underlying driving factors across both biological and physical applications.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
We thank Yedid Hoshen for early discussions. This work was supported by the Israeli Council for Higher Education Ph.D. fellowship (N.M)), the Center for Interdisciplinary Data Science Research at the Hebrew University of Jerusalem (N.M),, the Zuckerman Postdoctoral Program (M.R.). Azrieli Foundation Early Career Faculty Fellowship, Alon Fellowship, and the European Union (ERC, DecodeSC, 101040660) (M.N.). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.
10

---

Published as a conference paper at ICLR 2024
REFERENCES
Zin Arai, William Kalies, Hiroshi Kokubu, Konstantin Mischaikow, Hiroe Oka, and Pawel Pilarczyk. A database schema for the analysis of global dynamics of multparaamterr systems. SIAM Journal on Applied Dynamical Systems, 8(3):757-789, 2009.
Aimée Bastidas-Poncee Sophie Tritschler, Leander Dony, Katharina Scheibner, Marta Tarquis-Medina Ciro Salinno, Silvia Schirge, Ingo Burtscher, Anika Böttcher, Fabian Theis, et al. Comprehensive single cell mrna profiling reveals detailed roadmap for pancreatic endocrinogenesis. Development, 146(12):dev 178449, 2019.
Volker Bergen, Marius Lange, Stefan Peidli, F Alexander Wolf, and Fabian J Theis. Generalizing maa velocity to transient cell states through dynamical modeling. Nature biotechnology. 38(12): 1418-1414, 2020.
Volker Bergen, Ruslan A Soldatov, Peter V Kharchenko, and Fabian J Theis. Rna velocity-curees challenges and future perspectives. Molecular systems biology, 17(8)c 10282, 2021.
Petar Bevanda, Johannes Kirmayr, Stefan Sosnowski, and Sandra Hirche. Learning the koopman cigendccoompositonn: A diffcomorphic approach. In 2022 American Control Conference (ACC), pp. 2736-2741. IEEE, 2022.
Yusuf Billah and Robert H Scanlan. Resonance, tacoma narrows bridge failure, and undergraduate physics textbooks. American Journal of Physics, 59(2):118 1991.
Jason J Bramburger, Steven L Brunton, and J Nathan Kutz. Deep learning of conjugate mappings. Physica D: Nonlinear Phenomena, 427:133008, 2021.
Kevin S Brown and James P Sethna. Statistical mechanical approaches to models with many poorly known parameters. Phys. Rev. E Stat. Nonlin. Soft Matter Phys.. 68(2 Pt 1):021904, 2003.
Thomas M Bury, Daniel Dylewsky, Chris Bauch, Madhur Anand, Leon Glass, Alvin Shrier, and Gil Bub. Predicting discrete-time bifurcations with deep learning. arXiv preprint arXiv:2303.096: 2023.
Jack Carr. Applications of centre manifold theory, volume 35. Springer Science & Business Media, 2012.
Guoning Chen, Konstantin Mischaikow, Robert S Laramee, Pawel Pilarczyk, and Eugene Zhang. Vector field editing and periodic orbit extraction using morse decomposition. IEEE Transactions on Visualization and Computer Graphics, 13(4):769 -850 2007a.
Guoning Chen, Konstantin Mischaikow, Robert S Laramee, Pawel Pilarczyk, and Eugene Zhang. Vector field editing and periodic orbit extraction using morse decomposition. JEEE Transactions on Visualization and Computer Graphics, 13(4):769-785, 2007b.
Charles C Conley. Isolated invariant sets and the Morse index. American Mathematical Soc., 1978.
Lu Deng, Wei Wang, and Yang Yu. State-of-the-an art review on the causes and mechanisms of bridge collapse. Journal of Performance of Constructed Facilities, 30(2):0401 2016.
Grigorios Dimitriadis. Introduction to nonlinear aevoelasticity. John Wiley & Sons, 2017.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019.
Michael B Elowitz and Stanislas Leibler. A synthetic oscillatory network of transcriptional regulators. Nature, 403(6767) 2000.
Amin Ghadami and Bogdan I Epureanu. Data-driven prediction in dynamical systems: recent
developments. Philosophical Transactions of the Royal Society A, (,330222))2202122133 2022.
Gennady Gorin, Meichen Fang. Tara Chari, and Lior Pachter. Rna velocity unraveled. PLOS
Computational Biology, 1010992, 2022.

---

Published as a conference paper at ICLR 2024
James Helman and Lanbertus Hesselink. Representation and display of vector field topology in fluid flow data sets. Computer, 22(08):27-36, 1989.
James L Helman and Lambertus Hesselink. Visualizing vector field topology in fluid flows. LEEE Computer Graphics and Applications, 11(3):36-4,, 1991.
David Hilbert. Mathematical problems. In Mathematics: People- Problems- Results, pp. 273-278.
Chapman and Hal//CRC, 2019.
U Iben and C Wagner. Taylor mapping method for solving and learning of dynamic processes. Inverse Problems in Science and Engineering, 29(13):3190-33 2021.
Jonathan Karin, Yonathan Bornfeld, and Mor Nitzan, scPrisma infers, filters and enhances topological signals in single-cell data using spectral template matching. Nature Biotechnology, pp. 1-10, 2023.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics- machine learning. Nature Reviews Physics, 3(6))422-44) 440, 2021.
Farshad Khadivar, Ilaria Lauzana, and Aude Billard. Learning dynamical systems with bifurcations. Robotics and Autonomous Systems, 136:103706 2020.
Mikail Khona and 1la R Fiete. Attractor and integrator networks in the brain. Nature Reviews Neuroscience, 23(12):744-766, 2022.
Y. Kuznetsov. Elements of Applied Bifurcation Theorem. Springer-Verte April 1998.
Jose A Perea and John Harer. Sliding windows and persistence: An application of topological methods to signal analysis. Foundations of Computational Mathematics, 15:799- 838, 2015.
Laurent Potvin-Trottin Nathan D Lord, Glenn Vinnicombe, and Johan Paulsson. Synchronous
long-term oscillations in synthetic gene circuit. Nature, 538(7626):51-5517, 2016.
Zhilin Qu, Gang Hu, Alan Garfinkel, and James N Weiss. Nonlinear and stochastic dynamics in the heart. Physics reports, 543(2):61-162,, 2014.
Katherine N Quinn, Michael Abbott, Mark Transtrum, Benjamin B Machta, and James P Sethna. Information geometry for multiparameter models: New perspectives the origin of simplicity. arXiv preprint arXiv:2111.07176 2021.
Nalini Ravishanker and Renjie Chen. An introduction to persistent homology for time series. Wiley Interdisciplinary Reviews: Computational Statistics, e3(3) 2021.
Matthew Ricci, Noa Moriel, Zoe Piran, and Mor Nitzan. Phase2vec: Dynamical systems embedding with a physics-informed convolutional network. arXiv preprint arXiv:22/2.03887, 2022.
M Sáez, J Briscoe, and DA Rand. Dynamical landscapes of cell fate decisions. Interface focus, 12(4): 20220002, 2022.
Marten Scheffer, Jordi Bascompte, William A Brock, Vietor Brovkin, Stephen R Carpenter, Vasilis Dakos, Hermann Held, Egbert H Van Nes, Max Rietkerk, and George Sugihara. Early-warning signals for critical transitions. Nature, 461(7260):33-59 2009.
Daniel Schwabe, Sara Formichetti, Jan Philipp Junker, Martin Falcke, and Nikolaus Rajewsky, The transcriptome dynamics of single cells during the cell cycle. Molecular systems biology, 16(11); c9946, 2020.
Joseph D Skufca and Erik M Bollt. A concept of homeomorphic defect for defining mostly conjugate dynamical systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 18(1):013118, 2008.
Steve Smale. Mathematical problems for the next century. Mathematics; frontiers and perspectives, pp. 271-294, 2000.
James Johnston Stoker. Nonlinear vibrations in mechanical and electrical systems. Interscience Publishers, 1950.
12

---

Published as a conference paper at ICLR 2024
Steven Strogatz. Nonlinear dynamics and chaos: With applications to physies, biology. chemistry, and engineering. Perseus Books Publishing, 1994.
Itay Tirosh, Benjamin Izar, Sanjay M Prakadan, Marc H Wadsworth, Daniel Treacy, John Trombetta, Asaf Rotem, Christopher Rodman, Christine Lian, George Murphy, al. Dissecting the multi- cellular ecosystem of metastatic melanoma by single-cell ma-seeq Science, (52(6282):189-193 2016.
Mark Transtrum and Peng Qiu. Model reduction by manifold boundaries. Phys. Rev. Lett., 113(9): 098701, 2014.
Sarah Tymochko, Elizabeth Munch, and Firas A Khasawneh. Using zigzag persistent homology to detect hopf bifurcations in dynamical systems. Algorithms, 13(11):278, 2020.
Noelia Vallez, Gloria Bueno, Oscar Deniz, and Saul Blanco. Diffeomorphic transforms for data augmentation of highly variable shape and texture objects. Computer Methods and Programs in Biomedicine, 219:106775, 2022.
Anael Verdugo. Hopf bifurcation analysis of the repressilator model. American Journal of Compula- tional Mathematics, 8(02))137-152, 2018.
Conrad Hal Waddington. The strategy of the genes. Routledge, 2014.
Wen-Xu Wang, Rui Yang, Ying-Cheng Lai, Vassilios Kovanis, and Celso Grebogi. Predicting catastrophes nonlinear dynamical systems by compressive sensing. Physical review letters, 106 (15):154101, 2011.
Wen-Xu Wang, Ying-Cheng Lai, and Celso Grebogi. Data based identification and prediction of nonlinear and complex dynamical systems. Physics Reports, 644:1-76, 2016.
Terry A Weisshaar. Acroclasticity, an introduction to fundamental problems with an historical perspective, examples and homework problems, 2012.
Kurt Wiesenfeld. Noisy precursors of nonlinear instabilities. Journal of Statistical Physics, 38: 1071-1097, 1985.
C. Wissel. A universal law of the characteristic return time near thresholds. Oecologia, 65:101-107, 1984.
Lingyun Xiong and Alan Garfinkel. Are physiological oscillations physiological? The Journal of Physiology, 2023.
Han Zhang, Ian Goodfellow, Dimitris Metaxas, and Augustus Odena. Self-attention generative adversarial networks. In International conference on machine learning, pp. 7354-7363, PMLR, 2019.
Jian-Gang Zhang, Xian-Feng Li, Yan-Dong Chu, Jian-Ning Yu, and Ving-Xiang Chang. Hopf bifurcations, Iyapunov exponents and control of chaos for class of centrifugal flywheel governor system. Chaos, Solitons & Fractals, 39(5):2150-25) 2009.
13

---

Published as a conference paper at ICLR 2024
7 APPENDIX
7.1 DEFINITIONS
7.1.1
DYNAMICAL SYSTEMS INVARIANT SETS
This paper concerned with dynamical systems which can identified (potentially implicitly) with woo-dimensional systems of ordinary differentiall equations. For more general treatment, see Strogatz (1994); Kuznetsov (1998).
Let (z, 3) E R and
dz =
dr dy
f(,tt)
(1)
dt
dt
dt
z(0) 20
be initial value problem for assumed to be smooth. Definition (Continuous-time real dynamical system). The map & : R.2 x [0, - IR? given by
Φ<z0: T) =(0) + [' f(=(t).t) dt
(2)
is called a (continuous time real) dynamical system or flow. In this context, R2 is called the phase space, z(t) a particular state of at and f(2(t), t) is called the vector field associated to 6. For fixed 2g, the map 020- TR?2 TR.2 given by particular solution to Eq. is called trajectory through 20) = 20-
Hereafter, we assume f(2(1),)) f(=(t), in which case we call ø autonomous. Otherwise, would be called non-autonomous.
The study of dynamical systems characterized by its focus on long-term behavior of such systems. In our case, two long-term behaviors are of interest. First, if there exists 2 E R2 such that f(z") O. then we say = a fixed point of Φ) particular, o(=", T) = for all T.
The second long-term behavior of interest in our context is the limit cycle.
Definition (Limit cycle). Suppose there exist 2q and such that 2z (T) = 20. Then, the set
7 = E 31 s.t. (t) =
is closed, simple curve (simple since Eq. I always has unique solution). for each E E.yyhh there exists an e such that the ball (z(z) (z) of radius e around z is intersected by no other closed, simple curve, then we say isolated. A trajectory of & forming isolated, simple, closed curve, 7, is called a limit cycle.
Zero and one-dimensional structures like fixed points and limit cycles are examples of "invariant sets" and serve to distinguish different types of long-term behaviors in dynamical systems. Further distinctions can be made by considering the stability of these structures: informally, how the system behaves when it is slightly perturbed away from fixed point or limit cycle. In the following definition, set S will denote either a singleton fixed point or the curve in R coinciding with a limit cycle. For any let
d(z,S) inf - 2'112 Vz' E S
(4)
Definition 3 (Asymptotic stability). Let S be invariant given dynamical system, 6. If there exists a > such that when d(z,S) < 5 and d( t), - 0, then say S is (asymptotically) stable. Otherwise, we say is (asymptotically) unstable.
defined above with respect
Informally, S is stable when small perturbations away from S eventually converge onto S.
7.1.2 EQUIVALENCE OF DYNAMICAL. SYSTEMS
Fixed points and limit cycles are invariant sets not just because a solution starts in or enters the set it stays within it, but also because those sets preserved under the action of certain transformations
14

---

Published as a conference paper at ICLR 2024
to the whole system. Informally, these transformations continuously and smoothly warp the phase space in way that preserves fixed points, limit cycles and their stability properties. Systems which differ only by such transformation homeomorphism is continuous bijection between two topological spaces with a continuous inverse. diffeomorphism, similarly, is differentiable bijection with differentiable inverse. We define several notions of equivalence:
considered equivalent. More specifically, recall that
Definition 4. Let 6(x, t) and t) be dynamical systems (flows) on I2 Then,
1. Topological conjugacy. If there exists homemorphism, h RR? - IR?, such that (h(y),,) t) = h()(y,t))
(5)
for all y, then ④ and called "topologically conjugate."
2. Topological equivalence. Let O denote trajectory of Φ. Suppose there exists a homeomor- phism IR2 IR2 such that
h(O +)) (hoot :tE R] (6(h(v). =
(6)
and for each there is å > O such that if0 < |s| < t < ó and there exists an s such that (h((y),s) (y t), then > In other words, h maps orbits orbits and preserves the order of time. In this case, we say & and ab are "topologically equivalent".
3. Dynamical equivalence. Suppose the flows & and ab have associated vector fields f and g.
If there exists homeomorphism h I R² such that
dy = g 5
then we say that the flows are " "dynamically equivalent."
of these three notions, then there exists one-to- one correspondence between their invariant sets which may or may not preserve their stability properties. Note that throughout this paper, we use diffeomorphisms, which stronger notion than homeomorphism. In that sense, there are at least as many invariances implied in the diffeomorphic case as the homeomorphic one. In particular, for the notion of "dynamical equivalence", the one we use throughout this paper for its simplicity, we have the following claim:
If flows equivalent according
Claim:
Let and Y be two topological spaces and let that is continuous, bijective, and has a continuous inverse. Let à f(x) be dynamical system on X and define ý = g(y) = [(h (w)). Then, (1) f has fixed point, then so does y and (2) if f has periodic orbit, then does g.
Y be a homeomorphism, a function
Proof. If f has a fixed point at on then f(xo) = 0. In that case, g has fixed point at y0 = h(x))) since
(
f(xo)
0
Next, let (t) be solution to f which forms a periodic orbit, y, passing through xp =0)) y is simple, closed curve and therefore is its image, h(T),, since f is homeomorphism. However, we must still show that ((()) is periodic orbit of Since is periodic orbit of f. there exists 7 such that
x(T) x(0) + S f(x(t)) de
= x0))
15

---

Published as a conference paper at ICLR 2024
so that To (x(()) dt 0. Note that rrt) h hy((t) for all t. Thus,
v(r) u(0) + J s(v(1)) di
- v(0) + ! f(h s(h="(((t)) dt
- v( + L. f(a(t))
(0),
which shows that h(y) a periodic orbit of g with the same period, T.
Thus, two ODEs f(z and g(y) the plane "equivalent" (f g via h 1) if there exists homeomorphism h E H such that g(y) = ((hh((w) (*). The relation ~ an equivalence relation partitioning the space of ODEs into dynamical equivalence classes. In particular, this relation is symmetric. For, if holds, then g((()) = f(h- (h(x))) = f(x). Le. g ~ via
Since all flows within dynamical equivalence class share the same features that we care about in this study, if suffices often to work with a single example which has a simple functional form. In that case, refer that example as prototypical system.
7.1.3 BIFURCATIONS
Next, suppose that the ODEs da - (=; 0) associated family of flows is parameterized by 0 E Q which we assume to be one dimensional for simplicity. Then each homeomorphism partitions the family into dynamical equivalence classes corresponding to partition of e. We have the following definition:
Definition 5. Let (= E DEOCR,
(R} be parameterized familyy of ODEs corresponding family of flows dg. Suppose is a homeomorphism defining an equivalence relation on F. If there exists 8" such that, for all e >> 0, there are 0 such that ||0 0* < with fe + fas, then we say there is bifurcation ar
Informally, a bifurcation occurs when smooth changes in system parameter result in qualitative changes in dynamics. In this study, we focus classifying vector fields according to whether they belong to one of two equivalence classes family of flows exhibiting a supercritical Hopf bifurcation (definition in main text).
7.2 METHODOLOGY
7.2.1 AUGMENTATIONS
Topological data augmentation carried out by warping input vector fields with bounded, mono- tonic aational-quadraan splines proposed Durkan et al. (2019). The spline mapping h zd + yd for d-dimensional input, where equals the number of dimensions (in our case, d 2) is based on set of hyperparameters, including a boundary parameter B and the number of bins K. The mapping function h transforms each coordinate within the boundary range [-B, B] using a spline of function.
K +++ 1 knots with each bin in between following a random monotonic aational-quadren Coordinates outside the range retain the identity function (i.c., unchanged coordinates), see Durkan et al. (2019), Figure 1. In our experiments, we set B to 4 and use 5 bins (K = 5). Given these hyperparameters, exact parameters defining the splines (width, height and derivatives the knots) are retrieved from normalizing flow. Since we use random initialization values, this is equivalent to sampling these parameters directly, without using the normalizing flow architecture. Given a random diffeomorphism, denoted by h, and the system dynamics represented by f. we invert h to obtain coordinates X = h (Y) from the 64-by-64 grid coordinates Y. We then apply the dynamics f to the transformed coordinates X and consider these as the dynamics grid coordinates Y, that is g(Y) = f(). For our augmentations, we adopt an implementation by Andrej Karpathy available at https: / / thub. com /a rewttyaaahhhee http:/yytttoheoco werecch--rrraa norma 1 zizzz miing-2 flows
16

---

Published as a conference paper at ICLR 2024
After applying the augmentation, verified that the fixed point remained in-frame. This is proxy to determine if the attractor is still within our field of view. We refer to the dataset obtained by transforming the simple oscillator ("SO") with this augmentation the "Augmented SO" dataset.
7.2.2 ARCHITECTURE AND TRAINING
For vector field X = (æ, 3) of dimensions [2, 64, 64], we first compute its angular representation by replacing each vector with arctan (j/i), which returns a signed representation of the angle in radians.
Our architecture is composed convolutional attention feature extraction module based on the SAGAN discriminator (Zhang et al., 2019), coupled with multi-layered perception (Fig. A5).
The feature extraction module consists of four convolutional blocks. each utilizing kernel size of 3 x 3 and stride of 2 2. The number of channels in each layer is progressively increased to 64, 1 128, 256, and 512. As in SAGAN, we applied spectral normalization and LeakyReL.U non-linearities to all convolutional blocks and incorporated the self-attention mechanism into the last two convolutional layers.
The output of the convolutional layers was then fed into single-layer MLP (multi-layer perceptron) to map the activations to a 10-dimensional space, which we refer to as the "feature activations". The MLP consists of a hidden layer with 64 dimensions, ReLU activation, and dropout rate of 0.9. Using linear mapping, we predicted the logits for both point and periodic attractors from the feature activations.
We formulated the problem as a multi-class classification task, as there can simultaneous presence of point and cycle attractors. We used binary cross-entropy loss.
For training and evaluation, we prepared a dataset comprising 10,000 samples for training and 1,000 samples for testing. We employed the ADAM optimizer with learning rate of 1 10-1 and trained the model for 20 epochs, which we confirmed was enough time to the training data. For inference, we still perform dropout and average the output of 10 evaluations for each sample.
All experiments were carried out using pytorch v. 1.12 using an NVIDIA RTX 2080 GPU, taking 4-0 seconds per training experiment.
lolet
Sirlle ossilator vector de
Augmented
Angles
Corw SpectralNorm
Com Spectraloomm LeakyRelu
Carr Spect-alMorm Leakyficuu Attention
Corw SpectratNorrs Leacytelk Attention
Unear
Linaar
Linear
Rela
LeokyRell
Dropout
Figure Framework. Our woint-periodes classifier begins by augmenting the vector field converting it to an angular representation. Data is then passed through convolutional-ateniiio layers and an MLP classifier.
7.3 SYNTHETIC DATASETS
All input vector fields were computed by measuring the continuous system an evenly-spaced grid of 64 x 64 points. In Table A3 the governing equations of each system are described in terms of either polar coordinates (r, 8) cartesian coordinates (æ1, 222) in addition to the range of phase space examined, the parameter ranges we sample from, and bifurcation boundaries described Strogatz (1994). Examples of pre- and post-Hopf instances plotted in Fig. A6.
17

---

Published as a conference paper at ICLR 2024
Table A3: Descriptions of systems undergoing Hopf bifurcation. (*) Selkov periodic attractor condition: (1/2(1 2a vI 8a)) < (1/2(1 2a + vT -T 1/
Narte
Equation
Phase space
Paromeeerr ranges
Perbal altracie condition
Simple Ouciltttc
Ple
11
#20
pay
21.
M.
BZ Rusciion
80200
(2 A
Selkoe
Poytt 10.3
Saroid
Figure A6: Pre- and post- post-Hopf examples. Vector fields visualizations as streamplots of random instances of each system pre- and post-Hopf bifurcation (shown top and bottom panels respectively).
To assess the parametric proximity each system its bifurcation boundary, we calculate the parameter distance the minimum L2 distance from the bifurcation curve. In the left portion of Fig. A7, we present these signed distances (negative, or positive, distances indicate pre-Hopf, or post-Hopf, systems respectively) for the nonlinear bifurcation curves of the Selkov system. While parameters fully determine the bifurcation behavior, inferring this behavior from vector fields and trajectories of systems near bifurcations can be misleading. We showcase here instances of the Selkov systems, each depicted by its vector field and trajectories starting near and away from the fixed point. We see that for high values, despite having point attractor, trajectories exhibit periodic-like like behaviors, which can be mislcadingly similar to the behavior of a periodic attractor.
Seltiov
Point, 0.04,
Point, 08,
Point, 0.04, b
Point:
1.0
0.1
0.0
-0.2
0.04
0.08
Figure A7: Trajectories of pre-Hopf Selkov instances. examine instances Selkov system. show the parameter space (left), its true boundary (black line) and the four samples marked by 'x'. A score (color) is assigned to each system based on how close it is to a bifurcation (minimal distance from the bifurcation curve), and whether it pre- post-Hopf bifurcation (of negative positive scores, respectively). Colors scaled between value. On the right, we show vector fields with flow F' T") where zg is away (blue), or near (red) the fixed point. We set T = 100 and interpolation is over timesteps of 0.1.
maximum absolute value (red suggesting a cycle) and negation (blue) with white zero
We tested proposed method noisy data which was subjected to additive zero-centered Gaussian noise (μt = 0) and scaling o 0.1 (see Fig. A8 and Table 1) and on increasing scaling of for Augmented SO (see Table 2).
18

---

Published as a conference paper at ICLR 2024
LIITATE PLSMOOTM
Lamara Segmoid
Vin dur PM
BZ. Seuetior
Saftew
Figure A8: Perturbed pre- and post- Hopf examples. Gaussian noise of 0.1 added to vector fields of pre- and post-Hopl examples (shown in top and bottom panels respectively) of each system.
7.4 BASELINES
The advantage of learning topology-aware representation is emphasized by comparison to other across-system representations, namely: (1) "Critical Points", (2)"L-yapnoo"", (3)Phaseevee"" (4) "Autoencoder", and (5) Parameters"".
Critical Points. Since we focus here on Hopf bifurcations, a point and a periodic attractor can be distinguished by the identity of the fixed point, being an attractor for the former and repeller for the latter. Helman & Hesselink (1989; 1991) devised an algorithm for detecting and characterizing such critical points from an autoencoder representation. later implemented in https: : ////tttthb/// com/ zaman13 VVec or-cctor-Fae Tod-Topolggyyy 2D. Running this method in practice, we find that for a point attractor, a single point is often detected but for periodic attractor, the algorithm mistakenly outputs multiple coordinates as fixed points besides its repeller (mean number of points for point attractor is 1.00 and for periodic attractor 1.50). We therefore adjust the rule to classify system as periodic attractor if some identified critical point is of repelling type. Using this measure, we find that the classification of the synthetic systems is nearly perfect (see Table A4), but breaks immediately with noise; with Gaussian noise of a = 0.1 added to the Augmented SO dataset, accuracy drops to 59%.
Lyapunov. Given time-series data, the maximal Lyapunov exponent, Å1. captures the greatest separation of infinitesimally close trajectories across an n-dimensional space. Systems with A < are considered in state of fixed point, whereas periodic states have À = 0 (Stoker, 1950). Zhang ct al. (2009) uses this measure in the context of Hopf bifurcations towards the design of rotational machines. To evaluate this measure, we integrate each system from random point at an interval of 0.1 up to time T = 100 based on the vector field and compute the maximal Lyapunov exponent, λ of the trajectory using the python library https: :/ /pypi... orq/ pro ject holds. We find the optimal threshold of A₁ for classification based on Augmented SO ROC curve and use across all systems (see Table A4).
Phase2vec. For Phase2vec we follow the tutorial of data generation and training as described in Ricci al. (2022) and implemented in https: ////it thub. com/no ni ann ab/ phwssevec.. That is, we generate 10K train samples of sparse polynomial equations and train the phase2vec framework consisting of a convolutional encoder mapped by single linear layer into a d 100-dimensional embedding space and consequently decoded via another linear layer into predicted parameter values which are multiplied by the polynomial library to reconstruct the vector fields. We then train linear classifier of pre- and post- post-Hopf classes using the same Augmented SO dataset as used in our current framework and select for optimized hyperparameters (see Table A4).
For Autoencoder we deploy an autoencoder architecture consisting of an encoder of three convolutional blocks (kernel size of 3 x 3, a stride of 2 2, and 128 channels each layer), followed by a linear layer into 10 latent dimensions, and a decoder mirroring the encoder architecture. The autoencoder is trained on the same data used in phase2vec with ADAM optimizer with learning rate of x 10 for 20 epochs. We then train linear classifier of pre- mrof and post-Hopf classes using the same Augmented SO dataset as used in our current framework and select for optimized hyperparameters (see Table A4).
Autoencoder.
19

---

Published as a conference paper at ICLR 2024
Parameters. For Parameters, vector of 20 coefficients, (Cu, C1; ...,.. Cy. du, d1, dy), is fitted using
least squares to reconstruct the two-dimensional vector field (i: ÿ) with (q ylc), q3(2, y d)) where Ø3 g/))) chose polynomial representation of degree 3 following earlier work (Ihen & Wagner (2021)) which showed this value was sufficient for dynamical systems similar those we investigate here. We then train& a linear classifier of pre-Hopf and post-Hopf classes using the same Augmented SO dataset as used in our current framework and select for optimized hyperparameters (see Table A4).
co + c C2l + czz² + cay² + c52y + cg." d + cTy² + cgz y + Caryy. We
Table A4:
Test accuracy of baselines and their hyperparameter tuning, default: leamning rate=le-4, num epochs 20.
S0
Aupmented SO
Superenticall Hapf
Lienand Poly
Lionsed Signoid
Vin der Pel
BZ Resction
Sel'or
Avcrago
Dur Model
0,97=60.1
0.93+0001
0.98-0001
0.85= 14
0.92+0.11
0.4+0 14
0.54+0..0
0:65+003
0.87
leaming rale- 3e
0,9510009
0,9110.09
0.96..0.00
0,81
0.9510.09
0.9110.1
0.53110.12
0.6100.77
0.85
leaming 1e-5
0.97=:0.01
0.96+00.02
099= 01
0.58+001
0.7920.13
0.6440114
0.53+10.0
0.53+007
10
0:57=:0:01
0.92:0.31
0998=0001
0.88-0.112
0.8620013
033360115
0.74+:1.13
0.632002
= 40
0.96::0.11
0,91-00.01
0.98:0.01
0.86
093-00.06
0.5500.14
0.81+00.11
0.670006
Criccal Poims
0.33
1.00
0.70
0.92
[yspunov
0.85
0.71
0.89
0.89
0.87
0.84
0.78
0.68
0.81
Phoselves
0.8910:06
0,8610.02
0.1100002
0.5110012
0.68.10.0
0510000
05100
0.491000
0.64
leaming rate 1e-3
0,77-0.00
0:5==0022
0.61::0.16
0.740.14
8.7650177
0.51
0.49:002
0.72
lominng ratn le-5
0.85=-001
0.73+0.01
0.9-0001
0.4+0000
a 0.56+0007
048+00
0.51+0004
0.4++00
0.63
mum qpeche- - 10
0.80
0.82+0.00
091-0.02
0.4+000
0.46-0.00
048+00
0.500.0
0.4++00
0.62
mmm epochs- 40
0,9210.04
0,87
0.910011
0:4910.07
0,510005
0.5100.66
0.510.0
0.491000
0.65
Automnender
0,93-0000
0.86+0.01
0.960.011
0.65-0.15
092+0008
035+0.11
0.83+0113
0.5+000
0.81
leaming rate 1e:3
0,9420000
0,9210..01
0.0000001
0.63
0,9610004
0810.17
0.611..15
066101
0.81
leaming rate = 1e-5
0.912:004
0.65::0.00
09::00011
0:47:00:22
08:00
0.5420.11
0.53:10.11
0.2:00009
0.67
epochs 10
0,930003
0,8210.02
015500001
0.52::0.12
088110112
0.610.17
0.7110.34
0.491000
spocha 40
0.95m00.02
0.000.01
096
0,53-0.01
09240.08
0.+00.3
0.00-0133
0.5+0009
Parameters
0,53=0:05
0.50+0.01
0:56-0006
0.52-0.11
0.53-6.1
0.54+-012
0.53+011
025200077
leaming rate- 3
0.510.0
0,5710.01
0,47100
0,4810.0
0.10000
0.69100.
leaming rate- 1e-5
0.5=0:08
0.5+0:0
0.5e0.15
0.51=0.12
0.49=0.18
0.50.33
0.51+0012
0.51:0088
mmm epochs = 10
0:53:0:07
0.5+0.03
0:5-0007
0.51-0.12 12
0.51-0.2
0.535:0155
0:49::0.07
0.::0044
0.51
num spoche = 40
0.54-0:07
0.56+-0.0
0.57=0.08
0.52+6.1
0.55+0.14
0.55+0013
0.54+-0.11
0.5+0.05
0.54
7.5 ABLATIONS OF OUR MODEL
To evaluate the impact of our architecture, we compare the results with ablated versions of our full model, namely: (1) "No Attention" - where we exclude the two self-attention layers, resulting in network that purely feedforward, (2) "From Vectors" our model but trained on raw vector fields
instead of transformed angles as inputs, i.e. with two channels in the first convolutional layer, (3) "No Augmentations" our model but trained on the unaugmented SO data, and (4) CENN-baseline". a model including all ablations, see Table A5,
Table A5: Test accuracy on ablated models.
Augmented 5O
Supercitical Hopf
Liénard Poly
Liénard Sigmoid
Van Pol
Resction
Sel'kow
Average
Model
0.974-0:01
0.93=00011
0.98=00011
0.85e:
0.92=0.1
0.4400.14
01S4:-1.99
0,65+0.03
No Atention
0.97s0.01 0]
0.92s0.01
0.98a0.01
06340.11 13
0.66800 12
0.5t.007
0.8410.99
0.6250.05
0.77
From Vectoes
0.57200.12
0.9220011
0.982/001
0.6211.12
0.9620003
0,9329.05
062290.15
0.5220.05
0.82
No Augmentation
0L99=110
0.4500011
0.99:6.11
0.94+1U0
0.96sfLB
0.88-001
0.540.01
0.52+0.05
0.7%
CNN baseline
0 0.96.0.01
0.55.0001
0.99:0001
0.5910.13
0,7760.14
0.5120.1
0.4550.09
0.4910.0
0.68
We observe good performance when training on classical SO for systems that resemble the SO in having their fixed point located persistently at the center. In fact, as we increase the boundary parameter B of the neural spline augmentations, which directly affects the location distribution of the fixed point, we observe the change in test accuracy for two systems: the BZ reaction system, which exhibits a peripheral fixed point (with Φ ranging between 0.6 and 3.8 within the |0, 10j limits): and the Liénard Polynomial system, which has centrally located fixed point, see Fig. A9 and Table A6.
Furthermore, we examine the rules learned by the attention layers. We sum the values of the first attention mask and expand it to the original grid dimensions. Through this analysis, we observe that the attention is often sparse, focusing on the fixed point or its immediate neighbors, as depicted in Fig. A10. This observation raises intriguing follow-up questions like whether higher accuracy correlates with specific patterns in the attention masks.
Data size and resolution We assayed the average accuracy function of data size, showing that a substantially smaller train data size (2k instead of 10k) is sufficient in most cases, see Table A7. In the case of higher sampling resolution in space, our method is directly applicable. Indeed, for
20

---

Published as a conference paper at ICLR 2024
B
2
3
4
5
Reaction
0.42
0.65
0.79
0.88
U.71
Liénard Poly
0.92
0.90
0.86
0.79
0.89
Figure Table A6: Effect B. boundary parameter in Neural Spline Flows, on fixed point coordinate distribution (see Figure, dashed black lines mark the train data boundarics) and on test accuracy for BZ Reaction and Liénard Polynomial (see Table).
cyde
cycle
DOITT
pohi
jain.
Figure A10: Attention learned. Mask of the first attention layer is diluted back to the original grid dimensions for random constant) point shown grayscale (standardized by minimum and maximum values black and white respectively).
the current resolution (64x64), training was quick, 40s, while training twice the resolution resulted training time of 130s and the same test accuracy (87% average accuracy).
Table A7: Average test accuracy with different train sizes.
Train data size (k)
2
4
8
16
Average accuracy
0.84
0.87
0.89
Other architectural choices. We observe robust performance upon editing our network architecture, Table A8.
21

---

Published as a conference paper at ICLR 2024
Table A8: Test accuracy of alternative representations varying neural network architecture parameters.
SD
Augmented SO
Supersiiical Hopf
Limnard Poly
Lienan Sigmoid
Vin der Pol
BZ Reaction
Self knw
Avenge
Our Model
0:572001
0:38.0011
0.85211.14
0.92sll11
0.545).11 14
0.87
drapeut = 0.5
9K 01
0:930.:01
0.98e0.01
0.8440.1 17
0.9300.06
0,800
0.76+0. 12
0.68-0.06
0.86
kemel size 5
0.5720001
0.9220:11
0.98.0.01
0.9228).05
0.9220.00
0.8220.07
0.72+0.111
0.60005
0.87
latest dimeusion
0.98
0.33:0.01
0.9920.01
0388+) 12
0.3541)112
0.79ell.11
0.5=0004
0.87
latent dimension
0.9740.01
0:93:0.01
0.98.0.01
0.410..1 15
0.9130.09
0.8740.05
0,8240.1 11
0.65:0.05
0.87
channels
0.58:2001
0:5220.01
0,9920001
0.5629.113
0.9120.05
0.8720.09
0.76gt1.11
0.67:0.66
0,9720001
0.92:0.01
0,9820.01
0.3820.12
0.9260.1
0.3620.14
0.7960.11
0.6660.04
7.6 RECOVERING BIFURCATION BOUNDARIES
In Fig. All we plot for each system model's confidence function of its parametric proximity to its nearest bifurcation point. In Fig. A12 we depict the bifurcation diagrams for the alternative trained baselines.
Augmentad 0.54)
Supercritical [0.50, 0.42)
Linnard
Lienard Sigmald (0.87, 0.27)
Reartion (0.21, 0.68)
Selloov (0.17.
00
Distance from bifurcation
Low confidence near bifurcation boundaries. For each system, we plot mean (line) and standard deviation (shaded region) of the confidence over the 50 training runs as a function of the parametric distance from the bifurcation boundary for pre-Hopf (blue) and post- post-Hopf (red) systems with the respective Spearman's 's correlations noted in the title.
Figure A11:
Auguercad
Suneerrilea Hopf
Lirrard Paly
Lierard Signoidd
B2 Rieetion
-5.8
H0.8
+0.5
0.0
0.
Augmented
Supercntical
Unnand Poly
Lieaad Sigmald
Reaction
0.8
08
DE
-0.8
0.4
03
0A
044
055
04
0.5
0.0
03
0.5
0.
58
0.5
05
08
0.04
0.28
50
Augmented 50
Superentical Hop
Uecard Faiy
Ueeard Sigmood
BZ Seartonn
Sel'aov
10
5.0
09
0.6
6.4
16
30
-2.8
0.0
0.0
(c)
Figure A12: Inference of bifurcation boundary for alternative representations. As Fig. 2, each system plotted at coordinates of its simulated parameter values and colored according to its average cycle prediction across 50 training re-runs based on (a) Phase2vec, (b) Autoencoder, and (c) Parameters. The true parametric bifurcation boundary is plotted in black. Colors scaled between perfect cycle classification (red) and perfect point classification (blue) with white evenly split classification.
7.7 CO-OCCURRENCE OF ATTRACTORS
There are systems where both point and periodic attractors are present, for example in Suberitical Hopf bifurcation (Strogatz, 1994):
22

---

Published as a conference paper at ICLR 2024
f = por + p rb,
i + br
with point attractor where μ < -0.25, point and periodic attractor where -0.25 < la < 0, and periodic attractor where 1 0. The average (across runs) point logit our framework outputs for each regime is, respectively: 1.32, -6.72, -1.03.
7.8 REPRESSILATOR SYSTEM
Gene regulation of the repressilator system is formulated as set of six equations, or three sets of the following two equation (Elowitz & Leibler, 2000):
n4 - m, +
+
1 P3
1
p; -p;; m,)
where m,, or Pi describe the mRNA, or protein concentration, of gene which is one of Lacl, TetR and À phage cl genes. The mRNA expression gene is inhibited by the protein levels of gene j for the corresponding gene pairs of = LacI, TetR, cI and j cI, Lacl, TetR. ag is the rate of leaky mRNA expression (that is, transcription taking place despite being entirely repressed), a is transcription rate without repression, is the ratio of protein and mRNA degradation rates, and is the Hill coefficient.
We set the parameters ag 0.2, 2. and consider trajectories of cells of varying transcription rate a E (0, 30), and of varying ratio of mRNA and protein degradation 3 E (0, 10). A supercritical Hopf bifurcation occurs across both boundaries,
3.43 - 442--A-
- 8
AV9A
24AA 48
B1
4A+8
4A+8
3.4²
- 8
Av
24A 48
82
4A+8
4A
where A -amp(a-) (3) (1+pr and p
+ (Verdugo, 2018).
1+
We plot these boundaries in Fig. 3 and Fig. A13. The discontinuity between B1 and B2 in the plots is due to our discretized computation of the boundary; we take a grid of values and solve for their upper (31) and lower (82) boundaries.
To generate a vector field we simulate a trajectory of 50 timesteps sampled at 0.1 intervals starting at (m LacI PLacl, MT PTetR. Mel, P.1) (2.11, 2.28, 1.57, 1.71, 1.07, 1.14), see examples in Fig. A13. We sample 100 cells from this trajectory. add Gaussian noise a = 0.5 and evaluate the six-dimensional vectors of these cell states. In Fig. A14a show examples of these cell vectors as they are projected onto the TetR-Lacl protein plane together with their theoretical pre-Hopf or post-Hopf categorization. We then interpolate these vectors 64x64 grid (Fig. A14b), and follow-up with angle transformation (Fig. A14c).
7.9 PANCREAS DATASET
obtain the data, compute high-dimensional gene expression velocities per cell and their projection into UMAP space, and each cell's cell-cycle score described in sc Velo (Bergen et al., 2020), RNA velocity basics tutorial. Since cell types present distinct dynamics (Fig. A15), in order to enrich our dataset, partition cells into batches of at least 50 cells of common cell type (Fig. A16a). We then generate samples for our classifier by interpolating their sparse two-dimensional velocities to a 64-by-64 grid (Fig. A A16b), which are later transformed into angle representations (Fig. A A16e). From
23

---

Published as a conference paper at ICLR 2024
Repressilator
Periodic,
12.
-5
0
-10
Figure A13: Examples of repressilator trajectories. examine four cellular trajectories repressilator system. We show the parameter space (left), its true boundary (black line) and the four samples marked as 'x'. score (color) is assigned to each system hased on how close it is bifurcation (minimal distance from the bifurcation eurve), and whether is pre- post-Hopf bifurcation (of negative positive scores, respectively). Colors are scaled between the maximum absolute value (red suggesting cycle) and negation (blue) with white at zero value. On the right, we show complete cellular trajectories with T 50 and intervals of 0.1.
these, predict using classifier whether cells demonstrate cyclic dynamics, where in this case, such dynamics are characteristic of the cell cycle. To test our evaluation, we aggregate individual cell cycle scores into population score the fraction of cells that undergoing S-phase, where score threshold of 0.4 is employed (Fig. A16d).
24

---

Published as a conference paper at ICLR 2024
Periodic
Point
Periodic
Per odi.
(a)
Periodic
Point
Parledic
Per arddi
MERSIN VNI
Penboic
Poitt
Penodic
Perodi
Pertodic
Perlod::
Periodic
(c)
Figure A14: Repressilator simulated samples. We show here the data generation process for the repressilator system eight random instances (distinct a, values). We simulate noised trajectories of 100 cells and project their velocities onto the TetR and Lacl protein plane (a). From these, we interpolate vector values to 64-by-64 grid (b) and image their corresponding angles (c). Titles indicate the theoretical pre- and post Hopt-Hopf classification of each sample.
25

---

Published as a conference paper at ICLR 2024
Ductal
Ngn3 low EP
Ngn3 high EP
Prevendocring
Alpha
Epsilon
Figure A15: Cell type velocities, Depicting the velocitics UMAP representation for cach cell type separately.
26

---

Published as a conference paper at ICLR 2024
Prev aredee
Dustel
Alpra
Dacial
Honn
(a)
Non3 high
44hha
Ducaal
Ducal
(b)
Ihe endooro
Ductal
Apph
Ducal
Non'
Pre endiorrne
Ductal
Aphha
Ductal
acon
Non3 Nigh
Ducal
(d)
Figure A16: Pancreas data samples. We randomly split cells of common cell type into groups of 50 or more cells. Based on their velocities (a), we compute interpolated vectors (b) and their angles (c). We compute proxy for the cell-cycle signal in each sample. which we consider the fraction of cells undergoing S-phase according their S-score value (d).
27


Reviews:
[{"review_id": "OEfQDC8KDL", "reviewer": "Reviewer_pimK", "summary": "The paper proposes a method for classifying dynamical regimes in a data driven deep learning framework, combining a self-attention GAN for feature extraction and training an MLP for classification. \nTheir main goal is to extract invariant topological properties from dynamical systems around bifurcations, and use them to classify dynamical regimes. To this end they leverage a data augmentation method, based on the principle of topological/dynamical equivalence. The first pre-train their model on variants of a simple oscillator system featuring a supercritical Hopf bifurcations.\nThey apply their pre-trained model to classify two regimes (pre-Hopf and post-Hopf) in several two-dimensional, autonomous systems (both synthetic and experimental data), featuring the same bifurcation. To this end they show that they can extend knowledge of dynamical classes (distinguishing between fixed point and cycles) across different datasets, and determine bifurcation boundaries and closeness to bifurcations.", "questions": "Your examples indicate that augmentation actually decreases performance for the simpler baseline models, which goes a little counter to the argument of your paper. Why do you think this is the case?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The ideas in the paper are clearly presented. The idea of detecting bifurcations in dynamical systems is timely and important for applications. Their framework for data augmentation, based on the principle of dynamical/topological equivalence, is interesting, and the idea of transfer learning across different dynamical systems that share common topological properties seems like an important application. I also appreciate the application to experimental cell data, and the ablation studies in the supplement, investigating the effect of data augmentation and different modules, e.g. the attention modules, on outcomes.", "weaknesses": "Related Work/Dynamical Systems Theory:\n\nOne of my main concerns is that there is little reference to the rich mathematical body of work in dynamical systems theory anywhere in the paper. The concept of \"dynamical equivalence\" proposed in the supplement is well known in the literature (see e.g. Kutznetsov, \"Elements of Applied Bifurcation Theory), as are the proofs given in 6.1. in the context of the concepts of topological conjugacy&topological equivalence. However, both terms are only mentioned once in passing in the whole paper, and as far as I could see are not referenced anywhere. The term \"topological invariants\", while used in the abstract and title, is also never formally introduced nor explained in the rest of the paper.\n\nMy second main concern is that I don't understand your reasoning why you \"evaluate as baselines existing vector field representations, none of which explicitly encourage topological invariance\". Topological features (such as fixed points, cycles, stability of fixed points etc.) can also be numerically approximated directly from the flow fields, especially for simple 2-D models you investigate here where the vector field is already present (i.e. fixed points are simply locations where the flow field is zero). There are other topological properties (e.g. persistent homologies, Betti numbers) that one could also extract to aid with classification.  At least using this as a baseline for comparison would seem crucial in the context of this paper? If I'm not completely mistaken this should at least work better than random guessing, as the comparison methods on the Selkov model effectively do. Fitting polynomial/other representations to the vector fields that a priori are not tailored to extract a useful representation for topological invariants seems perhaps unsurprisingly not so effective?\nIt could also be interesting in this context to investigate to what extent what your model extracts can be intepreted, and how it relates to other topological invariants.\n\nRestriction to 2D\n\nThe restriction of all experiments to 2-D is also quite limiting. The authors indeed mention \"technical challenges\" in their conclusion for extensions to higher dimensions, but there are deeper theoretical reasons for these challenges, e.g. the lack of structural stability in higher dimensions, and more generally the scaling of your approach with the dimensionality of the dynamical system. Since this is a substantial limiting factor for extending your method to experimental settings, this could be made more explicit in the discussion.\n\n\nIn summary, in case I did not fundamentally misunderstand something about your approach and the general problem setting, it is unclear to me why you would need such a sophisticated machine learning architecture to extract features that can in the situation you investigate could be estimated from the flow directly? Since this is a crucial point to my mind point I am voting for rejection at this point, but am willing to adjust my score if these concerns are adressed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "iPrOZgynN6", "reviewer": "Reviewer_YHJc", "summary": "This paper introduces a method to identify the existence of fixed points and limit cycles in dynamical systems based on their vector fields in order to able to detect supercritical Hopf bifurcation boundaries. This method relies on a convolutional feed-forward network with attention and data augmentation that relies on warping. This method is then applied to some toy systems and single-cell data to identify different regimes in parameter space and differentiation transitions respectively.", "questions": "Why was a learning rate of $10^{-4}$ chosen for training? Why 20 epochs? How was it \"confirmed\" that that \"was enough to fit the training data\"? Why were 20 coefficients used for the Fitted Parameter Coefficients?\n\n\nWhat is the discontinuity of the theoretical boundary in Figure 3 and Figure A12 relating to?\n\nHow are cells partitioned into batches of at least 50 cells? How is the size furthermore determined for each batch? Does the partitioning have any consequences for the resulting model?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is to our knowledge the first method to perform topological invariant feature learning using deep learning.\nThe method has a high performance to detect bifurcation boundaries in various dynamical systems and even on real world data. Especially the application to detection of cell cycle score is a surprising contribution.\nFurthermore, the effectiveness of the approach is demonstrated against a couple of other methods highlighting its applicability.\nThe paper is well written and present with a good flow to explain clearly the contributions. The particular contributions of the different parts of the network architecture and training method provides useful insights.  Robustness to relatively high levels of noise seem to further suggest the usefulness of this method for real world applications.\nIdentifying bifurcations is an important problem in dynamical systems theory and this paper shows some promising results to tackle it.", "weaknesses": "\\paragraph{Minor comments}\n\"In the left portion of Fig. ??\" should have been \"In the left portion of Fig. A7\"?\n\nThe reference \"ML Cartwbight. Balthazar van der Pol, 1960.\" seems erroneous.\n\n\\paragraph{Performance on systems that have both features}\nIt would greatly benefit the demonstration of the usefulness of the method, if there would be an assessment of the performance on systems that have both a fixed point and a limit cycle.\n\n\n\n\\paragraph{Figure 2}\nIt is surprising that the correlation coefficient for the BZ reaction is higher than Supercritical Hopf, while the boundary seems quite off.\n\n\\paragraph{Use of different training set sizes for the different methods}\nThe comparison of the method proposed in the paper to phase2vec and the Vector Field Learned Representation seem unfair considering that the models have been trained on a different number of training samples. Why is your model trained on 10K training samples why the others are trained on only 7.5K. This could explain why phase2vec is performing worse. \n\nIt would further clarify the differences between the methods if you could furthermore show how the inference of the bifurcation boundaries look like for the other methods.\n\nFinally, how do the other methods perform with the addition of noise?\n\n\\paragraph{Use of the chosen hyperparameters}\nThe particular choice for the used hyperparameters is insufficient and makes the comparison to the other methods less convincing. See also the Questions. Showing the dependence of the performance of the methods on the different parameters would show how sensitive the method is to hyperparameter tuning. \nOr stating that the best performing hyperparameters were chosen for each method would make the comparison better.\n\n\n% Significance: \n\\paragraph{Comparison  to other methods}\nTo fully assess the contribution of this work a more extensive set of method to compare to should be considered.  \nFirst, of all, Lyapunov exponents could be used to track bifurcations and the existence of limit cycles (Witte, 2014). This method work for any dimension in principle, not just 2.\n\nFurthermore, the Conley-Morse graph (Arai, 2009) (has more information, works in a more general setting (other dimensions than 2)) should be considered.\nIs this method faster? More accurate than constructing the Conley-Morse graph or the method proposed in the paper?\n\nFinally, one could track down bifurcations through continuation algorithms (Veltz, 2020), see\n\\url{https://docs.juliahub.com/BifurcationKit/I1INQ/0.1.4/detectionBifurcation/}.\n\n\nThe performance of these other reliable methods would give a better idea of the usefulness of the proposed method.\n\n\n\nDe Witte, V., Govaerts, W., Kuznetsov, Y.A. and Meijer, H.G.E., 2014. Analysis of bifurcations of limit cycles with Lyapunov exponents and numerical normal forms. Physica D: Nonlinear Phenomena, 269, pp.126-141.\n\nVeltz, R., 2020. BifurcationKit. jl (Doctoral dissertation, Inria Sophia-Antipolis).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vSuj86tclV", "reviewer": "Reviewer_YdjJ", "summary": "Authors propose a data-driven, physically-informed deep-learning framework for classifying dynamical regimes and characterizing bifurcation boundaries based on the extraction of topologically invariant features. Authors further demonstrate the method's use in analyzing real data, recovering distinct proliferation and differentiation dynamics along pancreatic endocrinogenesis trajectory in gene expression space based on single-cell data.", "questions": "See above.", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "1. The experiments are extensive. \n2. The related works are well organized.", "weaknesses": "1. The paper is written poorly and hard to read. I cannot find the problem definition. What's the meaning of prototypical system? \n2. The main method is based on convolutional attention method and data augmentation, which is widely used in dynamical system modeling.\n3. The methodology part needs to be reorganized. Authors should divide that into different parts, which can be related to the contributions. Now, it only shows data augmentation and an existing network to me.   \n4. More SOTA methods should be compared. For example, there are extensive RNA velocity methods [1]. \n\n[1] Deep generative modeling of transcriptional dynamics for RNA velocity analysis in single cells, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "LyrbziQNoC", "reviewer": "Reviewer_sMRU", "summary": "In this paper, the authors present a deep-learning framework designed for the classification of dynamical regimes and the characterization of bifurcation boundaries. Their approach is rooted in the extraction of topologically invariant features. The effectiveness of this method is convincingly demonstrated through its application to real-world data. Notably, the framework excels in the analysis of such data, successfully delineating distinct proliferation and differentiation dynamics along the pancreatic endocrinogenesis trajectory within the gene expression space.", "questions": "In addition to the points mentioned in the 'Weaknesses' section, several questions arise:\n\n- Regarding the model setup, it would be valuable to understand why the current configuration is deemed the most suitable for this problem. An ablation study, exploring architectural variations such as layers with different sizes, would provide insights into the model's design choices.\n- Are there specific hyperparameters within this method that require tuning, and if so, what is the recommended approach for hyperparameter optimization?\n- A dedicated 'Limitations' section would be a useful addition to provide a comprehensive perspective on the constraints and challenges of the proposed approach.\"", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 4, "strengths": "- The paper introduces a novel approach to topologically invariant feature learning in dynamical systems, utilizing warped vector field data as augmentations to generalize across systems with similar dynamics in their topological representations within phase space. The method seems novel to me.\n- The paper demonstrates the effectiveness of the proposed method through its application to real data, successfully recovering distinct proliferation and differentiation dynamics within the gene expression space of the pancreatic endocrinogenesis trajectory. This suggests a broad range of applications across various fields, including biology, physics, and beyond.\n- The paper holds promise for making a substantial impact in the field of machine learning for dynamical systems, offering a new approach to addressing the challenge of predicting impending catastrophes in diverse real-world systems.\n- The paper also suggests promising future research directions, including expanding the scope of invariance types and enhancing equivalence notions.", "weaknesses": "- The methodology section appears disorganized. The addition of more subtitles would greatly enhance readability. Moreover, the frequent transition between the main content and the appendix disrupts the reading experience.\n- The model setup details being located in the appendix can make it less accessible, as it's separated from the primary content.\n- While the wording and sentences in the experiments section are well-constructed, there's room for improvement in the organization of tables and figures. Consider enlarging figures and addressing font size issues, particularly in tables. For instance, Table 3's column headers could be better aligned with their respective columns.\n- The process for reproducing the results remains unclear. It would be beneficial to include a dedicated \"Reproducibility Statement\" section, providing readers with the necessary guidance to replicate the findings.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)