A scalable generative model for dynamical system reconstruction from neuroimaging data

Authors: eric  volkmann
alena  brändle
daniel  durstewitz
georgia  koppe
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 09:25:52

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.

Main paper: 
 introduction: 
 Models and theories based on dynamical systems (DS) concepts have a long tradition in computational neuroscience in accounting for physiological phenomena and computational processes of the brain [74, 53, 33, 23] . Constructing such models from first principles (biophysics) is time-consuming and hard, and utilizing them to 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 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 diagnosing or predicting brain dysfunction and clinical characteristics based on DS features [66] , or for designing personalized therapies [54] .In computational neuroscience, multiple approaches to infer large-scale brain dynamics have been advanced over the past decades (e.g., [17, 46, 58, 32] ). Pioneering work introduced latent linear DS, including Gaussian process latent variable models [64, 81] , 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 to modeling a larger range of dynamical phenomena by combining several linear (or affine) DS, with a switching mechanism that selects the active system at each moment [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 on mean field neural population models, including neural mass [76, 5] 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, 58] . TVB incorporates biologically realistic brain connectivity into neural field models to generate simulations of large-scale brain activity. Alternatively, Dynamic Causal Modeling (DCM) describes a 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., [16, 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 (fMRI), only recently a shift in focus has introduced data-driven, deep-learning 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, 63] space have been proposed, partly using available structural information and hierarchical inference approaches [63] .A variety of Deep Learning (DL)-based models for approximating the generative dynamical processes underlying observed time series has been put forward in recent years ( [8, 69, 49, 40, 60, 42, 6, 71, 31, 13, 80] ; see also [24] for an overview). These include approaches which approximate a system's vector field, e.g., through a library of basis functions and penalized regression as in SINDy [8] , or through deep neural networks and neural Ordinary Differential Equations (ODE) [14, 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., [40, 29, 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 in 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-term 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 [51, 1, 68, 2] , by guiding the training process through optimally chosen control signals -modern variations of classical teacher forcing (TF) -that prevent exploding gradients [47, 6, 31, 7, 39] .Chaotic dynamics in particular, as typical for neural systems (e.g., [67, 22, 36, 26] ), poses a severe problem here as trajectories and hence loss gradients inevitably diverge due to the presence of a positive Lyapunov exponent [47] . Recent amendments of TF protocols, including sparse TF (STF; [47] ) and generalized TF (GTF; [31] ), keep trajectories and gradients in check by 'weakly synchronizing' them with the observed signals.Nonlinear SSMs distinguish between an underlying latent process that governs the dynamics of a system'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 latent space. This inversion is fairly straightforward when ones assumes that the current measurement depends solely on the present latent state, i.e, when x t = f (z t ), but not if it depends on a history of states {z t , z t-1 , . . . , z t-τ }. In practice, unfortunately, this assumption is often violated due to a signal's filtering properties. For instance, blood oxygenation level dependent (BOLD) signals, as assessed via fMRI, are only an indirect measurement of neural activity, with the hemodynamic response function (hrf ) broadly smearing out the signal across time. Each measurement is therefore a function of a history of past neural (latent) states whose dynamics we wish to infer [28, 78] . 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.Here we rectify this issue by developing a 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 [47, 31] . 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, as 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 on short time series accurately predict a system'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. 
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 [6, 31] . 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; [31] ). In the cshPLRNN, the temporal evolution of the (latent) system state z t ∈ R M is expressed asEQUATIONwhereϕ(•) = max(0, •) is an elementwise ReLU activation function, W 1 ∈ R M ×L , W 2 ∈ R L×M are connection weights, A ∈ R M ×Mis a diagonal matrix of autoregressive weights, and h 1 ∈ R M and h 2 ∈ R L are bias vectors. Its trajectories {z t } ∈ R M ×T will be bounded if the absolute eigenvalues of A are smaller than 1 [31] . The Markov property of the latent model is crucial to ensure it formally constitutes a DS with complete state space [50, 65] . Finally, Equation 1 can easily be extended to incorporate the effect of external inputs, such as experimental stimuli, by adding Cs t (with s t ∈ R K representing an input vector and C ∈ R M ×K its effect on the latent dynamics).However, here we consider input-free data from resting state experiments.In [31] , the latent state z t at each time point is assumed to be related to the actual observation x t ∈ R N by a linear (Gaussian) decoder modelEQUATIONreferred to simply as 'standard SSM' in the following. Here, B ∈ R N ×M is a matrix of regression weights, and η t describes Gaussian observation noise with diagonal covariance matrix Γ. A conventional mean squared error (MSE) type loss functionL = t L t = t ∥x t -x t ∥ 2 2between the generated (predicted) {x t } and the observed {x t } 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] .A fundamental issue in training such systems by SGD is the well-known 'exploding-and-vanishing 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 [26] ). 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-hyperbolic) solutions. A theoretically well founded approach that does not limit a system's dynamical expressivity, which we will adopt here, is GTF, 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 zt is computed as a linear interpolation between the PLRNN generated state z t = PLRNN(z t-1 ) and a data-inferred state d t that serves as a control signal [19] , i.e., zt := (1 -α) • z t + α • d t , α ∈ [0,(3) There is a theoretically optimal choice for α that can be approximated concurrently whilst training through a specifically designed annealing protocol [31] , but more simply α may just be determined by grid search (as done here). Control signals d t are obtained by inverting the decoder model (Equation 2). Since in general M ̸ = N , the matrix inverse of B ∈ R N ×M does not exist and is approximated by the Moore-Penrose (pseudo-) inverse B + :EQUATION) To keep the gradients on track, the interpolation is performed at each time step before applying the cshPLRNN mapping (Equation 1). These control signals are turned off during actual data generation by the model (i.e., in a test phase), where it runs completely autonomously.Decoder model for convolved signals GTF (and similar techniques like STF; [47] ) are powerful state-of-the-art (SOTA) tools for controlling gradients, especially in the context of DSR. However, they require a (pseudo-)invertible observation model for producing adequate control signals. Empirically, there are many situations where this requirement is not met. For instance, in BOLD time series the observed signal is a highly filtered and strongly smoothed version of the underlying neuronal process that we would like to recover [9, 10, 28] . This fairly complex hemodynamic process is often modeled by the hrf [28] .A decoder model that relates the neuronal processes given as latent time series {z t } to measured BOLD time series {x t } may be formulated as in [40] ,EQUATIONwith regression coefficient matrices B ∈ R N ×M and J ∈ R N ×P , nuisance variables r t ∈ R P (such as movement or respiratory artifacts) and a Gaussian observation noise term η t (with usually diagonal covariance Γ ∈ R N ×N ). Here, * denotes the convolution operation and z is a history of states z t-τ :t , the length of which depends on the observed sampling rate, commonly referred to as time of repetition (TR). The discrete hrf sequence is computed by evaluating the canonical hrf at the observed TR [78] . We will denote the hrf response for a given TR by hrf T R .By incorporating the hrf into the observation model, we disentangle the neural state and its dynamics -the processes of interest -from the neurovascular mechanics (or any filtering at 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 xt ) do not simply depend on the current state z t , but -due to the convolution -on a set of states across several previous time steps. We can thus not compute the control signal d t through straightforward decoder inversion anymore, but require a new type of inversion algorithm.Wiener deconvolution Following [78] , 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.4 for further details. Given an observed noisy signal {x t }, composed of the signal of interest {z t } convolved with a known impulse response hrf plus some noise term η t (distribution unknown, Wiener is optimal for Gaussian distribution),x t = (hrf * z) t + η t , (6) the Wiener deconvolution provides the estimate ẑt of the unknown signal z t through least-MSE estimation. Defining Conv -1 (•, hrf ) as the Wiener deconvolution operator, we can writeEQUATIONInversion of BOLD decoder model Using the notation introduced above, we can write the inversion to obtain control signals asEQUATIONwhere B + (x t -J r t ) is the time series that needs to be deconvolved. Note that this approach is quite general and we can simply exchange the hrf with alternative functions if we want to account for filtering in the original signal. As stated, since B and J are matrices of learnable parameters updated during training, we would need to perform this deconvolution step at 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 8asEQUATIONWith x deconv t and r deconv t denoting the respective deconvolved time series, this can be written asEQUATION)This now is computationally much more efficient, as we need to perform the deconvolution only once prior to 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 1 and Equation 5) trained with GTF and SGD as 'convSSM'. The full inversion algorithm is provided in Algorithm 1 with additional information given in Appx. A.6. Key components of SGD+GTF training are illustrated in Figure 1 . 
performance measures: 
 In DS theory in general, and DSR more specifically, we are 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 to short-term n-step ahead PEs, P E n , we 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 error (PSE), D P SE , in the following [47] , and 2. the Kullback Leibler divergence between the empirical and model-generated trajectories across state space, D stsp , measuring the overlap in attractor geometries [40] .To obtain a reference value for D stsp , we further included two references in which we assessed D stsp when all mass is centered on the expectation value (similar to a 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 a 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, see Appx. A.7.A Training: Gradient flowz t z t+1 z t-1 z t-τ B Training:Teacher forcing x deconv t-1 r deconv t-1 d t-1 z t-1 zt-1 z t α 1 - α P L R N N ∂zt ∂ zt-1 ∂ zt-1 ∂zt-1 ∂zt ∂zt-1 C (hrf * z) t xt r t L t x t hrf * ∂L ∂ xt ∂ xt ∂(hrf * z) t ∂(hrf * z) t ∂zt ∂(hrf * z) t ∂zt-1 ∂(hrf * z) t ∂zt-τ ∂zt+1 ∂zt ∂zt ∂zt-1 ∂zt-1 ∂zt-τAs a well-established and popular benchmark for a chaotic system, we first performed numerical experiments with the famous Lorenz63 system. The Lorenz63 is a 3-dimensional system introduced in [45] to describe atmospheric convection, and exhibits chaotic behaviour in the chosen regime (see Appx. B.1). To mimic BOLD observations, we generated 100 standardized chaotic Lorenz63 trajectories, convolved them with hrf T R functions at different sampling rates TR ∈ {0.2s, 0.5s, 1.2s} (Figure 2B ), and added Gaussian noise with standard deviation σ ∈ {0.01, 0.1}. This resulted in 6 benchmark settings with different levels of signal degradation by convolution and noise. Each data set was divided into a training and a test set of T = 5 • 10 4 time steps each.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 via SGD and no GTF, the standard SSM trained via SGD+GTF, and MINDy, a recently published method for DSR in fMRI [62] . convSSM and standard SSMs were trained with the shPLRNN, with M = 3, L = 50, and α = .1 (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 a variety of simulated and real-world data sets and is considered to be a SOTA model in the field [6, 31] . The performance measures D stsp , D P SE , and P E 20 were assessed on the test sets after training for 1, 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 on unconvolved, noiseless trajectories from the Lorenz63 system.The convSSM significantly outperformed all other methods, including the standard SSM in almost all cases, with the performance gap increasing with decreasing TR (see Appx. Table 2 for performance, and Figure 2A for example reconstructions, providing an intuition on how to interpret D stsp ). 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 empirical data is how well they scale with model size and convolution filter length. To assess this, we collected trajectories of length T = 10 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 convSSMs 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, λ max . λ max measures how quickly trajectories starting from nearby points in a system's state space converge or diverge with time. If λ max > 0, trajectories will exponentially diverge and the system, if bounded, will exhibit chaos. We show that we can successfully recover λ max (known for the Lorenz63 system; Figure 2C ) even from models trained on these just short 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 to accurately recover the underlying neural latent dynamics, and b) our measures (D stsp , D P SE ) 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, 11] for simulating whole brain (neural) activity. 100 data sets of length T = 10, 000 were simulated from this model using neurolib [12] , sampled at 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 s to mimic an experimentally realistic scenario (see Appx. B.2 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 fMRI 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 5, 001-10, 000 of the full simulation set) the ground truth (GT) test set (which would not be accessible in a real experiment). DSR performance was assessed on both the observed {x t } and latent {z t } time series, evaluated for a) the short pseudo-empirical test set of length 125, b) the full pseudo-empirical time series (i.e., of length 500), and c) the GT test set of length 5000 (which also assesses dynamics on the limit set and does not contain transients anymore).Figure 2D shows histograms over D stsp (on the full pseudo-empirical time series) for the convSSM, the standard SSM, and the benchmark conditions. The convSSM significantly outperformed the standard SSM in latent space (rank-sum test Z = 11.50, p ≤ .001), demonstrating an improved 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 fMRI. Moreover, the proposed performance measures (D stsp , D P SE ) 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 dynamics space (Figure 2F and Appx. Figure 5 ). Second, DSR assessed on the pseudo-empirical time series (either 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 2E and Appx. Figure 5 ).We finally tested convSSM on empirical data, for which we chose the LEMON study ('Leipzig Study for Mind-Body-Emotion Interactions') as a publicly available data set. This data set was collected at the Max-Planck-Institute 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 train = 489) and test (T test = 163) set, respectively. Data from participants with non-stable variance were discarded (i.e., non-stationary 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), α = .1, and hidden dimension L = 50 (where L and M refer to the dimensions of the connectivity matrices W 1 , W 2 in the cshPLRNN, Equation 1). Latent dimension and α were determined via grid search, by inferring systems using a 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 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 D ststp , D P SE , and P E n , we also assessed the trained models' maximum Lyapunov exponents, λ max , analyzed how reliably these can be inferred, and whether they distinguish between subjects. Note that obtaining an 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-generated trajectories, where here we used an algorithm proposed in [73] ).The DSR results are shown in Table 1 . We obtained successful reconstructions on average with a mean D stsp of 2.73, better than all other models (Figure 3A and D ). Interestingly, most recovered systems were characterized by a positive maximal Lyapunov exponent λ max (Figure 3E ), indicating the presence of chaotic attractors in these data (consistent with previous observations, [40, 36, 41] ). Moreover, λ max 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, p < .001.Table 1 : DSR measures evaluated for the convSSM, standard SSM, convSSM trained without GTF, as well as MINDy [62] , rSLDS [44] and LFADS [14]  
conclusions: 
 Methods for producing generative models of the underlying dynamics from time series observations is a rapidly expanding research field [24] . Current SOTA models for this purpose rely on control theoretically motivated training techniques like STF [47] and GTF [31] , but these require some means to generate from the actual observations a TF signal in the model'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 2+ 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 is 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 at large scale to understand inter-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 [31] . 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 on convolution filters like the hrf , alternatively we can parameterize the filter weights within the observation model, making them learnable through BPTT, with filter length either as a 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 an 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. [7] 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 9 to generate a TF signal which steers a 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 frameworks 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 mechanisms and may thus not be as suited to address specific questions relating to pharmacological or receptor mechanisms beyond functional-dynamical implications. Moreover, currently open questions are how to best deal with non-stationarity in the data, how to efficiently combine data from many subjects, and how trained models generalize to out-of-domain data. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary understand the results?Answer: [Yes]Justification: All experimental settings are included either in the main manuscript or in the Appx. subsection B.4, as referenced in 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 to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
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 deviations in all tables, and use statistical tests to back up 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 error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• 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 or the standard error of the mean. Justification: We include the system specifications of the machines which were used to perform all experiments, see Figure 4 and subsection B.4.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of workers or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as 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 it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [No] Justification: This paper presents a novel algorithm in the context of data-driven dynamical systems reconstruction, with applications to 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 in 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 no 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 (e.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 not tied to particular applications, let alone deployments. However, if there is a 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 as 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 or 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 ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [No] Justification: We do not release data and we believe our models do not have a high risk for misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 to 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a 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 our experimental section.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 an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [No] Justification: In this work, we performed data analysis on an anonymized publicly available dataset of human subjects. The data acquisition was carried out in accordance with the Declaration of Helsinki and the data acquisition protocol was approved by the ethics committee at the medical faculty of the University of Leipzig (reference number 154/13ff) [4] . An ethics approval in Germany requires study participants to 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 in this manuscript.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on 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 as the institution conducting the review. 
 
 Appendix: Code for the convSSM is available at https://github.com/humml-lab/GTF-ConvSSM.A Further methodological details X deconv : deconvolved time series as (N × T ) matrix Initialize X deconv := zeros (N, T ) for i = 1 to N do x t := X [i, :] σ, zt := VISUSHRINK(x t , ψ) if σ < σmin then σ := σmin end if Compute Fourier transforms F{•} and expectation values of spectral densities E |F {•} | 2 X k := F {x t } HRF k := F {hrf t } N k := E |F {η t } | 2 S k := E |F {z t } | 2W k := HRF * k •S k |HRF k | 2 •S k +N k Zk := W k • X kThe simplest form of the PLRNN is given byEQUATIONwhere z t ∈ R M is a latent state vector at time t, Φ(•) = max(0, •) is the ReLU activation function, W ∈ R M ×M is an off-diagonal matrix of connection weights and A ∈ R M ×M a diagonal matrix containing the autoregressive weights [21] .Neurobiologically motivated, the entries of the latent state z i,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 a certain firing threshold.By adding one hidden layer, we obtain the so-called shallow PLRNNz t = Az t-1 + W 1 ϕ (W 2 z t-1 + h 2 ) + h 1 (12) with z t ∈ R M latent states at time t, ϕ(•) = ReLU(•), diagonal matrix A ∈ R M ×M , rectangular connectivity matrices W 1 ∈ R M ×L and W 2 ∈ R L×M and thresholds h 1 ∈ R M and h 2 ∈ R L .The Wiener deconvolution filter is typically described in the frequency domain, and, for the setting in Equation 6 and Equation 7, is given byEQUATIONand returns the estimate ẑt asEQUATION• where time series denoted with capital letters correspond to the Fourier transformation of time series denoted by lowercase letters, i.e.,{X k } = F ({x t }), • W k is the Wiener filter, • S k = E[|Z k | 2 ]is the mean power spectral density of the original signal z t ,• N k = E[|H k | 2 ]is the mean power spectral density of the noise η t , • the superscript * denotes complex conjugation, and • F -1 is the inverse Fourier transform.The noise spectrum N k is typically unknown in practice, but can be reliably estimated based on the median estimator on the finest scale wavelet coefficients of x t . As approximation to the power spectrum of the original signal, we use the denoised signal xt which we obtain by applying the VISUSHRINK algorithm [18] , Algorithm 2, to the observed signal x t . This approximation works well in practice in absence of knowledge about the true underlying signal.Input: time series data {x t } of length N ; analyzing wavelet ψ 1. Apply the discrete wavelet transformation (DWT) to the input data {x t } to obtain the wavelet coefficients Θ 0EQUATION2. Calculate the MAD of the fine scale coefficients Θ 0 t , the estimate of the noise level σ and the universal estimator λ U . Let Θ = median(Θ 0 t ) be the median of the fine scale coefficients and N the length of the time seriesEQUATIONEQUATION3. Apply hard thresholding to the fine scale wavelet coefficientsEQUATION4. Apply the inverse DWT to the tresholded coefficients Θ 1 t to obtain an estimate of the denoised signalEQUATIONOutput: estimate of the noise level σ; estimate of denoised signal {x t }The discrete convolution of two time series {f t } and {g t } in its general definition is given byEQUATIONwhich in case of the decoder model given in Equation 5translates toEQUATIONdue to the finite non-zero value of τ , i.e., the maximum time difference at which a previous state z t-τ still influences the current state z t according to the hrf . Furthermore, due to causality, only past states are 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 a decoder model without noiseEQUATIONDue to the hrf being the same for all latent dimensions (i.e., a scalar and not a function), we interchange the matrix multiplication with B and the convolution with the hrfEQUATIONEQUATIONConsequentially, we divide the observations into two partsEQUATIONEQUATIONThe second part, Equation 28, is of the same form as Equation 2and permits the same inversion. Therefore, by computing x deconv t from {x t } once, we can perform GTF as in the standard SSM.To include the nuisance artifacts {r t }, one has to also swap the order of the matrix multiplication and the convolution in the full decoder model (Equation 5). This poses a problem since the {r t } are not convolved. Our solution is deconvolving the {r t } time series as wellEQUATION)EQUATIONWe use this as a simple mathematical trick to incorporate the artifacts into the convolutionEQUATIONFinally, we obtain the relation between the deconvolved time series and the latent states asEQUATION)EQUATIONIn order to deal with numerical instabilities, additional hyperparameters were introduced. A too low noise level σ determined by the VISUSHRINK Algorithm 2 can lead to high frequency artifacts, which can be dealt with by defining a lower noise level boundary σmin . 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 x deconv t may occur. With the hyperparameters cut_l, 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 is to be defined relative to the length of the hrf ).If the maximum Lyapunov exponent λ max of a dynamical system is larger than 0, a necessary condition for chaos, nearby trajectories will diverge exponentially. This limits the applicability of n-step ahead 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 a system's initial conditions and yet capture the most relevant dynamical properties. Here, we use two established measures [40, 47] 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 as average over 100 trajectories, generated by perturbing the initial state with a small Gaussian noise term (µ = 0, σ = .01).The n-step prediction error is given byEQUATIONi.e. , the mean squared error between ground truth data {x t } and n-step ahead predictions of the model {x t }.Given an observed N -dimensional time series {x t } of length T and a time series {x t } with the same dimension/length generated by a model, D stsp measures the geometrical overlap of orbits in state space [40] .For low dimensional systems, N ≤ 6, the state space is segregated into k N bins where k is the number of bins per dimension. Each bin is given an index i and we count the number of times n i the time series visited bin i. The relative frequency of visits is then obtained by dividing by time series lengthEQUATION)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-Leibler divergence of these empirical distributions can then be computed to assess the overlap of both systems in their state space geometry,EQUATION)where p i are the relative frequencies of the observed and q i of the predicted (generated) time series. Note that D stsp is not a metric in the mathematical sense, but a divergence that assesses the (dis-)agreement of probability distributions.The complexity of this binning approach scales exponentially with the observation dimension N and thus becomes intractable for larger N . To compute D stsp in higher-dimensional systems, [6] use Gaussian mixture models (GMMs) with centers (means) x t and diagonal covarianceΣ = diag(σ 2 , • • • , σ 2 ), where σ is a hyperparameter. The GMMs along the trajectory points are given byEQUATIONFollowing [30] , the Kullback-Leibler divergence of the two GMMs can be computed using a Monte-Carlo sampling approachD stsp ≈ 1 K K i=1 log f obs (y i ) f gen (y i ) = 1 K K i=1 log 1/T T t=1 N (y i ; x t ; Σ) 1/T T t=1 N (y i ; xt ; Σ) (38)where K is the number of samples drawn, f obs is the distribution of the observed time series, and f gen is the distribution of the generated time series. The binning and GMM-based measures correlate strongly in low dimensions (whereas determining the correlation in high dimensions is challenging for the stated reasons).The state space measure introduced above discards all temporal structure in the data. To include temporal information in the model evaluation, we compare the power spectra of observed and generated time series. For each dimension i ∈ {1, • • • , N }, the scalar time series {x i,t } is converted into a power spectrum density (PSD) {S i,k }. The components are computed from the Fourier transform of {x i,t }EQUATION)The agreement in power spectra is then evaluated in terms of the Hellinger Distance (HD) as suggested by [47] . Since HD is a probability measure, the computed power spectra have to be normalized byEQUATION)Using the power spectra of the observed time series p i,k = Si,k ({x i,t }) and the generated time series q i,k = Si,k ({x i,t }), the HD is then assessed asEQUATIONThe power spectrum error (D P SE ) is computed by averaging the HDs across all N dimensions of the observed system:EQUATIONBefore analysis, the power spectra usually have to be smoothed with a Gaussian kernel of width σ to reduce noise. σ can thus be considered a hyperparameter in the evaluation process and was set to σ = 1 in the current work as in [6] .The Lorenz63 introduced in [45] was designed to describe atmospheric convection based on three dynamic variables. A two-dimensional fluid layer is uniformly warmed from below and cooled from above. It is a continuous-time dynamical system given by the following set of differential equations:EQUATIONEQUATION)EQUATION)where x 1 is proportional to the rate of convection, x 2 to the horizontal temperature variation, and x 3 to the vertical temperature variation. The constants σ, ρ and β are system parameters proportional to the Prandtl number, Rayleigh number, and certain physical dimensions of the fluid-layer. For chaotic behavior, σ = 10, ρ = 28 and β = 8 3 are typical settings. These settings produce the so-called "butterfly attractor" characteristic of the Lorenz system. For each data set, we drew random initial conditions x 0 ∼ N (0, 1 3×3 ) and simulated the system using the DynamicalSystems.jl Julia package, ( [15] ). 10 5 time steps were saved and used as data set (1 : 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 hrf T R function. Gaussian white noise drawn from N (0, σ1 3×3 ) was added to all data points.Results on the Lorenz63 system are in Table 2 . 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-fire neurons (AdEx) with non-linear synaptic currents [3] .We used neurolib to create neural activity generated by an ALN model [12] . In neurolib, the firing rate of the excitatory subpopulation of every brain area is used to simulate the BOLD signal via the Balloon-Windkessel model (for formula see [11] ). As an alternative, we implemented the BOLD decoder model (Equation 5) as model linking the latent states z t (i.e., the latent neural activity corresponding to the excitatory firing rates) to the observed BOLD signal. In order to create interesting dynamics, certain values were altered from the authors' default settings, 'sigma_ou' = 0 and 'b' = 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.neurolib produces simulated (latent) neural activity with a sampling rate of 0.1ms. To stay in a comparable regime with the fMRI and Lorenz time series, we chose a sampling rate of 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 a 30 point finite impulse response (FIR) filter with Hamming window. Furthermore, the neural activity time series was smoothed with a Gaussian kernel with standard deviation σ = 1 and length 5, and then standardized.Outliers in DSR measures (Table 3 ) were removed using the interquartile range (IQR) method. The IQR method considers values as outliers if they are 1.5 IQRs above the third (Q 3 ) or below the first (Q 1 ) quantile (whereIQR = Q 3 -Q 1 ).In Table 3 and corresponding Figure 2D , we present the quantitative comparison between the standard SSM, the convSSM, the convSSM without GTF, and two benchmark conditions. Figure 5 depicts correlations between measures computed on different subsets of the ALN dataset. 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, we assessed the moving average of the variance over time (window size w = 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 https://ftp.gwdg.de/pub/misc/MPI-Leipzig_ Mind-Brain-Body-LEMON/. MINDy: For MINDy [62] we used the implementation at https://github.com/singhmf/MINDy. We trained models for each subject with the settings provided by the authors for fMRI data. To obtain trajectories for calculating P E 10 , D P SE , and D stsp , we used the provided deconvolution function to obtain an 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. LFADS: For LFADS [48] , we used the lfads-torch implementation at https://github.com/ arsedler9/lfads-torch, which provides an 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 as reconstruction target) to a Gaussian, and changed the start and stop learning rates from lr start = 4 • 10 -3 , lr end = 1 • 10 -5 , to lr start = 4 • 10 -4 , lr end = 1 • 10 -6 , which improved the fit to our 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. rSLDS: For rSLDS [44] , we used the implementation at https://github.com/lindermanlab/ ssm. We trained the rSLDS with the Laplace-EM method with the Structured Mean-Field Posterior, as recommended by the authors, with diagonal_Gaussian dynamics and Gaussian_id emissions. For a fair comparison, we used the same number of latent dimensions as for the observations (same as for the cshPLRNN). We determined the rSLDS training parameters α = 0.9 and K = 2 via grid search over α ∈ [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] and K ∈ [ [2, 15] ] by inferring systems using a 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 x and z and the test data as prefix input for the model.To illustrate the difference in "history dependence" for latent time series {z t } and convolved counterparts {x t }, an illustrative example was created: A latent time series {z t } was produced stochastically Justification: We discuss the limitations of our approach at the end of section 4. Guidelines:• The answer NA means that the paper has no limitation 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 in 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 on 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 a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but URLs data and code is permitted.Justification: The dataset we use containing human subjects was fully anonimized and all subjects gave consent for the the data to 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 the University of Leipzig (reference number 154/13-ff), see [4] for details.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: * From my understanding, there is actually no new model but different ways of training a particular model, although they are called ConvSSM, LinSSM, etc.
* 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.
* 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.
questions: questions: /
limitations: limitations: /
paper_summary: paper_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.
strength_weakness: strength_weakness: The work has several weaknesses and limitations detailed below. I’m open to discussion about any of them.

- 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.
- 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.
- 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.
- 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.).
- Some of the results go against the main motivation and claims of the paper. Detailed comments are included below.

**Unknown observation filter**

- 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.
- _“hrf with alternative functions if we want to account for filtering in the original signal”_

Related 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)?

