A scalable generative model for dynamical system reconstruction from neuroimaging data
Authors: Eric Volkmann
Alena Brändle
Daniel Durstewitz
Georgia Koppe
Venue: NeurIPS 2024 poster (2024)
Keywords: Dynamical Systems Reconstruction
Recurrent Neural Networks
Nonlinear Dynamics
Neuroscience
fMRI
Primary area: neuroscience_and_cognitive_science
TL;DR: We build on the recent success of control techniques for training SSMs for dynamical systems reconstruction (DSR), and propose a scalable DSR algorithm for empirical situations in which we deal with convolved observations, such as fMRI time series.

Abstract: Data-driven inference of the generative dynamics underlying a set of observed time series is of growing interest in machine learning and the natural sciences. In neuroscience, such methods promise to alleviate the need to handcraft models based on biophysical principles and allow to automatize the inference of inter-individual differences in brain dynamics. 
Recent breakthroughs in training techniques for state space models (SSMs) specifically geared toward dynamical systems (DS) reconstruction (DSR) enable to recover the underlying system including its geometrical (attractor) and long-term statistical invariants from even short time series. These techniques are based on control-theoretic ideas, like modern variants of teacher forcing (TF), to ensure stable loss gradient propagation while training. 
However, as it currently stands, these techniques are not directly applicable to data modalities where current observations depend on an entire history of previous states due to a signal’s filtering properties, as common in neuroscience (and physiology more generally). 
Prominent examples are the blood oxygenation level dependent (BOLD) signal in functional magnetic resonance imaging (fMRI) or Ca$^{2+}$ imaging data. 
Such types of signals render the SSM's decoder model non-invertible, a requirement for previous TF-based methods.
Here, exploiting the recent success of control techniques for training SSMs, we propose a novel algorithm that solves this problem and scales exceptionally well with model dimensionality and filter length. We demonstrate its efficiency in reconstructing dynamical systems, including their state space geometry and long-term temporal properties, from just short BOLD time series.

A scalable generative model for dynamical system reconstruction from neuroimaging data
Eric Volkmann 1.2.+, Alena Brändle Daniel Durstewitz! 1.34.f Georgia Koppe
Department of Theoretical Neuroscience, Central Institute of Mental Health (CIMH),
Medical Faculty Mannheim, Heidelberg University, Mannheim, Germany
IIntitutee for Machine Learning. Johannes Kepler University, Linz, Austria
Interdisciplinary Center for Scientific Computing, Heidelberg University, Heidelberg, Germany Facultyy of Physics and Astronomy, Heidelberg University, Heidelberg, Germany
Hector Institute for Al in Psychiatry & Dept. for Psychiatry and Psychotherapy, CIMH
Faculty of Mathematics and Computer Science, Heidelberg University, Heidelberg, Germany s, These authors contributed equally
Abstract
Data-driven inference of the generative dynamics underlying a set of observed time series is of growing interest in machine learning and the natural sciences. In neuroscience, such methods promise to alleviate the need to handcraft models based biophysical principles and allow to automatize the inference of inter-individual diffcrences in brain dynamics. Recent breakthroughs in training techniques for state space models (SSMs) specifically geared toward dynamical systems (DS) reconstruction (DSR) enable to recover the underlying system including its geo- metrical (attractor) and long-term statistical invariants from even short time series. These techniques are based on control-theoretic ideas, like modern variants of teacher forcing (TF), ensure stable loss gradient propagation while training. However, as it currently stands, these techniques are not directly applicable to data modalities where current observations depend on an entire history of previous states due to signal's 's filtering properties, as common in neuroscience (and physiology more generally). Prominent examples are the blood oxygenation level dependent (BOLD) signal in functional magnetic resonance imaging (IMRI) or Ca²+ imaging data. Such types of signals render the SSM's decoder model non-invertible, a requirement for previous TF- methods. Here, exploiting the recent success of control techniques for training SSMs, we propose a novel algorithm that solves this problem and scales exceptionally well with model dimensionality and filter length. We demonstrate efficiency reconstructing dynamical systems, including their state space geometry and long-term temporal properties, from just short BOLD time series.
1 Introduction
Models and theories based on dynamical systems (DS) concepts have a long tradition in computational neuroscience accounting for physiological phenomena and computational processes of the brain 774 Constructing such models from first principles (biophysics) is time-consuming and hard, and utilizing them account for inter-individual differences in brain dynamics, when model settings need to be personalized, is even more challenging. Yet. constructing valid models of the brain's 's functional dynamics is immensely important, not only for understanding the neurocomputational basis of inter-individual differences in cognitive and emotional style (23), but also when aiming at
53 33 23
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

