Tutorial 22: Source estimation

Authors: Francois Tadel, Elizabeth Bock, Rey R Ramirez, John C Mosher, Richard M Leahy, Sylvain Baillet

You have in your database a forward model that explains how the cortical sources determine the values on the sensors. This is useful for simulations, but what we need next is to solve the inverse problem: how to estimate the sources when we have the recordings. This tutorial introduces the tools available in Brainstorm for solving this inverse problem.

WARNING: The new interface presented here does not include all the options yet. The mixed head models are not supported: to use them, use the old interface (menu "Compute sources" instead of "Compute sources 2016", described in the ?old tutorials).

Ill-posed problem

Our goal is to estimate the activity of the thousands of dipoles described by our forward model. However we only have a few hundred spatial measurements as input (the number of sensors). This inverse problem is ill-posed, meaning there are an infinite number of source activity patterns that could generate exactly the same sensor topography. Inverting the forward model directly is impossible, unless we add some strong priors to our model.

Wikipedia says: "Inverse problems are some of the most important and well-studied mathematical problems in science and mathematics because they tell us about parameters that we cannot directly observe. They have wide application in optics, radar, acoustics, communication theory, signal processing, medical imaging, computer vision, geophysics, oceanography, astronomy, remote sensing, natural language processing, machine learning, nondestructive testing, and many other fields."

Many solutions to the inverse problem have been proposed in the literature, based on different assumptions on the way the brain works and depending on the amount of information we already have on the effects we are studying. Among the many methods available, in Brainstorm, we present three general approaches to the inverse problem that represent the most widely used methods in MEG/EEG source imaging: minimum-norm solutions, beamformers, and dipole modeling.

These approaches have the advantage of being implemented in an efficient linear form: the activity of the sources is a linear recombination of the MEG/EEG recordings, such that it is possible to solve the inverse problem by applying a linear kernel (in the form of a matrix that multiples the spatial data at each point in time) which is easily stored. Subsequent data manipulation and source visualization is then much simpler, as are comparisons among these techniques.

Below we first describe the minimum norm imaging approach and its options, followed by the beamformer and dipole modeling, both of which are actually quite similar and only use a subset of the options available in the minimum norm approach.

Source estimation options

Before we start estimating the sources for the recordings available in our database, let's start with an overview of the options available. This section focuses on the options for the minimum norm estimates. The other methods are described in advanced sections at the end of this page.



Selecting one of the three methods automatically switches the input screen to reveal the options available for each. The three inverse methods and their required inputs can be summarized as follows.

Minimum norm

LCMV beamformer

Measure [add link]

The minimum norm estimate computed by BrainStorm represents a measure of the current found in each point of the source grid (either volume or surface). As discussed elsewhere in the [rl1] , the units are strictly kept in A-m, i.e. we do not normalize by area (yielding A/m, i.e. a surface density) or volume (yielding A/m^2, i.e. a volume density); nonetheless, it is common to refer these units as a "source density" or "current density" maps when displayed directly.

More commonly, however, current density maps are normalized. The value of the estimated current density is normalized at each source location by a function of either the noise or data covariance. Practically, this normalization has the effect of compensating for the effect of depth dependent sensitivity and resolution of both EEG and MEG. Current density maps tend to preferentially place source activity in superficial regions of cortex, and resolution drops markedly with sources in deeper sulci. Normalization tends to reduce these effects as nicely shown by ref. We have implemented the two most common normalization methods: dSPM and sLORETA.

Source model: Dipole orientations

At each point in the source grid, the current dipole may point arbitrarily in three directions. In this section of the options, we describe alternatives for constraining orientation:


We automatically detect and display the sensors found in your head model. In the example above, two types of magnetometers are found, gradiometers and magnetometers. You can select one or all of the sensors found in your model, such as MEG and EEG.

However, cross-modality calculations are quite dependent on the accuracy by which you have provided adequate covariance calculations and consistency of the head models across sensor types. As of Fall of 2016, we have also elected to NOT account for cross-covariances between different sensor types, since regularization and stability of cross-modalities is quite involved. For multiple sensor types, the recommendation is that you try each individually and then combined, to test for discordance.

Computing sources for an average

Using the above selections, we now discuss explicit directions on how to compute and visualize.

Display: Cortex surface

