Bayesian analysis of the neuromagnetic inverse problem with ℓ^p-norm priors

Abstract

Magnetoencephalography (MEG) allows millisecond-scale non-invasive measurement of magnetic fields generated by neural currents in the brain. However, localization of the underlying current sources is ambiguous due to the so-called inverse problem. The most widely used source localization methods (i.e., minimum-norm and minimum-current estimates (MNE and MCE) and equivalent current dipole (ECD) fitting) require ad hoc determination of the cortical current distribution (ℓ²-, ℓ¹-norm priors and point-sized dipolar, respectively). In this article, we perform a Bayesian analysis of the MEG inverse problem with ℓ^p-norm priors for the current sources. This way, we circumvent the arbitrary choice between ℓ¹- and ℓ²-norm prior, which is instead rendered automatically based on the data. By obtaining numerical samples from the joint posterior probability distribution of the source current parameters and model hyperparameters (such as the ℓ^p-norm order p) using Markov chain Monte Carlo (MCMC) methods, we calculated the spatial inverse estimates as expectation values of the source current parameters integrated over the hyperparameters. Real MEG data and simulated (known) source currents with realistic MRI-based cortical geometry and 306-channel MEG sensor array were used. While the proposed model is sensitive to source space discretization size and computationally rather heavy, it is mathematically straightforward, thus allowing incorporation of, for instance, a priori functional magnetic resonance imaging (fMRI) information.

Introduction

Magnetoencephalography (MEG) allows non-invasive measurement of the magnetic fields generated by neural activity of the living brain (e.g., Baillet et al., 2001, Hämäläinen et al., 1993, Vrba and Robinson, 2001). Along with clinical applications, MEG is used in studies of basic sensory (auditory, visual, and somatosensory) processes as well as cognitive functions. Time resolution of this method is excellent (∼milliseconds), but in order to locate the underlying source currents accurately on the basis of MEG data, one needs to solve the so-called electromagnetic inverse problem, which does not have a unique solution (Sarvas, 1987). Therefore, additional constraints are needed to select the most feasible estimate from the multitude of possible solutions.

A traditional approach to the MEG inverse problem is to employ the equivalent current dipole (ECD) model, which relies on the assumption that the extents of the activated areas are small enough to be adequately modeled with dipolar point-like sources. Using fully automatic or manually guided, often partly heuristic, fitting methods, the model giving best fit to the measured data is obtained. A downside is that the number and locations of the source dipoles need to be known to a certain extent (although, see Mosher et al., 1992). This is a problem especially when complex cognitive brain functions are studied.

Other widely used methods employ distributed source current estimates (e.g., Hämäläinen et al., 1993, Pascual-Marqui, 2002, Uutela et al., 1999). In the well-known minimum-norm (Dale and Sereno, 1993, Dale et al., 2000, Hämäläinen and Ilmoniemi, 1984, Hauk, 2004) and minimum-current (Uutela et al., 1999) estimates (MNE and MCE), extra information is embedded to the model as mathematical ℓ²- and ℓ¹-norm constraints on the source currents, respectively. Specifically, the least squares error function is combined with an additional penalty term consisting of a weighted norm of the current distribution. Unlike dipole fitting, the exact number and approximate locations of the sources do not need to be known in advance. However, the resulting estimate may be quite diffuse, especially in the case of the minimum-norm estimate and, therefore, it may be equally difficult to discern the number of distinct activated areas in practice.

In Bayesian interpretation, MNE and MCE correspond to ℓ²- and ℓ¹-norm priors for the source currents with a Gaussian likelihood for the measurements (Uutela et al., 1999). The use of predefined values 1 or 2 for the ℓ^p-norm order p is somewhat arbitrary as it leads to prior-wise feasible inverse models even though any value between 1 and 2 could be used. The ℓ²-norm prior produces overly smooth and widely spread estimates whereas ℓ¹-norm estimates might be too focal. The choice of p is subject to uncertainty, hence p should be treated as an unknown variable utilizing Bayesian inference, which has lately gained popularity in solving the electromagnetic inverse problem (e.g., Baillet and Garnero, 1997, Phillips et al., 1997, Schmidt et al., 1999). Markov chain Monte Carlo (MCMC) methods have become popular in this methodology due to rapid expansion of computing resources (e.g., Kincses et al., 2003, Schmidt et al., 1999).

In this paper, we perform a Bayesian analysis of the MEG inverse problem with ℓ^p-norm priors, using MCMC methods and simulated source currents with a realistic MRI-based forward head model. Furthermore, we apply the model on a set of real MEG measurement data. The purpose of this study is to focus on the Bayesian interpretation of the problem, determine an optimal source space discretization size when the discretized points are assumed independent of each other, and to determine whether there is enough information in the data to clarify which ℓ^p-norm prior should be used. We specifically hypothesize that there is no single value for p that would be optimal for all cases, but instead the value depends on the grid discretization size and also on the underlying source configuration and, therefore, it should be inferred from the data rather than determined ad hoc.