We thank the referee for the enthusiastic support of our work and the thoughtful comments!

Weaknesses

W1 (description of training method): Thank you for raising this point, this was indeed described a bit too briefly. We have now thoroughly rewritten sect. 3 and expanded on these points; at the same time streamlined this sect. and moved details which are less essential (and may have disrupted the flow) to the Appx.. We also included a new figure (Fig. 8) to further clarify this aspect, and amended the model overview in Fig.1.

In essence, the trainable parameters of the overall framework consist of (a) the group-level matrices $\\mathbf{P}\_{\\mathbf{A}}$, $\\mathbf{P}\_{\\mathbf{W}_1}$ etc. (as listed on lines 170-174), (b) the subject-specific features $\\mathbf{l}\_j, j=1...S$, and (c) scaling parameters $\\Sigma\^{(j)}$ in eq. 5. These are all trained simultaneously, end-to-end, on the MSE loss eq. 5. While training, GTF replaces subject-specific latent states in eq. 4 by an optimal linear mixture of model-forward-iterated states $\\mathbf{z}\_{t-1}$ and data inferred states $\\mathbf{z}\_{t-1}=\\mathbf{x}\_{t-1}$ as described in A.1.1. As each of the matrices in eq. 4 is a linear combination of group-level matrices $\\mathbf{P}$ and subject-specific features $\\mathbf{l}_j$ (i.e., could be rewritten in terms of those as on lines 172-174), qua chain rule GTF’s gradient control transfers directly to these matrices (see first pg. of sect. 3.2 and new Fig. 8).

W2 (alternative models): We had initially chosen the shPLRNN as our base model based on its SOTA performance compared to many other architectures, including SINDy (Hess et al. [1]). However, we completely acknowledge this point and addressed it by repeating the experiments in sect. 4.1 with vanilla RNNs (still trained by GTF) and LSTMs (as e.g. used in Vlachas et al. [2]) as DSR models, see Appx. A.1.6. New Table 3 confirms that the basic pattern of results stays the same.

W.r.t. PINNs and SINDy in particular, however, note these are different from our approach (and from neural ODEs) in that they both require physical priors about the functional form of the underlying equations, which for complex systems like the EEG or fMRI data considered here is a bit difficult to specify. In fact, it has been shown previously that SINDy fails when its function library doesn’t match the ground truth (Göring et al. [4]) or in empirical scenarios where such priors are difficult to obtain (Hess et al. [1]; with related problems reported for PINNs in e.g. Krishnapriyan et al. [5]; we further observed that SINDy also runs into severe numerical issues on higher-dimensional datasets like those used here, with models practically always diverging).

In general, a physical prior should also be sufficiently broad to accommodate all systems trained on, not just one particular system, a point on which we added a note to the Limitations sect.. One such more generic prior could be strong dissipativity in a class of DS. This assumption is used in Markov Neural Operators [4], and so we tested MNOs as an architecture with a physics-inspired regularization loss within our hierarchisation framework (see new Appx. A.5.1); new Fig. 17 essentially reproduces the observations from Fig. 3 with the MNO, likewise obtaining an about linear relation between the inferred feature value and the true control parameter.

Questions

Q1 (Eq. 5): See our response to W1 above: The $\Sigma^{(j)}$ are trained simultaneously with all other model parameters, as we now clarified in the first pg. of revised sect. 3.2.

Q2 (how training works): We hope this is now sufficiently addressed by our response to W1 above, by the textual changes and additions we have done in this context (rewritten sect. 3.2), and the updates to Fig. 1 and the new Fig. 8.

References

[1] F. Hess et al. “Generalized Teacher Forcing for Learning Chaotic Dynamics,” in Proceedings of the 40th International Conference on Machine Learning

[2] Vlachas et al., “Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 474, no. 2213, p. 20170844, May 2018

[3] Brenner et al., “Almost-Linear RNNs Yield Highly Interpretable Symbolic Codes in Dynamical Systems Reconstruction,” Oct. 18, 2024, arXiv: arXiv:2410.14240

[4] N. A. Göring et al. “Out-of-Domain Generalization in Dynamical Systems Reconstruction,” in Proceedings of the 41st International Conference on Machine Learning

[5] A. S. Krishnapriyan et al., “Characterizing possible failure modes in physics-informed neural networks,” in Advances in Neural Information Processing Systems, vol. 34, 2021.

We thank the referee for the general appreciation of our work, and the constructive comments.

Weaknesses

