A scalable generative model for dynamical system reconstruction from neuroimaging data

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.

Paper

Similar papers

Peer review

Reviewer qvr37/10 · confidence 3/52024-06-14

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.

Strengths

1. The proposed method (convSSM trained via SGD+GTF) can be efficiently scales with model size and convolution length 2. The model selection strategy based on short empirical time series are practically useful. 3. 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. 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

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?

Rating

7

Confidence

3

Soundness

4

Presentation

3

Contribution

3

Limitations

Using AR(1) for system equation and simultaneously making current ovservation explicitly depending on the whole latent history may make inference difficult.

Reviewer 6CX87/10 · confidence 3/52024-07-09

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.

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

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

Rating

7

Confidence

3

Soundness

4

Presentation

4

Contribution

3

Limitations

The authors have adequately addressed the limitations of their work

Reviewer 9P9G5/10 · confidence 3/52024-07-10

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.

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. The writing quality is good, which makes it easy to follow the arguments and contributions. Most 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. - 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

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?

Rating

5

Confidence

3

Soundness

3

Presentation

3

Contribution

2

Limitations

See the weaknesses section.

Reviewer oqqk5/10 · confidence 3/52024-07-11

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.

Strengths

* The two techniques seems effective for training an SSM model and dynamical system reconstruction. * Experimental results on both synthetic and real-world datasets show that ConvSSM with GTF is better than others. * 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. * 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

/

Rating

5

Confidence

3

Soundness

2

Presentation

2

Contribution

2

Limitations

/

Reviewer qvr32024-08-09

Thanks for the authors for their detailed explanations. Most my concerns are resolved, and hence I raised my rating with 1 more score.

Authorsrebuttal2024-08-09

Response to comment

We thank the reviewer for engaging with our rebuttal, and are happy to hear we could address the open points!

Authorsrebuttal2024-08-11

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.

Reviewer 6CX82024-08-12

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.

Reviewer 9P9G2024-08-12

Updated review

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.

Authorsrebuttal2024-08-13

Response to updated review

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.

Authorsrebuttal2024-08-13

Response to updated review

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

Reviewer 9P9G2024-08-13

Reply

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

Authorsrebuttal2024-08-13

Response to reply

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!

Reviewer 9P9G2024-08-13

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.

Reviewer oqqk2024-08-12

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.

Program Chairsdecision2024-09-25

Decision

Accept (poster)

© 2026 NYSGPT2525 LLC