- 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_.

**Results against the claims/motivation**

- _“\[40\] proved that for chaotic systems gradient-based training techniques for RNNs will inevitably lead to diverging loss gradients.”_

This 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).

- 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?

**Fig. 2**

- Fig. 2E,F are labeled incorrectly?
- Unclear what’s shown in Fig. 2A.
- What are the takeaways from this figure?
- In Fig. 2C, can you include $\\lambda_{max}$ distributions from other models (specifically convNoGTF model)?

**Long-term performance measures**

- 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?
- 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?
- 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?
- 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).
- 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?
- 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?

**More comments**

_“This is not automatically given for standard RNN”_

No 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.

_“guiding the training process through optimally chosen control signals – also referred to as teacher forcing (TF) signals”_

At 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.

_“Chaotic systems in particular, as typically encountered in neuroscience (e.g., \[53, 21, 32\]), pose a severe problem here”_

Another paper considered the theoretical aspects of training RNNs on chaotic data \[1\], please cite and discuss it.

_“but not if it depends on a history of states”_

The assumption of history dependence is a very common assumption made in most latent variable models.

\[1\] Engelken, Rainer, Fred Wolf, and Larry F. Abbott. "Lyapunov spectra of chaotic recurrent neural networks." Physical Review Research 5.4 (2023): 043044.