W1 (clarity): We reworked major parts of the text, in particular completely rewrote the Methods sect. 3.1 & 3.2, in an attempt to streamline it and improve clarity (reducing sentence length, rearranging material, adding explanation, revising Fig. 1 & adding new Fig. 8; please compare revised to previous version). This is important, as we feel there were indeed some major misunderstandings, as explained in more detail in our responses below.

W2 (novelty): Our model is not primarily a time series (TS) model, but a dynamical systems reconstruction (DSR) approach. It is important to note that, while TS modeling and DSR are certainly related, they are not the same and have partly different objectives. DSR is by now a burgeoning field in its own right in scientific ML (as reviewed e.g. in [1,2]), a trend acknowledged by ICLR by including ‘dynamical systems’ as a topic term fairly recently. Unlike conventional TS modeling, in DSR the central goal is to learn from data a generative model of the underlying system dynamics, which means the model must also express similar long-term statistics (as assessed by, e.g., agreement in auto-corr functions or power spectra) and a similar state space geometry with the same attractor states (as assessed via our geometrical criterion $D_{stsp}$). This is crucially important for scientific purposes, where we seek a model that reflects the ground truth in its dynamical mechanisms. But it might be less relevant for conventional TS analysis, where we may be content with a model that produces good ahead predictions or downstream classification results.

This crucial difference we have now illustrated in new Fig. 28, where we compared the shPLRNN trained with GTF to various popular TS models, including N-HITS, temporal CNNs, and Transformers. While all models produce reasonable short-term forecasts (up to the prediction horizon given by the exponential divergence of trajectories in chaotic systems), only the shPLRNN trained by GTF produces sensible reconstructions of attractor geometry and long-term statistics.

To the best of our knowledge, ours is the first comprehensive hierarchical framework for DSR that allows for integration from many different systems, living in different dynamical regimes (e.g. with chaotic vs. periodic dynamics), into the same underlying latent model. The only somewhat related frameworks (LEADS and Coda) were profoundly out-performed by our construction, and often were not able to actually reconstruct the systems in question (see Table 1 & Fig. 16). More importantly, we report an important connection between subject/system-specific feature vectors and crucial control parameters of the ground truth (GT) system. This semantic structuring of the feature space with respect to GT parameters that control transitions between different dynamical regimes (now backed up by further results on simultaneous observations from completely different systems, new Fig. 27), is, we think, extremely valuable for scientific/ medical interpretation. And, although these scientific applications were our prime objective, we believe it may also have strong implications for downstream classification tasks, usually a domain of classical TS models, as shown in Fig. 7.

W3 (practicality): We would argue that our model is far less complex, as, for instance, the Transformer architecture (Vaswani et al., [3]) with its many different mechanisms and components (temporal encoding, embeddings, encoder/ decoder/ MLP layers, self-/cross-attention etc.). In essence, there are parameters at only two levels, the group-level (common to all subjects) and the subject/system level (the subject-specific feature vectors). These are related to the NN parameters that instantiate the subject-level models (which all share the same architecture, eq. 4), through a simple matrix multiplication and thus just one step of the chain rule, computed by autograd [4]. Hopefully some of these aspects are clearer now after our comprehensive rewrite of sect. 3.

The framework’s feasibility for end-to-end training on DSR tasks we believe we have thoroughly demonstrated in Figs. 2-7, 23-26, using a variety of different benchmarks and challenging real-world datasets (see also Table 1).

I am sorry that the authors find my comments not "honestly helpful." However, I am not sure I appreciate the the authors' tone.

My critique is precise. The authors had the opportunity to provide a rebuttal to improve their score (which I was open to considering), but their response failed due to its excessive length, inconsistency, and logical gaps. Unless these issues are resolved, it is not productive for me to provide specific examples or citations to guide the authors. Repeated inquiries from the authors did not and can not change this fact. It is the authors' responsibility to address these problems and clarify their work. It is not the reviewer's duty to guide the authors' research.

My only reason for posting my previous response was that the authors asked for clarification again. It is weird that the authors are now accusing me of "drowning the discussion." I find this behavior highly unprofessional. For this, I am lowering my score to a strong reject.

We thank the referee for the very supportive and encouraging feedback!

Critical Weaknesses

W1 (Related Work on DSR): Thank you for pointing out these references, which we now have integrated into our discussion of related work.

