Let's do the time-warp-attend: Learning topological invariants of dynamical systems

We thank the reviewer for the helpful and practical comments. We have adjusted the manuscript accordingly, specifically:

1. We have clarified the text and divided the Methodology section into subjects. To improve clarity and flow, we moved some of the model details from the Appendix to the Methodology and kept the more technical details in the Appendix. We framed and centered tables’ content and adjusted font size where plausible.

2. We now provide an anonymous link to our preliminary repository as a footnote at the end of the Introduction section.

3. We added results for varying hyperparameters (learning rate and training epochs) and architectural choices (dropout rate, kernel size, MLP latent dimension size, number of convolutional channels, and number of convolutional layers), see Tables A4 and A8. We found no significant effect on performance for these hyperparameters and ablations, except for a reduced performance with learning rate = 1e-5.

4. We added a dedicated limitations section, touching upon the challenges of combinatorial annotations and scaling to higher-dimensional systems.

We thank the reviewer for their remarks.

Consequently, we have extensively revised the introduction and methods to clarify our problem setting. In brief, our aim is to learn a limit cycle detector which is robust both to nuisance features (cycle shape, frequency, etc) and noise. We achieve this with the topological data augmentation method described in the manuscript, which ensures that those nuisances are ignored in favor of the underlying topological structure which defines the bifurcation classes. Our system is trained on augmented versions of a simple system (the “prototype system”, in this case the simple oscillator system whose equations can be found in “Dynamical equivalence and topological augmentation” subpart in Methodology) which concisely captures the main features of these classes. We have adjusted the text to emphasize these details.

Further, in the Related Works and in the Results sections, we now argue and demonstrate that existing topological methods are either very sensitive to noise, or require system-specific prior knowledge, like the explicit form of the governing equation or the true bifurcation parameter.

To our knowledge, current data-driven methods for bifurcation detection deal largely with the problem of predicting the onset of a bifurcation from time series data measured from a single system. Our work is the first to tackle a much more general problem: a data-driven approach to automatically delineate dynamical classes across a huge number of real and synthetic systems varying both in their functional and parametric form. And, while the reviewer is correct that convolutional architectures are not uncommon in data-driven dynamical systems, they are certainly rare in this novel problem setting.

Our paper makes several key contributions in service of this goal. First, it introduces a novel approach for learning topologically invariant features (further clarified in the revised text) in dynamical systems by leveraging warped vector field data. Second, it achieves invariant classification of complex synthetic data. Third, the paper demonstrates the use of classifier confidence to effectively identify bifurcation boundaries. Last, the proposed method is successfully applied to analyze biological systems, specifically uncovering distinct dynamics in the pancreatic differentiation trajectory using single-cell gene expression data where alternative baselines fail (see updated Results section 4.4).

Last, veloVI (Gayoso et al. 2023), the method the reviewer mentioned, is a generative framework of RNA dynamics. It is an alternative to the upstream steps prior to our method. In our analysis we use the current convention in the field, applying standard preprocessing steps (e.g. filtering genes, normalizing expression, etc.) and extraction of velocities with scVelo. For specifically analyzing single-cell RNA velocity data, the coupling of veloVI and our framework in the future can be very interesting. It would require introducing confidence across phase space but will also better mitigate its partial coverage.

Given the interdisciplinary nature of our framework, which spans dynamical systems, topology, and machine learning, we aimed to offer comprehensive background information and diverse references to cater to a broad readership. We believe that our manuscript and the works it cites provide adequate background for a machine-learning audience. Indeed, we believe that the manuscript is also much clearer than the initial draft, as has been recognized by other reviewers. In the final draft, we will reproduce standard mathematical definitions in a glossary in the appendix.

Thanks for your response. Now the paper is till challenging for many machine learning researchers to understand due to the use of specific technical terms like "cycle detector" and "prototype system." To enhance clarity and accessibility, it would be beneficial for the authors to provide definitions for each technical term they use. This approach can make the paper more comprehensible to a broader audience, including those who may not be familiar with these specific concepts.

