Tutorial 22: Source estimation
Authors: Francois Tadel, Elizabeth Bock, Rey R Ramirez, John C Mosher, Richard Leahy, Sylvain Baillet
You have in your database a forward model matrix that explains how the cortical sources determine the values on the sensors. This is useful for simulations, but what we really need is to build the inverse information: how to estimate the sources when we have the recordings. This tutorials introduce the tools available in Brainstorm for solving this inverse problem.
- Ill-posed problem
- Minimum norm vs. beamformers
- Computing sources for a single data file
- Sources visualization
- Computing sources for multiple data files
- Minimum norm options
- Additional discussions on the forum
- Inverse model
- Explore the sources
- Rey on sLORETA
- Issues with dSPM average
- On the hard drive
- Additional documentation
Our goal is now to estimate the activity of the 45,000 dipoles described by our forward model. However we only have a few tens or hundreds of variables to estimate this activity (the number of sensors). This inverse problem is ill-posed, there is an infinity of combinations of source activity that can generate exactly the same sensor topography. Inverting the forward problem directly is impossible, unless we add some strong priors in 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 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, two classes of inverse models have been widely used in MEG/EEG source imaging in the past few years: minimum-norm solutions and beamformers.
Both approaches have the advantage of being linear: the activity of the sources is a linear recombination of the MEG/EEG recordings. It is possible to solve the inverse problem independently of the recordings, making the data manipulation a lot easier and faster.
Minimum norm vs. beamformers
The main difference between those two approaches is the type of inputs they need. They both need an estimation of the noise at the level of the sensors (noise covariance matrix), but the beamforming solutions additionally need an accurate representation of the effect we are trying to localize in the brain (data covariance matrix).
Provided that we know at what latencies to look at, we can compute a correct data covariance matrix and may obtain a better spatial accuracy with a beamformer. However, in many cases we don't exactly know what we are looking at, the risks of misinterpretation of beamforming results are high. Brainstorm tends to favor minimum norm solutions, which have the advantage of needing less manual tuning for getting acceptable results.
Both approaches are available in Brainstorm, so you can use the one which is the most adapted to your recordings or to your own personal expertise. Only the minimum norm estimates will be described in this tutorial, but the other solutions work exactly in the same way.
Computing sources for a single data file
With this window you can select the method you want to use to estimate the cortical currents, and the sensors you are going to use for this estimation. The default "Normal mode" only let you edit the following options:
Comment: This field contains what is going to be displayed in the database explorer.
Method: Please select wMNE. The other methods dSPM and sLORETA are also based on wMNE. They may give better and/or smoother results depending on the cases.
Sensors type: Modalities that are used for the reconstruction. Here we only have one type of MEG sensors (axial gradiometers), so nothing to change.
Expert mode: Show other options we do not care about right now.
- Click on Run.
It is displayed inside the recordings file ERF, because it is related to this file only.
- Meaning of that weird filename: "MN" stands for "Minimum Norm", and "Constr" stands for "Constrained orientation" of the dipoles (the estimated dipoles orientations are constrained to be normal to the cortex).
You can have a look to the corresponding matrix file (right-click > File > View file contents). You would find all the options of forward and inverse modeling, and only one interesting field : ImagingKernel, which contains the inversion kernel. It is a [nVertices x nChannels] matrix that has to be multiplied with the recordings matrix in order to get the activity for each source at all the time samples.
The minimum norm solution being a linear operation (the time series for each source is a linear combination of all the time series recorded by the sensors), we make this economy of saving only this linear operator instead of the full source matrix (nVertices x nTime)
Do the same for the Left / ERF file
There are two main ways to display the sources: on the cortex surface and on the MRI slices.
Sources on cortex surface
Double-click on recordings Right / ERF, to display the time series (always nice to have a time reference).
Double-click on sources Right / ERF / MN: MEG.
Equivalent to right-click > Cortical activations > Display on cortex.
Then you can manipulate the sources display exactly the same way as the surfaces and the 2D/3D recordings figures: rotation, zoom, Surface tab(smoothing, sulci, resection...), colormap, sensors, predefined orientations (keys from 0 to 7)...
Three new controls are available in the Surfaces tab, in panel Data options:
Amplitude: Only the sources that have a value superior than a given percentage of the colorbar maximum are displayed.
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 sources on the cortex.
Take a few minutes to understand what this threshold value represents.
The colorbar maximum depends on the way you configured your Sources colormap. In case the colormap is NOT normalized to current time frame, and the maximum is NOT set to a specific value, the colorbar maximum should be around 68 pA.m.
- On the screenshot above, the threshold value was set to 35%. It means that only the sources that had a value over 0.35*68 = 23.8 pA.m were visible.
- If you set the threshold to 0%, you display all the sources values on the cortex surface; and as most of the sources have values close to 0, the brain is mainly blue.
- The figure on the right shows the most active area of the cortex 46ms after an electric stimulation of the right thumb. As expected, it is localized in the left hemisphere, in the middle of post central gyrus (projection of the right hand in the primary somatosensory cortex).
Sources on MRI (3D)
Close all the figures (Close all button). Open the time series view for Right / ERF.
Right-click on Right / ERF / MN: MEG > Cortical activations > Display on MRI (3D).
- A new menu is available in the popup menu of this figure: MRI Display
MIP Anatomy: for each slice, display the maximum value over all the slices instead of the original value in the structural MRI (fig 1)
MIP Functional: same thing but with the layer of functional values (fig 2)
This view can be used to lots of different types of contact sheets: in time or in space, for each orientation. You can try all the menus. Example: Right-click on the figure > Snapshot > Volume contact sheet: axial:
Sources on MRI (MRI Viewer)
Right-click on Right / ERF / MN: MEG > Cortical activations > Display on MRI (MRI Viewer).
Minimum norm values are not only positive
You should pay attention to a property of the current amplitudes that are given by the wMNE method: they can be positive of negative, and they oscillate around zero. It's not easy to figure out what is the exact meaning of a negative value respect with a positive value, and most of the time we are only interested in knowing what is activated at what time, and therefore we look only at the absolute values of the sources.
In some other cases, mainly when doing frequency analysis, we need to pay attention to the sign of these values. Because we cannot do a frequency decomposition of the absolute values of the sources, we need to keep the sign all along our processes.
Display again the sources for Right / ERF on the cortex surface (double-click on the source file), and uncheck the Absolute option for the colormap "Sources" (right-click on the figure > Colormap Sources > Absolute values). Decrease the threshold to observe the pattern of alternance between positive and negative values on the surface. Then double click on the colorbar to reset it to its default.
Minimum norm units
For information about the units used to represent the minimum norm source activation, pA.m (pico Ampere meter), please refer to the following forum post:
Computing sources for multiple data files
The sources file we are observing was computed as an inversion kernel. It means that we can apply it to any similar recordings file (same subject, same run, same positions of sensors). But in our TutorialCTF database, the MN: MEG node only appears in the the ERF file, not in the Std one. What is it necessary to share an inversion kernel between different recordings ?
Compute another source estimation: but instead of clicking on the Compute sources from the ERF recordings popup menu (which would mean that you only want sources for this particular recordings file), get this menu from the Right condition. This means that you want the inversion model to be applied to all the data in the condition.
- Select "Minimum Norm Imaging", click on Run.
The actual inversion kernel you have just computed (1), contains the same information as the one from the previous section (Computing sources for a single data file). Note that you cannot do anything with this file: if you right-click on it, you can see that there are no Display menu for it.
Two links (2) that allow you apply this inversion kernel to the data files available in this condition (ERF and Std). In their popup menus, there are all the display options introduced in the previous section.
- These links are not saved as files but as specific strings in the database: "link|kernel_file|data_file". This means that to represent them, one should load the shared kernel, load the recordings, and multiply them.
The sources for the Std file do not have any meaning, do not even try to open it. It was just to illustrate the way a kernel is shared
Double-click on both sources files available for Right / ERF (link and non-link), and verify at many different times that the cortical maps are exactly the same in both cases.
You can estimate the sources for many subjects or conditions at once, as it was explained for the head models in previous tutorial: the Compute sources menu is available on all the subjects and conditions popup menus.
- Delete the shared kernel (1), we don't need redundant and confusing information for the next steps. Observe that both links disappear at the same time.
Minimum norm options
Let's introduce briefly the other options offered for the source estimation. Right-click again on Right / ERF > Compute sources. Click on "Expert mode", you see more options appearing in the window. If you click on Run, you have access to all the options of the wMNE algorithm.
wMNE: Whitened and depth-weigthed linear L2-minimum norm estimates algorithm inspired from Matti Hamalainen's MNE software. Introduced in Brainstorm 3.1 by Rey Ramirez. For a full description of this method, please refer to the MNE manual, section 6, "The current estimates". Download MNE manual here.
dSPM: Noise-normalized estimate (dynamical Statistical Parametric Mapping [Dale, 2000]). Its computation is based on the wMNE solution.
Basically, the dSPM value at each location is equal to the wMNE value divided by the projection of the estimated noise covariance matrix onto each source point. After whitening, the operational noise covariance matrix is by definition the identity matrix, and hence the projection of the noise is equal to the L2 norm of the row vector of the wMNE inverse operator (in the case of fixed dipole orientations). So, dSPM is what you get when the rows of the wMNE inverse operator all have unit norm (i.e., they all point in different directions but lie in a unit hyper-sphere).
sLORETA: Noise-normalized estimate using the sLORETA approach (standardized LOw Resolution brain Electromagnetic TomogrAphy [Pasqual-Marqui, 2002]). sLORETA solutions have in general a smaller location bias than either the expected current (wMNE) or the dSPM.
Full results: Saves in one big matrix the values of all the sources (15000) for all the time samples (375). The size in memory of such a matrix is about 45Mb for 300ms of recordings. This is still reasonable, so you may use this option in this case. But if you need to process longer recordings, you may have some "Out of memory" errors in Matlab, or fill your hard drive quickly.
Kernel only: Saves only the Inversion kernel, a matrix that describes how to compute the sources when you know the values at the sensor level. So its size is: number of sources (15000) x number of sensors (151). This is possible because these minimum norm methods are linear methods.
- To get the sources time series, you just need to multiply this kernel by the MEG recordings.
Full results = Inversion kernel * Recordings
- The size of this matrix is about 18Mb. In this case, the difference is not very important because we only process 375 time samples. But this inversion kernel is independent from the recordings length, so you can easily scale its computation to much longer recordings.
Probably "Kernel only", as it is faster and produces smaller files.
- All the following operations in Brainstorm will be exactly the same either way. Each time you access the sources values, the program has to do the multiplication Kernel * Recordings, but this is done in a totally transparent way.
The only reason that would make you choose the "Full results" options would be any interest in having the full matrix in one file, in case you want to process the sources values by yourself (filtering, statistics, display...).
Constrained: We consider that at each vertex of the cortex surface, there is only one dipole, and that its orientation is the normal to the cortex surface at this point.
- The size of the inverse operator is [nVertices x nChannel].
- 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. But it's hard to know if we can really rely on it at this level, for this algorithm: Is it the case everywhere in the cortex? Are we supposed to use the inner (grey matter/white matter) or the outer (grey/CSF) surface of the cortex? Can we really be that precise in terms of MRI/MEG registration? These are questions that do not have final answers yet.
- A technical advantage of this method, it produces one value per vertex instead of three. As a consequence 1) the output size is smaller, 2) it's faster to compute and display, and 3) the results are much more intuitive to display because we don't have to think about how to combine three values in one on a cortical map in a 3D figure
Unconstrained: At each vertex of the cortex surface, we define a base of three dipoles with orthogonal directions, and then we estimate the sources for the three orientations independently.The size of the inverse operator is [3*nVertices x nChannel].
Loose: A version of the "unconstrained" method that integrates a weak orientation constraint. It generates an inverse operator that is the same size as the unconstrained one, but that emphasizes the importance of the sources that have an orientation that is close to the normal to the cortex. The value associated with this option set how "loose" should be the orientation constrain (recommended values in MNE are between 0.1 and 0.6, --loose option). This is the default in MNE software.
Truncated:An SVD of the gain matrix for each source point is used to remove the dipole component with least variance, which for the single sphere head model, corresponds to the radial silent component. Only useful when combining a spherical head model and sLORETA.
Default: At the present time, the default in Brainstorm is the "constrained" option, but we will probably switch to the unconstrained model soon. Because 1) we are not exactly sure that the orientation constraint is 100% correct, 2) the unconstrained sources look smoother and nicer, 3) the computation and storage capacities of the average computer increased a lot since the 1990s, so we can now afford to multiply by three the size of all the data. But right now, there are still some issues to fix in the processing pipeline of the unconstrained sources.
Signal to noise ratio (SNR): An estimate of the amplitude SNR of the recordings, as defined in MNE (--snr option in MNE), used to compute the regularization parameter (lambda2 = 1/SNR2). The default value is SNR = 3. Automatic selection of the regularization parameter is currently not supported.
PCA Whitening: Parameter introduced by Rey Ramirez. For more information, see the code of bst_wmne function.
Noise covariance matrix
Full noise covariance: Use the full noise covariance matrix available in Brainstorm database. If the noise covariance file previously computed in is a diagonal matrix (as it is the case in this tutorial), this value is ignored, and the "diagonal noise covariance" option is used instead.
Diagonal noise covariance: Discard the off-diagonal elements of the noise covariance matrix (assuming heteroscedastic uncorrelated noise). Corresponds in MNE to the --diagnoise option.
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).
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).
Do it yourself
Typically, the only options you want to change is the type of estimator that is computed (wMNE, dSPM, sLORETA) and the orientation of the dipoles (constrained, unconstrained). Don't try to play with the other options if you are not really comfortable with what they represent. Now, compute and compare the following combinations of options for Right / ERF:
wMNE, Kernel only, Constrained (already computed)
Check the contents of the file (right-click > File > View file contents): The inverse operator is saved in the field ImagingKernel [nVertices x nSensors]
wMNE, Full results, Constrained:
Check the contents of the file: The full results matrix is saved in the field ImageGridAmp [nVertices x nTime]
- Open solutions #1 and #2 at the same time and check visually at different time points that the results are exactly the same
wMNE, Kernel only, Unconstrained:
- Check the contents of the file: The inverse operator contains now the information for three dipoles (three orientations) per vertex, [3*nVertices, nSensors]
- Open solutions #1 and #3 at the same time, and observe that the unconstrained solution is much smoother
dSPM, Kernel only, Unconstrained
sLORETA, Kernel only, Unconstrained
- Open solutions #3 (wMNE), #4 (dSPM) and #5 (sLORETA), all unconstrained.
- Notice that the units are different: wMNE values are in pAm, dSPM and sLORETA are in arbitrary units (never try to compare these values to anything but the exact same type of inverse solution)
- Observe around 46ms the respective behavior of these three solutions:
- 3) wMNE tends to highlight the top of the gyri and the superficial sources,
- 4) dSPM tends to correct that behavior and may give higher values in deeper areas,
- 5) sLORETA produces a very smooth solution 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.
wMNE: constrained (kernel and full), and unconstrained
Now delete all these files when you're done, and keep only the initial solution: wMNE, Constrained.
Additional discussions on the forum
Imaging resolution kernels: http://neuroimage.usc.edu/forums/showthread.php?1298
Spatial smoothing of sources: http://neuroimage.usc.edu/forums/showthread.php?1409
Units for dSPM and sLORETA: http://neuroimage.usc.edu/forums/showthread.php?1535
Sign of the MNE values: http://neuroimage.usc.edu/forums/showthread.php?1649#post7014
Combining magneto+gradiometers: http://neuroimage.usc.edu/forums/showthread.php?1900
Residual ocular artifacts: http://neuroimage.usc.edu/forums/showthread.php?1272
From continuous tutorials
Right-click on (Common files), on the head model or on the subject node, and select "Compute sources". A shared inversion kernel is created in (Common files); a link node is now visible for each recordings file, single trials and averages. For more information about what these links mean and the operations performed to display them, please refer to the tutorial #8 "Source estimation".
Explore the sources
Right-click on the sources for the left average > Cortical activations > Display on cortex, or simply double click on the file. Go to the main response peak at t = 30ms, and increase the amplitude threshold to 100%. You see a strong activity around the right primary somatosensory cortex, but there are still lots of brain areas that are shown in plain red (value >= 100% maximum of the colorbar ~= 280 pA.m).
The colorbar maximum is set to an estimation of the maximum source amplitude over the time. This estimation is done by finding the time instant with the highest global field power on the sensors (green trace GFP), estimating the sources for this time only, and then taking the maximum source value at this time point. It is a very fast estimate, but not very reliable; we use it because calculating the full source matrix (all the time points x all the sources) just for finding the maximum value would be too long.
In this case, the real maximum is probably higher than what is used by default. To redefine the colorbar maximum: right-click on the 3D figure > Colormap: sources > Set colorbar max value. Set the maximum to 480 pA.m, or any other value that would lead to see just one very focal point on the brain at 30ms.
Go back to the first small peak at t=16ms, and lower the threshold to 10%. Then do what you did with at the sensor level: follow the information processing in the brain until 100ms, millisecond by millisecond, adapting the threshold and the camera position when needed:
16 ms (top-left): First response, primary somatosensory cortex (S1 right)
30 ms (top-right): S1 right
60 ms (bottom-left): Secondary somatosensory cortex (S2 right)
70 ms (bottom-right): Activity ipsilateral to the stimulus (S2 left + S2 right)
Select the two imported folders at once, right-click > Compute sources
Select dSPM and keep all the default options.
- Then you are asked to confirm the list of bad channels to use in the source estimation for each run. Just leave the defaults, which are the channels that we set as bad earlier.
- One inverse operator is created in each condition, with one link per data file.
For more information: Source estimation tutorial.
Average in source space
- Now we have the source maps available for all the trials, we average them in source space.
Select the folders for Run01 and Run02 and the [Process sources] button on the left.
Run process "Average > Average files": Select "By trial group (subject average)"
- Double-click on the source averages to display them (standard=top, deviant=bottom).
Note that opening the source maps can be very long because of the online filters. Check in the Filter tab, you probably still have a 100Hz low-pass filter applied for the visualization. In the case of averaged source maps, the 15000 source signals are filtered on the fly when you load a source file. This can take a significant amount of time. You may consider unchecking this option if the display is too slow on your computer.
Standard: (Right-click on the 3D figures > Snapshot > Time contact sheet)
Rey on sLORETA
Yes in sLORETA the noise covariance is not used at all for the standardization process. It can be used modeling correlated noise and whitening, but that is optional.
I have noticed that a lot of folks are confused about this and I have seen many statements in papers spreading this awful confusion. The sLORETA is standardized by the resolution matrix (diagonal for dipole orientations constraints, or block diagonals for free orientations).
That is why sLORETA has zero localization bias for ALL point-spread functions, and why I always prefer sLORETA over dSPM, MNE, or any beamformer. This is all in the math .... but ..... just so that you know Fas Hsuan Lin's paper comparing sLORETA with dSPM, and MNE has a big mistake, the assumed source covariance matrix is not the identity matrix, and that violates the beauty of the math and results in non-zero localization bias. That's why in Brainstorm the prior source covariance matrix used for sLORETA automatically uses no depth bias compensation (identity matrix). sLORETA accomplishes depth bias compensation via the resolution matrix, NOT via the prior source covariance matrix. Trying to use a depth exponent of 0.7 or 0.8 like we do for MNE and dSPM will mess up sLORETA.You will not find this in a paper, but I checked it all out many years ago. This is critical.
Explain sLORETA units (see email exchanges from Feb 2015)
Issues with dSPM average
Average(dSPM) is NOT equal to dSPM(Average).
There is no problem for the MNE and sLORETA solutions, because the scaling of the noise covariance matrix doesn't impact the results. wMNE(Data, NoiseCov) = wMNE(Data, NoiseCov / N) So when we average we get: Average(wMNE(Trials, NoiseCov)) = wMNE(Average, NoiseCov) = wMNE(Average, NoiseCov / N) But for dSPM we have: dSPM(Data, NoiseCov) = dSPM(Data, NoiseCov / N) ./ sqrt(N) So when we average we get: Average(dSPM(Trials, NoiseCov)) = dSPM(Average, NoiseCov) = dSPM(Average, NoiseCov / N) ./ sqrt(N)
Rey: "dSPM is really a source mapping of SNR, not of activity. Hence, it's not all so surprising that the single trial SNR maps are smaller...
"Rey: "Perhaps, dSPM should be used only for averaged data (i.e., ERF, ERP), at least until it's all figure out. In a way, dSPM is just MNE followed by the noise normalization. Thus, you could do all the single trial processing with the MNE algorithm, and only do the noise normalization when needed (e.g., after averaging or on single trials only if they are not going to be averaged)."
On the hard drive
Document file tags
Document file structure