W2 (Related Work on Hierarchical Modeling): Thank you, we now mention these in the 2nd pg. of our related work sect.. LFADS, however, in our perception is not strictly designed for DSR but more for posterior inference. Volkmann et al. [5], for instance, show LFADS is not capable of reconstructing dynamical invariants, like attractor geometries, from fMRI data. Combining data from multiple sessions is also not quite the same as reconstructing dynamics from many different systems and dynamical regimes: While LFADS allows for session-specific en-/decoders to accommodate, e.g., different neuron numbers, as far as we can see it does not have subject-specific parameters for the recursion model that would allow for reproducing subject-specific regimes, but the generative model and dynamics are assumed to be the same across sessions. The other paper [4] seems very recent and currently does not provide any code which would allow for a direct comparison (to us it appears it is under review for this same conference). Hence, both methods are difficult to directly compare with (we did, however, include a number of other new comparisons in the revision, see list of additional experiments in general response).

W3 (dependence on number of subjects): Yes, fair point! We now systematically tested this, both as a function of the number of subjects and the length of the time series, see new pg. ‘Scaling and robustness’ at end of sect. 4.1 and new Fig. 13. The referee’s intuition is right, as the number of subjects $S$ increases, performance gains are initially strong and then start to level off logarithmically. Importantly, however, multi-subject information always helps to improve performance, and very profoundly so when the time series are rather short (see also Fig. 12).

W4 (overlap between subjects required): Thank you for bringing this up. First, please note that in the experiments on the cardiovascular, the EEG or the fMRI data, indeed many different subjects, conditions, or brain regions were combined, with, in part, largely different dynamics (e.g. epileptic activity vs. normal EEG), showing that this is not only possible, but rather a strength of our approach. We would argue that these differences in dynamical regime (e.g. chaotic vs. periodic) that our model is capable of integrating (as also systematically explored in the simulated benchmarks), is the more crucial aspect, rather than the different systems (subjects) the time series come from.

To further address the referee’s point, we have now added to our revision new experiments where the hierarchical model simultaneously reconstructed dynamics from three completely different systems (Lorenz, Rössler and Chua system). We found that still only low-dimensional ($\leq 4$) feature vectors were required to faithfully reproduce dynamics from all three systems and the different dynamical regimes, with a semantically clearly interpretable feature space, see new Fig. 27 and end of revised pg. on ‘Multidimensional parameter spaces’ (sect. 4.2).

Of course, at some point, if dynamics from vastly different systems evolving on completely different time scales are to be combined (like the referee’s example of neuropixels vs. annual climate data), our approach would likely need to assign rather high-dimensional feature vectors to still accommodate these differences. On the other hand, from a practical perspective, this may also yield a recipe for how to use our framework to decide how dynamically related different systems are: If they occupy largely orthogonal regions in the model’s state and feature spaces, this would indicate they are dynamically rather unrelated. We have included these thoughts in our revised Discussion.

W5 (number of hyper-parameters): Not really, probably our presentation has evoked a bit of a wrong perception here: On top of the usual parameters one would need to specify in any deep learning framework (like model size L), only the size of the subject-specific feature vectors needs to be provided, as now clarified in the new pg. at the end of sect. 3.2 (see also Appx. A.1.11 and Table 5 for more details). The GTF rate $\alpha$ can in fact also be determined automatically (called ‘adaptive GTF’, see Hess et al. [6]), but for the range of datasets we explored, simply adopting default settings from [6] worked (i.e. w/o further grid search). In fact, all-in-all, there were only two hyper-parameters which needed tuning.

First, I sincerely apologize for my late reply. I read the authors' comments early on and wanted to ensure I had the proper time to thoroughly review their additions and explanations. Please know that the delay does not reflect a lack of interest nor a lack of appreciation for your work and rebuttal efforts.

Regarding the authors' responses:

W2: I understand and agree with the comments about LFADS and [4]. While I still think LFADS is relevant to mention, I recognize the differences. Also, the [4] I referenced seems to be from May.

W3: I appreciate the additional experiments, as they clarify many aspects!

W4: I think the exact clarification you provided in your response to W4 would be a valuable addition to the discussion in your paper.

M3: While I am sure including the figure in the appendix was not a light decision, I maintain that it is crucial to have it in the main text, as you have now done.

Q6: I understand the authors' point about non-stationarity. However, in real-world applications, particularly with biomedical data, non-stationarity is a critical consideration. Equation 4 is piecewise linear, implying an assumption that the main non-stationarity arises from abrupt system behavior changes. While this is a common assumption in recent works (e.g., [1, 2]), it may not align well with how biological data are typically assumed to evolve. I recognize that this is not the central focus of the paper and want to clarify that I do not expect the authors to address this experimentally or modify their dynamic equation model within the scope of this version of the paper. However, given the importance of non-stationarity, particularly in biomedical contexts, I believe it would be valuable to discuss this in the paper’s discussion section or acknowledge it as a limitation or a direction for future work.