We appreciate the important and constructive comments raised by the reviewer. We have made the following adjustments to the manuscript:

1. Minor comments - These are now corrected in the manuscript.

2. We added a new experiment in the appendix regarding systems that have mixed limit cycle / fixed point behavior. Note that we use a supervised, multi-class setup so that the classifier can predict both the existence of a limit cycle and a fixed point. Empirically, we find that our method learns that the point and periodic classes are exclusive. For emphasis, we tested our trained model on samples of a subcritical Hopf bifurcation (see Appendix 7.6 having regimes of point attractor, point and periodic attractors, and periodic attractor). We find that, even though our system has never seen this bifurcation, it can detect the monotonic relationship between class score and the order parameter. In the future, we aim to extend to combinatorial annotations of invariant sets. We report this as a limitation and in the discussion converse about how we intend to address this in follow-up work.

3. Figure 2 - Pearson correlation assesses linear relationships. Since the confidence of the BZ reaction oscillator gradually increases with the distance from the bifurcation, its Pearson correlation is very high, as opposed to systems like Supercritical Hopf whose confidence quickly plateaus. We now plot the confidence as a function of the distance from the bifurcation boundary in Figure A11 alongside Spearman’s correlation to quantify the monotonic relationship.

4. Use of different training set sizes for the different methods - We adopted the reviewer’s suggestion of training baselines on train data of equal size of 10K. This had little effect on the performance of the baselines (see Table A4).

5. Effect of noise and bifurcation boundary diagrams across methods - In Tables 1 and 2 we now benchmark the performance of baselines with noise and their bifurcation boundary diagrams in Figure A12.

6. Use of the chosen hyperparameters - In Table A4 we report test accuracies for Phase2vec, Autoencoder and Parameters across learning rates and epoch iterations. Some hyperparameters improved the average accuracy, e.g. for Phase2vec, changing the learning rate from 1e-4 to 1e-3 increased accuracy from 66% to 72%, however, the ranking of methods across systems is maintained. For all learning rates except 1e-5, our method achieved consistently high accuracies, even when varying architectural parameters, see Table A8.

7. Comparison to other methods - To motivate our approach to learning and generalization, we have added new numerical experiments which compare to methods making use of data topology, namely, via (1) heuristics based on critical points, and (2) Lyapunov exponent thresholding, see Appendix 7.3 for details and Tables A4 and 1 for performance on noiseless and noisy data. We observe that the Lyapunov exponent method is more sensitive to noise than our approach (see Table 2). The Conley Morse Graph Database proposed in (Arai et al. 2009) assumes knowledge of the governing equations, so we did not consider it. The Conlety-Morse graph computation from vector representation developed in (Chen et al. 2007; Chen et al. 2008) was not accessible as the software only works in Windows 10. Finally, numerical continuation algorithms operate by iteratively constructing intersecting solution curves in conjoined phase-parameter space. In most practical settings, we do not have access to parameter space and must contend with a dataset of isolated measurements, for example, time series or vector fields. We do not know how those measurements relate to one another in parameter space, and we certainly lack the ability to parametrically interpolate between those samples in a continuous way. The practical data scientific question then becomes ”how does one detect invariant sets in individual examples?”, which is what our paper strives to address.

8. As we now note in the revised manuscript, we chose polynomial representation of degree 3 following the earlier work of (Iben and Wagner 2021) which showed this value was sufficient for dynamical systems similar to those we investigate here.

9. The discontinuity of the theoretical boundary of the repressilator system is due to our discretized computation of the boundary. We take a grid of alpha values and solve for their upper and lower boundaries. We mention this now in Appendix 7.7.

10. For constructing the pancreas dataset, we first partition cells according to their cell type. Given n cells and minimal size s=50, we generate k=floor(n/s) labels, bin the cells into k buckets and then permute the labels (see notebooks/pancreas.ipynb in our repository) . Having observed the 2D grids that emerge with the current batch size (see Appendix Fig A13) we were satisfied since we observed a clear distinction between the dynamic regimes and were able to generate a sizable cohort of 69 samples.

