|Deletions are marked like this.||Additions are marked like this.|
|Line 60:||Line 60:|
|* '''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 [[http://www.nmr.mgh.harvard.edu/meg/manuals/MNE-manual-2.7.pdf|MNE manual]], section 6, "The current estimates". Units: picoamper-meter (pA-m).||* '''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 [[https://mne.tools/mne-c-manual/MNE-manual-2.7.3.pdf|MNE manual]], section 6, "The current estimates". Units: picoamper-meter (pA-m).|
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
Source estimation options
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).
- 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.
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]
Size of the inverse operator: [Nvertices x Nchannels].
Size of the inverse operator: [3*Nvertices x Nchannel].
Size of the inverse operator: [3*Nvertices x Nchannels].
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.
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?
Display: MRI Viewer
Right-click on the sources for the deviant average > 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 sources for the deviant average > 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
In Run#01, right-click on the average response for the deviant stim > Compute sources .
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.
In Run#01, right-click on the average recordings for the deviant stim > Compute sources .
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.
dSPM: Tends to correct this behavior and may give higher values in deeper areas. The values obtained are unitless and similar to Z-scores, therefore they are easier to interpret. They are by default not scaled with the number of averages. To obtain correctly scaled dSPM values, one has to call the process "Sources > Scale averaged dSPM", as explained in the advanced section Averaging normalized values.
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
Right-click on the head model or the folder for Run#01 > Compute sources .
Select: Minimum norm imaging, Current density map, Constrained: Normal to cortex
Because we did not request to compute an inverse model for a specific block of recordings, it computed a shared inverse model. If you right-click on this new file, you get a warning message: "Inversion kernel". It does not contain any source map, but only the inverse operator that will allow us to convert the recordings into source maps.
Averaging in source space
Computing the average
First compute the same source model for the the second acquisition run.
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.
Run process "Average > Average files":
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.
The file comments say "2 files" because they were computed from two averages each (one from each run), but the number of corresponding trials is correctly updated in the file structure.
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").
Note that opening the source maps can be very long because of the filters for visualization. Check in the Filter tab, you may have a filter applied with the option "Filter all results" selected. In the case of averaged source maps, the 15,000 source signals are filtered on the fly when you load a source file. This filtering of the full source files can take a significant amount of time, consider unchecking this option if the display is too slow on your computer.
- Clear the Process1 list, then drag and drop the new averages in it.
Run process "Pre-process > Band-pass filter": [0,40] Hz
Epochs are too short: Look at the filter response, the expected transient duration is at least 78ms. The first and last 78ms of the average should be discarded after filtering. However, doing this would get rid of almost all the 100ms baseline, which we need for normalization. As mentioned here, we should have been importing longer epochs in order to filter and normalize the averages properly.
- 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
Averaging normalized source maps within a single subject requires more attention than averaging current density maps. Since averaging reduces variance, the resulting source maps will have a different statistical distribution than the nominal distribution of the individual maps.
For example, averaging z-score normalized maps will result in maps with variance less than 1. The same holds true for dSPM maps. Assuming independent samples, the variance of an average of N maps drops by 1/N. For this reason, it is generally recommended to select the "Weighted average" option in the ‘Average files’ process when averaging trials or source maps (which performs mean(x) = (N1*sum(x1(i)) + N2*sum(x2(i)) + …)/ (N1+N2+…) ) in order to keep track of the number of samples and the actual variance of averaged statistical maps.
- An averaged dSPM map has variance equal to 1/N (and thus standard deviation equal to 1/sqrt(N)). Therefore one could multiply the averaged dSPM map by sqrt(N) in order to maintain variance 1 under the null hypothesis. In previous versions of Brainstorm, this was done automatically when visualizing the files, and when averaging source maps with the option "Adjust normalized source maps for SNR increase". To simplify the interface and make the interpretation of maps more intuitive and consistent with other cases (min-norm, z-scored), we now dropped this option. Thus averaging dSPM maps now results in maps with variance less than 1, and is consistent with handling min-norm, z-scored and sLORETA maps.
Ajusting an averaged dSPM file by this sqrt(N) factor is still possible manually, eg. in order to visualize cortical maps that can be interpreted as Z values. Select the average dSPM files in Process1 and run process "Sources > Scale averaged dSPM". This should be used only for visualization and interpretation, scaled dSPM should never be averaged or used for any other statistical analysis.
- The same SNR issues arise while averaging Z-scores: the average of the Z-scores is lower than the Z-score of the average.
Average the current density maps, then normalize.
- This normalization is not based on the SNR of signal, but rather on the spatial smoothness of the maps. Managing these maps is similar to min-norm maps: the variance of the individual maps is not explicitly modeled or known analytically.
- As in other cases, sLORETA(Average(trials)) = Average(sLORETA(trials)), and this relationship is guaranteed to hold with averaging uneven number of samples when using the option "Weighted average".
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
Right-click on the deviant average for Run#01 > Compute sources .
Click on the button [Show details] to bring up all the advanced minimum norm options.
Noise covariance regularization [TODO]
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_2018.m" for notes and details.
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.
Data covariance regularization
Advanced options: Dipole modeling
Noise covariance regularization
Similarly, use "median eigenvalue".
Combining MEG+EEG for source estimation
Magnetoencephalography and EEG sensor data can be processed jointly to produce combined source estimates. Joint processing presents unique challenges because EEG and MEG use head models that exhibit differing sensitivities to modeling errors, which can in turn lead to inconsistencies between EEG and MEG with respect to the (common) source model. In practice joint processing is relatively rare (Baillet et al., 1999). However, these data are complementary, which means that joint processing can potentially yield insights that cannot be seen with either modality alone.
For example, in the evoked responses in the data set used here, the first peak over the occipital areas is observed in MEG (90 ms) slightly before EEG (110 ms). This delay is too large to be caused by acquisition imprecisions. This indicates that we are not capturing the same brain processes with the two modalities, possibly because the orientation and type of activity in the underlying cortical sources is different.
MEG and EEG have different sensitivities to source orientation and depth. Given the challenges of joint processing, our advice is to first look at the source reconstructions for the two modalities separately before trying to use any type of fusion technique.
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
ImageGridAmp: [Nsources x Ntime] Full source time series, in Amper.meter. If this field is defined, ImagingKernel must be empty. If you want this field to be set instead of ImagingKernel, make sure you select the advanced option Full results when estimating the sources.
Time: [1 x Ntime] Time values for each sample recorded in F, in seconds.
Function: Type of values currently saved in the file: 'mn', 'mnp', 'dspm', 'dspm2018', 'dspm2018sc', 'sloreta', 'lcmv', 'lcmvp', 'lcmvnai', 'lcmvpow', 'gls', 'glsp', 'glsfit', 'glschi', 'zscore', 'ersd'...
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.
Leff: Effective number of averages. For averaged files, number of trials that were used to compute this file. For source files that are attached to a data file, we use the Leff field from the 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
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: 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: Debugging weird sLORETA results: https://neuroimage.usc.edu/forums/t/dont-trust-the-source-power-spectrum-results/21265