PaperlayerView exact revision on arXiv

Source-linked AI summary

Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI

Georgia Koppe, Hazem Toutounji, Peter Kirsch, Stefanie Lis, Daniel Durstewitz

arXiv:1902.07186v2cs.LGq-bio.NCstat.ML

TL;DR

Fitting observations alone does not guarantee that a latent model captures the underlying system dynamics. The paper introduces a generative PLRNN state-space model with stepwise training and a state-space KL measure; on the van der Pol system, its trajectories follow the limit cycle while LDS trajectories decay or diverge.

Abstract

from arXiv · show

A major tenet in theoretical neuroscience is that cognitive and behavioral processes are ultimately implemented in terms of the neural system dynamics. Accordingly, a major aim for the analysis of neurophysiological measurements should lie in the identification of the computational dynamics underlying task processing. Here we advance a state space model (SSM) based on generative piecewise-linear recurrent neural networks (PLRNN) to assess dynamics from neuroimaging data. In contrast to many other nonlinear time series models which have been proposed for reconstructing latent dynamics, our model is easily interpretable in neural terms, amenable to systematic dynamical systems analysis of the resulting set of equations, and can straightforwardly be transformed into an equivalent continuous-time dynamical system. The major contributions of this paper are the introduction of a new observation model suitable for functional magnetic resonance imaging (fMRI) coupled to the latent PLRNN, an efficient stepwise training procedure that forces the latent model to capture the 'true' underlying dynamics rather than just fitting (or predicting) the observations, and of an empirical measure based on the Kullback-Leibler divergence to evaluate from empirical time series how well this goal of approximating the underlying dynamics has been achieved. We validate and illustrate the power of our approach on simulated 'ground-truth' dynamical (benchmark) systems as well as on actual experimental fMRI time series, and demonstrate that the latent dynamics harbors task-related nonlinear structure that a linear dynamical model fails to capture. Given that fMRI is one of the most common techniques for measuring brain activity non-invasively in human subjects, this approach may provide a novel step toward analyzing aberrant (nonlinear) dynamics for clinical assessment or neuroscientific research.

Abstract

The paper develops a generative PLRNN state space model to infer governing dynamics from neural measurements, addressing limits of linear models that cannot reproduce many nonlinear phenomena. Its stepwise training, state-space KL measure, and fMRI observation model target recovery of latent dynamics rather than observation fit alone.

Introduction

  • Linear state space models can approximate latent trajectories but cannot adequately model dynamics such as multistability, complex limit cycles, or chaos, motivating nonlinear dynamical models.RNNs can approximate the flow of other dynamical systems and yield neuronally interpretable equations.Source: Introduction · p. 3

Author Summary

  • The approach combines a stepwise training procedure, a state-space KL measure, and an fMRI observation model to push latent dynamics toward the underlying system rather than merely fitting observations.The measure compares the prior generative dynamics with inferred latent states; the fMRI decoder accounts for hemodynamic filtering.Source: Introduction · pp. 3–4

PLRNN-based state space model (PLRNN-SSM)

  • A PLRNN models latent dynamics using diagonal autoregression weights, off-diagonal recurrent connections, and an element-wise ReLU transfer function.The model also includes external inputs and Gaussian process noise; its parameters can be interpreted in neural-rate terms.Source: PLRNN-based state space model… · p. 4

Observation model for BOLD time series

The BOLD observation model links latent PLRNN states to fMRI measurements through HRF convolution and nuisance predictors. The authors also use stepwise training and state-space KL comparisons to assess recovery of underlying dynamics.

Observation model for BOLD time series

  • The BOLD observation model generates measurements from HRF-convolved latent states using a linear-Gaussian model, with nuisance predictors accounting for recording artifacts such as movement.Stimuli act through latent neuronal dynamics, whereas nuisance variables directly affect the observations.Source: Observation model for BOLD time series · p. 5

Stepwise initialization and training protocol

  • The stepwise protocol shifts the fitting burden from the observation model to the latent dynamics, eventually fixing observation parameters and increasing the latent model's temporal-consistency contribution.Ordinary EM can favor observation fit over a good representation of the underlying dynamics.Source: Stepwise initialization and training… · p. 6

Evaluation of training protocol

  • Annealing significantly lowered normalized KL divergence and raised expected joint log-likelihood versus random initialization, and produced more estimates that reconstructed the chaotic attractor.Across training steps, KL_x decreased while expected joint log-likelihood increased; the two measures were only moderately correlated (r=-27, p<.001).Source: Evaluation of training protocol · pp. 8–9

Reconstruction of benchmark dynamical systems