Why does it look so noisy?

The source maps look very noisy and discontinuous, they show a lot of disconnected patches. This is due to the orientation constraint we imposed on the dipoles orientations. Each value on the cortex should be interpreted as a vector, oriented perpendicular to the surface. Because of the brain’s circumvolutions, neighboring sources can have significantly different orientations, which also causes the forward model response to change quickly with position. As a result, the orientation-constrained minimum norm solution can produce solutions that vary rapidly with position on the cortex resulting in the noisy and disjointed appearance.

It is therefore important not to always interpret disconnected colored patches as independent sources. You cannot expect highly spatial resolution with this technique (~5-10mm at best). Most of the time, a cluster of disconnected source patches in the same neighborhood that show the same evolution in time can be interpreted as "there is some significant activity around here, but with some uncertainty as to its precise location".

To get more continuous maps for visualization or publication purposes, you can either smooth the values explicitly on the surface (process "Sources > Spatial smoothing") or use unconstrained source models.

For data exploration, orientation-constrained solutions may be a good enough representation of brain activity, mostly because it is fast and efficient. You can often get a better feeling of the underlying brain activity patterns by making short interactive movies: click on the figure, then hold the left or right arrows of your keyboard.

Activity patterns will also look sharper when we compute dSPM or sLORETA normalized measures (later in this tutorial). In most of the screen captures in the following sections, the contrast of the figures has been enhanced for illustration purposes. Don't worry if it looks a lot less colorful on your screen. Of course, ultimately statistical analysis of these maps is required to make scientific inferences from your data.

Display: MRI Viewer

Display: MRI 3D

Sign of constrained maps

You should pay attention to the sign of the current amplitudes that are given by the minimum norm method: they can be positive or negative and they oscillate around zero. Display the sources on the surface, set the amplitude threshold to 0%, then configure the colormap to show relative values (uncheck the "Absolute values" option), you will see those typical stripes of positive and negative values around the sulci. Double-click on the colorbar after testing this to reset the colormap.

This pattern is due to the orientation constraint imposed on the dipoles. On both sides of a sulcus, we have defined dipoles that are very close to each other, but with opposite orientations. If we have a pattern of activity on one side of a suclus that can be modeled as a current dipole (green arrow), the limited spatial resolution of the minimum norm model will blur this source using the dipoles that are available in the head model (red and blue arrows). Because of the dipoles’ orientations, the minimum norm images produces positive values (red arrows) on one side of the sulcus and negative on the other side (blue arrows).

When displaying the cortical maps at one time point, we are usually not interested in the sign of the minimum norm values but rather by their amplitude. This is why we always display them by default with the colormap option "absolute values" selected.

However, we cannot simply discard the sign of these values because we need these for other types of analysis, typically time-frequency decompositions and connectivity analysis. For estimating frequency measures on the source maps it is essential that we retain the sign of the time course at each location so that the correct oscillatory frequencies are identified.

Unconstrained orientations

In cases where the orientation constraint imposed on the dipole orientations produces implausible results, it is possible to relax it partially (option "loose constraints") or completely (option "unconstrained"). This produces a vector (3 component) current source at each location which can complicate interpretation, but avoids some of the noisy and discontinous features in the current map that are often seen in the constrained maps. Unconstrained solutions are particularly appropriate when using the MNI template instead of the subject's anatomy, or when studying deeper or non-cortical brain regions for which the normal to the cortical surface obtained with FreeSurfer or BrainSuite is unlikely to match any physiological reality.

In terms of data representation, the option "unconstrained" and "loose constraints" are very similar. Instead of using one dipole at each cortical location, a base of three orthogonal dipoles is used.
Here we will only illustrate the fully unconstrained case.

Source map normalization

The current density values returned by the minimum norm method have a few problems:

Normalizing the current density maps with respect to a reference level (estimated from noise recordings, pre-stimulus baseline or resting state recordings) can help with all these issues at the same time. In the case of dSPM and sLORETA, the normalizations are computed as part of the inverse routine and based on noise and data covariances, respectively. While dSPM does produce a Z-score map, we also provide an explicit Z-score normalization that offers the user more flexibility in defining a baseline period over which BrainStorm computes the standard deviation for normalization.