\[2\] Brenner, Manuel, Georgia Koppe, and Daniel Durstewitz. "Multimodal teacher forcing for reconstructing nonlinear dynamical systems." arXiv preprint arXiv:2212.07892 (2022).

\[3\] Schimel, Marine, et al. "iLQR-VAE: control-based learning of input-driven dynamics with applications to neural data." bioRxiv (2021): 2021-10.
questions: questions: The effect of nuisance variables is considered to be linear, why is this a good assumption?

The 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)?

Comparisons 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: limitations: See the weaknesses section.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- 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.
questions: 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?
- How does the standard SSM technique perform on the LEMON dataset? Does the convSSM method provide better insights?
limitations: limitations: The authors have adequately addressed the limitations of their work
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

However, 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
questions: 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: limitations: Using AR(1) for system equation and simultaneously making current ovservation explicitly depending on the whole latent history may make inference difficult.

Comments: title: Response to reply
Reply
Response to updated review
Response to updated review
Updated review
Response to comment
comment: Great, thanks for running the analysis very quickly. Indeed the conditional independence test would be more appropriate and interesting, but the current result is already consistent with what you'd like to show. The large mutual information could indeed be driven by the consecutive time points and Markov property. I wonder if this could also be used as a strategy to pick the best hyperparameters or learn the filter in case the optimal filter needs some tuning (i.e. find the filter or find the hyperparameters that result in the highest conditional independence of $p(z_t|z_{t-1}),p(z_t|z_{t-1:t-\tau})$). Anyway, I think my comments are mostly addressed and I look forward to the final version of the paper.
As the reviewer suggested, we computed the MI between the deconvolved states $z_t$ and $z_{t-1:t-\tau}$ and between observables $x_t$ and $x_{t-1:t-\tau}$ on ten time series of the Lorenz63, using the suggested package. We indeed find that the MI for the latent process is significantly lower (with a mean of $20.59 \pm 0.01$ for the MI evaluated on the states, and $26.10 \pm 0.92$ on the observables, $p<.001$). However, given that we are yet unfamiliar with this package, we do not wish to overstate these results. Also, based on the arguments we had provided in our previous response, we still believe the more appropriate test is the one on conditional independence, on which we will continue to work.



