Tutorial 22: Source estimation
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. (For backward compatibility, see the old tutorials)
- Ill-posed problem
- Source estimation options
- Computing sources for an average
- Display: Cortex surface
- Why does it look so noisy?
- Display: MRI Viewer
- Display: MRI 3D
- Sign of constrained maps
- Unconstrained orientations
- Source map normalization
- Delete your experiments
- Computing sources for single trials
- Averaging in source space
- Note for beginners
- Averaging normalized values
- Display: Contact sheets and movies
- Model evaluation
- Advanced options: Minimum norm
- Advanced options: LCMV beamformer
- Advanced options: Dipole modeling
- Combining MEG+EEG for source estimation
- On the hard drive
- Additional documentation
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.
Minimum norm imaging
Estimates the sources as the solution to a linear imaging problem, that can be interpreted in various ways (Tikhonov regularization, MAP estimation). The method finds a cortical current source density image that approximately fits the data when mapped through the forward model. The "illposedness" is dealt with by introducing a regularizer or prior in the form of a source covariance that favors solutions that are of minimum energy (or L2 norm).
Min norm requires specification of a noise and a source covariance matrix. Users can estimate a noise covariance matrix directly from recordings (for example, using pre-stim recordings in event related studies) or simply assume a white-noise identify matrix covariance as described below.The source covariance prior is generated from the options discussed in detail below.
- In contrast to the LCMV beamformer, in which the data covariance is estimated directly from the data, for minimum norm the data covariance is determine by the choice of source and data covariances and the forward model.
- Linearly constrained minimum variance (LCMV) beamformers compute an estimate of source activity at each location through spatial filtering. The spatial data are linearly combined with weights (the spatial filter) chosen separately for each location to ensure that the strength of a dipolar source at that location is correctly estimated (assuming a perfect head model).
- The remaining degrees of freedom in selecting the weights are used to minimize the total output power. This has the effect of suppressing contributions of sources from other locations to the estimated signal at the location of interest.
- It should be noted, however, that correlation between sources can at times lead to partial or full signal cancellation, and the method can be sensitive to accuracy of the head model.
LCMV beamformers require specification of the data covariance matrix, which is assumed to include contributions from background noise and the brain signals of interest. In practice, the data covariance is estimated directly from the recordings. A linear kernel (matrix) is formed from this data covariance matrix and the forward model. This kernel defines the spatial filters applied at each location. Multiplying by the data produces an output beamformer scanning image. These images can either be used directly, as is common practice with LCMV methods, or the largest peak(s) can be fit with a dipolar model at every time instance. See section below on LCMV Beamformer Modeling.
Dipole modeling [TODO]
- In some sense this is the simplest model: we fit a single current dipole at each point in time to the data. We do this by computing a linear kernel (similar to the min norm and LCMV methods) which when multiplied by the data produces a dipole scanning image whose strongest peak represents the most likely location of a dipolar source.
As with LMCV, the dipole scanning images can be viewed directly, or the single best dipole fit (location and orientation) computed, as described in (LINK ?). More details.
- Still under much debate, even among our Brainstorm team. In cases where sources are expected to be focal (e.g. interictal spikes in epileptic patients, or early components of sensory evoked responses) the single dipole can be precise in terms of localization. For cases where sources are expected to be distributed, the min norm method makes the least restrictive source assumptions. LCMV beamformers fall somewhere between these two cases.
- One advantage of Brainstorm is that all three approaches can be easily run and compared. If the results are concordant among all three techniques, then our underlying assumptions of source modeling, head modeling, and data statistics are confirmed. If the results are disparate, then a more in depth study is needed to understand the consequences of our assumptions and therefore which technique may be preferred. The next several sections discuss in detail the options associated with the "mininum norm imaging" method.
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 on the user forum, 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 (LINK ?). We have implemented the two most common normalization methods: dSPM and sLORETA.
Current density map: Produces a "depth-weighted" linear L2-minimum norm estimate current density using the method also implemented in Matti Hamalainen's MNE software. For a full description of this method, please refer to the MNE manual, section 6, "The current estimates". Units: picoamper-meter (pA-m).
dSPM: Implements dynamical Statistical Parametric Mapping (Dale, 2000). The MNE is computed as above. The noise covariance and linear inverse kernel are then used to also compute estimates of noise variance at each location in the current density map. The MNE current density map is normalized by the square root (standard deviation) of these variance estimates. As a result dSPM gives a z-score statistical map. Units: unitless "z".
sLORETA: Standardized LOw Resolution brain Electromagnetic TomogrAphy (Pasqual-Marqui, 2002). As with dSPM, the MNE current density map is normalized at each point. While dSPM computes the normalization based on the noise covariance, sLORETA replaces the noise covariance with the theoretical data covariance, as is assumed in the minimum norm estimation. The theoretical data covariance is the noise covariance plus the theoretical signal covariance. As discussed in (Pasqual-Marqui 2002), this theoretical data covariance simplifies sLORETA to an alternative form that results in a "resolution" kernel (eq.(17) of (Pasqual-Marqui 2002). (We note that the theoretical data covariance is not the experimental data covariance estimated directly from the data, as is used in beamformers). Units: unitless.
Recommended option: Discussed in the section Source map normalization below.
Source model: Dipole orientations [TODO]
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:
Constrained: Normal to cortex: Only for "surface" grids. At each grid point, we model only one dipole, oriented normally to the cortical surface. This is based on the anatomical observation that in the cortex, the pyramidal neurons are mainly organized in macro-columns that are perpendicular to the cortex surface.
Size of the inverse operator: [Nvertices x Nchannels].
Loose: Only for "surface" grids. As introduced by (LINK ?), at each point in the surface grid the dipole direction is constrained to be normal to the local cortical surface. Two additional elemental dipoles are also allowed, in the two directions tangential to the cortical surface. As contrasted with "unconstrained," these two tangential elemental dipoles are constrained to have an amplitude that is a fraction of the normal dipole, recommended to be between 0.1 and 0.6. Thus the dipole is only "loosely" constrained to be normal to the local cortical surface.
Size of the inverse operator: [3*Nvertices x Nchannel].
Unconstrained: Either "surface" or "volume" grids. At each grid point, we leave undefined the assumed orientation of the source, such that three "elemental" dipoles are needed to model the source. In Brainstorm, our elemental dipoles are in the x, y, and z ("Cartesian") directions, as compared to other software that may employ polar coordinates. Thus for "N" vertices, we are calculating the estimate for "3*N" elemental dipoles.
Size of the inverse operator: [3*Nvertices x Nchannels].
Recommended option: The constrained options use one dipole per grid point instead of three, therefore the source files are smaller, faster to compute and display, and more intuitive to process because we don't have to think about recombining the three values into one. On the other hand, in the cases where its physiological assumptions are not verified, typically when using an MNI template instead of the anatomy of the subject, the normal orientation constraint may fail to represent certain activity patterns. Unconstrained models can help in those cases. See further discussion on constrained vs unconstrained solutions below in section Why does it looks so noisy.
We automatically detect and display the sensors found in your head model. In the example above, only one type of sensors is found ("MEG"). 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 Spring of 2018, 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.
In Run#01, right-click on the average response for the deviant stim > Compute sources .
Select the options: Minimum norm imaging, Current density map, Constrained: Normal to cortex.
- The other menu "Compute sources" launches the interface that was used previously in Brainstorm. We are going to keep maintaining the two implementations in parallel for a while for compatibility and cross-validation purposes.
The result of the computation is displayed as a dependent file of the deviant average because it is related only to this file. In the file comment, "MN" stands for minimum norm and "Constr" stands for "Constrained: normal orientation".
Display: Cortex surface
Right-click on the sources for the deviant average > Cortical activations > Display on cortex.
Double-click on the recordings for the deviant average to have a time reference.
In the filter tab, add a low-pass filter at 40Hz.
- Change the current time (click on the time series figure or use the keyboard arrows) and note it updates the source maps in the 3D figure. You can also use all the menus and shortcuts introduced in the anatomy tutorial (like setting the view with the keys from 0 to 6).
- You can edit the display properties in the Surface tab:
Min size: Hide all the small activated regions, ie. the connected color patches that contain a number of vertices smaller than this "min size" value.
Transparency: Change the transparency of the source activity on the cortex surface.
Take a few minutes to understand what the amplitude threshold represents.
On the screen capture below, the threshold value is set to 20%. It means that only the sources that have a value over 0.20*150 = 30 pA.m are visible.
The threshold level is indicated in the colorbar with a horizontal white line.
At the first response peak (91ms), the sources with high amplitudes are located around the primary auditory cortex, bilaterally, which is what we are expecting for an auditory stimulation.
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 high 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
Right-click on the sources for the deviant average > Cortical activations > Display on MRI (MRI Viewer).
This figure shows the sources computed on the cortical surface and re-interpolated in the MRI volume. If you set the amplitude threshold to 0%, you would see the thin layer of cortex in which the dipoles where estimated.
- You can configure this figure with the following options:
MIP Functional: Same as for MIP Anatomy, but with the layer of functional values.
Smooth level: The sources values can be smoothed after being re-interpolated in the volume. Right-click on the figure to define the size of the smoothing kernel (in number of slices).
Amplitude threshold: In the Surface tab of the Brainstorm window.
Current time: At the top-right of the Brainstorm window (or use the time series figure).
Display: MRI 3D
Right-click on the sources for the deviant average > Cortical activations > Display on MRI (3D).
This view was also introduced in the tutorials about MRI and surface visualization.
Right-click and move your mouse to move the slices (or use the Resect panel of the Surface tab).
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.
Select the options: Minimum norm imaging, Current density map, Unconstrained.
S = sqrt(Sx2 + Sy2 + Sz2)
Source map normalization
The current density values returned by the minimum norm method have a few problems:
- The values tend to be higher at the surface of the brain (close to the sensors).
- The maps are sometimes patchy and difficult to read.
Select successively the two normalization options: dSPM, sLORETA, (constrained).
Double-click on all of them to compare them (screen capture at 143ms):
Current density maps: Tends to highlight the top of the gyri and the superficial sources.
In Process1: Select the constrained current density maps (file MN: MEG(Constr)).
Do not select "Use absolute values": We want the sign of the current values.
Double-click on the new normalized file to display it on the cortex (file with the "| zscore" tag).
Use non-normalized current density maps for:
- Computing shared kernels applied to single trials.
- Averaging files across MEG runs.
- Computing time-frequency decompositions or connectivity measures on the single trials.
Use normalized maps (dSPM, sLORETA, Z-score) for:
- Estimating the sources for an average response.
- Exploring visually the average response (ERP/ERF) at the source level.
- Normalizing the subject averages before a group analysis.
- Avoid averaging normalized maps (or computing any additional statistics)
- Recommended normalization approach:
Delete your experiments
Select all the source files you computed until now and delete them.
Computing sources for single trials
Select: Minimum norm imaging, Current density map, Constrained: Normal to cortex
Averaging in source space
Computing the average
In Run#02, right-click on the head model or the folder > Compute sources .
Select: Minimum norm imaging, Current density map, Constrained: Normal to cortex
The two following approaches are equivalent:
- Averaging the sources of all the individual trials across runs,
- Averaging the sources for the sensor averages that we already computed for each run.
Select "By trial group (subject average)" to average together files with similar names.
Select "Arithmetic average" function.
Check "Weighted average" to account for the different numbers of trials in both runs.
Right-click on each of them > File > View file contents, and check the Leff field:
78 trials for the deviant condition, 378 trials for the standard condition.
Open the sensor-level averages as a time reference.
Use the predefined view "Left, Right" for looking at the two sides at the same time (shortcut: "7").
- Clear the Process1 list, then drag and drop the new averages in it.
Run process "Pre-process > Band-pass filter": [0,40] Hz
- In Process1, select the two filtered averages.
Run process "Standardize > Baseline normalization", baseline=[-100,-1.7]ms, Z-score.
Four new files are accessible in the database: two filtered and two filtered+normalized.
Double-click on the source averages to display them (deviant=top, standard=bottom).
Delete the non-normalized filtered files, we will not use them in the following tutorials.
Note for beginners
Everything below is advanced documentation, you can skip it for now.
Averaging normalized values
Average the current density maps, then normalize.
Display: Contact sheets and movies
Contact sheet: Right-click on any figure > Snapshot > Time contact sheet: Figure
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).
Open side-by-side the original and simulated MEG recordings for the same condition:
Advanced options: Minimum norm
Click on the button [Show details] to bring up all the advanced minimum norm options.
Noise covariance regularization [TODO]
Regularization parameter [TODO]
Signal-to-noise ratio: Use SNR of 3 as the classic recommendation, as discussed above.
- Full results [15000x361] = Inverse kernel [15000x274] * Recordings [274x361]
Advanced options: 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. On top of the noise covariance matrix, you need to estimate a data covariance matrix in order to enable the option "LCMV beamformer" in the interface. Note that pre-whitening with the noise covariance matrix has not yet been implemented for the LCMV beamformer, and only the data covariance is used in the current version. Thus, until noise pre-whitening is implemented, select the "No noise modeling (identity matrix)" option in the contextual menu for the noise covariance to enable the option "LCMV beamformer" in the "Compute sources" interface.
Data covariance regularization
Advanced options: Dipole modeling
Noise covariance regularization
Similarly, use "median eigenvalue".
Combining MEG+EEG for source estimation
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
Time: [1 x Ntime] Time values for each sample recorded in F, in seconds.
HeadModelType: Type of source space used for this head model ('surface', 'volume', 'mixed').
HeadModelFile: Relative path to the head model used to compute the sources.
SurfaceFile: Relative path to the cortex surface file related with this head model.
Atlas: Used only by the process "Sources > Downsample to atlas".
GridAtlas: Atlas "Source model" used with mixed source models.
GoodChannel: [1 x Nchannels] Array of channel indices used to estimate the sources.
Comment: String displayed in the database explorer to represent this file.
History: Operations performed on the file since it was create (menu "View file history").
Whitener: Noise covariance whitener computed in bst_inverse_linear_2018.m
DataWhitener: Data covariance whitener computed in bst_inverse_linear_2018.m
Std: For averaged files, number of trials that were used to compute this file.
DisplayUnits: String, force the display of this file using a specific type of units.
ChannelFlag: [Nchannels x 1] Copy of the ChannelFlag field from the original data file.
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:
It contains the full time series (ImageGridAmp) instead of an inverse operator (ImagingKernel).
The Z-score process updated the field Function ('mn' => 'zscore')
Minimum norm: Baillet S, Mosher JC, Leahy RM
Electromagnetic brain mapping, IEEE SP MAG 2001.
Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron 2000 Apr, 26(1):55-67
Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details, Methods Find Exp Clin Pharmacol 2002, 24 Suppl D:5-12
Westner BU, Dalal SS, Gramfort A, Litvak V, Mosher JC, Oostenveld R, Schoffelen JM
A unified view on beamformers for M/EEG source reconstruction,
NeuroImage. 2022 Feb 1;246:118789
Tutorial: Volume source estimation
Tutorial: Deep cerebral structures
Tutorial: Computing and displaying dipoles
Tutorial: Dipole fitting with FieldTrip
Tutorial: Maximum Entropy on the Mean (MEM)
Forum: Minimum norm units (pA.m): http://neuroimage.usc.edu/forums/showthread.php?1246
Forum: Imaging resolution kernels: http://neuroimage.usc.edu/forums/showthread.php?1298
Forum: Spatial smoothing: http://neuroimage.usc.edu/forums/showthread.php?1409
Forum: Units for dSPM/sLORETA: http://neuroimage.usc.edu/forums/showthread.php?1535
Forum: EEG reference: http://neuroimage.usc.edu/forums/showthread.php?1525#post6718
Forum: Sign of the MNE values: http://neuroimage.usc.edu/forums/showthread.php?1649
Forum: Combine MEG+EEG: https://neuroimage.usc.edu/forums/t/combining-eeg-and-meg-for-source-analysis/1684/4
Forum: Combine mag+gradiometers: http://neuroimage.usc.edu/forums/showthread.php?1900
Forum: Combine EEG+fMRI: http://neuroimage.usc.edu/forums/showthread.php?2679
Forum: Residual ocular artifacts: http://neuroimage.usc.edu/forums/showthread.php?1272
Forum: Dipole fitting: http://neuroimage.usc.edu/forums/showthread.php?2400
Forum: Simulate recordings from sources: http://neuroimage.usc.edu/forums/showthread.php?2563
Forum: Simulate recordings from simulated signals: https://neuroimage.usc.edu/forums/t/simulate-scalp-recording/2421/3
Forum: Pre-whitening: https://neuroimage.usc.edu/forums/t/10459
Forum: LCMV beamformer and noise covariance: https://neuroimage.usc.edu/forums/t/30943
Forum: Debugging weird sLORETA results: https://neuroimage.usc.edu/forums/t/21265
Forum: Subset of sensors: https://neuroimage.usc.edu/forums/t/26496