On experimental fMRI, PLRNN models better captured generative dynamics than linear state-space models and represented task-related structure. Benchmark tests showed that the method could recover nonlinear attractors, while the KL measure offered a useful—but imperfect—indicator of reconstruction quality.

Reconstruction of benchmark dynamical systems

Reconstruction of experimental data

  • The latent-space KL proxy correlated with observation-space KL on benchmarks and lower normalized latent KL accompanied better agreement between inferred and generated state trajectories.The authors caution that time-series agreement alone does not necessarily establish that the underlying dynamics were captured.Source: Reconstruction of experimental data · pp. 13–14
  • After about three prediction steps, PLRNN-BOLD-SSM retained better observation-prediction accuracy than LDS-BOLD-SSM, whose errors and latent-state discrepancies rose sharply.LDS-generated trajectories often converged to fixed points, whereas trained PLRNNs often reproduced oscillatory activity in the data.Source: Reconstruction of experimental data · pp. 14–15
  • For M>2, PLRNN prior-generated trajectories significantly reduced task-classification error, including without stimulus inputs, while LDS-generated trajectories were near chance without inputs.Models with inputs enabled almost perfect classification with as few as M=8 states; LDS performance was also worse on inferred trajectories.Source: Reconstruction of experimental data · pp. 16–17, Figure 7 · p. 16
  • Trained PLRNNs produced stable limit cycles, chaotic attractors, and multistability in fMRI fits; their oscillations also reflected task-phase and task-sequence periodicities.In the CMT, better performance was associated with fewer stable and more unstable objects (F(1,24)=5.277, p=.031); the authors call for further study of this link.Source: Reconstruction of experimental data · pp. 18–19, Figure 10 · p. 19

Discussion

Neuroscience theories often describe computation through nonlinear dynamical systems, motivating tools that recover governing equations as well as state-space representations. The paper proposes a generative PLRNN state-space model and evaluates its ability to capture dynamics in empirical fMRI data.

Discussion

  • Neuroscience theories describe computation through attractor states, transitions, and other properties of nonlinear dynamical systems.Source: Discussion · p. 19
  • The proposed PLRNN state-space model is trained to reproduce underlying dynamics, with a state-space KL divergence assessing reconstruction when ground truth is unavailable.Source: Discussion · p. 20
  • In a working-memory paradigm, the approach recovered nonlinear dynamics and trajectory flows from human fMRI recordings related to task structure and behavioral performance.Source: Discussion · p. 20

Choice of model formalism and nonlinearity

The discussion presents PLRNNs as expressive, analyzable models whose generative dynamics can capture nonlinear structure, and describes an HRF-aware observation model for fMRI. It contrasts this framework with alternative methods and notes broader applications to filtered recordings.

Choice of model formalism and nonlinearity

  • The PLRNN is chosen for its dynamical-system expressiveness, while ReLU supports efficient inference and approximates input-output functions of some cortical cell types.Source: Choice of model formalism and… · p. 20
  • Task-specific coupling matrices are possible, and latent dimensionality may need to exceed the number of observed channels rather than reduce it.Source: Choice of model formalism and… · p. 21

Comparison to other approaches for identifying dynamical systems

  • Delay embedding does not yield governing equations, while linear latent dynamics generally cannot reconstruct nonlinear system dynamics; earlier RNN approaches emphasized inferred states over systematic generative analysis.Source: Comparison to other approaches for… · p. 21
  • PLRNN equations have a neural-rate interpretation, admit continuous-time translation, and allow analytical fixed-point and stability analysis.Source: Comparison to other approaches for… · p. 22
  • Eight latent states and one 1,000-sample time series were sufficient for the PLRNN-SSM to rebuild the chaotic Lorenz attractor.Some other methods required hundreds of hidden states and millions of samples for the van der Pol or Lorenz attractors.Source: Comparison to other approaches for… · p. 22

Applications in fMRI research and beyond

  • The fMRI observation model accounts for HRF filtering by allowing current observations to depend on a sequence of previous states.Source: Applications in fMRI research and beyond · p. 22
  • fMRI remains valuable for clinical research despite coarser temporal resolution than MEG or EEG, and the framework could help test dynamical hypotheses about psychiatric conditions.The cited proposal about altered attractor basins in schizophrenia still lacks direct experimental evaluation.Source: Applications in fMRI research and beyond · pp. 22–23
  • The framework can accommodate other modalities with low- or band-pass filtering by replacing the observation model.The paper names calcium imaging, voltage-sensitive dyes, EEG, MEG, single-unit data, mobile sensors, ecological momentary assessments, and electronic health records.Source: Applications in fMRI research and beyond · p. 23

Open issues and outlook