We will add the information on RC2 and Q2 to the Appendix, as requested. 



We thank you again for your detailed engagement with our paper!
The new discussion on extensions looks great, thanks for sharing it here. Regarding UF1, [1] can help pick the right estimator (they also have an accompanying code package). I think it'd be a very compelling result to show that the filter is indeed helping with the Markovian dynamics, which would strengthen the paper's motivation. It'd be great if you could rephrase what you have above in response to RC2 and Q2 and add it to the supplementary, this information could help practitioners find the right hyperparameters. 

I look forward to the updated manuscript with this information and comparisons against latent variable models.

[1] Czyż, Paweł, et al. "Beyond normal: On the evaluation of mutual information estimators." Advances in Neural Information Processing Systems 36 (2024).
Thank you.
RC2) Sure. Most hyperparameters were adopted from previous experiments (see Appx., Table 4). For the LEMON data, we additionally determined the latent dimensions (and $\alpha$) via grid search, by inferring systems using a subset of the data and assessing the performance measures ($D_{stsp}$, $D_{PSE}$) on the held-out set (see also [2]). The latent dim. of the ALN dataset was then set to match that of LEMON. For the Lorenz63, we selected a hidden dim. in which the no-filter model performed well [5], to have a strong baseline comparison. All other hidden dims. were aligned with this setting. We will provide additional details to clarify.

