Tutorial 22: Source estimation [TODO]
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 is to build the inverse information: how to estimate the sources when we have the recordings. This tutorials introduces the tools available in Brainstorm for solving this inverse problem.
- Ill-posed problem
- Source estimation options [TODO]
- 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
- Averaging normalized values
- Display: Contact sheets and movies
- Model evaluation
- Advanced options [TODO]
- Equations [TODO]
- On the hard drive
- Additional documentation [TODO]
Our goal is to estimate the activity of the thousands of dipoles described by our forward model. However we only have a few hundred variables in input (the number of sensors). This inverse problem is ill-posed, meaning there is an infinite number of combinations of source activity patterns that can generate exactly the same sensor topography. Inverting the forward model directly is impossible, unless we add some strong priors to our model.
Many solutions 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 hundreds of methods available, in Brainstorm, we initially present three general approaches to the inverse models widely used 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 as a linear kernel 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 found in the minimum norm.
Source estimation options [TODO]
Selecting one of the three Methods automatically switches the input screen automatically to reveal the options available for each. Above are the options revealed for the mininum norm imaging. First we brieflly describe the Methods options:
Minimum norm: Bayesian Estimation of the sources, assumes a noise and signal prior.
User may synthesize the noise covariance as the "identity matrix", or estimate a noise covariance matrix from recordings. The source covariance prior is generated from the options discussed in detail below. The signal covariance prior (i.e. the source prior as "propagated" to the sensor array) is generated from the options selected. The data covariance (the sum of the signal and noise covariance) is therefore synthesized.
LCMV beamformer: "Linearly constrained minimum variance," requires a data covariance.
The data covariance matrix is assumed to also contain the signal(s) of interest. In practice, the data covariance is estimated directly from the recordings. A linear kernel is formed from this data covariance matrix, after which the beamformer scanning image can be stepped through the data, and the best peak can be fit with a dipolar model at every time instance, as described later.
Dipole modeling: Single equivalent current dipole, requires a noise covariance.
The noise covariance matrix is assumed to NOT contain the signal of interest. In practice, this noise covariance matrix is calculated directly from recordings. In MEG, these recordings are often "empty room." In EEG, no such "empty room" is possible, so "quiet periods" of ongoing brain data are selected, as can also be selected in MEG. A linear kernel is formed from this noise covariance matrix, after which a scanning image (analogous to the beamformer) can be viewed, and the best dipole location and orientation can be found for every time instance.
Recommended option: Still under much debate, even among our Brainstorm team. 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 indepth 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 as revealed by selecting the "mininum norm imaing" Method above.
As directly formed and computed, the minimum norm estimate produces a measure of the current found in each point of the source grid (either volume or surface). As discussed elsewhere in the forum (link to Mosher's discuss on units of the minimum norm), the units are strictly kept in A-m, i.e. we do not attempt to divide out the area (yielding A/m, i.e. a surface density) or volume (yielding A/m^2, i.e. a volume density) of each source point; nonetheless, it is common to refer these units as a "source density" or "current density" map when displayed directly.
More commonly, however, users often prefer to "normalize" each source point in order to compensate for the rapid falloff in signal intensity for deeper dipoles. Because of their popularity, we provide directly in the Measure box two forms of standardization, the dSPM and sLORETA. Summarizing the Measures selection:
Current density map: The Minimum Norm Estimate (MNE)
"Whitened" (i.e. noise covariance matrix is accounted for) and "depth-weighted" (if user selected the option) linear L2-minimum norm estimates algorithm inspired from 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: dynamical Statistical Parametric Mapping [Dale, 2000]
The MNE is computed as above, but then we use exactly the same linear kernel with given noise covariance matrix to essentially "z-score" each point in the current density map. Units: unitless "z".
sLORETA: standardized LOw Resolution brain Electromagnetic TomogrAphy [Pasqual-Marqui, 2002])
As originally presented (2002), sLORETA also uses the same MNE solution from above, then we use exactly the same linear kernel with the MNE data covariance matrix (instead of noise covariance) to again weight each point in the current density map. Subsequent theoretical developments show an alternative form of a "resolution" kernel that may be calculated instead. We use the "resolution" form, but the interested researcher is invited to review the details of our coding.
Recommended option: Discussed in the section "Source map normalization" below.
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 prescribe if the inverse method will constrain the orientation, leave it unconstrained, or if we impose a "soft" constraint to the 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 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 Lin (cite), 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 only modeled 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 a MNI template instead of the anatomy of the subject, the normal orientation constraint may fail representing certain activity patterns. Unconstrained models can help in those cases.
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. 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.
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 negate the effects of depth weighting. At each point in the source grid, the deeper points 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. Details on depth weighting may be found in the MNE manual. If using dSPM or sLORETA, little difference in using depth weighting should be noted.
Noise covariance regularization [TODO]
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 cause the eigenspectrum of the covariance matrix to be somewhat deficient. 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.
Regularize noise covariance: The L2 matrix norm is defined as the largest eigenvalue of its eigenspectrum. This option adds to the covariance matrix a diagonal matrix whos entries are a fraction of the matrix norm. The default is 0.1, such that covariance matrix is stabilized by adding to it an identity matrix that is scaled to 10% of the largest eigenvalue.
Median eigenvalue: The eigenspectrum of MEG data can often span many decades, due to highly colored spatial noise, but this broad spectrum is generally confined to the first several modes only. Thus the L2 norm is many times greater than the majority of the eigenvalues, and it is difficult to prescribe a conventional regularization parameter. Instability in the inverse is dominated by defects found in the smallest eigenvalues. This approach stabilizes the eigenspectrum by replicating the median (middle) eigenvalue for the remainder of the small eigenvalues.
Diagonal noise covariance: Deficiencies in the eigenspectrum often arise from numerical inter-dependencies found among the channels, particularly in covariance matrices computed from relatively short sequences of data. One common method of stabilization is to simply take the diagonal of the covariance matrix and zero-out the cross-dependencies. Each channel is therefore modeled as independent of the other channels. The eigenspectrum is now simply the (sorted) diagonal values.
No covariance regularization: We simply use the noise covariance matrix as computed or provided by the user.
Automatic shrinkage: Stabilization method of Ledoit and Wolf (2004), still under testing in the Brainstorm environment. Basically tries to estimate a good tradeoff between no regularization and diagonal regularization, using a "shrinkage" factor. See Brainstorm code "bst_inverse_linear_2016.m" for notes and details.
Recommendation: This author (Mosher) votes for the median eigenvalue as being quite generally effective.
Regularization parameter [TODO]
Output mode [TODO]
Computing sources for an average
Using the above selections, we now discuss explicit directions on how to compute and visualize.
Select the options: Minimum norm imaging, Current density map, Constrained: Normal to cortex.
- The other menu "Compute sources" brings 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.
Display: Cortex surface
Right-click on the sources for the deviant average > Cortical activations > Display on cortex.
In the filter tab, add a low-pass filter at 40Hz.
- You can edit the display properties in the Surface tab:
Transparency: Change the transparency of the source activity on the cortex surface.
Take a few minutes to understand what the amplitude threshold represents.
The threshold level is indicated in the colorbar with a horizontal white line.
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 has to be interpreted as a vector, oriented perpendicular to the surface. Because of the brain circumvolutions, all the sources have different orientations, two adjacent sources have very little chance to have the same orientation in this model, therefore the minimum norm method may attribute completely different values to them. This causes all these gaps we see here.
Visually, you should not always interpret disconnected colored patches as independent sources. You cannot expect a very spatial resolution with this technique (~5-10mm). 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 we don't know where exactly".
For data exploration, this is a good enough representation of the brain activity, mostly because it is fast and efficient. You can 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 normalized measures (later in this tutorial). In most of the screen captures in this 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.
Display: MRI Viewer
Right-click on the source file > Cortical activations > Display on MRI (MRI Viewer).
- You can configure this figure with the following options:
MIP Functional: Same as for MIP Anatomy, but with the layer of functional values.
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 source file > Cortical activations > Display on MRI (3D).
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 would 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 assimilated to an electric dipole (green arrow), the minimum norm model will try to explain it with the dipoles that are available in the head model (red and blue arrows). Because of the dipoles orientations, it translates into 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 by 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 them for other types of analysis, typically time-frequency decompositions and connectivity analysis. For estimating frequency measures on the source maps, we need to keep the oscillations around zero.
In the cases where the orientation constraint imposed on the dipoles orientations looks too strong, it is possible to relax it partially (option "loose constraints") or completely (option "unconstrained"). This is typically something to consider when using a MNI template instead of the subject's anatomy, or when studying deeper or non-cortical brain regions for which the normal to the FreeSurfer cortex surface is unlikely to match any physiological reality.
Here we will only illustrate the fully unconstrained case.
Select the options: Minimum norm imaging, Current density map, Unconstrained.
S = sqrt(Sx2 + Sy2 + Sz2)
The maps we observe here look a lot smoother than the constrained sources we computed earlier. This can be explained by the fact that there is no sharp discontinuity between two adjacent points of the grid, while the normal to the surface between two nearby points can be very different.
Source map normalization
The current density values returned by the minimum norm method have a few problems:
Their amplitude is therefore difficult to interpret directly.
- 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.
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. Some normalizations can be computed independently from the recordings, and added to the linear inverse operator (dSPM or sLORETA). Another way of proceeding is to divide the current density maps by the standard deviation estimated over a baseline (Z-score).
The normalization options do not change the temporal dynamics of your results, they are just different ways for looking at the same minimum norm maps. If you look at the time series associated with one given source, it would be exactly the same for all the normalizations, except for a scaling factor. Only the relative weights change between the sources, and these weights do not change over time.
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.
sLORETA: Produces smoother maps where all the potentially activated area of the brain (given to the low spatial resolution of the source localization with MEG/EEG) is shown as connected, regardless of the depth of the sources. However, the units are difficult to interpret.
The Z-transformation converts the current density values to a score of deviation from a baseline. We define a baseline period in our file (in this case, the pre-stimulus baseline) and compute the average and standard deviation for this segment. Then for every time point we subtract the baseline average and divide by the baseline standard deviation. Z = (Data - μ) / σ
- This measure tells how much a value deviates from the baseline average, in number of times the standard deviation. This is done independently for each source, so the sources with a low variability during baseline will be more salient in the cortical maps post-stimulus.
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).
You can see that the cortical maps obtained in this way are very similar to the other normalization approaches, especially with the dSPM maps. The units are different but the global observations are the same.
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.
- Recommended normalization approach:
- It is difficult to declare that one normalization technique is better than another. They have different advantages and may be used in different cases. Ideally, they should all converge to similar observations. If you obtain results with one method that you cannot reproduce with the others, you should question your findings.
- dSPM and sLORETA are linear measures and can expressed as imaging kernels, therefore they are a lot easier to manipulate in Brainstorm. sLORETA maps can be smoother but its units are difficult to understand. dSPM units are much easier to understand and interpret.
- Z-normalized current density maps are also easy to interpret. They represent explicitly a "deviation from experimental baseline" while dSPM indicates a "deviation from noise" (as represented by the noise covariance matrix).
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.
It is a lot faster because it needs to read 4 files (one average file per run and per condition) instead of 456 files (total number of good trials in the two runs).
Select "By trial group (subject average)" to average together files with similar names.
Select "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 nAvg 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.
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.
- From the expression of the dSPM, we know exactly what factor should be applied to compensate for this effect. When computing the average of dSPM or other normalized values, we have to also multiply the average with the square root of the number of files averaged together.
dSPM(Average(trials)) = sqrt(Ntrials) * Average(dSPM(trials))
Average the current density maps, then normalize.
- sLORETA(Average(trials)) = Average(sLORETA(trials))
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 [TODO]
Click on the button [Show details] to bring up all the advanced minimum norm options.
Depth weighting [TODO]
The minimum-norm estimates have a bias towards superficial currents. This tendency can be alleviated by adjusting these parameters. To understand how to set these parameters, please refer to the MNE manual (options --depth, --weightexp and --weightlimit).
Noise covariance regularization [TODO]
Explain the influence of this parameter.
Regularize noise covariance: Regularize the noise-covariance matrix by the given amount for each type of sensor individually (value is restricted to the range 0...1). For more information, please refer to the MNE manual, section 6.2.4 (options --magreg, --gradreg and --eegreg).
Diagonal noise covariance:
No covariance regularization:
Regularization parameter [TODO]
Explain the influence of this parameter.
Signal-to-noise ratio: Fixed SNR...
RMS source amplitude:
- Full results [15000x361] = Inverse kernel [15000x274] * Recordings [274x361]
LCMV beamformer [TODO]
Explain the method briefly.
Neural activity index:
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_2016.m
DataWhitener: Data covariance whitener computed in bst_inverse_linear_2016.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')
Additional documentation [TODO]
Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron 2000 Apr, 26(1):55-67
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 mag+gradiometers: http://neuroimage.usc.edu/forums/showthread.php?1900
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