State-space solutions to the dynamic magnetoencephalography inverse problem using high performance computing

Determining the magnitude and location of neural sources within the brain that are responsible for generating magnetoencephalography (MEG) signals measured on the surface of the head is a challenging problem in functional neuroimaging. The number of potential sources within the brain exceeds by an order of magnitude the number of recording sites. As a consequence, the estimates for the magnitude and location of the neural sources will be ill-conditioned because of the underdetermined nature of the problem. One well-known technique designed to address this imbalance is the minimum norm estimator (MNE). This approach imposes an $L^2$ regularization constraint that serves to stabilize and condition the source parameter estimates. However, these classes of regularizer are static in time and do not consider the temporal constraints inherent to the biophysics of the MEG experiment. In this paper we propose a dynamic state-space model that accounts for both spatial and temporal correlations within and across candidate intracortical sources. In our model, the observation model is derived from the steady-state solution to Maxwell’s equations while the latent model representing neural dynamics is given by a random walk process.

💡 Research Summary

The paper addresses the notoriously ill‑posed inverse problem of magnetoencephalography (MEG): estimating the three‑dimensional locations and amplitudes of neuronal current sources from magnetic field measurements recorded at a few hundred sensors on the scalp. Classical approaches such as the Minimum‑Norm Estimator (MNE) impose an L2 regularization term that stabilizes the solution but treat the regularization matrix as static, ignoring the intrinsic temporal continuity of neural activity. To overcome this limitation, the authors formulate a full‑rank linear Gaussian state‑space model that incorporates both spatial and temporal correlations of the sources.

The observation equation is Yk = H Jk + εk, where Yk ∈ ℝS is the sensor vector at time k, H ∈ ℝS×3P is the lead‑field matrix derived from a realistic head model (MRI‑based geometry and conductivity), Jk ∈ ℝ3P stacks the three‑component current dipoles for P candidate cortical locations, and εk is zero‑mean Gaussian sensor noise with covariance Σε. The latent dynamics are modeled as a random walk: Jk = Jk‑1 + ωk, with ωk ∼ N(0,Q). The process covariance Q encodes spatial smoothness (e.g., a Laplacian prior linking neighboring dipoles) and determines the strength of temporal smoothing.

Given the linear‑Gaussian structure, optimal Bayesian inference reduces to Kalman filtering (forward pass) and Fixed‑Interval Smoothing (backward pass). The Kalman filter recursively updates the posterior mean and covariance of Jk using all measurements up to time k; the smoother then refines each state estimate by incorporating future data, yielding the minimum‑variance estimate for every time point. This framework naturally yields time‑varying error covariances, unlike MNE’s fixed regularization, and propagates information across the entire trial.

A major computational obstacle is the dimensionality: typical MEG setups have S≈300 sensors and P≈10 000 cortical locations, leading to state vectors of size 30 000 and covariance matrices of size 30 000 × 30 000. Direct dense matrix operations are infeasible on a single workstation. The authors therefore implement the Kalman filter and smoother on the NSF Teragrid supercomputing platform, exploiting MPI‑based parallelism. They partition the lead‑field matrix and covariance updates into blocks distributed across compute nodes, use Cholesky factorizations that respect the symmetry and positive‑definiteness of Q and Σε, and avoid explicit storage of full covariance matrices whenever possible. This high‑performance implementation makes it practical to process whole‑trial MEG recordings (thousands of time points) within a few hours.

Two empirical evaluations are presented. First, synthetic data with known source locations and added Gaussian noise demonstrate that the state‑space approach recovers both the true positions and amplitudes far more accurately than MNE, which exhibits a bias toward superficial sources and underestimates amplitudes. Second, real somatosensory MEG data from human subjects are analyzed. The state‑space estimates reveal a coherent, temporally smooth activation pattern in the primary somatosensory cortex (S1) that aligns with known neurophysiology, whereas MNE produces noisier, less spatially precise maps. Importantly, the state‑space method automatically adapts the regularization strength via the estimated Q and Σε, eliminating the need for manual tuning of a scalar λ.

The authors conclude that (i) incorporating spatiotemporal priors through a linear Gaussian state‑space model substantially improves MEG source localization, (ii) Kalman filtering and fixed‑interval smoothing can be scaled to the high dimensionality of realistic MEG problems using modern HPC resources, and (iii) the approach outperforms static regularization methods in both simulated and real data. They suggest that this framework could be extended to real‑time brain‑computer interfaces, clinical diagnostics such as epilepsy focus localization, and multimodal integration with EEG or fMRI, where exploiting temporal dynamics is essential.