PM3) Yes exactly. The PE does not assess the generative long-term statistics and, especially when the system is chaotic, will not be informative about how well the underlying attractor has been reconstructed. 
The purpose of Fig. 2E, however, was to demonstrate that assessing $D_{stsp}$ on short sequences, as common in empirical studies, still provides a good indicator for the actual long-term performance, i.e. as assessed from long sequences. This is essential for using $D_{stsp}$ to evaluate reconstruction performance on relatively short empirical data.

Q2) We apologize for misunderstanding your question. Piecewise linear models have demonstrated good reconstruction performance with far less than 50 pieces [12]. However, increasing the number of pieces enhances the model's flexibility and improves performance. How many pieces are required also depends on the complexity of the underlying attractor, and whether one is interested in capturing the full geometry (as assessed via $D_{stsp}$) or just general topological features. Since our primary goal was to compare the filter and non-filter options, ensuring that performance differences were due to the model's efficacy in well-performing regimes [5], rather than being limited by too few dimensions, we opted for higher dim. here. Further, SLDS models differ from our approach in that they allow different parameter matrices for each regime ("piece"), while here we ensure continuity in the states and the switching in Jacobians is solely regulated through the ReLU function, i.e. matrices {A, W1, W2} are constant across state space. 


Q3) We promise to include such comparisons in our revised manuscript. We already started runs with SLDS, but since there is only 1 day left, we are skeptical we will be able to finish this by tomorrow. We note that previous experience with LFADS did not yield comparable performance, as LFADS often tended to converge to fixed points in the long run (i.e., did not capture the underlying attractor and long-term statistics).
W1-3,Q1) This is a very good idea, thank you for the great suggestions! We add here two paragraphs that we will integrate into the Discussion:

“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. The GTF training framework would remain unchanged as the control signal and the latent state update (eqn. (3)) are not affected by such modifications [5]. 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 on convolution filters like the HRF, alternatively we can parameterize the filter weights within the observation model, making them learnable through BPTT, with filter length either as a 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 [12]. 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 an initial stage would have the advantage that subsequent training would enjoy the same speed benefits as in our suggested method.

Finally, we would like to highlight that our framework could be adapted to accommodate noise in the latent process. For example, in Brenner et al. [4] the GTF procedure has been modified to work in the context of stochastic DSR models using variational inference. The key idea lies in introducing a (multimodal) variational autoencoder that takes the observed variables as input and maps them into the DSR model’s latent space, thereby generating the control signal required for GTF. In a similar fashion, we could replace the MVAE with the inversion in eqn. (10), thereby providing a TF signal to be used to steer a probabilistic latent DS model, i.e. its distributional mean, via eqn. (3), and use the reparameterization trick [17,18] for BP in the latent space. However, although probabilistic frameworks are appealing, ‘deterministic’ BPTT has previously been shown to be (at least) comparable in terms of DSR performance, even for noisy observations and/or processes [2], such that the benefits for DSR would need to be further examined.”


17) Rezende, D. J., Mohamed, S., and Wierstra, D. Stochastic Backpropagation and Approximate Inference in Deep Generative Models. In Proceedings of the 31st International Conference on Machine Learning, 2014. 