The normalization options do not change the temporal dynamics of your results when considering a single location but they do alter the relative scaling of each point in the min norm map. If you look at the time series associated with one given source, it will be exactly the same for all normalizations, except for a scaling factor. Only the relative weights change between the sources, and these weights do not change over time.



Typical recommendations

Delete your experiments

Computing sources for single trials

Because the minimum norm model is linear, we can compute an inverse model independently from the recordings and apply it on the recordings when needed. We will now illustrate how to compute a shared inverse model for all the imported epochs.

Averaging in source space

Computing the average

Visualization filters

Low-pass filter

Z-score normalization


Averaging normalized values

Averaging normalized source maps within a single subject requires more attention than averaging current density maps. The amplitude of the normalized measures increase with the SNR of the signal, the higher the SNR the higher the normalized score. For instance, the average of the dSPM for the single trials is lower than the dSPM of the averaged trials (by a factor of sqrt(N), where N is the number of trials).





Display: Contact sheets and movies

A good way to represent what is happening in time is to generate contact sheets or videos. Right-click on any figure and go to the menu Snapshot to check out all the possible options. For a nicer result, take some time to adjust the size of the figure, the amplitude threshold and the colormap options (hiding the colorbar can be a good option for contact sheets).

A time stamp is added to the captured figure. The size of the text font is fixed, so if you want it to be readable in the contact sheet, you should make you figure very small before starting the capture. The screen captures below where produced with the colormap "hot".


Model evaluation

One way to evaluate the accuracy of the source reconstruction if to simulate recordings using the estimated source maps. This is done simply by multiplying the source time series with the forward model:
MEG_simulated [Nmeg x Ntime] = Forward_model [Nmeg x Nsources] * MN_sources [Nsources x Ntime]
Then you can compare visually the original MEG recordings with the simulated ones. More formally, you can compute an error measure from the residuals (recordings - simulated).

To simulate MEG recordings from a minimum norm source model, right-click on the source file, then select the menu "Model evaluation > Simulate recordings".

Open side-by-side the original and simulated MEG recordings for the same condition:


Advanced minimum norm options

Right-click on the deviant average for Run#01 > Compute sources [2016].
Click on the button [Show details] to bring up all the advanced minimum norm options.


Depth weighting

Briefly, the use of various depth weightings was far more debated in the 1990s, before the introduction of MNE normalization via dSPM, sLORETA, and other "z-scoring" methods, which mostly cancel the effects of depth weighting (put another way, after normalization min norm results tend to look quite similar whether depth weighting is used or not). By modifing the source covariance model at each point in the source grid, deeper sources are "boosted" to increase their signal strength relative to the shallower dipoles; otherwise, the resulting MNE current density maps are too dominated by the shallower sources. If using dSPM or sLORETA, little difference in using depth weighting should be noted. To understand how to set these parameters, please refer to the MNE manual. (options --depth, --weightexp and --weightlimit).

Noise covariance regularization

MNE and dipole modeling are best done with an accurate model of the noise covariance, which is generally computed from experimental data. As such, these estimates are themselves prone to errors that arise from relatively too few data points, weak sensors, and strange data dependencies that can cause the eigenspectrum of the covariance matrix to be illconditioned (i.e. a large eigenvalue spread or matrix condition number). In Brainstorm, we provide several means to "stabilize" or "regularize" the noise covariance matrix, so that source estimation calculations are more robust to small errors.

Regularization parameter

In MNE, as mentioned above in the comparisons among Methods, the data covariance matrix is essentially synthesized by adding the noise covariance matrix to a modeled signal covariance matrix. The signal covariance matrix is generated by passing the source prior through the forward (head) model. The source prior is in turn prescribed by the source model orientation and the depth weighting. A final regularization parameter, however, determines how much weight the signal model should be given relative to the noise model, i.e. the "signal to noise ratio" (SNR). In Brainstorm, as of Fall 2016, we follow the definition of SNR as first defined in the original MNE software of Hamalainen. The signal covariance matrix is "whitened" by the noise covariance matrix, such that the whitened eigenspectrum has elements in terms of SNR (power). We find the mean of this spectrum, then take the square root to yield the average SNR (amplitude). The default in MNE and in Brainstorm is "3", i.e. the average SNR (power) is 9.

Output mode

