A scalable generative model for dynamical system reconstruction from neuroimaging data

TL;DR

This work proposes a novel algorithm that is efficiency in reconstructing dynamical systems, including their state space geometry and long-term temporal properties, from just short BOLD time series and scales exceptionally well with model dimensionality and filter length.

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.

Figures (7)

• Figure 1: Schematic of training protocol and gradient flow. A: Before training, observations $\Bqty{x_t}$ and nuisance artifacts $\Bqty{r_t}$ are deconvolved. B: The deconvolved time series are used to generate a forcing signal $d_{t-1}$ which is used for guiding cshPLRNN training. C: Latent states $z_{t-\tau:t}$ and nuisance artifacts $r_t$ are used to predict $\hat{x}_t$ through the decoder model. Gradients are computed on the squared error loss $\mathcal{L}_t$, propagated from the decoder model back to the latent states (blue), and from the latent DS model backwards in time (orange).
• Figure 2: Validations on Lorenz63 and ALN. A: Illustration of reconstruction performance as assessed by the geometrical agreement measure $D_{stsp}$. Average $D_{stsp}$ values for the convSSM were $D_{stsp}< 0.30$ at noise level $\sigma = .01$ and $D_{stsp} < 0.71$ at noise level $\sigma = .1$, indicating successful reconstructions in the majority of cases. B: Example trajectory from the Lorenz63 system in latent space (top) and observation space (convolved with $hrf_{0.2}$) (bottom). C: Probability density over maximal $\lambda_{max}$ values (orange) assessed on 1000 convSSMs trained on Lorenz63 time series of length $1000$ (example shown in right panel). Black line denotes the known $\lambda_{max}\approx 0.9056$ of the Lorenz system. D: Comparison of standard SSM ('standard'), convSSM ('conv'), and convSSM trained without generalized teacher forcing ('conv (NoGTF)') on the ALN data set. Histograms over $D_{stsp}$ assessed on the observed space (left panel) and latent space (right panel). E: $D_{stsp}$ for convSSM evaluated on the full pseudo-empirical time series of typical empirically available length ($T=500$; x-axis) vs. the long GT test set ($T=5,000$; y-axis). F: $D_{stsp}$ for convSSM evaluated on the observed time series (x-axis) vs. on the latent time series (y-axis).
• Figure 3: Results on empirical LEMON data set. A: Distribution over $D_{stsp}$ for $1020$ systems inferred with convSSM. B: Example of a good and C: poor reconstruction. D: Illustration of reconstruction performance as a function of $D_{stsp}$. E: Histogram over maximum Lyapunov exponents $\lambda_{max}$. F: Distribution over $\lambda_{max}$ for 5 selected participants ($n=100$ systems with $10$ trajectories each). G: Within- as compared to between-subject variance in $\lambda_{max}$ distribution after filtering models by DS performance measures (selecting the 20 best by $D_{stsp}$ and 10 best by $D_{PSE}$).
• Figure 4: Training duration per epoch (y-axis) in seconds for different TRs (A), hidden dimensions L (B), latent dimensions M (C), time series length (D), and observation dimensions (E). Mean, standard error (SEM) and linear curve fits (gray dashed lines) are displayed. The per-epoch-runtime increases approximately linearly with dimensions $L$, $M$, and $N$; explained variance $R_L^2=0.989$, $R_M^2=0.993$, and $R_N^2=0.996$ for linear regressions with predictors '$L$', '$M$', and '$N$', respectively. Experiments were performed on a standard notebook with Intel i5-8250U 1,60 G CPU and 8GB RAM.
• Figure 5: A: Agreement in DSR measures assessed on the observed (x-axis) vs. latent (y-axis) space of the short pseudo-empirical test set (top) and the full pseudo-empirical time series (bottom). Correlations between $D_{stsp}$ (left), $D_{PSE}$ (middle), and $PE_{10}$ (right) are displayed, respectively. B: Top: Agreement in DSR measures assessed on the pseudo-empirical test set (short) vs. GT test set (long). Bottom: Same for full pseudo-empirical time series (short) vs. GT test set (long). Correlations between $D_{stsp}$ (left), $D_{PSE}$ (middle), and $PE_{10}$ (right) are displayed, respectively. C: Correlations between DSR measures between pseudo-empirical test set and full pseudo-empirical time series, same order as in B.