18) Kingma, D. P. and Welling, M. Auto-Encoding Variational Bayes. In Proceedings of the 2nd International Conference on Learning Representations, 2014.


UF1) In a deterministic DS observed within the full state space, the complete future of a trajectory is fully determined by its initial condition, so we would not expect the MI to decay to 0 with $\tau \rightarrow \infty$ (and in fact, we checked this and it does not), but we would expect that the previous state contains all information about the current state if the deconvolution worked as intended. So - in our minds - the crucial question would be whether we have $p(z_t|z_{t-1:t-\tau})=p(z_t|z_{t-1})$ in latent space, but $p(x_t|x_{t-1:t-\tau}) \neq p(x_t|x_{t-1})$ in observation space (i.e., confirming the Markov property in latent but not observation space). These high-dimensional multiple state probabilities are, however, hard to access. For the revision, we will look into different ways we can approximate them, to illustrate this property.
Sorry for the late response. I think the response from authors resolved my questions. And I have no further questions. In this case, I have raised my score from 4 to 5.
Thank you for replying to all questions and comments. In particular, the addition of new experiments and results are very helpful for a more in-depth understanding of the paper. I will slightly increase my score and would appreciate it if the authors include the following.

**W1-3,Q1)** Can you include these two in the discussion? Specifically, can you describe in the discussion part of the paper how to extend the method to learnable filters (W1), how to extend to variational schemes for models with noise in the latent space (W2), how to extend GTF to arbitrary choices of the architecture in addition to PLRNN (W3), and how to extend to nonlinear nuisance model (Q1)? This could open up new applications of the method and make the paper more accessible to a wider audience. Would you include the updated discussion here for a review?