I am still reading and reviewing the new experiments and additions and may provide additional comments by the end of the discussion period. Thank you again for the responses.

[1] Zhang, Y., & Saxena, S. Inference of Neural Dynamics Using Switching Recurrent Neural Networks. In The Thirty-eighth Annual Conference on Neural Information Processing Systems.

[2] Karniol-Tambour, O., Zoltowski, D. M., Diamanti, E. M., Pinto, L., Brody, C. D., Tank, D. W., & Pillow, J. W. Modeling state-dependent communication between brain regions with switching nonlinear dynamical systems. In The Twelfth International Conference on Learning Representations.

[3] Mohammadi, Z., Ashwood, Z. C., International Brain Laboratory, & Pillow, J. W. (2024). Identifying the factors governing internal state switches during nonstationary sensory decision-making. bioRxiv, 2024-02.

We thank the referee very much for the valuable comments and appreciation of our work!

Weaknesses

W1 (applic. to other models): This is a fair point. We chose the shPLRNN because it’s a SOTA in the field [2], but now have repeated some of the experiments with vanilla RNN (trained by GTF as in [2]), LSTM (as in Vlachas et al. [4]), and Markov Neural Operators (Li et al. [5]), to highlight the generality of our approach; see reworked sect 3.1 (below eq. 4), Appx. A 1.6, Appx. A.5.1, new Table 3, and new Fig. 17.

W2 (contributions): Thank you for acknowledging that while, indeed, our approach built on previous work, it required many components and a specific construction for the group- and subject-level parameters to work as successfully as presented here. We view this and the demonstration of the efficacy of the approach as one major contribution. Another one is the observation that the subject feature spaces had a clear semantic structure with an interpretable relation between features and crucial control parameters of the underlying system that has not been reported before to our knowledge, but is important for scientific or medical applications.

Questions

Q1 (scaling with $T_{max}$): This is a good question, which we had addressed in Appx. A.1.8, Fig. 12 (formerly Fig. 10), where we investigated relative performance (hierarchical/baseline) as a function of time series length for the 64-subject Lorenz-63 dataset. As this figure (and new Fig. 13) show, increases in performance are essentially obtained for all time series lengths, but are particularly profound for $T_{max} \leq 500$. Admittedly this point was a bit lost in the main text, and we made it more prominent now in a new pg. on ‘Scaling and robustness’ at the end of sect. 4.1. In addition, we now also studied more systematically how performance depends on both, number of subjects and time series length, see new Fig. 13.

Q2 (SINDy & PINNs): Both SINDy and PINNs are different from our approach in the sense that they both require specific prior knowledge about the functional form of the presumably underlying ODE system. In several previous publications where SINDy was pitched against other approaches, including the shPLRNN used here as a basic reconstruction model (see [2]), SINDy performed comparably as long as the pre-specified function library matches the ground truth (GT) system (i.e., contains the same terms, like multinomials of correct order for the Lorenz-63). But it breaks down as soon as the function terms describing the GT are lacking from its library [1]. In particular, SINDy failed on essentially all real-world datasets checked in [2], as for those sufficient prior knowledge is often not available (and using larger libraries is computationally infeasible, as we now explicitly confirmed for all three real-world applications considered in our paper, where even just numerically it became impossible to train SINDy on the complete set of (high-dimensional) observations and trained models always diverged). This is a crucial point of course, as at the end of the day we are primarily interested in those real-world applications where sufficient priors may not be available.

For PINNs, similarly, prior physical knowledge about the functional form of the presumably underlying ODE system is required (implying somewhat similar issues as with SINDy, e.g. Krishnapriyan et al. [3]). In general, the prior should also be sufficiently broad to accommodate all systems trained on, not just one particular system. One such more generic prior could be strong dissipativity in a class of DS. This assumption is used in Markov Neural Operators [4], and so we now tested MNOs as an architecture with a physics-inspired regularization loss within our hierarchisation framework (Appx. A.5.1); new Fig. 17 essentially reproduces the observations from Fig. 3 with the MNO, likewise obtaining an about linear relation between the inferred feature value and the true control parameter.

Thus, incorporation of physical priors as in PINNs is fully compatible with our framework (and our new experiments suggest won’t change the results). But it may require more detailed study in its own right, as we now also commented on at the end of our Conclusions sect. (in Limitations).