The authors identify unresolved questions in both training and evaluation: the stepwise procedure is heuristic, and empirical data lack ground truth for assessing recovered dynamics. They motivate evaluating freely generated trajectories and state-space distributions, while noting practical limits for fMRI.

Open issues and outlook

  • The stepwise training heuristic appears to improve recovery over a bootstrap sample of standard-trained models, but needs systematic derivation and testing of alternatives and explicit optimization criteria.The authors note that the criteria could directly use equation 10.Source: Open issues and outlook · p. 23
  • For empirical data without ground truth, how to quantify confidence in reconstructed equations of motion remains open; time-series prediction errors are not clearly suitable for complex or chaotic dynamics.The authors note that standard squared-error or likelihood measures can be unhelpful for nonlinear systems with complex or chaotic dynamics.Source: Open issues and outlook · p. 23
  • Because small initial-condition or noise differences can quickly decorrelate trajectories, the authors compare state-space distributions and propose evaluating newly generated trajectories rather than short-horizon predictions.Their rationale is that dynamical equivalence concerns similar long-term behavior, not literal time-series overlap.Source: Open issues and outlook · p. 24
  • fMRI’s high-dimensional, short data series make the observation-space distribution difficult to estimate, so the original comparison measure is not directly applicable.The authors instead compare observation-constrained and freely generated latent-state distributions, approximating the former with a Gaussian mixture to enable an analytical KL approximation.Source: Open issues and outlook · p. 24

Model specification and inference

The model is inferred with EM, using a global Gaussian approximation to the latent-state posterior and piecewise optimization for the E-step. The M-step updates model parameters from posterior expectations, with a global procedure suited to the temporal dependencies introduced by the fMRI hemodynamic response function.

Model specification and inference

  • EM maximizes an evidence lower bound on the observation log-likelihood, and the bound becomes exact when the proposal latent-state density matches the true posterior.The E-step aims to determine or approximate that posterior; the M-step updates parameters with the proposal held fixed.Source: Model specification and inference · pp. 24–25
  • Because explicit integration over the high-dimensional mixture of piecewise Gaussians is infeasible at large T and M, the E-step uses a MAP estimate and a local Gaussian approximation.The state covariance is approximated from the inverse negative Hessian around the MAP estimate.Source: Model specification and inference · pp. 25–26
  • The E-step optimization is piecewise quadratic in the latent states and activation-pattern indices, with constraints enforcing consistency between each state’s sign and its index.A Newton-type iteration alternates between solving for states with fixed activation indices and flipping inconsistent indices until convergence; a quadratic-programming step may be added.Source: Model specification and inference · pp. 25–26
  • The global E-step supports fMRI’s HRF-induced dependencies across multiple past latent states, which prevent the posterior from factorizing as required by common unscented or extended Kalman filters.The authors state that efficient block-tridiagonal inversion still scales linearly in T.Source: Model specification and inference · p. 26
  • In the M-step, fixing the posterior approximation reduces parameter updates to linear regression because the ReLU nonlinearities have already been resolved in the required expectations.The observation-model and latent-model parameters can be solved separately; latent parameters are efficiently estimated row-wise given the off-diagonal structure of W.Source: Model specification and inference · p. 26
  • For the observation-model update, the derivation concatenates the observation and nuisance matrices and uses sums of cross-products involving HRF-convolved states.The E- and M-steps are iterated until the expected joint log-likelihood converges.Source: Model specification and inference · pp. 26–27

Stepwise model training procedure

The stepwise training scheme shifts fitting responsibility from the observation model toward the latent dynamics, aiming to make the latent model capture the data-generating system. It starts with a linear latent initialization and progressively tightens temporal consistency.

Stepwise model training procedure

  • Training begins with a linear latent model and mild latent-process constraints, then gradually decreases the ratio of Σ to Γ to enforce temporal consistency.The scheme may later fix the observation-model parameters, placing the full burden of observation fitting on the latent model.Source: Stepwise model training procedure · p. 27

Algorithm-1

Algorithm-1 initializes with a linear dynamical system under Σ = I and includes an eigenvalue constraint that biases models toward stability.

Algorithm-1

  • A maximum absolute eigenvalue below 1 for A + W biases the model toward stable dynamics.Source: Algorithm-1 · p. 28
  • The procedure initializes by fixing Σ = I and running a linear dynamical-system state-space model.Source: Algorithm-1 · p. 28

Reconstruction of benchmark dynamical systems

The benchmark evaluation tests whether trained PLRNNs reproduce the underlying Lorenz and van der Pol dynamics when run freely, rather than merely fitting observed trajectories. It compares annealed training with random-initialized EM and evaluates spatial distribution, stability, and dynamical properties.