We thank the reviewer for the valuable insights and constructive feedback. We have made significant revisions to the manuscript to extend the discussion on related works and underscore the motivation behind our approach over and above existing methods which make use of data topology. As we now demonstrate in several new numerical experiments, these earlier methods fare significantly worse in classifying noisy vector field data, both synthetic and real, as they are either tricked by spurious signals from fixed points or struggle to learn a single decision boundary across all systems. In contrast, our method learns features which are both robust and universal across systems, making our method more viable in practical settings than previous techniques.

We address the main points raised by the reviewer:

1. Following the reviewer’s suggestion, we have extended the Methodology section to better contextualize our method in the dynamical systems literature. In particular, in the subsection “Topological equivalence in dynamical systems”, we discuss the literature cited by the reviewer (Kutznetsov and the standard notions of topological conjugacy and equivalence) and precisely what we mean by the term topological “invariant”. In the subsequent subsection, “Dynamical equivalence and topological augmentation”, we elaborate on our augmentation scheme in light of this literature. We are now particularly careful to distinguish the term “dynamically equivalent” from similar notions in the literature while highlighting its role in preserving limit cycles and fixed points, which is what we prove in Appendix 7.1.1.

2. We thank the reviewer for recommending new baselines to justify our approach. We emphasize first that our aim is eventually to classify real data, and our reasoning is that generalization based on learned neural features will be more robust to the noisy conditions inherent to real data than existing measures based on direct, hardwired measurements of flow. For example, the reviewer suggested using a heuristic based on fixed points might obviate the need for learning altogether. We implemented such a heuristic and found that it achieves nearly perfect classification accuracy for noiseless synthetic systems but entirely fails when even a little noise is added ( see Appendix 7.3 for method details and Tables A4 and 1 for performance). Similarly, classification accuracy using this heuristic for the single-cell pancreas dataset is only 54%. Persistent homology, as well as other dynamical measures like Lyapunov exponents (Stoker 1992; Zhang et al. 2009), rely on time-series data, which is not always available, e.g., in single-cell RNA velocity data. We have added a benchmark that generates trajectories from the vector field and computes its Lyapunov exponent. This method, too, was more sensitive to noise than our approach (see Table 2). Concerning persistent homology, we currently followed (Berwald et al. 2013) in using the two longest bars of degree 1 persistence barcode for classification, but observed low accuracy on the train data (63%), hence, we did not include it in the attached submission but are currently pursuing its integration. We are highly motivated to further explore the rules learned by the model, to compare these to existing theory, and to leverage such understanding for extensions of the framework.

3. In a new section on the limitations of our approach, we address the restriction to two dimensions along with other constraints. We would first note that preprocessing by dimensionality reduction can mitigate difficulties from high dimensions, as we use for the analysis of the six-dimensional repressilator, or for the single-cell data with over 2K dimensions. What’s more, low-dimensional dynamics are still an active area of study, as our experiments in single-cell trajectories show. Indeed, superficially high-dimensional dynamics in biology are often intrinsically low-dimensional and suited to our approach (see for instance Xiong and Garfinkel 2023). Exponential scaling with dimensions is inherent to the vector field format, although we imagine several ways to circumvent this problem. Some recent works have used graph neural networks or transformers to learn compact representations of phase space from sparse data. We will examine these extensions in future work.

4. We did not entirely follow the reviewer’s observation: "Your examples indicate that augmentation actually decreases performance for the simpler baseline models, which goes a little counter to the argument of your paper." In general, the augmentations mitigate the overfitting to simple systems, as discussed in Appendix 7.4. Although we train on Augmented SO data, performance on unseen, simpler test systems, like SO, achieve higher accuracy across frameworks, as opposed to the trained Augmented SO system, because the augmented data is very complex, being asymmetric, off-center, and sometimes with limit cycles partially out of frame.