**UF1)** A simple check is to look at the mutual information between x_t, x_{t-1:t-\tau} and compare it to the mutual information between z_t, z_{t-1:t-\tau}. If these two are largely different, it shows that deconvolution has successfully removed the history dependence.

**UF3, RC1)** The new experiments address these two.

**RC2)** Can you add some details on model selection to the discussion? Specifically, how did you set the latent dimension for specific datasets? This would be important for practitioners who’d be interested in applying this model to other datasets.

**PM3)** Apologies for the unreferenced question. This question is mostly addressed by Fig. R3,4,6 but just to double-check; from what I understand after reading the rebuttal the authors are arguing against the use of PE (either 1-step or 10-step) as it doesn’t represent “long-term” behavior of the dynamical systems and instead suggest using D_{stsp} or D_{PSE}. Is the point of Fig. 2E to say that increasing the horizon of prediction doesn’t help with better quantification of DSR?

**Q2)** Previous switching linear approaches have shown successful reconstruction using 2 pieces. I guess my question is why do we need the number of states to be this large for a good reconstruction?

**Q3)** At least a comparison against one of the switching models (e.g. SLDS or rSLDS) and a nonlinear model (e.g. LFADS) would be very helpful just to make your point.
Thank you. Much appreciated.
I thank the authors for their detailed rebuttal and for providing additional material in such a short time frame. I fell most of my concerns are now resolved and I have raised my score (1 point) to reflect the addition of further analyses.
Dear reviewer, 

thank you once more for your time and effort in reviewing our work. We would kindly like to ask if we were able to address your concerns, and whether there are any remaining questions we can clarify.
Dear reviewer, 

thank you once more for your time and effort in reviewing our work. We would kindly like to ask if our rebuttal adequately addressed your questions and concerns, and whether there are any remaining questions we can clarify.
Dear reviewer, 

thank you again for your time and effort in reviewing our manuscript. We would kindly like to ask if our rebuttal adequately addressed your major concerns, or if there are any remaining questions we can clarify. 

Due to the strict character limit, our responses were partly brief. We are happy to elaborate on points that may yet be unclear.
We thank the reviewer for engaging with our rebuttal, and are happy to hear we could address the open points!
Thanks for the authors for their detailed explanations. Most my concerns are resolved, and hence I raised my rating with 1 more score.