Reconstruction of benchmark dynamical systems

  • The evaluation used Lorenz and van der Pol systems, with 100 ground-truth trajectories of length 1000 and PLRNN models spanning 5–20 latent states.Annealed stepwise training was compared with random-initialized EM for models with 8, 10, 12, and 14 latent states.Source: Reconstruction of benchmark dynamical… · p. 28
  • The free-running PLRNN was evaluated by comparing its state-space distribution with the ground-truth distribution, assessing whether the reconstructed attractor geometry was captured.The divergence is computed across state space rather than time, avoiding trajectory mismatch expected even between runs of the same chaotic system.Source: Reconstruction of benchmark dynamical… · pp. 28–29
  • Unstable estimates were excluded from the normalized KL and likelihood comparisons; strong likelihood outliers occurred only in about 3.8% of random-initialized fits and were removed.About 14% of unstable estimates were removed.Source: Reconstruction of benchmark dynamical… · p. 29
  • For Lorenz, reconstruction was also assessed by deviation in maximal Lyapunov exponents; for the non-chaotic van der Pol oscillator, power-spectrum agreement was assessed by average Pearson correlation.The Lyapunov exponent was estimated from the initial growth of distance between nearby trajectories.Source: Reconstruction of benchmark dynamical… · p. 29

Reconstruction of dynamical systems from experimental data

The experimental evaluation applies the PLRNN-BOLD-SSM to fMRI time series from working-memory and control tasks. Because the true dynamics are unknown, it uses agreement between data-constrained latent states and the model’s unconstrained generative prior as a reconstruction proxy.

Reconstruction of dynamical systems from experimental data

  • The fMRI dataset comprised 26 healthy participants, with 20 regional time series per participant selected from regions relevant to the n-back task.The paradigm included two working-memory tasks and a 0-back choice-reaction control task.Source: Reconstruction of dynamical systems… · p. 30
  • The PLRNN-BOLD-SSM and LDS-BOLD-SSM used the regional observations and six movement nuisance predictors; a direct PLRNN fit tied each latent state one-to-one to an observation.The direct-observation model set M=N and constrained B to be diagonal.Source: Reconstruction of dynamical systems… · p. 30
  • Because the ground-truth dynamics are unavailable for fMRI, reconstruction is assessed by the state-space KL divergence between the inferred posterior and the data-unconstrained latent prior.The proxy assumes these distributions should match closely when the observations represent the underlying dynamics well.Source: Reconstruction of dynamical systems… · p. 31
  • To approximate the latent-state KL, the inferred posterior and generative prior were represented as Gaussian mixtures across trajectory times, with Monte Carlo or variational approximation used for evaluation.The variational approximation was reported to agree almost exactly with Monte Carlo results and was suggested for high-dimensional latent spaces.Source: Reconstruction of dynamical systems… · pp. 31–32

S1 Text. Model specification and inference.

The supplementary text details PLRNN-BOLD-SSM inference and parameter estimation, including the HRF observation convolution and closed-form latent-parameter updates. It also reports latent-state settings for benchmark systems and checks the KL measure’s sensitivity to bin size.

S1 Text. Model specification and inference

  • The observation model convolves latent states with the hemodynamic response function, whose vector length depends on the time-series temporal resolution.The HRF convolution is expressed as a weighted sum of latent states over the HRF window.Source: S1 Text. Model specification and… · p. 38
  • The E-step estimates latent states by maximizing the log-joint distribution, treating the posterior as characterized by its first two moments and its MAP estimate as an approximation to the mean.The active states are defined by z_t>0 in the binary activation vector.Source: S1 Text. Model specification and… · p. 38
  • The M-step solves for latent transition, recurrent, bias, and input parameters together using state, activation, constant, and input predictors.Because the parameter matrices are constrained, the solution is computed one row at a time using only the defined parameters.Source: S1 Text. Model specification and… · pp. 39–40
  • M=14 was about optimal for van der Pol and M≈16 for Lorenz, based on dependence of normalized KL on latent-state count.These settings are reported in S1 Fig.Source: S1 Text. Model specification and… · p. 40
  • The likelihood landscape can form piecewise-Gaussian regions that cluster near zero states, producing an elevated plateau of high-likelihood solutions around the MAP estimate.The illustrated landscape varies a single latent state across two consecutive time steps.Source: S1 Text. Model specification and… · p. 41
  • KL estimates across bin widths were nearly monotonically related, preserving rank information; Δx=1 gave the largest spread and was selected for the main evaluation.The comparisons used Lorenz-system bins of Δx=1, 0.5, 0.2, and 0.1.Source: S1 Text. Model specification and… · p. 41

Related papers

Loading 1902.07186v2…