diagnosing or predicting brain dysfunction and clinical characteristics based on DS features [66], or for designing personalized therapies [54).
1.1 Dynamical models in neuroscience
In computational neuroscience, multiple approaches to infer large-scale brain dynamics have been advanced over the past decades (c.g.. [77 46
32 ). Pioneering work introduced latent linear DS, including Gaussian process latent variable models [64 [1 and extensions like Switching Linear DS (SLDS) [27]. Whereas linear DS models are very limited in their dynamical repertoire, SLDS offer a more flexible approach modeling larger range of dynamical phenomena by combining several linear (or affine) DS, with switching mechanism that selects the active system at each moment
58
44 43 (27] These approaches have become common tools for inferring and visualizing neural trajectories within low-dimensional state spaces.
Arguably the most popular class of generative models in whole brain simulations relies mean field neural population models, including neural mass [76 [ and neural field models [35] [5] These model the moments of the activity of cortical areas by averaging over properties (like firing rates) of neural (sub)populations and are often biophysically motivated 17). Many popular large-scale mean
field modeling approaches are implemented in The Virtual Brain (TVB) environment [56] 588. TVB
incorporates biologically realistic brain connectivity into neural field models to generate simulations of large-scale brain activity. Alternatively, Dynamic Causal Modeling (DCM) describes set of more statistically motivated and mostly linear techniques, primarily for the purpose of inferring the effective connectivity between brain regions based on invasive or non-invasive brain recordings (e.g.. [66 37 46)). While many of these models may account for aspects of the dynamics, like patterns of functional connectivity and their task modulation, they are not, strictly, dynamical systems reconstruction (DSR) tools, entailing that they may miss important dynamical phenomena by being constrained through the biological assumptions and simplifications imposed.
In the analysis of functional magnetic resonance imaging (fMRD, only recently a shift in focus has introduced data-driven, deep-learnigg based methods for inferring generative models of system dynamics [40] 62. 63]. Models that implement dynamics either directly on the observed [62]. or within an underlying latent [40. 633 space have been proposed, partly using available structural information and hierarchical inference approaches [63
1.2
Dynamical systems reconstruction (DSR)
variety of Deep Learning (DL)-based models for approximating the generative dynamical processes
underlying observed time series has been put forward in recent years (18 69 49
40
60
42
6
71
31
80
see also 24 for an overview). These include approaches which approximate a system's
vector field, c.g., through a library of basis functions and penalized regression in SINDy g
or through deep neural networks and neural Ordinary Differential Equations (ODE) [4
57 39
Alternatively, methods that approximate the associated flow (solution) operator directly have been suggested, often employing state space model (SSM)-type architectures which distinguish between an observation process and a latent process commonly instantiated through recurrent neural networks (RNNs; e.g., [4] 2 [70 60 (7).. Recurrent SLDS (rSLDS) 44], an extension of SLDS, and Latent Factor Analysis via Dynamical Systems (LFADS)
48 also fall into this category.
In DSR we ask for models that are generative the sense that after training - they provide an executable surrogate model of the observed system, from which we can simulate samples that agree with their empirical counterparts in topological, geometrical and temporal characteristics (in contrast to (24. This required agreement in long- temporal and geometric properties is not automatically guaranteed for standard training of common RNNs or neural ODEs, which may yield good short-term predictions but may fail to recover the full system dynamics [34 24 52 Recent breakthroughs in
data-driven DSR build on insights from the field of chaos control and synchronization [5T 68 2]
by guiding the training process through optimally chosen control signals modern variations of
classical teacher forcing (TF) - that prevent exploding gradients [47 31 39..
Chaotic dynamics in particular, as typical for neural systems (e.g., 67
22
36 26
poses severe
problem here as trajectories and hence loss gradients inevitably diverge due to the presence of positive Lyapunov exponent 47). Recent amendments of TF protocols, including sparse TF

---

(STF; 477 and generalized TF (GTF; [MI, keep trajectories and gradients in check by 'weakly synchronizing' them with the observed signals.
1.3 Specific contributions
Nonlinear SSMs distinguish between an underlying latent process that governs the dynamics of system's 's state, and an observation process (referred to as observation or decoder model), that links the system states to the actual measurements [20]. Invertible (or pseudo-invertible) decoder models play a crucial role in control-theoretic approaches, like STF or GTF, for training SSMs, in order to project observations into the model's 's latent space. This inversion is fairly straightforward when ones assumes that the current measurement depends solely on the present latent state, i.,, when I, f++). but not if it depends this assumption is often violated due to a signal's is filtering properties. For instance, blood oxygenation level dependent (BOLD) signals, assessed via fMRI, are only indirect measurement of neural activity, with the hemodynamic response function (hr f) broadly smearing out the signal across time. Each measurement is therefore function of history of past neural (latent) states whose dynamics we wish to infer 288 [8). Similar challenges arise in calcium imaging when spike times are to be inferred 79 72 or, in fact, any other observation process where the actual measurement is a filtering of the process of interest.
history of states {ze, Rt-1...., 2t-3). J. In practice, unfortunately,
Here we rectify this issue by developing particularly efficient SSM approach which works for measurements that depend on longer histories of latent states, yet allows to take advantage of recent powerful training strategies for DSR (7 B1J. In particular, our contributions are threefold: First, we create and validate a novel SSM-based DSR algorithm for observation models which involve convolutions with latent state series, and demonstrate its scalability with SSM size, well as convolution filter length. Second, we introduce an evaluation scheme for selecting DSR models on short empirical time series, by demonstrating that the used DSR measures assessed short time series accurately predict a system's 's long-term temporal and geometric properties. This is of high practical relevance, as in many empirically relevant scenarios, like fMRI, we only have access to comparatively short time series. Finally, we show that the proposed models can reliably extract key DS features that, moreover, differentiate between subjects.
2 Convolution SSM model (convSSM)
2.1 Latent DSR model
We start by defining our generative model used for DSR, a variant of a piecewise linear RNN (PLRNN). Its specific architecture has first been suggested in (21) and later expanded to increase the PLRNN's expressivity 3I]. While the present approach is generic and independent of the specific DSR architecture, here we use the so-called shallow PLRNN (shPLRNN; Appx. A.3) and the clipped shallow PLRNN (cshPLRNN; BI)). In the cshPLRNN, the temporal evolution of the (latent) system state 2, E R is expressed
(I) where ((-) = ma. is an elementwise ReLU activation function, W1 E MRL W2 E RLXM are connection weights, A E RMM is a diagonal matrix of autoregressive weights, and h E M and h2 E R2 are bias vectors. Its trajectories {z+} E MxT will be bounded if the absolute eigenvalues of A are smaller than [31]. The Markov property of the latent model is crucial to ensure formally constitutes a DS with complete state space [50] 65). Finally, Equation caan easily be extended to incorporate the effect of external inputs, such as experimental stimuli, by adding Cs (with 8 E R representing an input vector and C E R Mx x K its effect on the latent dynamics). However, here we consider input-free data from resting state experiments.
2, A + W1 16 (W24-1 + f2) - ó (W2 + h1
2.2 Teacher forcing for invertible decoder models
In 31J. the latent state 2; at each time point is assumed to be related to the actual observation
It R.V by linear (Gaussian) decoder model
Ft Bz + The Tpe IO,
(2)
3

---

referred to simply as 'standard SSM' in the following. Here, B E RNX M is a matrix of regression weights, and 7i describes Gaussian observation noise with diagonal covariance matrix r. A con- ventional mean squared error (MSE) type loss function L >, L; = E, l| - = 1l between the generated (predicted) } and the observed { sequence is commonly used to optimize parameters by stochastic gradient descent (SGD) with GTF [31. Regularization terms to enforce a structure in latent space that helps to map slowly evolving processes may further be added to this loss [60].
fundamental issue in training such systems by SGD is the gradients' problem (EVGP). preventing systems from capturing relevant time scales in the data. In fact, [47] proved that for chaotic systems gradient-based training techniques for RNNs will inevitably lead to diverging loss gradients (see also (2)).. Successful DSR algorithms need to address this problem. Based on this connection between chaos and diverging gradients, Engelken 25 suggested regularizing the system's Lyapunov spectrum, thereby also biasing the dynamics toward certain (non-hy moon-hyperbolic) solutions. A theoretically well founded approach that does not limit a system's 's dynamical expressivity, which we will adopt here, is GTE, proposed in [31] GTF is designed to keep model generated trajectories on track and, theoretically, can completely abolish the EVGP without constraining model expressivity. The main idea is that during training the latent state & is computed as a linear interpolation between the PLRNN generated state 2,- = PLRRN(22-1) and a data-inferred state d, that serves control signal
'exploding and-vanishing
[9]] i...
== (1 z ta de, a a 10
There is theoretically optimal choice for a that can be approximated concurrentyy whilst training through a specifically designed annealing protocol [3I]. but more simply may just be determined by grid byigrdd search (as done here). Control signals de are obtained by inverting the decoder model (Equation 2). Since in general M * N, the matrix inverse of B E IR N does not exist and is approximated by the Moore-Penrose (pseudo-) inverse B":
(4) To keep the gradients on track, the interpolation is performed at each time step before applying the cshPLRNN mapping (Equation I). These control signals are turned off during actual data generation by the model (i.e., in a test phase), where it runs completely autonomously.
de d B+
2.3 Teacher forcing for decoder models with signal convolution
Decoder model for convolved signals GTF (and similar techniques like STF; (47)) are powerful state-of-the-ar (SOTA) tools for controlling gradients. especially in the context of DSR. However, they require a there are many situations where this requirement not met. For instance, in BOLD time series the observed signal is highly filtered and strongly smoothed version of the underlying neuronal process that we would like to recover [9 28) This fairly complex hemodynamic process often modeled by the hrf (28).
observation model for producing adequate control signals. Empirically,
28).
A decoder model that relates the neuronal processes given as latent time series {z+) measured BOLD time series (1x) may be formulated as in [40].
(5) with regression coefficient matrices B E NxM and J E RNPP, nuisance variables r, E R (such as movement or respiratory artifacts) and Gaussian observation noise term 7) (with usually diagonal covariance r E R~XNN). Here, denotes the convolution operation and is a history of states zt- the length of which depends on the observed sampling rate, commonly referred to as time of repetition (TR). The discrete hrr ʃ sequence is computed by evaluating the canonical hr ʃ at the observed TR 788. We will denote the hr f response for a given TR by hr fIR-
Z3 B f * 2)t)+ + Jr. Te, 7~ NO,
By incorporating the fr into the observation model, we disentangle the neural state and its dynamics the processes of interest from the neurovascular mechanics (or any filtering the level of observed signals). We thereby eliminate the history dependence present in the observations, and thus help unfolding trajectories in latent space and satisfying the uniqueness assumptions required in reconstructing dynamical systems [50] (see also Appx. Figure 7). However, Equation 5) poses a major complication for applying TF techniques, as observations (and model outputs ax)) do not simply depend on the current state 2t, but - due to the convolution on set of states across several previous time steps. We can thus not compute the control signal de through straightforward decoder inversion anymore, but require new type of inversion algorithm.

---

Wiener deconvolution Following 778 we use a Wiener filter [75] to invert Equation 5 We briefly introduce this approach here in the context of our specific problem, and refer to Appx. A.44 for further details. Given an observed noisy signal (*t) ), composed of the signal of interest { convolved with a known impulse response hr plus some noise term T).. (distribution unknown, Wiener optimal for
Gaussian distribution).
Iy = f 2), + Tj.
(6)
the Wiener deconvolution provides the estimate & of the unknown signal & through least-MSE estimation. Defining Conv (,, as the Wiener deconvolution operator, we can write
{&,} Conv fr
Inversion of BOLD decoder model Using the notation introduced above, we can write the inversion to obtain control signals as
{de} Conv
- (x;
(8)
where [B"(x; (x,-Jr() is the time series that needs to be deconvolved. Note that this approach is quite general and simply exchange the fr with alternative functions if we want to account for filtering in the original signal. As stated, since B and are matrices of learnable parameters updated during training, would need to perform this deconvolution step every training epoch, which is computationally very costly. We therefore make use of the linearity of convolutions and separate the deconvolution step from the learnable parameters, rewriting Equation 8 as
{d+} B ( ([te),,r]]] hr Conv Conv-"((re)
With (ndccom)) and { [rdccom)] denoting the respective deconvolved time series, this can be written as de B (ndeccoy - J deconv)
(10)
This now is computationally much more efficient, as we need to perform the deconvolution only prior training. During training then, only the decoder model parameters need to be inferred to obtain the control signal. We will refer to the convolutional model for DSR (Equation I and Equation 5) trained with GTF and SGD as 'COnySSM'. The full inversion algorithm provided in Algorithm| [ with additional information given in Appx. A.6 Key components of SGD+GTF training are illustrated in Figure 1
3 Results
3.1 Performance measures
In DS theory in general, and DSR specifically, we mostly concerned with invariant properties of a system, such as attractor geometry and long-term temporal statistics (24). In chaotic systems in particular, in which trajectories diverge exponentially fast with time, the mean squared prediction
error (PE) is a useful statistic only on relatively short time scales [77] 40 60) (see Appx. Figure 6). To evaluate our model's performance, in addition short-term n-step ahead PEs, PE,. assess the following two established performance measures to capture the temporal and geometrical structure:
1. The deviation in power spectra between the (smoothed) empirical and model-generated power spectra, assessed in terms of the Hellinger distance and referred to as power spectrum enror (PSE), DPSE, in the following [47]. and
2. the Kullback Leibler divergence between the empirical and model-generated trajectories across state space, Dttap, measuring the overlap in attractor geometries 40)
To obtain a reference value for Datap, we further included two references which we assessed Datsp when all mass is centered on the expectation value (similar to fixed point solution), and when the state space is populated by points drawn from a Gaussian with mean and variance equal to the data (similar to fixed point solution plus measurement noise).
For comparability with experimental data, we evaluated performance on comparatively short time series obtained from the adaptive linear-nonlinear (ALN) cascade model and the LEMON data set. In these cases, performance metrics were assessed on 100 generated trajectories per model and then averaged. For more details, Appx. A.7
5

---

A Pro-Traiuing Deconvolution
B
Training: Teacher lorcing
Training: Gradient flow
Dir-1 diger
8thefao)
p-
#z-1
Arf
(hrf xx)
21
Dorome
Bit Bie-1
"aihe fmals 2014
PLIAR
21
2E-1
Figure 1: Schematic of training protocol and gradient flow. A: Before training, observations {x*) and nuisance artifacts {re} are deconvolved. B: The deconvolved time series are used to generate foreing signal d;-1 which used for guiding cshPLRNN training. C: Latent states 2t-rit and nuisance artifacts Fr are used to predict 2n through the decoder model. Gradients are computed on the squared error loss L., propagated from the decoder model back the latent states (blue), and from the latent DS model backwards in time (orange).
3.2 convSSM validation & scalability on Lorenz63
As well-established and popular benchmark for a chaotic system, we first performed numerical experiments with the famous Lorenz63 system. The Lorenz63 is 3-dimensional system introduced in [45] to describe atmospheric convection, and exhibits chaotic behaviour in the chosen regime (see Appx. B.T. To mimic BOLD observations, we generated 100 standardized chaotic Lorenz63 trajectories, convolved them with hr frR functions at different sampling rates TR E: {0.2s, 0.5s, 1.2s} {0.01, 0.1}. This resulted benchmark settings with different levels of signal degradation by convolution and noise. Each data set was divided into training and a test set T = 10 time steps each.
(Figure 23), and added Gaussian noise with standard deviation
We trained 100 models on each of these 6 data sets. The following models were compared: the ConvSSM trained via SGD+GTF, the convSSM trained SGD and no GTF, the standard SSM trained via SGD+GTF, and MINDy, recently published method for DSR in fMRI [62]. convSSM and standard SSMs were trained with the shPLRNN, with M = 3, 50, and a
(see Appx.
Table 4 for all additional hyperparameters). The hidden dimension was selected such that the standard SSM (no-hrf model) performed well [31]. We emphasize that the standard SSM has already been extensively benchmarked on variety of simulated and real-world data sets and is considered to be SOTA model in the field 6 31. The performance measures D Dstsp- DYSE, and PE were assessed on the test sets after training for 000 epochs. We used the same hyperparameters for all networks (aside from TR) to show that performance increases can be solely attributed to the improved decoder model. Hyperparameters were chosen such that the shPLRNN achieved near perfect performance unconvolved, noiseless trajectories from the Lorenz63 system.
The CoMvSSM significantly outperformed all other methods, including the standard SSM in almost all cases, with the performance gap increasing with decreasing TR (see Appx. [Talle 2 for performance, Figure 2 for example reconstructions, providing an intuition on how to interpret Dstsp). The more heavily the signal was degraded by the convolution filter, the larger was the performance gap in favor of the convSSM.
An important consideration especially for large-scale applications of such models to empir- ical data is how well they scale with model size and convolution filter length. To as- sess this, we collected trajectories of length T = 10a from the chaotic Lorenz63 sys- tem, and studied training epoch times
function of convSSM latent dimension M (3, 10, 50, 100,500), hidden dimension I. = (10, 50, 100, 500, 1000). TR (..2, 1.2,3). time scries length T (500, 1000, 5000, 10000, 50000, 100000} and observation dimension

---

D 0.37
D
2.73
D Aop
6.45
D 9.65
D 13:49
A
2/2
B
C
conYSSM conVSSM (NoGTF)
-rree
probability density
1.0
250
300
750
1000
0.8
0.4
0.2
0.0
-1
3
0
250
500
750
1000
Amt
time
D
standard white noise
conv (NoGTF) fixed point
conv E
20
F 20
0,81
couns
15
0.5
0.6
0.3
long
0.4
0.2
0.3
0.2
0.1
5
0.1
0.99***
0.0
0.0
10.0 D stop (obaerved)
5.0
20 D stop (latent)
10 short
15
20
10 observed
15
20
Validations on Lorenz63 and ALN. A: Illustration of reconstruction performance as assessed by the geometrical agreement measure Dstsp- Average Dstsp values for the convSSM were D Dstap 0.30 at noise level a .01 and Datsp 0.71 noise level .1, indicating successful reconstructions in the majority of cases. B: Example trajectory from the Lorenz63 system in latent space (top) and observation space (convolved with hr fu.2) (bottom). C: Probability density over maximal Amaz values (orange) assessed on 1000 convSSMs trained on Lorenz63 time series of length 1000 (example shown in right panel). Black line denotes the known Amaz As 0.9056 of the Lorenz system. D: Comparison of standard SSM ('standard'), ) convSSM ('conv"), and convSSM trained without generalized teacher forcing ('conv (NoGTFF)) on the ALN data set. Histograms over Dstsp assessed on the observed space (left panel) and latent space (right panel). E: Datsp for convSSM evaluated on the full pscudo-empirical time series of typical empirically available length (T = 500; x-axis) vs. the long GT test set (T' 5, 000; y-axis). F: Datep for convSSM evaluated on the observed time series (x-axis) vs. on the latent time series (y-axis).
Figure 2:
N = {10, 30, 50, 100, 500, 1000). Shorter/longer TRs directly implicate longer/shorter convolution filters since the filters assume a constant time interval. Results are displayed in Appx. Figure 4/.. The runtime per epoch did not significantly depend on TR, which means that time series convolved with long impulse response functions can trained in times comparable short ones. The per-epoch- runtime increases approximately linearly with dimensions L,, M, and N (Appx. Figure 4 B,C.E), implying that models can be scaled up efficiently.
Finally, empirical data is often short, yet we want to reliably infer DS features that characterize the underlying dynamics. To demonstrate that our model can robustly reconstruct dynamics based on short time series, we inferred 1000 conv SSMs on n 100 convolved Lorenz63 time series (TR 0.5) of length T 1000 only (see Figure 2C). We then assessed the degree of chaoticity in the recovered trajectories by examining the trained models' maximum Lyapunov exponents, Amaz- Amaz measures how quickly trajectories starting from nearby points a system's 's state space converge or diverge with time. If Amar 0, trajectories will exponentially diverge and the system, if bounded, will exhibit chaos. We show that we can successfully recover Xmaz (known for the Lorenz63 system; Figure 2C) even from models trained these just short series.
7

---

3.3 Validating performance measures on short time series
In empirical situations, we do not have access to the latent dynamics of the true system, of course, but we still rely on our reconstruction measures evaluated on the observed signals to yield results valid for the (unobserved) latent space. It is therefore a practically very relevant question whether a) the convSSM trained on such short time series would be able accurately recover the underlying neural latent dynamics, and b) our measures (D stop; DPSE) evaluated on such short time series, and directly on the observations, would yield results similar to what would be expected if much longer time series and access to the ground truth latent space were available.
To tackle these questions, we used a more realistic simulation model, the ALN model 3 for simulating whole brain (neural) activity. 100 data sets of length T -= 10, 000 were simulated from this model using neurolib [12, sampled 0.1 ms, and filtered through [Equation 5 (with TR 0.1 ms) to compute the corresponding BOLD time series. Subsequently, these time series were downsampled to a TR of 0.5 mimic an experimentally realistic scenario (see Appx. [B. for details). Most hyperparameters were adopted from previous experiments (see Appx. Table 4). For the latent dimension, we chose M 16 to match the dimensions of the empirical LEMON data set, see subsection 3.4]
To mimic real IMRI experiments, we then pretended that only the first 500 time points are available for model estimation (called 'pseudo-empirical here to distinguish it from the actual empirical LEMON data set). We trained 10 convSSM models on the first 375 time steps of each of these virtual experiments, treating the left out 125 time points as pseudo-empirical test set and call the last 5, 000 time points of the entire trajectory (i.e., time steps 001-100, 000 of the full simulation set) the ground truth (GT) test set (which would not be accessible in real experiment). DSR performance was assessed on both the observed } and latent { time series, evaluated for a) the short pseudo-empiriee 500), and c) the GT test set of length 5000 (which also assesses dynamics on the limit set and does not contain transients anymore).
test set of length 125, b) the full pseudo-empiric time series (i.e., of length
Figure 2D shows histograms over Dstop (on the full pseudo-empirical time series) for the convSSM, the standard SSM, and the benchmark conditions. The ConvSSM significantly outperformed the improved
standard SSM in latent space (rank-sum test Z 11.50, p < .001), demonstrating recovery of the ground truth DS and indicating that the deconvolution acts as an inductive bias that forces the model to learn a latent space structured in agreement with our biophysical understanding of IMRI. Moreover, the proposed performance measures (Detsp DSSE) successfully discriminated between good and poor reconstructions even on these short time series more typical for empirical data: For one, evaluating DSR on the observations was consistent with evaluating DSR directly on the latent dynamies space (Figure 2F and Appx. Figure 5). Second, DSR assessed on the useudo-empirice time series (cither full or only test set) was strongly correlated with performance assessed on the long GT test set (which, again, in empirical situations we do not have, Figure 2 and Appx. Figure 5p.
3.4 Application to experimental fMRI data
We finally tested convSSM on empirical data, for which we chose the LEMON study ( ('Leipzig Study for Hind-Body-Emoole Interactions') ') as publicly available data set. This data set was collected at the Daa--Planc-Innsttuuee Leipzig 4 and consists of 227 healthy participants, each of whom completed a battery of tests, including a 15min 30s resting state fMRI (rsfMRI) session sampled at TR = 1.4s (thus comprising T 652 time points). We used the preprocessed rsfMRI data sets as provided, and selected 16 regions from which we extracted a subset of the available time series. These were subsequently smoothed, band-pass filtered, and standardized as in [40]. The time series were split 3 1 into training (Térain 489) and test (Trest = 163) set, respectively. Data from participants with non-stable variance were discarded (i.c., non ntt-stati data, see Appx. B.3 for details), leaving N = 51 participants for analysis.
We trained 20 models on data of each participant with latent dimension M = N 16 (i.e., equal to the observation dimension), a = .1, and hidden dimension L = 50 (where L and M refer to the dimensions of the connectivity matrices W,, W₂ in the cshPLRNN, Equation 1). Latent dimension and a were determined via grid search, by inferring systems using subset of the data and assessing the performance on the held-out set (6. Otherwise the same hyperparameters as used in (31 for
EEG data were applied (see Appx. Table 4 for all details). In addition to the model comparisons
8

---

A
B
C
Floed poins
D Dasp 104
D Dag 1611
bornaklied pounte
200
BOLD signal
60
HOLD
40
200
200
400
6
10
time
time
stop
LO1
D stop
2.02
3.03
D stop 5.29
D
7.0
2
3
Z1
2
3/3
Z
Z2
normalized counts m
9
25
0.05
20
0.00
count
6
15
0.05
10
-..00
3
0.15
0.20
0
...2
0.1
010028
010050
010091
010194
010311
1.0
1.5
Amax
partoriponn (
Participants
Figure 3: Results on empirical LEMON data set. A: Distribution over Dstsp for 1020 systems inferred
with convSSM. B: Example of good and poor reconstruction. D: Illustration of reconstruction performance as function of D stsp. E: Histogram over maximum Lyapunov exponents Amaz. F: Distribution over Amaz for 5 selected participants (n 100 systems with 10 trajectories each). G: Within- compared to between-subject variance in Amaz distribution after filtering models by DS performance measures (selecting the 20 best by Dstsp and 10 best by DPSF).
discussed earlier, we also compared our method to the performance of rSLDS 44 and LFADS 48 (see Appx. C.2 for details). On top of Datetp. DPSE, and PE,. we also assessed the trained models' maximum Lyapunov exponents. Amaz. analyzed how reliably these can be inferred, and whether they distinguish between subjects. Note that obtaining estimate of the Lyapunov exponent is an advantage of the generative model, as empirical time series are often too short to compute it reliably. Also, since we do have access to the cshPLRNN's Jacobians, the computations can be performed analytically (although practically we need to evaluate these along model-genneatee trajectories, where here we used algorithm proposed in 73
The DSR results are shown in Table We obtained successful reconstructions average with mean Decap of 2373, better than all other models (Figure 3/4 and D). Interestingly, most recovered systems were characterized by a positive maximal Lyapunov exponent Xmaz (Figure 3E). indicating the presence of chaotic attractors in these data (consistent with previous observations, 40 36 4T).. Moreover, Amaz values could be inferred reliably (Figure 3F), and differentiated between individuals, as indicated by lower within- as compared to between-subject variation (Figure 3G, T(50) -11.53, .001.
Table 1: DSR measures evaluated for the convSSM. standard SSM, convSSM trained without GTF, well as MINDy (62). rSLDS [44] and LFADS [14), trained on the LEMON dataset. Model runs were excluded if the 1-step PE > on the training data.
metric
CoM
standurd SSM
No GTF 0)
MINDy
ISLDS
LEADS
Noise process
Fixed point
Dacep
2.73 1.09
2.77 0.93
3.77 1.22
6.79 1 1.92
15.51 10.02
3.24 1.17
4.62 0.91
5.27 1.27
Dyer
0.14 0.03
0.15 0.08
0.34 0.11
0.27 4 0.06
0.24 0.03
0.43 0.09
0.76 0.02
1D-step PE
1.78 & 0.38
2.00 0.44
1.43 0.61
1.97 4 0.31
1.78 à 0.42
2.45 0.55

---

4 Conclusions
Methods for producing generative models of the underlying dynamics from time series observations
is rapidly expanding research field [24]. Current SOTA models for this purpose rely on control theoretically motivated training techniques like STF [7 and GTF [31, but these require some means to generate from the actual observations TF signal in the model's 's latent space for guiding trajectory and gradient flows. This becomes difficult if the current observation depends on a whole series of latent states, as common if the actual measurements are some filtering of an underlying process of interest, such as in fMRI or Ca imaging. Here we provide a novel technique that efficiently deals with this problem, exploiting linearity of Wiener deconvolution. A hallmark of our technique is that it efficiently scales with model size and convolution length.
Another major contribution of this work is to numerically demonstrate that the short time series obtained in typical fMRI experiments are actually sufficient for proper model selection according to established DSR performance measures, and that these can indeed be properly evaluated in observation space and do not require access to the unobserved dynamics/ latent space. This is of major empirical relevance for many scientific scenarios, beyond fMRI, in which time series sampling costly or restricted for technical reasons. Finally, using our DSR technique, we showed that experimental fMRI signals mostly exhibit properties of chaotic oscillators (consistent with [36)). and that these can be reliably inferred and differ between subjects. Taken together, these contributions pave the way for deploying data-driven fMRI DSR models large scale to understand inter- inder-individual differences in brain dynamics and explore the predictive value of nonlinear DS features for cognitive or clinical assessment.
We emphasize that the proposed framework is highly flexible due to its modular structure, and may be easily adapted to meet diverse requirements. First, the latent model can be replaced with any other differentiable and recursive dynamical model, such as e.g. LSTMs [59] The GTF training framework would remain unchanged as the control signal and the latent state update (Equation 3) are not affected by such modifications [I] Likewise, the observation model can easily be adapted to account for nonlinear effects of nuisance covariates, e.g. through basis expansions in these variables, or through learnable but regularized MLPs. While our model was designed as a scalable method to integrate biological prior knowledge convolution filters like the fr alternatively we can parameterize the filter weights within the observation model, making them learnable through BPTT, with filter length either hyperparameter, or by imposing a regularization that truncates filter length by driving coefficients to zero. To prevent conflicts between filter adjustment and latent model, a viable strategy may be stage- wise learning as suggested in 40). Once the filter is adjusted, one may reduce the learning rate on the observation model, or even fix its parameters, to prioritize learning of the dynamics. Fixing the filter parameters after initial stage would have the advantage that subsequent training would enjoy the same speed benefits as in our suggested method.
We furthermore highlight that our framework could be adapted to accommodate noise in the latent process. For example, in Brenner et al. 77 the GTF procedure has been modified to work in the context of stochastic DSR models using variational inference. Instead of the multimodal encoder model in Brenner et al. (7]. one may use the inversion in Equation too generate a TF signal which steers probabilistic latent DS model, i.e. controls its distributional mean, via Equation 3 and using the reparameterization trick (55 38 for BPTT in latent space. However, although probabilistic frame- works are appealing, 'deterministic' BPTT has previously been shown to be (at least) comparable in terms of DSR performance, even for clearly noisy observations and latent processes 6 such that the benefits for DSR would need to be further examined.
Limitations Data-driven approaches such as the one proposed here lack detailed biophysical mech- anisms and may thus not be suited to address specific questions relating to pharmacological or receptor mechanisms beyond aunttionaldyyaamical implications. Moreover, currently open questions how to best deal with non-stationarity in the data, how to efficiently combine data from many subjects, and how trained models generalize to oft-oi-domain data.
Software and Data
Code for the ConvSSM is available at btttpa///githbb. httt::://tthabb.ccm//mmm---eaa
10

---

Acknowledgements
This work was supported by the German Research Foundation (DFG) within the collaborative research center TRR 265, subproject BO8, granted to GK. TRR 265 subproject A06 granted to DD and GK, Germany's s Excellence Strategy EXC 2181/1 -309900948 (STRUCTURES), and the Hector Il foundation.
References
[1] H. D. Abarbanel. Predicting the future: Completing models of observed complex systems. Springer, 2013.
[2] H. Abarbanel. The statistical physics of data assimilation and muchine learning. Cambridge University Press, 2022.
[3] M. Augustin, J. Ladenbauer, Baumann, and K. Obermayer. Low-dimensional spike rate models derived from networks of adaptive integrate ontegrate-and-fin and--ire neurons: Comparison and imple- mentation. PLoS Computational Biology, 1305)ee100544,, 2017.
[4] A. Bahayan, M. Erbey, D. Kumral, J. D. Reinelt, A. M. Reiter, J. Röbbig, H. L.. Schaare, M. Uhlig, A. Anwander, P.-... Bazin, et al. mind-brain-body dataset of MRI, EEG, cognition, emotion, and peripheral physiology in young and old adults. Scientific Data, 6(1):1-21, 2019.
[5] M. Breakspear. Dynamic models of large-scale brain activity. Nature Neuroscience, 20(3): 340-352, 2017.
[6] M. Brenner, F. Hess, J. M. Mikhaeil, L. F. Bereska, Z. Monfared, P.-C. Kuo, and D. Durstewitz. Tractable dendritic RNNs for reconstructing nonlinear dynamical systems. In International Conference Machine Learning, volume 162, pages 2292-2320, PMLR, 2022.
[7] M. Brenner, G. Koppe, and D. Durstewitz. Multimodal teacher forcing for reconstructing nonlinear dynamical systems. AAAI Workshops, 2023.
[8] S. L. Brunton, J. L. Proctor, and J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences USA, 1 113(15):3932-3937 2016.
[9] R. Buckner. Event-related fMRI and the hemodynamie response. Human Brain Mapping, (5-6):373-377, 1998.
[10] R. Buxton, Uludağ, J. Dubowitz, and T. Liu. Modeling the hemodynamic response to brain activation. Neurolmage, 23:S220-S223 S233, 2004.
[11] C. Cakan and K. Obermayer. Biophysically grounded mean-field models of neural populations under electrical stimulation. PLoS Compulational Biology, 16(4)x 2020.
[12] C. Cakan, N. Jajcay, and K. Obermayer. neurolib: A simulation framework for whole-brait brain neural mass modeling. Cognitive Computation, ((4)::132-1152) 2023.
[13] J. Chen and K. Wu. Deep-OSG: Deep learning of operators in semigroup. Journal of Computa- lional Physics, 493:1 12498, 2023.
[14] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud. Neural ordinary differential equations. In Advances in Neural Informatinn Processing Systems, volume 31, pages 6571-6588 2018.
[15] G. Datseris. DynamicalSystems.j.. Julia software library for chaos and nonlinear dynamics. Journal of Open Source Software, 3(23):598, 2018.
[16] O. David and K. Friston. A neural mass model for MEG/EEG:: Coupling and neuronal dynamics. Neurolmage, 20(3); 2003.

---

[17] G. Deco, V. K. Jirsa, P. A. Robinson, M. Breakspear, and K. Friston. The dynamic brain: From spiking neurons neural masses and cortical fields. PLoS Computational Biology. 4(8): e1000092, 2008.
[18] D. L. Donoho and I, M. Johnstone. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 8 81(3):425-44 -55, 1994.
[19] K. Doya. Bifurcations in the learning of recurrent neural networks. International Symposium on Circuits and Systems, volume 6, pages 2777-2780. IEEE, 1992.
[20] J. Durbin and S. J. Koopman. Time series analysis by state space methods, volume 38. OUP Oxford, 2012.
[21] D. Durstewitz. state space approach for piecewise-linear recurrent neural networks for identifying computational dynamics from neural measurements. PLoS Computational Biology, 1005544 2017.
[22] D. Durstewitz and Gabriel. Dynamical basis of irregular spiking in NMDA-driven prefrontal cortex neurons. Cerebral Cortex, 17(4):894-900 2007.
[23] D. Durstewitz, Q. Huys, and G. Koppe. Psychiatric illnesses disorders of network dynamics. Biological Psychiatry: Cognitive Neuroscience and Neuroimaging, 6(9):865-876 2021.
[24] D. Durstewitz, G. Koppe, and M. Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, pages 1-18, 2023.
[25] R. Engelken. Gradient flossing: Improving gradient descent through dynamic control of Jacobians. In Advances in Neural Information Processing Systems, volume 36, pages 10412- 10439, 2023.
[26] R. Engelken, F. Wolf, and F. Abbott. Lyapunov spectra of chaotic recurrent neural networks. Physical Review Research, 5(4):043044, 2023.
[27] Z. Ghahramani. Variational learning for switching state-space models. Neural Computation, 12 (4):831-844 2000.
[28] R. Henson and K. Friston. Chapter 14 Convolution Models for fMRI. In Statistical Parametric Mapping, pages 178-192. Academic Press, London, 2007.
[29] D. Hernandez, A. K. Moretti, Z. Wei, S. Saxena, J. Cunningham, and L. Paninski. Nonlinear evolution via -petially-deeendent linear dynamics for electrophysiology and calcium data. Neurons, Behavior, Data analysis, and Theory, 3(3), 2020.
[30] R. Hershey and P. A. Olsen. Approximating the Kullback Leibler divergence between Gaussian mixture models. In International Conference on Acousties, Speech and Signal Processing, volume 4, pages 317-320. IEEE, 2007.
[31] F. Hess, 7. Monfared, M. Brenner, and D. Durstewitz. Generalized teacher forcing for learning chaotic dynamics. International Conference on Machine Learning. volume 202 page 13017-13049. PMLR, 2023.
[32] R. Hindriks, M. H. Adhikari, Y. Murayama, M. Ganzetti, D. Mantini, N. K. Logothetis, and G. Deco. Can sliding-window correlations reveal dynamic functional connectivity in resting- state fMRI? Neurolmage, 127:242-256, 2016.
[33] E. M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT press, 2006.
[34] R. Jiang, Lu, E. Orlova, and R. Willett Training neural operators to preserve invariant mea- sures of chaotic attractors. In Advances in Neural Information Processing Systems, volume 36, pages 27645-27669, 2023.
[35] V. K. Jirsa and H. Haken. Field theory of electromagnetic brain activity. Physical Review Letters, 77(5):960, 1996.
12

---

[36] S. Keilholz, E. Maltbie, X. Zhang. B. Youseñi, W.-J. Pan, N. Xu, M. Nezafati, T. J. LaGrow. and Y. Guo. Relationship between basic properties of BOLD fluctuations and calculated metries of complexity in the human connectome project. Frontiers in Neuroscience, 14:550923, 2020.
[37] S. J. Kiebel, M. I. Garrido, R. J. Moran, and K. J. Friston. Dynamie causal modelling for EEG and MEG. Cognitive Neurodymamics, 2:121-136, 2008.
[38] D. Kingma and M. Welling. Auto-encoding variational bayes. In International Conference Learning Representations. ICLR, 2014.
[39] J.-H. Ko, H. Koh, N. Park, and W. Jhe. Homotopy-base based training of NeuralODEs for accurate dynamics discovery. In Advances in Neural Information Processing Systems, volume 36, pages 64725-64752, 2023. 2023.
[40] G. Koppe, H. Toutounji, Kirsch, S. Lis, and D. Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI. PLoS Computa- tional Biology. 15(8):e 1007263, 2019.
[41] D. Kramer, P. L. Bommer, C. Tombolini, G. Koppe, and D. Durstewitz. Reconstructing nonlinear dynamical systems from multi-modal time series. In International Conference on Machine Learning, volume 162, page 11613-11633. PMLR, 2022.
[42] F. Lejarza and M. Baldea. Data-driven discovery of the governing equations of dynamical systems via moving horizon optimization. Scientific Reports, 12(1):11836 2022.
[43] S. Linderman, M. Johnson, A. Miller, R. Adams, D. Blei, and L. Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial Intelligence and Statistics, pages 914-922. PMLR, 2017.
[44] S. W. Linderman, A. C. Miller, R. P. Adams, D. M. Blei, L. Paninski, and M. J. Johnson. Recurrent switching linear dynamical systems. In Advances in Neural Information Processing Systems, volume 33, pages 14867-14878, 2020.
[45] E. N. Lorenz. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2):130 130-141, 1963.
[46] A. Marreiros, S. Kiebel, and K. Friston. Dynamic causal modelling for fMRI: A two-state model. Neurolmage, (9(1):69--278, 2008.
[47] Mikhaeil, Z. Monfared, and D. Durstewitz. On the difficulty of learning chaotic dynamics with RNNs. In Advances in Neural Information Processing Systems, volume 35, pages 11297-11312, 2022.
[48] C. Pandarinath, D. J. O'ea, J. Collins, R. Jozefowicz, S. D. Stavisky, J. C. Kao, E. M. Trautmann, M. T. Kaufman, S. I. Ryu, L. R. Hochberg, al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10):00-881, 2018.
[49] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott. Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Physical Review Letters, 120(2): 024102, 2018.
[50] L. Perko. Differential equations and dynamical systems, volume 7. Springer Science & Business Media, 2013.
[51] A. Pikovsky, M. Rosenblum, and Kurths. Synchronization: universal concept in nonlinear sciences. Cambridge Nonlinear Science Series. Cambridge University Press, 2001.
[52] J. A. Platt, S. G. Penny, T. A. Smith, T.-C. Chen, and H. D. Abarbanel. Constraining chaos: En- forcing dynamical invariants in the training of reservoir computers. Chaos: An Interdisciplinary Journal of Nonlinear Science,
2023.
[53] M. I. Rabinovich, Varona, A. Selverston, and H. D. Abarbanel. Dynamical principles in neuroscience. Reviews of Modern Physics, 78(4):1213, 2006.
13

---

[54] J. P. of major depression: the role of glutamate dysfunction on finull-frontte network dynamics. Cerebral Cortex, 27(1):660-679 2017.
A. Roxin, H. S. Mayberg, and A. Compte. A computational model
[55] D. J. Rezende, S. Mohamed, and D. Wierstra. Stochastic backpropagation and approximate inference in deep generative models. In International Conference on Machine Learning. volume 32, pages 1278-1286. PMLR, 2014.
[56] Ritter, M. Schirner, A. R. McIntosh, and V. K. Jirsa. The virtual brain integrates computational modeling and multimodal neuroimaging. Brain Connectivity, 322):22 2013.
[57] N. A. Roy, J. H. Bak, A. Akrami, C. Brody, and J. W. Pillow. Efficient inference for time-varying behavior during learning. In Advances Neural Information Processing systems, volume 31, 2018.
pages
[58] P. Sanz-Leon, S. Knock, M. Woodman, L. Domide, J. Mersmann, A. McIntosh, and V. Jirsa. The Virtual Brain: simulator of primate brain network dynamics. Frontiers in Neuroinformatics, 10, 2013.
[59] J. Schmidhuber, S. Hochreiter, al. Long short-term memory. Neural Compulation, 9()):
1735-1780, 1997.
[60] D. Schmidt, G. Koppe, Z. Monfared. M. Beutelspacher. and D. Durstewitz. Identifying nonlinear dynamical systems with multiple time scales and long-range dependencies. In International Conference on Learning Representations, 2021. 2021.
[61] A. R. Sedler and C. Pandarinath. Ifads-torch: A modular and extensible implementation of latent factor analysis via dynamical systems. arXiv preprint TrXii:2309.0/22/0 2023.
[62] M. Singh, T. Braver, M. Cole, and S. Ching. Estimation and validation of individualized dynamic brain models with resting state fMRI. Neurolmage, 221:117046, 2020.
[63] V. Sip, M. Hashemi, T. Diekscheid, K. Amunts, S. Petkoski, and V. Jirsa. Characterization of regional differences in resting-state fMRI with a data-driven network model of brain dynamics. Science Advances, 9(1 a(14), 2023.
[64] A. C. Smith and E. N. Brown. Estimating state-space model from point process observations. Neural computation, 1555):655-91, 2003.
[65] S. H. Strogatz. Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering. CRC press, 2018
[66] J. Thome, R. Steinbach, J. Grosskreutz, D. Durstewitz, and G. Koppe. Classification of amyotrophic lateral sclerosis by brain volume, connectivity, and network dynamics. Human Brain Mapping, 43(2):681 119999,,1669,,, 2022.
[67] C. van Vreeswijk and H. Sompolinsky. Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science, 274(5293):1724 1724 1726, 1996.
[68] Verzelli, C. Alippi, and L. Livi. Learn to synchronize, synchronize learn. Chaos: Interdisciplinary Journal of Nonlinear Science, 31(8), 2021.
[69] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos. 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, 474(2213):201705) 70844, 2018.
[70] P. R. Vlachas. J. Pathak, B. R. Hunt, T. P. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos. Backpropagation algorithms and reservoir computing recurrent neural networks for the forecasting of complex spatiotemporal dynamics. Neural Networks, 126:191-217, 2020.
[71] P. R. Vlachas, G. Arampatzis, C. Uhler, and P. Koumoutsakos. Multiscale simulations of complex systems by learning their effective dynamics. Nature Machine Intelligence, 4(4): 359-366, 2022.
14

---

[72] J. T. Vogelstein, A. M. Packer, T. A. Machado, T. Sippy, B. Babadi. R. Yuste, and L. Paninski. Fast nonnegative deconvolution for spike train inference from population calcium imaging. Journal of Neurophysiology, 104(6)3391-3774, 2010.
[73] R. Vogt, M. Puelma Touzel, E. Shlizerman, and G. Lajoie. On Lyapunov exponents for RNNs: Understanding information propagation using dynamical systems tools. Frontiers in Applied Mathematics and Staristics, 8:818799, 2022.
[74] X.-J. Wang. Probabilistic decision making by slow reverberation in cortical circuits. Neuron, 36(5):955-968, 2002.
[75] N. Wiener. Extrapolation, interpolation, and smoothing of stationary time series: With engi- neering applications. MIT Press, 1949.
[76] H. R. Wilson and D. Cowan. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal, 12(1):1-24, 1972.
[77] S. N. Wood. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466 102-0):1002-1104, 2010.
[78] G.-R. Wu, N. Colenbier, S. Van Den Bossche, K. Clauw, A. Johri, M. Tandon, and D. Marinazzo. rsHRF; A toolbox for resting-state HRF estimation and deconvolution. Neurolmage, 244: 1 118591, 2021.
[79] E. Yaksi and R. W. Friedrich. Reconstruction of firing rate changes across neuronal populations by temporally deconvolved Ca imaging. Nature Methods, 3(5):377-383, 2006.
[80] L. Yang, X. Sun, Hamzi, H. Owhadi, and Xie. Learning dynamical systems from data: A simple pross-validation perspective, part V: Sparse kernel flows for 132 chaotic dynamical systems. arXiv preprint arXiv:2301.1032/. 2023.
[81] B. M. Yu, J. P. Cunningham, G. Santhanam, S. Ryu, K. Shenoy, and M. Sahani. Gaussian- process factor analysis for fow-dimensional single- analysis of neural population activity. In Advances in Neural Information Processing Systems, volume 21, pages 1881-1888, 2008.
15

---

A Further methodological details
A.1 Deconvolution algorithm in convSSM
Algorithm 1 Deconvolution in convSSM
Input:
X: measured time series as (N x 7) matrix
Parameters:
v analyzing wavelet
mmin: minimum noise level
culy, cul,: edge cutoffs
fr kernel of hemodynamic response function
Output:
Edeconei deconvolved time series as (N x T) matrix
Initialize Xeccooe zeros (N, T)
for i = 1 N do
Ir x]
o, 2t :- VISUSHRINK (I+, v)
if < minn then
σ := Min
end if
Compute Fourier transforms FC) and expectation values of spectral densities TE [ { 12]
Xu = F {
HRF: : F fhr
N₃ = =F 1
S = B
Compute and apply Wiener filter HRF
WE HRF S
7 WR Xx
Transform back to time domain
8. = F-
7
Remove signal edges
T, bbegin + cut, end] = NaN
T. beginn end cut,]] : NaN
Xxccono [i, = Zt
end for
16

---

A.2
Scalability
A
B
C
time per epoch [s]
1.5
3.0
10.0
1.0
5.0
0.5
1.0
1.2
3.0
500
1000
50
100
TR[s]
D
E
M
1.5
7.5
1.0
5.0
2.5
10
10 T
10
500
1000
Figure 4: Training duration per epoch (y-axis) in seconds for different TRs (A), hidden dimensions L (B), latent dimensions (C), time series length (D), and observation dimensions (B). Mean, standard error (SEM) and linear curve fits (gray dashed lines) displayed. The per-epoch-runtinng increases approximately linearly with dimensions L, M,, and N; explained variance R2 = 0.989, R3, = 0.993, and R2 0.996 for linear regressions with predictors LL, 'M", and 'N', respectively. Experiments were performed on a standard notebook with Intel i5-8250U 1,60 GHz CPU and 8GB RAM.
A.3 PLRNNs
The simplest form of the PLRNN is given by
24 = A W (2t-1) +
(1I)
where E RM is latent state vector time t, o(-) max(0. :) is the ReLU activation function, R MMM is an dff-diagonal matrix of connection weights and E R RMX a diagonal matrix containing the autoregressive weights 21.
Neurobiologically motivated, the entries of the latent state 21,t may be interpreted as membrane potentials, the diagonal elements in A as the neurons' individual membrane time constants, and the off-diagonal elements in W as synaptic connections between neurons. The ReLU activation emulates that neurons only start spiking above certain firing threshold.
By adding one hidden layer, obtain the so-called shallow PLRNN
4 = A + W1 (W2 h2) + h
(12)
with 2, E R latent states at time t, <0-) ReLU(-), diagonal matrix A E M rectangular connectivity matrices W, E M and W2 E kgxx and thresholds hi E RM and ho R
A.4
Wiener filter
The Wiener deconvolution filter typically described in the frequency domain, and, for the setting in Equation 6 and Equation 7 is given by
HRF Sk K HRF S + No.
(13)
and returns the estimate & as
(W. (Conv (Coov--((aee (fxe) hrr))),
(14)
where time series denoted with capital letters correspond to the Fourier transformation of time series denoted by lowercase letters, i.e., {X,) = F ((xx)),
Wx is the Wiener filter,
Sg Z [2] is the mean power spectral density of the original signal zt.
17

---

M H is the mean power spectral density of the noise the superscript denotes complex conjugation, and is the inverse Fourier transform.
The noise spectrum N₂ is typically unknown in practice, but can be reliably estimated based on the median estimator on the finest scale wavelet coefficients of z₁. As approximation the power spectrum of the original signal, we use the denoised signal ž, which we obtain by applying the VISUSHRINK algorithm [18], Algorithm 2| the observed signal zy. This approximation works well in practice in absence of knowledge about the true underlying signal.
Algorithm 2 VISUSHRINK algorithm
Input: time series data } of length N; analyzing wavelet ab
1. Apply the discrete wavelet transformation (DWT) to the input data {2z} to obtain the wavelet coefficients Og
+ (
(15)
2. Calculate the MAD of the fine scale coefficients O%, the estimate of the noise level and the universal estimator Ay. Let θ = median(O)) be the median of the fine scale coefficients
and N the length of the time series {2r²)
MAD median( - -
(16)
MAD
(17)
0.6745
λ = N)
(18)
3. Apply hard thresholding the fine scale wavelet coefficients
0, Op.
O Au O 2 Au
θ = 1
(19)
4. Apply the inverse DWT to the tresholded coefficients O to obtain an estimate of the denoised signal {
(a.) (0:)
(20)
Output: estimate of the noise level ē; estimate of denoised signal {
A.5 Additional details on BOLD observation model
The discrete convolution of two time series {fr} and {9i} in its general definition is given by
( * g):
= Σ
Ss 9t-s
(21)
which in case of the decoder model given in Equation 5 translates to
(hr * 2),
Σ ==0
hr 2t-2
(22)
due the finite non-zero value of i.e., the maximum time difference at which previous state 2;-, still influences the current state 2; according to the hrf. Furthermore, due to causality, only past states allowed to influence the current state.
Since the convolution operation as well as matrix multiplication are linear, the order of convolution and matrix multiplication in Equation 5 can be inverted to yield Equation 9 This can be seen more easily if we consider decoder model without noise
B(hr z)e
(23)
18

---

Due to the hrf being the same for all latent dimensions (i.e., a scalar and not function), we interchange the matrix multiplication with B and the convolution with the hrf
B(hrf 2)2 B
E
fr 2--
(24)
M
M Σ fs ==0
(B(rrf x) Σ
Bij
Σ
hr 2t-
Σ
Big 2t-
Σ ==0
hr fa
21-s);
1=1
8=0
j=1
(25)
B(hrf * 2), = (hrf (Bz))r-
(26)
Consequentially, we divide the observations into two parts
(hr x decouv)
+ {x+) Conv (f (1troom)) }, ..hr)
deconv)
(27)
B
(28)
The second part, |Equation 28 is of the same form as Equation 2 and permits the same inversion. Therefore, by computing { [ from {2+} once, we can perform GTF as in the standard SSM.
To include the nuisance artifacts (rr). one has to also swap the order the matrix multiplication and the convolution in the full decoder model (Equation 5). This poses problem since the {re} are not convolved. Our solution deconvolving the {Te) time series as well
Jr: [hr] (Jydbcom)) } doccoov Conv ((re) hr.).-
(29)
(30)
We use this as a simple mathematical trick to incorporate the artifacts into the convolution
Ft = Bhr f * z), + Jrt- = (hr f * Bz), + (hr a (J)deo))) (hr * (Bz + J decony
(31)
Finally, we obtain the relation between the deconvolved time series and the latent states as
deconv
Bz, Jr. deconv
(32)
2 B
(a (anco J..coo
(33)
A.6 Additional information on Algorithm 1
In order to deal with numerical instabilities, additional hyperparameters were introduced. A too low noise level ō determined by the VISUSHRINK [Algriithm 2) can lead to high frequency artifacts, which can be dealt with by defining a lower noise level boundary mmi- Although this is unlikely to occur in empirical (noisy) data, the lower noise level boundary helps to study artificial noise free data. Since the convolution treats the finite signal as periodic, artifacts at the boundaries of the computed deconvolved signal a, may occur. With the hyperparameters cut_ I. cut_r (corresponding to start and end of signal, respectively) one can therefore define absolute cutoff times, either by integer or by floating point values (if, e.g., a cutoff time be defined relative the length of the hr f).
A.7 Performance measures
If the maximum Lyapunov exponent Amax of dynamical system larger than 0, a necessary condition for chaos, nearby trajectories will diverge exponentially. This limits the applicability of n-step ahcad prediction errors (PEs), as conventionally used in machine learning, to evaluate model performance, as even small numerical errors will lead to exponentially growing PEs. In processes in which we can expect chaotic behavior (like neural recordings), we therefore need performance measures which are insensitive to system's 's initial conditions and yet capture the most relevant dynamical properties. Here, we use two established measures [40] to evaluate the DSR, the state space divergence, capturing geometric overlap of (ergodic) distributions in state space, and the power spectrum error, capturing agreement in long-term temporal properties.
On the ALN and LEMON benchmark, we computed these two measures average over 100 trajectories, generated by perturbing the initial state with a small Gaussian noise term (μ 01)..
19

---

A.7.1 Prediction Error PE
The n-step prediction error is given by
PE(n)
|eetn Betnll2:
(34)
N(T n)
i.c. the mean squared error between ground truth data {at} and n-step ahead predictions of the model {
A.7.2 State space divergence Dttap
Given an observed d-dimensional time series {at} of length T and a time series {āt} | with the same dimension/length generated by model, Detsp measures the geometrical overlap of orbits in state space [40.
For low dimensional systems, N < 6, the state space is segregated into KN bins where k is the number of bins per dimension. Each bin given index and we count the number of times n, the time series visited bin i. The relative frequency of visits then obtained by dividing by time series length T
V4g T
(355
Combining these frequencies across all bins in space results in a probability distribution which approximates the state space distribution (the occupation measure) of the underlying dynamical system. The Kullback-Leible divergence of these empirical distributions can then be computed to assess the overlap of both systems in their state space geometry,
Dstsp Σ
Pe log
Ps
(36)
where Pt the relative frequencies of the observed and of the predicted (generated) time series. Note that Dassp is not a metric in the mathematical sense, but divergence that assesses the (dis- Jagreement of probability distributions.
The complexity this binning approach scales exponentially with the observation dimension N and thus becomes intractable for larger N. To compute Dstep in higher-dimensional systems, 6 use Gaussian mixture models (GMMs) with centers (means) 2, and diagonal covariance Σ diag(o?,.. R where a is hyperparameter. The GMMs along the trajectory points given by
T
f(v)
Σ Ny Σ)
(37)
Following (30), the Kullback-Leible divergence of the two GMMs can be computed using a Monte- Carlo sampling approach
K log
(abo) Igen (y1) K
1/T N (yazg) E) 1/T N N E)
Dstap Rd
log
(38)
where K is the number of samples drawn, fobe is the distribution of the observed time series, and fgen is the distribution of the generated time series. The binning and GMM- -bsedd measures correlate strongly in low dimensions (whereas determining the correlation in high dimensions is challenging for the stated reasons).
A.7.3 Power spectrum error DrsE
The state space measure introduced above discards temporal structure the data. To include temporal information in the model evaluation, compare the power spectra of observed and
20

---

generated time series. For each dimension E {1, (1,... .N). the scalar time series {Is,e) is converted into power spectrum density (PSD) {S The components are computed from the Fourier transform of {3
Tr.kl T
g (
S
(39)
The agreement in power spectra is then evaluated terms of the Hellinger Distance (HD) as suggested by [47]. Since HD is probability measure, the computed power spectra have to be normalized by
S.ik Σラ=1 St.j
S
(40)
Using the power spectra of the observed time series Pi,k ((xr,e)) and the generated time series li.k SI, (L,(e the HD is then assessed as
DH
Σ
(41)
k=
The power spectrum error (DpsF) computed by averaging the HDs across all N dimensions the observed system:
N
DPSE
DHA
(42)
Before analysis, the power spectra usually have to be smoothed with a Gaussian kernel of width or to reduce noise. can thus be considered hyperparameter in the evaluation process and was set to in the current work as in [6.
B Details on dynamical systems benchmarks
B.1 Lorenz63 system
The Lorenz63 introduced in was designed describe atmospheric convection based on three dynamic variables. A two-dimensional fluid layer is uniformly warmed from below and cooled from above. It continuous-time dynamical system given by the following set of differential equations:
dr 0(I2 - x1)
(43)
dt
dx2
x (p - [a3) - 122
(44)
dt
drs
T1 ..33
(45)
dt
where is proportional to the rate of convection, to the vertical temperature variation. The constants d, and are system parameters proportional to the Prandtl number, Rayleigh number, and certain physical dimensions of the fluid-layer. For chaotic behavior, a 10, p = 28 and "butterfly attractor" characteristic of the Lorenz system. For each data set, we drew random initial conditions 2g ~ (00, 13x3) and simulated the system using the DynamicalSystems/ Julia package, (115). 10 time steps were saved and used data set (1 training vs. test split) after discarding the first 1, 000 time points to remove transients from the data. To create the benchmark data sets, these trajectories were then convolved with the hr f function. Gaussian white noise drawn from A/ (0, olaxa) was added to all data points.
the horizontal temperature variation, and 2g
are typical settings. These settings produce the so-called
Results the Lorenz63 system in Table 2
21

---

Table 2: Quantitative comparison between standard SSM and convSSM noisy Lorenz63 data
(N conwergee is the number of converged models).
Convolution
a
Obs. model
N Neoaverged
PEz0
Dstep
DPSE
fr
0,01
Standard
80
0.0013 + 0.0001
0.15 + 0.32
0.06 + 0.01
hrf12
0.01
Conv
82
0.0011 0.0001
0.16 0.46
0.06 0.02
fo.5
0.01
Standard
92
0.0117 = 0.0018
0,42 0.55
0.10
hr fo.5
0.01
Conv
90
0.0012 0.0001
0.09 = 0.02
0.06
hr fo.2
0.01
Standard
22
0.0418 + 0.0099
2.97 = 0.75
0.64 0.22
hr fo.2
0.01
Cony
92
0.0023 0.0001
0.29 0.47
0.14 0.03
hrf12
0.1
Standard
86
0.0387 = 0.0002
0.45 0.02
0.09 + 0.01
f1.2
0.1
Conv
86
0.0158 0.0001
0.46 = 0.17
fo.5
0.1
Standard
93
0.0538 0.0016
0.47
0.10
fr.5
0.1
Conv
85
0.0145 0.0001
0.44 0.27
0.09 0.01
Jor2
0.1
Standard
54
0.0893 = 0.0087
2.79 = 0.90
0.58 + 0.22
0.22
hr fo.2
0.1
Conv
84
0.0191 = 0.0002
0.70 = 0.08
0.24 = 0.02
B.2 ALN model
The adaptive linear-nonlinear (ALN) cascade model is a population model of spiking neural networks. The dynamical variables of the ALN model describe the average firing rate and other macroscopic variables of a randomly connected, delay-coupled network of excitatory and inhibitory adaptive exponential integrate-and and-fire neurons (AdEx) with non-linear synaptic currents [3]
We used neurolih to create neural activity generated by an ALN model [2 In neurolib, the firing rate of the excitatory subpopulation of every brain area is used to simulate the BOLD signal via the Balloon-Windkesse BOLD decoder model (Equation 5) as model linking the latent states 2t (i.e., the latent neural activity corresponding to the excitatory firing rates) to the observed BOLD signal. order to create interesting dynamics, certain values were altered from the authors' default settings, 'sigma_ou" = 0 and 5.0. Furthermore, only the first 16 dimensions of the structural connectivity matrix 'Cmat' and the delay matrix 'Dmat' (see explanatory notebook provided by [12] were used for comparability with the empirical LEMON data set.
model (for formula see (T]. As an alternative, we implemented the
neurolib produces simulated (latent) neural activity with sampling rate of 0.1ms. To stay in comparable regime with the fMRI and Lorenz time series, we chose sampling rate 0.5s for the simulated (observed) data. To achieve this without loss of critical information, the neuronal activity as well as the BOLD time series were decimated using 30 point finite impulse response (FIR) filter with Hamming window. Furthermore, the neural activity time series was smoothed with Gaussian kernel with standard deviation o and length 5, and then standardized.
Outliers in DSR measures TTblle 3) were removed using the interquartile range (IQR) method. The IQR method considers values outliers if they are 1.5 IQRs above the third (Qa) or below the first (Q1) quantile (where IQR = Q3 - Q1).
In Table 3 and corresponding Figure 2D, we present the quantitative comparison between the standard SSM, the CSM, the ConVSSM without GTF, and two benchmark conditions. |Figure 5) depicts correlations between measures computed on different subsets of the ALN dataset.
22

---

Table 3: DSR measures evaluated on the ALN data set for the convSSM, the standard SSM, and the COnvSSM trained without generalized teacher forcing by setting O. Measures were evaluated on the ground truth latent space and the noisy observation space on the different created test sets.
metric
ConvSSM
Standard SSM
No GTF (a 0)
White noise
Fixed point
full pseudo-
Dawpp
2.35 + 1.01
2.16 + 0.86
4.97 +: 1.48
3.49
5.15
empirical
DESE
0.28 0.04
0.3 0.04
0.41
0.15
0.79
observational
time series
10-step
0.23
4.14
1.11 0.29
pscudo-
Datep
4.49 + 1.49
3.82 + 1.19
5.37 + 1.98
3.68
4.81
empirical
DrSE
0.22 1:2::0000 0.02
0,22 + 0.03
0.29 + 0.1
0.76
test set
10-step PE
1.46 0.447
1.91 0.51
1.12 +: 0.41
Ground
Datep
2.32 1.53
3.13 1.7
5,19 1: 0.59
2.77
4.89
truth
DESE
0.37 0.11
0.43 +.1
0.54 + 0.23
0.81
test set
10-step PE
1.55 0.2
1.86 + 0.22
1.2 0.19
full pseudo-
Datop
6.17 1.75
7.26 1.53
15.03 2.4
9.84
15.2
empirical
DrSE
0.04
0.55 + 0.03
0.48 + 0.18
0.58
time series
10-step PE
0.66
3.31 + 0.48
2.76 - 0.43
Tatent
pseudo-
Dasp
9.25 3.08
9.98 2.99
15.67 4.4
11.71
15.27
empirical
DrsE
0.41
0.5 0.05
0.42 1.2 0.13
0.56
test set
10-step PR
4.47
3.75
2.85 + 0.73
Ground truth
Dasp
2.13
7.64 2.19
14.99
7.35
14.76
DEEE
0.5 0.08
0.62 0.05
0.56 +.25
0.58
10-step
0.41
3.65
2.81 0.26
23

---

A
D
par
PEn
10
latent
16
latent
latent
0.3
0.2
5
0.1
0.96
0.47
0.69
0
10
20
25
0.1
0.2
0.3 0.4
2
10
observed
observed
observod
D Pys8
0.6
5.0
latent
Iatent
latent
10
0.4
2.5
6
0.58
0
10
15
20
0.0
0.2
0.4
0.8
0.8
2.5
5.0
7.5
observed
obeerved
observed
B
Jtp
Pys
PE b
15
0.75
0.84
0.50
long
Inng
0.25
0
0.00
0
10
15
0.00
0.25
0.50
0.75
0
1
2
short
short
short
D
17
PE
1.0
0.8
long
long
0.6
lang
0.4
0.2
0.8
0
0.0
10
15
0.0
0.2
0.4
0.6
0.8
1.0
0
3
short
short
short
C
Dase
Dysx
PE
full empiricaal data set
full ampirical data set
Bull empirical data set
0.5
d
10
0.4
0.3
10
6
0.2
0.1
0.5
0.51
0.85
0.0
empirical
empirical
emppirical
Figure 5: A: Agreement in DSR measures assessed the observed (x-axis) latent (y-axis) space of the short pseudo-empirical test set (top) and the full pseudo-empirical time series (bottom). Correlations between D step (left), DESE (middle), and PE1o (right) are displayed, respectively. B: Top: Agreement in DSR measures assessed on the pseudo-empirical test set (short) vs. GT test set (long). Bottom: Same for full pseudo-empirical time series (short) vs. GT test set (long). Correlations between Dstsp (left), DrsE (middle), and P (right) are displayed, respectively. C: Correlations between DSR measures between pseudo-empirical test set and full pseudo-empirical time series, same order as in B.
24

---

B.3 LEMON data set
For the LEMON data set, we observed non- stationarity in the data that partly resulted in performance degradation. We therefore only included participants with nearly constant variance over time. To remove non-stationary data sets, assessed the moving average of the variance over time (window size 40 time steps). We then discarded data sets in which the variance changed with time (assessed by computing the correlation with time, with threshold set to r| .16).
The LEMON dataset can be found at httpss//ftpprrree.e////mmmsssssp Mind- BrwBBrrinnn--- Srain-Body-LLMON
de/
B.4 Hyperparameter settings for the different experiments
Table 4: Hyperparameter settings for the experiments conducted. 'Varies' means the respective
hyperparameter was varied in the experiment.
Hyperparameter
Lorenz benchmark
ALN benchmark
LEMON data set
latent_ dim
16
16
gaussian moiee_oieee_eeeeeeee level
0.05
0.05
0.05
optimizer
RADAM
RADAM
RADAM
start Ir
0,001
0.001
0.001
batch_sizz
16
16
model
shPLRNN
cshPLRNN
cshPLRNN
batches batches_per__eeee
50
50
observation_model
Identity
Identity
Regressor
ro__mdee__eguularzzttoon end_Ir
0.0001
0.0
0.0
le-06
le-06
le-06
device
cpu
cpu
cpu
gradient clipping norm
10.0
10.0
0.0
hidden_dii dim
50
Varies
requence_lengee
500
200
200
MAR _raio
0.0
0.0
0.0
obs jos__odel__reguu mrode__regularizztion
0.0
0.0
0.0
epochs
1,000
1,000
1,000
MAR_lambda
0.005
0.0
0.01
weak_tf_alp
0.1
0.2
Varies
min_conv_nois
1.0e-5
5e-06
1.0e-5
train_test_split
50,000
0.75
0.75
TR
Varies
0.5
1.4
cut
0.25
0.25
cut
0.5
0.5
All experiments
system with Xeon Gold 6248 CPU and 768 GB of RAM.
25

---

C Further Details
C.1 Illustration of reconstruction measures
A
B
1 0.007
Data (ground truth + noise)
Chaotic aftractor reconstruction
Limit cycle reconstruction
500
1000
1500
PE
0.011
22
1
500
1000
1500
Figure 6: A. Ground truth Lorenz trajectory sampled with noise (black), a good reconstruction with low Datp (orange) that accurately recovers the attractor, and poor reconstruction with high Dassp (green) that represents the attractor inaccurately, yielding an oscillatory (limit cycle) instead of a chaotic solution. B. Trajectories of systems in A. unfolded in time. The inaccurate reconstruction (top) achieves a lower prediction error (PE) than the accurate reconstruction (bottom), due to trajectory divergence in chaotic systems. This example illustrates that PEs are inadequate to capture the reconstruction of chaotic DS.
C.2 Comparison methods
MINDy: For MINDy (62] we used the implementation at htttts:://gtthub. cww himm/innghfffii /MINDy We trained models for each subject with the settings provided by the authors for IMRI data. To obtain trajectories for calculating PE1o- DESE, and Dxtsp- we used the provided deconvolution function to obtain initial condition in latent space. We iterated the model forward in time for the latent trajectory and then applied the authors' observation function to output a BOLD time series which is compared with the test data.
com/ grsedlee9/1ffe wawr999ffassooooo torch which provides LFADS re-implementation in the deep learning library Pytorch [61]. The default hyperparameters provided are optimized for neural spiking data. We changed the observation model (in the framework this is referred to reconstruction target) to Gaussian, and changed the start and stop learning rates from lr atstrtt = 4 10- 3 (rend = 10 to Itstart 4 10 (rend 1- 10 10-6 6, which improved the fit to data. To obtain trajectories, we iterated the trained models forward in time with initial conditions from the test dataset using the provided model.predict step function.
LFADS: For LFADS (48), we used the Ifads-torch implementation at
rSLDS: For rSLDS [44], we used the implementation at httttt::/
cwm htnm/lindermanlab.
Bam We trained the rSLDS with the Laplace-EM. method with the Structured Mean-Field Posterior, as recommended by the authors, with diagonal_ Gaussian diagonal_Gaussian dynamics and Gaussian emissions. used the same number of latent dimensions for the observations (same as for the cshPLRNN). We determined the rSL.DS training parameters a = 0.9 and K = via grid search over [0.0, 0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 0.9, 1.0) and K E [2, 15] by inferring systems using subset of the data and assessing the performance on the held-out set [6]. To generate trajectories for the model comparison, we first approximated the posterior of the test data and then sampled with model.sample using the approximated posteriors r and z and the test data as prefix input for the model.
For fair comparison,
C.3 History dependence
To illustrate the difference in "history dependence" for latent time series {2t} and convolved counter- parts {xt). illustrative example was created: A latent time series {2e) produced stochastically
26

---

A
Autocorrelation
1.0
1.0
Z
X
0.5
(residual)
0.5
X (residual)
(residual, model)
0.0
0.0
0.5
0
50
100
150
(1
50
100
150
time lag
time lag
Mututal
information
3
3
(residual)
(residual)
2-
(residual, model)
2
10
20
30
40
50
60
10
20
30
40
50
60
time lag
time lag
Figure 7: A: Full (solid) and residual (dashed) average (across dimensions) auto-correlation functions for the latent (left) and observed (right) time series. For the residual auto-correlation, the immediately preceding time step was regressed out. For the dotted curve, cshPLRNN( (x--)) instead 2 was regressed out. B: Same A for the mutual information as a function of time lag.
by using a cshPLRNN equipped with a process noise term cx, i.e.
Zt cshPLRNN((-11) + E,, & NN di))(
The cshPLRNN employed here was trained to become a surrogate model for the Lorenz63 system, and e, is Gaussian white noise with standard deviation o = 0.1 on each dimension. Note that there correlations between noise terms different time points. A time series of length T = 10th was simulated. The corresponding observation time series { was created by convolving { with the hemodynamic response function hr fo.s for TR = 0.5
We then computed auto-correlation functions for the actual and residual time series, where for the latter the linear effect of 2 on 2 (and, likewise, Iz- on Xτ)) was removed (similar to partial auto-correlation). As Figure 7A shows, the autocorrelation of the residual time series drops much faster, instantaneously at first, for the latent states (left) as compared to the observed/convolver variables on the right, illustrating the convolution effect is removed in the model's latent space. is not completely gone if only linear dependencies are removed if the model forwarded-iteree -tterated states (aMPLRNN(2(11) are regressed out instead, the auto-correlation immediately drops to zero for the latent states (dotted lines in Figure 7| left), it should by model definition. Likewise, the (nonlinear) mutual information (Figure 7 shows there are no temporal dependencies left the residual latent series, while still present in the residual observed series, 'empirically' confirming our approach does what it is supposed to do.
C.4 Acronyms
SSM: State space models
RNN: Recurrent neural network
DS: Dynamical systems
rSLDS: Recurrent SLDS
DSR: Dynamical systems reconstruction
LFADS: Latent Factor Analysis via Dynamical
TF: Teacher forcing
Systems
BOLD: Blood oxygenation level dependent
STF: Sparse TF
IMRI: Functional magnetic resonance imaging
GTF: Generalized TF
SLDS): Switching Linear
HRF: Hemodynamic response function
TVB: The Virtual Brain
PLRNN: Piecewise linear RNN
DCM: Dynamic Causal Modeling
shPLRNN: Shallow PLRNN
DL: Deep Learning
cshPLRNN: Clipped shallow PLRNN
ODE: Ordinary Differential Equation
MSE: Mean squared error
27

---

SGD: Stochastic gradient descent
Amax: Maximum Lyapunov exponent GT: Ground truth
EVGP: Exploding and- -ad-ianih gradients problem
LEMON: "Leipzig Study for Mind-Body-
SOTA: State of the art
Emotion Interactions"
TR: Time of repetition
rsfMR:: Resting state fMRI EEG: Electroen- cephalography
PE: Prediction error
DPSE: Hellinger distance/Power spectrum
LSTM: Long Short-Term Memory
ror
MLP: Multi layer perceptron
D sisp Kullback Leibler/ State space divergence
BPTT: Backpropagation through time
ALN: Adaptive nonlinrr cascade model
28

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract accurately states the scope of paper, the introduction and
validation generative model for dynamical systems reconstruction for BOLD time series.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the limitations of approach the end of secciio
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]

---

Justification:
Our paper does not include theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced
statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We provide the data preprocessing steps and hyperparameters used to train our models, the settings to compute the benchmark datasets with open source libraries and we use publicly available dataset. We provide code with an implementation of the convSSM framework.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing verifying the results.
5. Open access to data and code
30

---

reproducible, they Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: All data used in the paper is generated with open source libraries or taken from publicly available datasets, We provide public access our model implementations and training data.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.. cc/
publice/guides/CodeSuttiii ottp://ooooeeaaaisssooooooooooio Psilc for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/nipss.cc/pulli//uu.de htt///CooSSSSiiaaeeeeeee swtsssssssiiiiyyyy for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permirted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All experimental settings included either in the main manuscript or in the Appx. (ubsectioo B.4 as referenced the manuscript.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We include error bars in all figures, where appropriate, report standard devia- tions in all tables, and use statistical tests to back claims made in the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
31

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We include the system specifications of the machines which were used to perform all experiments, see sigure and subsection B.4
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: http:://e..i...eeee wwwuurips.ccpuullice.ooooc cc/publli//thicccuueeeelll..
Answer: [Yes]
The dataset we use containing human subjects was fully anonimized and all subjects gave consent for the the data be used. The data acquisition experiments were carried out in accordance with the Declaration of Helsinki and the experimental protocol was approved by the ethics committee at the medical faculty of University of Leipzig (reference number 154/13-ff), see 4 for details.
Justification:
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- cration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative socictal impacts of the work performed?
Answer: [No]
32

---

Justification: This paper presents novel algorithm in the context of data-driven dynamical systems reconstruction, with applications fields such as neuroscience and psychiatry. We believe the primary consequences of the intended applications, such as characterizing human cognition in health and disease, are positive. We therefore point to these applications. However, we cannot entirely rule out the exploitation of such methods for other, currently not known to us, potentially unethical, purposes that may arise the future.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (c.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out, For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [No]
Justification:
We do not release data and we believe models do not have a high risk for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We cite the libraries used to create our benchmarks and the fMRI dataset used in experimental section.
33

---

Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperssithcone com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re- ee-ackagedd both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Along with the manuscript, upload code for training the COnvSSM model, as well as evaluation and plotting routines.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [No]
Justification: Full experimental details on acquisition of the human publicly available LEMON data set given in 4 and the dataset is properly referenced in the paper. Here, we only include the information relevant to our analyses and replication. The data set consists of far more details, beyond the scope of the paper.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
34

---

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [No]
Justification: In this work, we performed data analysis anonymized publicly available dataset human subjects. The data acquisition was carried out accordance with the Declaration of Helsinki and the data acquisition protocol approved by the ethics committee at the medical faculty the University Leipzig (reference number 154/13- ff) 4 An ethics approval in Germany requires study participants be made aware of foreseeable risks. The data is specifically made available for research purposes. To the best of our knowledge, no additional risks were incurred by the analyses presented this manuscript.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such the institution conducting the review.
35


Reviews:
[{"review_id": "4GjBkwTLuR", "reviewer": "Reviewer_oqqk", "summary": "This paper introduces two techniques, pseudo-inverse and deconvolution, for dynamical system reconstruction. The two techniques are used to help the teacher forcing for the latent sequence $z_t$ so that the learning can be more efficient. Experimental results show the effectiveness of the proposed ConvSSM + GTF compared with alternative variants.", "questions": "/", "limitations": "/", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "* The two techniques seems effective for training an SSM model and dynamical system reconstruction.\n* Experimental results on both synthetic and real-world datasets show that ConvSSM with GTF is better than others.\n* This works might have broader impact to the field of computational neuroscience, since lots of models are based on dynamical system.", "weaknesses": "* From my understanding, there is actually no new model but different ways of training a particular model, although they are called ConvSSM, LinSSM, etc.\n* In most recent dynamical systems, the latent RNN procedure is not a deterministic progress, but latent sequence $z_t$ is usually treated as latent variable with some noise at each time steps, such as Gaussian. It is not clear whether this work can be generalized to them.\n* The presentation should be improved. It is quite hard, at least for me, to get the main idea of the model until I arrive at line 125. Lots of sentences in the abstract and introduction section seems a bit verbose and distracting.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Summary**\n\nWe thank the reviewer for taking the time to comment on and read our manuscript in detail.  We provide new results and figures in an uploaded PDF. A complete list of references can be found in the general rebuttal to all reviewers.\n\n**Weaknesses**\n\n[“*W1*”] \n\nWe see our major contributions as follows: \n\n1) We developed the first of its kind approach that makes efficient large-scale inference of dynamical systems (DS) models from empirical systems observed through convolved signals feasible. This is achieved by reformulating a SSM model such that it a) becomes amenable to recent, highly efficient control-theoretic training methods for DS reconstruction, and b) performing the deconvolution in a computationally very effective way. This leads to an algorithm which efficiently scales up to large amounts of data. In our minds this is of huge practical relevance for the field, as it enables to construct a large battery of single subject-level models in comparatively short time (and, after all, the tremendous success of some recent methods, e.g. structured SSMs like Mamba [1], lies less with the novelty of the ingredients per se, rather than with their high computational efficiency).\n\n2) Apart from methods development, we also, for the first time, explicitly demonstrate that training and DS assessment of the models is indeed valid even on such short time series as provided by fMRI, which we think is important for the community to know. We also show that the Lyapunov spectrum can be retrieved from the trained models, something that is not possible directly on the experimental data itself. \nWe believe that the novelty of our contribution should therefore be judged by how this whole framework and its validation could advance the field by opening up new possibilities, not just by how novel the underlying latent model is.\n\n[“*W2*”] \n\nThe choice to use conventional (deterministic) backpropagation through time (BPTT) as training algorithm is based on previous findings that, surprisingly, in the context of DS reconstruction BPTT outperforms probabilistic training algorithms like those based on variational inference (VI; [2]), even when the data are in fact highly noisy. Nonetheless, we can easily adapt the framework presented here to work within a VI training framework, as suggested in works by [2-4]. The crucial step is to make the convolutional filter amenable to generalized teacher forcing (GTF) which we have made here.\n\n[“*W3*”] \n\nThank you for pointing this out. We will rework the presentation to emphasize the crucial ideas behind the model.", "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": "gxlVJAZm6o", "reviewer": "Reviewer_9P9G", "summary": "The paper proposes a teacher forcing (TF) mechanism for a latent variable model where the dynamics evolve according to a (deterministic) piecewise linear RNN model and the observations are linear projections of the latent space convolved with a filter (in the case of BOLD signals, the form of that filter is known and it’s called hemodynamic response). In the presence of nuisance parameters (such as residual motion, etc.) a linear term is added to the observation to control for their effect and isolate them from dynamics. TF proves useful for robust training of the model and does a better job in recovering long-term and topological features of the chaotic attractors. Two examples in low-d and high-d synthetic data and one example in a real fMRI dataset (called LEMON) are shown.", "questions": "The effect of nuisance variables is considered to be linear, why is this a good assumption?\n\nThe latent space for the Lorenz model is 3 dimensional, why do you need L=50 to model this dataset (or even larger in the next experiment)?\n\nComparisons with many other latent variable models are not shown. Specifically, the neuroscience community has developed a suite of latent variable models with linear, piecewise linear, or nonlinear dynamic models and linear or nonlinear observation models. How do those models compare to the datasets presented here?", "limitations": "See the weaknesses section.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The literature on using control theoretic ideas for model fitting and training has been expanding. This work takes a natural step towards extending TF to latent variable models with deterministic dynamics.\n\nThe writing quality is good, which makes it easy to follow the arguments and contributions.\n\nMost time series models focus on short-term prediction and use MSE to quantify the errors. Although I’m not fully convinced (see the questions and weaknesses) developing metrics that assess the topological properties of attractors is an interesting idea.", "weaknesses": "The work has several weaknesses and limitations detailed below. I’m open to discussion about any of them.\n\n- As the title and narrative of the paper suggest, the work is mostly useful for the fMRI data where a well-established convolution filter exists. Extending this to arbitrary (and learnable) observation filters or nonlinear transformations sounds non-trivial.\n- The model does not incorporate noise in the dynamics space. Given the main motivation for developing the work is to analyze real-world fMRI data the common assumption of the existence of noisy dynamics is a sensible assumption to make and lack of noise in the dynamics can limit the expressivity of the model. Some control theoretic ideas for noisy dynamics are proposed before \\[2,3\\], a discussion on this would be helpful to orient readers.\n- I didn’t quite understand why the authors used a PLRNN in the latent space. To me, it sounds like any generic RNN should be able to take advantage of teacher forcing in the proposed framework. If this is true, please include examples with other generic RNNs (with tanh or other nonlinearities). Otherwise please explain why it’s critical to use PLRNN as the dynamics model.\n- The introduction could be improved. Some relevant models are not introduced or discussed. In my understanding, some of the mean-field theory is developed not to fit the data, but to gain theoretical insight into the population dynamics. In contrast, many latent variable models are developed in the field that are not discussed in the introduction (such as nonlinear LDS, switching linear LDS and its variants, gaussian process latent variable models, data-constrained RNNs, etc.).\n- Some of the results go against the main motivation and claims of the paper. Detailed comments are included below.\n\n**Unknown observation filter**\n\n- For the effect of _hrf_, deconvolution methods have been proposed previously. In your framework, you still run deconvolution (through using convolution in the generative model). It’s unclear whether the history-dependence still holds after controlling for the effect of _hrf_. This statistical dependency needs to be shown to motivate your work.\n- _“hrf with alternative functions if we want to account for filtering in the original signal”_\n\nRelated to this, what about the cases in which the _hrf_ is not known and an a priori observation model does not exist (e.g. a neural network maps the latents to the observations)?\n\n- In addition, the noise spectra are another unknown which is estimated using the VISUSHRIN algorithm (as discussed in the appendix). How are these choices robust to the misspecification? This robustness analysis is important given that it’s almost never the case that we can capture the true noise spectra or _hrf_.\n\n**Results against the claims/motivation**\n\n- _“\\[40\\] proved that for chaotic systems gradient-based training techniques for RNNs will inevitably lead to diverging loss gradients.”_\n\nThis is true under some assumptions, but for specific architectures, this might not hold. In order to motivate your work, it’s crucial to empirically show that the gradients diverge without TF. In fact, your results show that the model without TF does a pretty good job in learning the task and it even outperforms the models with TF on short-term performance measures. Given that the main motivation for the paper is avoiding exploding and vanishing gradients through teacher forcing, it's important to first show that it’s indeed a problem for this specific model (and datasets).\n\n- On the LEMON dataset the latent dimension is set to be equal to the observation dimension. Some parts of the motivation of the paper come from the low dimensionality of the latent space. If the best-performing models on real data are the cases in which the latent and observation dimensions are equal, then the traditional teacher-forcing methods are equally applicable. What do the authors think about this?\n\n**Fig. 2**\n\n- Fig. 2E,F are labeled incorrectly?\n- Unclear what’s shown in Fig. 2A.\n- What are the takeaways from this figure?\n- In Fig. 2C, can you include $\\\\lambda_{max}$ distributions from other models (specifically convNoGTF model)?\n\n**Long-term performance measures**\n\n- The precise definition of performance measures is in the appendix. Can you reference the appendix section in the main text for readers who want to learn more about the definitions of performance measures?\n- It’s important to see how long-term and short-term performance measures are correlated. These results are in the appendix. Can the authors include references in the main text?\n- In the figure, the authors show that short-term and long-term prediction errors are correlated. However, in the LEMON dataset, this trend does not seem to hold. How do the authors explain this?\n- The 10-step prediction error is always better for the models without GTF. This is very surprising. First, it shows that the GTF is not necessary for this model and the models still learn the task without suffering from vanishing or exploring gradients (which goes against the main motivation of the paper).\n- It looks like estimating $D_{PSE}$ in high dimensions is associated with challenges. This is perhaps confounded with dynamics and noise itself too. How trustworthy are these estimates? In other words, how much should we read into these performance measures?\n- On the LEMON dataset (based on Table 1) it looks like that model without GTF achieves lower prediction error but higher $D_{PSE}$ and $D_{stsp}$. What exactly does this mean? Given the correlations reported between long and short-term performance measures don’t we expect a successful model to outperform baseline models in both measures?\n\n**More comments**\n\n_“This is not automatically given for standard RNN”_\n\nNo model that I’m aware of has the capability of precisely predicting neural data if there’s no clear trial structure due to non-stationarity.\n\n_“guiding the training process through optimally chosen control signals – also referred to as teacher forcing (TF) signals”_\n\nAt least the definition that I’m familiar with does not fully coincide with this. In my understanding, there’s no notion of optimality in the control signals and they’re usually driven by a teacher model or the output feedback.\n\n_“Chaotic systems in particular, as typically encountered in neuroscience (e.g., \\[53, 21, 32\\]), pose a severe problem here”_\n\nAnother paper considered the theoretical aspects of training RNNs on chaotic data \\[1\\], please cite and discuss it.\n\n_“but not if it depends on a history of states”_\n\nThe assumption of history dependence is a very common assumption made in most latent variable models.\n\n\\[1\\] Engelken, Rainer, Fred Wolf, and Larry F. Abbott. \"Lyapunov spectra of chaotic recurrent neural networks.\" Physical Review Research 5.4 (2023): 043044.\n\n\\[2\\] Brenner, Manuel, Georgia Koppe, and Daniel Durstewitz. \"Multimodal teacher forcing for reconstructing nonlinear dynamical systems.\" arXiv preprint arXiv:2212.07892 (2022).\n\n\\[3\\] Schimel, Marine, et al. \"iLQR-VAE: control-based learning of input-driven dynamics with applications to neural data.\" bioRxiv (2021): 2021-10.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Summary**\n\nThank you for your detailed review! All refs. are found in the general rebuttal to all revs.\n\n\n**Weaknesses**\n\n[“*W1*”] \nThe conv. filter could in principle also be learnable by parameterizing its length and weights. Here our focus lay on incorporating biological prior knowledge on the HRF to more accurately reconstruct the dynamics of the underlying system (Figs. R1&R2).\n\n [“*W2*”] \nWe chose BPTT as it performs competitively even for highly noisy data [2], focusing on a strong and scalable approach. However,  our framework can easily be adapted to variational inference [2-4]. The key is making the conv. filter amenable to generalized teacher forcing (GTF). Will discuss.\n\n [“*W3*”] \nYes, as shown in [4], the architecture for dynamical systems reconstruction (DSR) is very flexible. Other act. func. also yield good performance (eg $D_{stsp}(ReLU)= 0.29 \\pm 0.47$; $D_{stsp}(tanh)=0.61 \\pm 1.61 $ on Lorenz). The PLRNN was chosen as it consistently yields best performance in low dims. [5], and is math. tractable in the sense that many of its topological properties can be determined semi-analytically [6,7]. \n\n [“*W4*”] \nWe will elaborate on neurosci. models like DCMs [8], SLDS [9], and gaussian process models [10], and their usages apart from DSR [6].    \n\n\n**Unknown observation filter (UF)**\n\n\n[“*UF1*”] \nWe are unsure we fully understand the question. The biophys. assumption is that the obs. signal is generated from the latent process via conv. By incorporating this prior, we achieve more accurate reconstructions (Figs. R1&2). Should we demonstrate that before deconv. $p(x_t|z_t) \\neq p(x_t|z_t…z_{t-\\tau})$, while after deconv. $p(x_t|z_t) = p(x_t|z_t…z_{t-\\tau})$? We are unaware of a model-independent method to show this?\n\n [“*UF2*”] \nSee *W1*.\n\n [“*UF3*”] \nIn Fig. R6, we now show noise level inference in VISUSHRINK for time series with varying noise levels and conv. filters. We also emphasize that the Lorenz exps. (sect. 3.2) demonstrate the model's robustness in inferring GT systems using *default* VISUSHRINK settings, indicating resilience to misspecification.\n\n**Results against the claims (RC)**\n\n[“*RC1*”] \nThe proof in [40] is not related to architecture, but follows from the properties of GD techniques and the chain rule (the same product series of Jacobians occurs in the def. of the max. Lyap. exponent as in the loss derivatives, causing the problem). So it is indeed quite general as long as GD-based training algos. are used. We now illustrate div. (Fig. R4A) and the problem of models that architecturally prevent div. (Fig. R4B).\n\n [“*RC2*”]  \nWe assume by trad. TF the ref. means that during training the obs. are provided as inputs? GTF is designed to achieve a fine balance between obs.-inferred and forward-iterated latent states that optimally controls trajectory and gradient flows. Traditional ‘ad-hoc’ methods do not work here [5, 11]. Also, trad. TF for DSR requires a mechanism to ensure consistent interpret. of model inputs during training and runtime; directly using obs. as inputs necessitates replacing them with predicted obs., causing the latent DS to lose its Markov property due to conv. Finally, the optimal model dims. are not always of equal dim., but depend on dataset and no. of obs. \n\n**Fig. 2**\n\nFig. 2E,F: Thanks, will be adapted. \n\nFig. 2A: Visualizes the reconstruct. performance to provide an intuition on the applied measueres. Will clarify. \n\nFig. 2C: Yes, Fig. R3.\n\n**Long-term performance measures (PM)**\n\n[“*PM1,2*”] \nWe are happy to reference Appx. defs. and figs. in the main text.\n\n [“*PM3*”] \nWe apologize for any confusion and assume the question pertains to Fig. 5 (Appx.). We compare predictions assessed on both short and longer data to show that our PMs can be accurately evaluated on short data. This is shown only on ALN data (since LEMON data is short). Fig. 5B bottom is most relevant for empirical eval, comparing PMs on 500 vs. 5000 data points. All 3 scores correlate with $r \\geq 0.72$, which we find satisfactory for assessing DSR. Will clarify.\n\n [“*PM4*”] \nThere may be a misunderstanding. Our goal with DSR is to achieve agreement in long-term behavior between recons. and true systems [6]. In chaotic systems, short-term forecasting measures like PEs *can be lower for worse reconstructions* [5-6,12], as poor models might capture mean trends or dominant osc. periods that *deviate less from the true signal on average* than accurate models that capture chaos but therefore cause exponential diverg. of trajectories [11] (Fig. R6). We now show that noGTF models frequently *inaccurately* converge to fixed points and osc. (Figs. R3&4).\n\n [“*PM5*”] \nFor $D_{PSE}$, dimensionality is irrelevant since it is performed dimension-wise. For $D_{stsp}$, efficient high-dim. approximations using GMMs are available and used here (see [5,12] for eval). Thus the main question was whether these measures are reliable on short time series.\n\n [“*PM6*”] \nSee *PM4*.\n\n\n**More comments (MC)**\n\n[“*MC1*”] \nAgreed, the point here is not precise prediction, but rather capturing long-term temporal and geometric structure [6]. Will clarify.\n\n[“*MC3*”] \nHappy to cite and discuss.\n\n[“*MC4*”] \nYes, true, but it still constitutes a problem for DSR. A proper DS model must be Markovian by def. to ensure the uniqueness of trajectories [13]. \n\n**Questions**\n\n\n[“*Q1*”] \nWe followed the established practice of modeling fMRI data with linear nuisance vars. but emphasize that our approach can easily accommodate nonlinear effects.\n\n [“*Q2*”] \nThe applied RNN is piecewise linear and more pieces are needed to approx. the nonlinearities of the true eqns.\n\n[“*Q3*”] \nWe specifically chose SOTA models for *DSR*. Linear lat. models cannot produce DS properties like limit cycles and chaos. Most neurosci. models focus on inferring connectivity params., with only a few addressing DSR [eg 14]. While open to comparisons, prior experiences makes us confident that these models are no match in DSR performance.", "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": "J46xUg8kSA", "reviewer": "Reviewer_6CX8", "summary": "The authors introduce a novel algorithm for dynamic system reconstruction (DSR) suited for systems where current observables depend on an entire history of previous states, which notably includes fMRI signals (BOLD signals) and calcium imaging, as both signal are filtered with a response function. The algorithm extends a previous class of methods (teacher forcing - TF) and employs a deconvolution pre-training step where the deconvolved signal is recovered via Wiener deconvolution, which is then used as in the standard TF paradigm. The pre-training step scales linearly with system size enabling efficient application to large datasets. Performance on synthetic dataset show better performance of the novel algorithm as opposed to standard TF.", "questions": "- In Figure 2D, leftmost panel, the conv model does not seem to improve on the quality of the observables. Does this have implications for the generative-mode of the inferred models, i.e. the conv-SSM do not offer better generative performances?\n- How does the standard SSM technique perform on the LEMON dataset? Does the convSSM method provide better insights?", "limitations": "The authors have adequately addressed the limitations of their work", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- *Originality & Significance:* The paper presents a novel algorithm (convSSM) that extends the applicability of DSR to systems where the observables depend on the history of the latent variables. This is particularly significant as this class encompasses relevant experimental settings in neuroscience such as BOLD signals in fMRI or calcium imaging. Moreover the presented algorithm shows improved performance as compared to the previous (sota) techniques.\n- *Quality & Clarity:* The paper is nicely crafted and written, with clearly organized sections. The experiments are designed with careful control to ensure meaningful comparisons.", "weaknesses": "- In the application of the DSR technique to the empirical LEMON dataset the comparison with the competing standard SSM technique is missing. It would have been helpful to identify whether the two techniques differed significantly in a real-world application and whether the novel algorithm offered better insights.\n- Minor point: In Figure 2B/C is not clear whether the plots of $x_1$ and $x_3$ belong to panel B or not as panel C sits in between. Panel labels have generally exaggerated font size.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Summary**\n\nWe thank the reviewer for taking the time to comment on and read our manuscript in detail, as well as for this supportive and positive evaluation. We provide new results and figures in an uploaded PDF.\n\n**Weaknesses**\n\n[“*W1*”] \n\nWe apologize for the error in the table labeling. The standard is indeed found in Table 1, where it is referred to as LinSSM (for linear SSM). We will correct this. \n\nOverall, the convSSM outperforms the standard SSM on the Lorenz63 benchmark and on the ALN data set in latent space. It performs comparably to the standard SSM on the LEMON data where the latent space is inaccessible. We therefore dissected performance contributions of the convSSM in detail by adding an additional (simple to visualize 2D) benchmark, the Van der Pol nonlinear oscillator (VdP), which we could specifically adjust to produce oscillations in a frequency range consistent with the empirical fMRI data. We show that the convSSM outperforms the standard SSM in latent space (where the VdP lives in; Fig. R1), even if we deconvolve the standard SSM in latent space. This demonstrates that the standard SSM, while a very powerful tool for dynamical system reconstruction when trained by GTF, is not able to reproduce the true underlying process without bias. We also show that this can result in biased estimates of maximum Lyapunov exponents (as now illustrated on the Lorenz system where the max. Lyapunov exponent is known, Fig. R2). \n\nIn a sense, compared to the standard SSM trained with GTF, the convSSM implements a biological prior, which enables it to often capture strongly low-pass filtered (convolved) processes more efficiently, and to provide a more accurate description of the underlying latent dynamics. In that sense, yes, we believe the novel algorithm offers better mechanistic insights, as the main findings with respect to the empirical data is that we predominantly find chaotic systems.\n\n[“*W2/Minor point*”] \n\nThank you for pointing this out. We will adjust Figure 2B/C to make it more clear, and reduce the font size in the panel labels.\n\n\n**Questions**\n\n[“*Q1*”] \n\nAs pointed out in the response to *W1*, the latent dynamics is what we want to analyze in terms of computational mechanisms, and therefore care about reproducing accurately. The standard SSM trained by GTF produces biased estimates of this system (as now further illustrated in Figs. R1& R2 in the uploaded PDF). Therefore, while yes, the prediction performance is comparable, the analysis of the generative mechanisms will be flawed. \n\n\n[“*Q2*”] \n\nPlease see our response to *W1*, where we address these questions in detail.", "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": "ct8JEKbyLN", "reviewer": "Reviewer_qvr3", "summary": "This paper proposes a SSM-based DSR algorithm, convolution SSM model (convSSM), to recover the underlying systems. To prevent exploding gradients, the model is trained by SGD with generalized teacher forcing (GTF), for both (pseudo-)invertible and signal convolved decoders. After validating and studying the proposed method by simulations (Lorenz63 and ALN), they apply it to a experimental fMRI dataset (LEMON). Overall, the paper provides a scalable model for dynamic system reconstruction, which can be helpful for both inference and prediction task.", "questions": "The paper is clearly written. But the major question is the same as weakness, why it is necessary put correlation modeling into the observation equation, and hence make it more difficult to inference?", "limitations": "Using AR(1) for system equation and simultaneously making current ovservation explicitly depending on the whole latent history may make inference difficult.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The proposed method (convSSM trained via SGD+GTF) can be efficiently scales with model size and convolution length\n2. The model selection strategy based on short empirical time series are practically useful.\n3. The proposed models can reliably extract key dynamical system features.", "weaknesses": "From the model and model inference perspective, the main contribution for the paper is training convSSM with GTF for signal convolution model, by exploiting linearity of Wiener deconvolution. The proposed method is mainly to solve the difficulty that the current observation depends on the whole series latent states, which is common in real application, especially in neuroscience. \n\nHowever, the evolution of latent states described by the system equation can capture some correlation to previous states, i.e., the current latent states depends on the history of previous states. The modeling of observation equation as a function of the whole latent history further captures the remaining correlation structure, which is missed by the AR(1) assumption in system equation. But doing this makes inference harder. It would great to show the necessity of complicated modeling on observation equation, rather than put all correlation modeling into latent space (such as using GP, which is clearer and usually easier for inference), based on some evaluation criter", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Summary**\n\nWe thank the reviewer for taking the time to comment on and read our manuscript in detail, as well as for this supportive and positive evaluation. We provide new results and figures in an uploaded PDF, and replies to the questions and comments below. A full list of references is listed in the general rebuttal to all reviewers.\n\n\n**Weaknesses**\n\n[“*W1*”] \n\nThank you for your detailed review. We would like to clarify a few points regarding our model, which may not have been fully articulated in the paper. First, the AR(1) component, representing the Markov property in our latent model, is crucial because it ensures that the model constitutes a true dynamical system. If the dynamics in state space were not Markovian, it would mean the state space is not complete (lacking dimensions/ variables) and trajectories not uniquely resolved [13, 15]. For a true state space, the future must be determined solely by the current position, and the same state cannot be associated with different futures. This is where the delay embedding theorems come in in partially observed systems [16].  \nAn AR(1) process is therefore essential to recover the underlying dynamical system, and be able to analyze its dynamics unambiguously, such as the maximum Lyapunov exponent (e.g., Fig. R2). \n\nSecond, our model is designed to reconstruct the generative dynamical system underlying neural signals measured through functional magnetic resonance imaging (fMRI). These neural signals are not directly accessible due to the delay in the delivery of nutrients and oxygen to cells, which results in a convolution with the hemodynamic response function (HRF).  To address this, we convolve the underlying neural state with the HRF, a standard practice in the fMRI research community. This approach allows us to accurately model the neural dynamics by separating them from the neurovascular response. \nBy incorporating the HRF into our observation model, we disentangle the neural state and its dynamics from the neurovascular mechanics, which are not the focus of interest. We will further clarify this aspect in the updated version of our paper.\n\n\n**Questions**\n\n[“*Q1*”] \n\nThe necessity of having an observation model that relates a history of states to the observed measurement stems from the actual biophysical properties of the blood-oxygenation-level-dependent (BOLD) signal measured using magnetic resonance imaging. As alluded to in response to *W1* above, when neurons in the brain are active they require oxygen. The fMRI signal measures changes in blood oxygenation in response to neural activity (the BOLD signal). To provide oxygen and nutrients, there is therefore a change in blood flow, however, it comes at a delay and with a characteristic shape, accounted for by the HRF. So the reason we have to convolve with the HRF is simply due to the biological properties that give rise to the data. We therefore agree that it makes inference more difficult, however, it is necessary to capture the underlying neural process in an unbiased manner (as now also illustrated in Figs. R1 & R2).\n\n**Limitations**\n\n[“*L1*”]\n\nAs mentioned earlier, while incorporating additional lags into the latent model might simplify inference from a function approximation perspective, it would not result in a true dynamical system (see also response to *W1*). Our primary interest lies in accurately recovering the underlying system dynamics, which requires preserving the Markov property inherent to a true dynamical system. This focus allows us to explore the system's fundamental dynamics, as highlighted in [6]. Our approach ensures that we maintain the integrity of the system recovery, which is crucial for our analysis and objectives.", "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:
**General response**

We thank all reviewers for their positive and supportive feedback, for taking the time to review our work, and for providing many helpful comments and suggestions, which we address in detail below. We have prepared a PDF file with additional material and results.

Specifically, we provide additional analyses that demonstrate once more the superiority of the introduced convSSM model over the standard SSM, both trained with generalized teacher forcing (GTF). Our results now illustrate that the convSSM produces less bias in dynamical systems reconstruction, resulting in improved performance measures (Fig. R1), as well as more accurate estimation of dynamical systems phenomena (Fig. R2).

We also explicitly show that training without GTF leads to higher gradient explosions (Fig. R4A), or alternatively, when gradient explosions are avoided through architectural adjustments, it results in more bias in model estimates due to the inability to accurately reconstruct chaotic phenomena (Fig. R4B & Fig. R3). Additionally, we illustrate the robustness of the VISHUSHRINK algorithm to misspecification (Fig. R5) and explain why prediction errors (e.g., $PE_{10}$) are not an appropriate performance measure for dynamical systems reconstruction (Fig. R6). Finally, we provide an update to Fig. 2C (Fig. R3).

We hope these additions address the reviewers' main questions and concerns. Due to the character limitations, we had to significantly condense our responses at times. We are happy to provide more detailed answers to any further questions upon request.



**References**

1) Gu & Dao. (2023). Mamba: Linear-time sequence modeling with selective state spaces. *arXiv preprint* arXiv:2312.00752

2) Brenner et al. Tractable dendritic RNNs for reconstructing nonlinear dynamical systems. In *Proc. 39th International Conference on Machine Learning* (eds. Chaudhuri, K. et al.) 2292–2320 (PMLR, 2022).

3) Kramer et al. Reconstructing nonlinear dynamical systems from multi-modal time series. In *Proc. 39th International Conference on Machine Learning* (eds Chaudhuri, K. et al.) 11613–11633 (PMLR, 2022).

4) Brenner et al. Multimodal teacher forcing for reconstructing nonlinear dynamical systems. In *Proc. 41st International Conference on Machine Learning* (2024).

5) Hess et al. Generalized teacher forcing for learning chaotic dynamics. In *Proc. 40th International Conference on Machine Learning* (eds Krause, A. et al.) 13017–13049 (PMLR, 2023). 

6) Durstewitz et al. (2023). Reconstructing computational system dynamics from neural data with recurrent neural networks. *Nature Reviews Neuroscience, 24(11)*, 693-710.

7) Eisenmann et al. (2024). Bifurcations and loss jumps in RNN training. Advances in *Neural Information Processing Systems2*, 36.

8) Friston et al. (2003). Dynamic causal modelling. *Neuroimage 19*, 1273–1302 

9) Ghahramani & Hinton. (2000). Variational learning for switching state-space models. *Neural Computation 12*, 831–864 

10) Yu et al.  (2009). Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. *Journal of Neurophysiology 102*, 614–635

11) Mikhaeil et al. On the difficulty of learning chaotic dynamics with RNNs. In *Proc. 35th Conference on Neural Information Processing Systems* (eds. Koyejo, S. et al.) (Curran Associates, Inc., 2022).

12) Koppe et al. (2019). Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI. *PLoS Computational Biology, 15(8)*, e1007263.

13) Perko, L (2013). Differential equations and dynamical systems (Vol. 7). Springer Science & Business Media.

14)  Singh et al.  (2020). Estimation and validation of individualized 460 dynamic brain models with resting state fmri. *Neuroimage, 221*:117046. 

15) Strogatz, SH. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (CRC, 2018).

16) Sauer et al. (1991). Embedology. *Journal of Statistical Physics 65*, 579–616

Meta-review: The reviewers were generally positive about the contribution, the proposed method was found well presented and sound. The reviewers ackonwledged the novelty and interest of the proposed approach to dynamical system reconstruction combining convolution state-space models with generalized teacher forcing. The experimental validation was also found appropriate. Specific questions were on the novelty and generality of the contribution, and on the deterministic nature of the dynamics considered in the latent space. It was also recommended to provide a comparison with respect to models based on simpler assumptions. 

The rebuttal provided further experimental assessment on the benefit of the proposed approach over a range of benchmarks (e.g. to assess the advantages of training with GTF as compared to SGD). The authors provided further arguments to position the contributions with respect to the state-of-the-art, in particular with respect to models implementing stochastic latent dynamics. The discussion phase was rich of interactions and allowed the authors to further clarify the contributions and provide additional experimental material. Notably, the discussion identified further investigation directions, for instance relating the Markov properties in latent and observations spaces, which were found of substantial interest. 

As a result, several reviewers raised the score to recommend acceptance to this work. Overall, this paper makes interesting contributions to the relevant domain of state-space modeling, and the proposed approach may have impact beyond the neuroimaging domain. For all these reason, the paper is accepted for presentation to NeurIPS 2024.
Decision: Accept (poster)