As mentioned above, these methods create a convenient linear imaging kernel that is "tall" in the number of elemental dipoles (one or three per grid point) and "wide" only in the number of sensors. At subsequent visualzation time, we efficiently multiply the kernel with the data matrix to compute the min norm images.

For some custom purposes, however, a user may find it convenient to pre-multiply the data matrix and generate the full source estimation matrix. This would only be recommended in small data sets, since the full results can become quite large.

LCMV beamformer

As mentioned in the introduction above, two other Methods can be selected for source estimation, a beamformer and dipole modeling. In this section, we review the options for the beamformer. The LMCV beamformer uses an estimated data covariance matrix to compute spatial filters that define the beamformer. It is important that this covariance matrix contains contributions from all sources that contribute to the data over the time window to be analyzed. We have found, has have others, that it is important use a relatively long time window for robust beamformer performance. This also helps to avoid numerical instability when inverting the data covariance, as required in computing the beamformer kernel matrix. Our recommendation for evoked responses to use a window that spans prestim through the end of the response of interest, with a minimum of 500ms total duration. For analysis of spontaneous data, the full time window to be analyzed should be used, with a minimum of at least 30 seconds of data.

By selecting "LCMV beamformer, the options automatically switch to the figure below (shown with details).



The Measure has switched to a single option, the "Pseudo Neural Activity Index," (PNAI), after the definition of the Neural Activity Index (NAI). We have modified Van Veen’s definition to rely strictly on the data covariance, without need for a separate noise covariance matrix, but the basic premise is the same as in dSPM, sLORETA, and other normalizations. Viewing the resulting "map," in an identical manner to that with MNE, dSPM, and sLORETA described above, reveals possibly multiple sources as peaks in the map. The PNAI scores analogously to z-scoring. The PNAI can then be dropped into the Process Box and the optimal dipole and orientations extracted for every time instance as described in where (advanced tutorial)?

Source Model: Dipole orientations

The definitions here are identical to the MNE above; however, we do not recommend you select "constrained" (although technically allowed). Choose only "unconstrained" and let the later "dipole scanning" process, which finds the best fitting dipole at each time point, optimize the orientation with respect to the data.


Same considerations as in MNE.

Data covariance regularization

Same definitions as in MNE, only applied to the data covariance matrix, rather than the noise covariance marix. Our recommendation is to use median eigenvalue.

Output mode

We recommend you use kernel only.

Dipole modeling

Dipole modeling fits a single dipole at each potential source location to produce a dipole scanning map. This map can be viewed as a indication of how well, and where, the dipole fits at each time point; however, we recommend using the subsequent best-dipole fitting routine to determine the final location and orientation of the dipole (one per time point). Please note that this function does not fit multiple simultaneous dipoles.

Although not widely recognized, dipole modeling and beamforming are more alike than they are different – when comparing the inverse operators required to compute the dipole scanning map (dipole modeling) and the beamformer output map (LCMV), we see that they differ only in that the former uses an inverse noise covariance matrix while the latter replaces this with the inverse of the data covariance.

Below is the options screen when dipole modeling is selected.


You will notice that "Measure" is now missing, but the resulting imaging kernel file is directly analogous to the PNAI result from LCMV beamforming. The user can display this scanning measure just as with the LCMV case, where again the normalization and uunits are a form of z-scoring.

Source model, Sensors, Noise covariance regularization, and output mode all follow from the explanations given above for the LCMV and for the MNE.

Again, the recommendation is to ALWAYS use "unconstrained" source modeling and let the Process "dipole scanning" optimize the orientation of the dipole for every time instance. Similarly, use "median eigenvalue" for the optimization.

The tutorial "MEG current phantom (Elekta)" demonstrates dipole modeling of 32 individual dipoles under realistic experimental noise conditions.

Equations [TODO]



On the hard drive

Constrained shared kernel

Right-click on a shared inverse file in the database explorer > File > View file contents.

Structure of the source files: results_*.mat

Mandatory fields:

Optional fields:

Full source maps

In Intra-subject, right-click on one of the normalized averages > File > View file contents.

This file has the same structure as a shared inverse kernel, with the following differences:

Filename tags

Useful functions

Additional documentation [TODO]

Articles [TODO]


Forum discussions

Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):

Tutorials/SourceEstimation (last edited 2016-10-17 15:05:27 by FrancoisTadel)