25699
Comment:

69242

Deletions are marked like this.  Additions are marked like this. 
Line 1:  Line 1: 
= Tutorial 15: Source estimation = ''Authors: Francois Tadel, John C Mosher, Richard Leahy, Sylvain Baillet'' Now 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 you have the recordings. Many solutions are described in the literature, some of them are implemented in Brainstorm, and only one is presented in this tutorial: the minimumnorm estimation. It is not really the most advanced solution, but it is one of the most robust. 
= Tutorial 22: Source estimation = ''Authors: Francois Tadel, Elizabeth Bock, Rey R Ramirez, John C Mosher, Richard M Leahy, Sylvain Baillet'' You have in your database a forward model that explains how the cortical sources determine the values on the sensors. This is useful for simulations, but what we need next is to solve the inverse problem: how to estimate the sources when we have the recordings. This tutorial introduces the tools available in Brainstorm for solving this inverse problem. (For backward compatibility, see the [[Tutorials/TutSourceEstimationold tutorials]]) 
Line 8:  Line 8: 
== Computing sources for a single data file == 1. Rightclick on ''Subject01 / Right / ERF'' > ''Compute sources''.<<BR>><<BR>> {{attachment:popupComputeSources.gif}}  {{attachment:panelComputeSources.gif}} 1. 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:<<BR>> * '''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. 1. A new file is available in the database explorer.<<BR>><<BR>> {{attachment:treeMinNorm.gif}} * 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 (rightclick > 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)<<BR>><<BR>> {{attachment:resultsMat.gif}} 1. Do the same for the ''Left / ERF'' file == Sources visualization == There are two main ways to display the sources: on the cortex surface and on the MRI slices. === Sources on cortex surface === 1. Doubleclick on recordings ''Right / ERF'', to display the time series (always nice to have a time reference). 1. Doubleclick on sources ''Right / ERF / MN: MEG''. <<BR>>Equivalent to rightclick > Cortical activations > Display on cortex. 1. Go to the main peak around 46ms (by clicking on the times series figure)<<BR>><<BR>> {{attachment:sources1.gif}} 1. 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)... 1. 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. 1. Take a few minutes to understand what this threshold value represents.<<BR>> * 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. * Move the slider and look for a threshold value that would give you a really focal source.The following figures represent the sources activations at t=46ms respectively with threshold at 0% and 90%.<<BR>><<BR>> {{attachment:threshold0.gif}} {{attachment:threshold90.gif}} * 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) === 1. Close all the figures (''Close all'' button). Open the time series view for Right / ERF. 1. Rightclick on Right / ERF / MN: MEG > Cortical activations > Display on MRI (3D). 1. This view was also introduced in the tutorial about MRI and surface visualization. Try to rotate, zoom, move the slices, move in time, change the threshold.<<BR>><<BR>> {{attachment:sources3D.gif}} {{attachment:popupFigMri.gif}} 1. 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) * '''Smooth level''': The sources values are smoothed after being reinterpolated in the volume. Those menus define the size of the smoothing kernel (fig2: smooth=2; fig3: smooth=5).<<BR>><<BR>> {{attachment:mriMipAnat.gifheight="143px",width="176px"}} {{attachment:mriMipFunc.gifheight="142px",width="175px"}} {{attachment:mriSmooth.gifheight="141px",width="173px"}} 1. 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: Rightclick on the figure > Snapshot > Volume contact sheet: axial: <<BR>><<BR>> {{attachment:popupSnapshot.gif}} {{attachment:contactAxial.gifheight="288px",width="322px"}} === Sources on MRI (MRI Viewer) === 1. Rightclick on Right / ERF / MN: MEG > Cortical activations > Display on MRI (MRI Viewer). 1. This view was also introduced in the tutorial about MRI and surface visualization. Try to move the slices (sliders, mouse wheel, click on the views), move in time, change the threshold.<<BR>><<BR>> {{attachment:sourcesMriViwer.gifheight="331px",width="359px"}} === 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 those 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 (doubleclick on the source file), and uncheck the Absolute option for the colormap "Sources" (rightclick 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. {{attachment:relValues.gif}} === 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: [[http://neuroimage.usc.edu/forums/showthread.php?1246DoubtaboutcurrentdensityunitspA.morpAm2http://neuroimage.usc.edu/...DoubtaboutcurrentdensityunitspA.morpAm2]] == 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 ? 1. 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. 1. Select "Minimum Norm Imaging", click on Run. 1. Three new nodes are available in the tree:<<BR>><<BR>> {{attachment:popupComputeMulti.gif}} {{attachment:treeMinNormMulti.gif}} 1. 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 rightclick 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. * Those links are not saved as files but as specific strings in the database: "linkkernel_filedata_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 1. Doubleclick on both sources files available for ''Right / ERF'' (link and nonlink), and verify at many different times that the cortical maps are exactly the same in both cases. 1. 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. 1. 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. Rightclick again on Right / ERF > Compute sources. Click on "Expert mode", you more options appearing in the window.If you click on Run, you have access to all the options of the wMNE algorithm. {{attachment:panelExpert.gif}} {{attachment:panelOptions.gif}} 
== Illposed problem == Our goal is to estimate the activity of the thousands of dipoles described by our forward model. However we only have a few hundred spatial measurements as input (the number of sensors). This inverse problem is illposed, 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 wellstudied 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." . {{http://neuroimage.usc.edu/brainstorm/Tutorials/HeadModel?action=AttachFile&do=get&target=forward_inverse.gifforward_inverse.gifclass="attachment"}} 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: '''minimumnorm 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. 
Line 96:  Line 25: 
* '''wMNE''': Whitened and depthweigthed linear L2minimum 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". [[http://www.nmr.mgh.harvard.edu/meg/manuals/MNEmanual2.7.pdfDownload MNE manual here]]. * '''dSPM''': Noisenormalized estimate (dynamical Statistical Parametric Mapping [Dale, 2000]). Its computation is based on the wMNE solution. <<BR>>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 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 hypersphere). * '''sLORETA''': Noisenormalized estimate using the sLORETA approach (standardized LOw Resolution brain Electromagnetic TomogrAphy [PasqualMarqui, 2002]). sLORETA solutions have in general a smaller location bias than either the expected current (wMNE) or the dSPM. 
<<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:minnorm_options.gif}} <<HTML(</A></div>)>> '''Minimum norm imaging''' * Estimates the sources as the solution to a linear imaging problem, that can be interpreted in various ways ([[https://en.wikipedia.org/wiki/Tikhonov_regularizationTikhonov regularization]], [[https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimationMAP 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 prestim recordings in event related studies) or simply assume a whitenoise identify matrix covariance as [[http://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Noise_covariance_regularizationdescribed 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. '''LCMV beamformer''' * 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. [[http://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Advanced_options:_LCMV_beamformerSee 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 ?'''). [[http://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Advanced_options:_Dipole_modelingMore details]]. __'''Recommended option'''__ * 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. === Measure [TODO] === <<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:minnorm_options_measure.gif}} <<HTML(</A></div>)>> 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 [[http://neuroimage.usc.edu/forums/showthread.php?1246DoubtaboutcurrentdensityunitspA.morpAm2on the user forum]], the units are strictly kept in Am, 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 "depthweighted" linear L2minimum 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/mnecmanual/MNEmanual2.7.3.pdfMNE manual]], section 6, "The current estimates". Units: picoampermeter (pAm). * '''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 zscore statistical map. Units: unitless "z". * '''sLORETA''': Standardized LOw Resolution brain Electromagnetic TomogrAphy (PasqualMarqui, 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 (PasqualMarqui 2002), this theoretical data covariance simplifies sLORETA to an alternative form that results in a "resolution" kernel (eq.(17) of (PasqualMarqui 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 [[http://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Source_map_normalizationSource 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: <<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:minnorm_options_orient.gif}} <<HTML(</A></div>)>> * '''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 macrocolumns that are perpendicular to the cortex surface.<<BR>>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.<<BR>>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. <<BR>>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 [[http://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Why_does_it_look_so_noisy.3FWhy does it looks so noisy]]. === Sensors === <<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:minnorm_options_sensors.gif}} <<HTML(</A></div>)>> 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, crossmodality 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 crosscovariances between different sensor types, since regularization and stability of crossmodalities 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, rightclick on the average response for the '''deviant''' stim > '''Compute sources [2018]'''.<<BR>>Select the options: '''Minimum norm''' imaging, '''Current density''' map, '''Constrained''': Normal to cortex. <<BR>><<BR>> {{attachment:minnorm_single.gifwidth="529",height="433"}} * 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 crossvalidation 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". <<BR>><<BR>> {{attachment:minnorm_single_tree.gif}} == Display: Cortex surface == * Rightclick on the sources for the deviant average > Cortical activations > '''Display on cortex'''.<<BR>><<BR>> {{attachment:minnorm_single_popup.gifwidth="380",height="167"}} * Doubleclick on the '''recordings '''for the deviant average to have a time reference. <<BR>>In the filter tab, add a '''lowpass filter at 40Hz'''.<<BR>><<BR>> {{attachment:display_cortex.gifwidth="482",height="163"}} * 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: * '''Amplitude''': Only the sources that have a value superior to 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 source activity on the cortex surface. * Take a few minutes to understand what the '''amplitude threshold''' represents. * The colorbar maximum depends on the way you configured your ''Sources ''colormap. If the option "Maximum: Global" is selected, the maximum should be around 150 pA.m. This value is a rough estimate of the maximum amplitude, and this default value is not always adapted to your figure. To edit the maximum value, use the colormap option "Maximum: Custom". * 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. <<BR>>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. <<BR>><<BR>> {{attachment:display_sliders.gifwidth="509",height="215"}} == 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 orientationconstrained 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 (~510mm 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, orientationconstrained 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 == * Rightclick on the sources for the deviant average > Cortical activations > '''Display on MRI (MRI Viewer)'''. * The MRI viewer was introduced in tutorials [[Tutorials/ImportAnatomy#2]] and [[Tutorials/ExploreAnatomy#3]]. <<BR>>Additionally you can change the current time and amplitude threshold from the Brainstorm window. * This figure shows the sources computed on the cortical surface and reinterpolated 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. <<BR>><<BR>> {{attachment:display_mriviewer.gifwidth="330",height="356"}} * You can configure this figure with the following options: * '''MIP Anatomy''': Checkbox in the MRI Viewer figure. For each slice, display the maximum value over all the slices instead of the original value in the structural MRI ("glass brain" view). * '''MIP Functional''': Same as for MIP Anatomy, but with the layer of functional values. * '''Smooth level''': The sources values can be smoothed after being reinterpolated in the volume. Rightclick 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 topright of the Brainstorm window (or use the time series figure). * {{attachment:display_smooth.gifwidth="363",height="356"}} == Display: MRI 3D == * Rightclick 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.<<BR>>Rightclick and move your mouse to move the slices (or use the Resect panel of the Surface tab). <<BR>><<BR>> {{attachment:display_mri3d.gifwidth="405",height="203"}} == 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. Doubleclick on the colorbar after testing this to reset the colormap. . {{attachment:display_negative.gifwidth="452",height="173"}} 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). . {{attachment:minnorm_sketch.gifwidth="467",height="155"}} 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 timefrequency decompositions and connectivity analysis. For estimating frequency measures on the source maps it is essential that we retain the sign of the time course at each location so that the correct oscillatory frequencies are identified. == Unconstrained orientations == In cases where the orientation constraint imposed on the dipole orientations produces implausible results, it is possible to relax it partially (option "loose constraints") or completely (option "unconstrained"). This produces a vector (3 component) current source at each location which can complicate interpretation, but avoids some of the noisy and discontinuous features in the current map that are often seen in the constrained maps. Unconstrained solutions are particularly appropriate when using the MNI template instead of the subject's anatomy, or when studying deeper or noncortical brain regions for which the normal to the cortical surface obtained with FreeSurfer or BrainSuite is unlikely to match any physiological reality. In terms of data representation, the option "unconstrained" and "loose constraints" are very similar. Instead of using one dipole at each cortical location, a base of three orthogonal dipoles is used. Here we will only illustrate the fully unconstrained case. * In Run#01, rightclick on the average response for the '''deviant''' stim > '''Compute sources [2018]'''.<<BR>>Select the options: '''Minimum norm''' imaging, '''Current density''' map, '''Unconstrained'''. * Doubleclick on the new source file for the deviant average, open the time series simultaneously. The two brain maps below represent the same file at 91ms, with different colormap options (absolute values on the left, relative values on the right). Explanations below. <<BR>><<BR>> {{attachment:minnorm_unconstr_all.gifwidth="652",height="413"}} * We have to be careful with the visual comparisons of constrained and unconstrained source maps displayed on the cortex surface, because they are very different types of data. In unconstrained source maps, we have '''three dipoles with orthogonal orientations at each cortex location''', therefore we cannot represent all the information at once. To display them as an activity map, Brainstorm computes the '''norm of the vectorial sum of the three orientations at each vertex'''. <<BR>>S = sqrt(Sx^2^ + Sy^2^ + Sz^2^) <<BR>><<BR>> {{attachment:minnorm_unconstr_sketch.gifwidth="476",height="158"}} * This explains that we only observe '''positive values''' (no blue values when the colormap is set to display positive and negative values): the norm displayed at each vertex is always positive. The underlying values along each orientation (x,y,z) can be positive or negative and oscillate around zero in time, but we cannot get access to this information with these static cortical maps. * 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 in the forward model between two adjacent points of the grid for a vector dipole represented in Cartesian coordinates while the normal to the surface for two nearby points can be very different, resulting in rapidly changing forward models for the constrained case. * '''Delete''' the unconstrained file, we will not explore this option in the introduction tutorials. You can refer to the tutorial [[http://neuroimage.usc.edu/brainstorm/Tutorials/EpilepsyEEG and epilepsy]] for an example of analysis using unconstrained sources. == Source map normalization == The current density values returned by the minimum norm method have a few problems: * They depend a lot on the SNR of the signal, which may vary significantly between subjects. 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, prestimulus baseline or resting state recordings) can help with all these issues at the same time. In the case of dSPM and sLORETA, the normalizations are computed as part of the inverse routine and based on noise and data covariances, respectively. While dSPM does produce a Zscore map, we also provide an explicit Zscore normalization that offers the user more flexibility in defining a baseline period over which Brainstorm computes the standard deviation for normalization. The normalization options do not change the temporal dynamics of your results when considering a single location but they do alter the relative scaling of each point in the min norm map. If you look at the time series associated with one given source, it will be exactly the same for all normalizations, except for a scaling factor. Only the relative weights change between the sources, and these weights do not change over time. ==== dSPM, sLORETA ==== * In Run#01, rightclick on the average recordings for the '''deviant''' stim > '''Compute sources [2018]'''.<<BR>>Select successively the two normalization options: dSPM, sLORETA, ('''constrained''').<<BR>><<BR>> {{attachment:minnorm_normfiles.gif}} * Doubleclick on all of them to compare them (screen capture at '''143ms'''): <<BR>><<BR>> {{attachment:minnorm_normalized.gifwidth="628",height="156"}} * '''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 Zscores, 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 [[https://neuroimage.usc.edu/brainstorm/Tutorials/SourceEstimation#Averaging_normalized_valuesAveraging normalized values]]. * '''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. The maps are unitless, but unlike dSPM cannot be interpreted as Zscores so are more difficult to interpret. ==== Zscore ==== * The '''Ztransformation''' converts the current density values to a score that represents the number of standard deviations with respect to a baseline period. We define a baseline period in our file (in this case, the prestimulus 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  <<HTML(μ)>>) / <<HTML(σ)>>''' * 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 location, so the sources with a low variability during baseline will be more salient in the cortical maps poststimulus. * In Process1: Select the constrained current density maps (file ''MN: MEG(Constr)''). * Run process "'''Standardize > Baseline normalization'''", '''[100,1.7]ms''', '''Zscore transformation''' <<BR>>Do not select "Use absolute values": We want the sign of the current values.<<BR>><<BR>> {{attachment:zscore_process.gifwidth="558",height="530"}} * Doubleclick on the new normalized file to display it on the cortex (file with the " zscore" tag). <<BR>><<BR>> {{attachment:zscore_cortex.gifwidth="657",height="153"}} * You can see that the cortical maps obtained in this way are '''very similar''' to the other normalization approaches, especially with the dSPM maps. * A value of 3 in this figure means: at this vertex, the value is 3 times higher than the standard deviation from zero during the baseline. If the values during the baseline follow a normal distribution '''N(<<HTML(μ)>>,<<HTML(σ<SUP>2</SUP>)>>)''', then the values we computed follow a N(0,1)='''Z distribution'''. We can get a level of significance from this well known distribution, for instance a value Z=1.96 corresponds to a pvalue of 0.05. These questions will be discussed in more details in the statistics tutorial. * The Znormalized source maps are '''not impacted by the''' '''visualization filters'''. If you open simultaneously the time series and all the files you have now (MN, dSPM, sLORETA, Zscore) and modify the options in the Filter tab, all the figures are updated except for the Zscore one. We can filter easily all the linear models (MN, dSPM, sLORETA), but we would lose the interesting properties of the Zvalues if we were filtering them (the values would not follow a Zdistribution anymore). * If the baseline and the active state are not in the same file, you can use the '''Process2 tab''': place the baseline in the left list (Files A) and the file to normalize in the right list (Files B). ==== Typical recommendations ==== * Use nonnormalized '''current density maps''' for: * Computing shared kernels applied to single trials. * Averaging files across MEG runs. * Computing timefrequency decompositions or connectivity measures on the single trials. * Use '''normalized maps''' (dSPM, sLORETA, Zscore) 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: * 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 and inferences. 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 easier to manipulate in Brainstorm. sLORETA maps can be smoother but they are difficult to interpret. dSPMs, as zscore maps, are much easier to understand and interpret. * Znormalized current density maps are also easy to interpret. They represent explicitly a "deviation from experimental baseline" as defined by the user. In contrast, dSPM indicates the deviation from the data that was used to define the noise covariance used in computing the min norm map. == Delete your experiments == * Select all the source files you computed until now and '''delete''' them. <<BR>><<BR>> {{attachment:delete_norm.gifwidth="474",height="172"}} == Computing sources for single trials == Because the minimum norm model is linear, we can compute an inverse model independently from the recordings and apply it on the recordings when needed. We will now illustrate how to compute a shared inverse model for all the imported epochs. * Rightclick on the '''head model''' or the '''folder '''for Run#01 > '''Compute sources [2018]'''.<<BR>>Select: '''Minimum norm''' imaging, '''Current density''' map, '''Constrained''': Normal to cortex<<BR>><<BR>> {{attachment:minnorm_shared_popup.gifwidth="500",height="273"}} * Because we did not request to compute an inverse model for a specific block of recordings, it computed a '''shared inverse model'''. If you rightclick 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.<<BR>><<BR>> {{attachment:minnorm_shared_kernel.gif}} * The database explorer now shows one '''source link''' to this inverse model for each block of recordings available in this folder, single trials and averages. These links are not real files saved on the hard drive, but you can use them exactly like the previous source files we calculated for the deviant average. If you load a link, Brainstorm loads the corresponding MEG recordings, loads the inverse kernel and multiplies the two on the fly before displaying it. This optimized approach saves a lot of computation time and lot of space on the hard drive.<<BR>><<BR>> {{attachment:minnorm_links.gifwidth="513",height="176"}} == Averaging in source space == ==== Computing the average ==== * First compute the same source model for the the second acquisition run.<<BR>>In Run#02, rightclick on the '''head model''' or the '''folder '''> '''Compute sources [2018]'''.<<BR>>Select: '''Minimum norm''' imaging, '''Current density''' map, '''Constrained''': Normal to cortex<<BR>><<BR>> {{attachment:minnorm_run02.gifwidth="245",height="242"}} * Now we have the source maps available for all the recordings, we can '''average them in source space''' across runs. This allows us to average MEG recordings that were recorded with different head positions (in this case Run#01 and Run#02 have different channel files so they could potentially have different head positions preventing the direct averaging at the sensor level). * Thanks to the linearity of the minimum norm model: MN(Average(trials)) = Average(MN(trials))<<BR>>The two following approaches are equivalent: 1. Averaging the sources of all the individual trials across runs, 1. Averaging the sources for the sensor averages that we already computed for each run. * We will use the second option: using the sources for the sensorlevel averages. 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). * Drag and drop to the Process1 tab the average recordings for '''Run01 '''and '''Run02''', then press the ['''Process sources'''] button on the left to select the source files instead of the MEG recordings. * Run process "'''Average > Average files'''": <<BR>>Select "'''By trial group (subject average)'''" to average together files with similar names. <<BR>>Select "'''Arithmetic average'''" function. <<BR>>Check "'''Weighted average'''" to account for the different numbers of trials in both runs.<<BR>><<BR>> {{attachment:average_process.gifwidth="526",height="565"}} * The two averages that are produced (one for each condition) are saved in the folder '''Intrasubject'''. This is where all the files computed using information from multiple folders within the same subject are saved. If you prefer to have them somewhere else, you can create new folders and move them there, just like you would do with a regular file explorer. <<BR>><<BR>> {{attachment:average_files.gifwidth="268",height="159"}} * 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. <<BR>>Rightclick on each of them > File > View file contents, and check the '''Leff''' field: <<BR>>78 trials for the deviant condition, 378 trials for the standard condition. * Doubleclick on the source averages to display them (deviant=top, standard=bottom). <<BR>>Open the sensorlevel averages as a time reference. <<BR>>Use the predefined view "Left, Right" for looking at the two sides at the same time (shortcut: "7").<<BR>><<BR>> {{attachment:average_display_mn.gifwidth="530",height="296"}} ==== Visualization filters ==== * 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. <<BR>><<BR>> {{attachment:filter_sources.gifwidth="207",height="163"}} * It was not a problem until now because the source files were saved in the compact form (Kernel*recordings) and the visualization filters were applied on the recordings, then projected to the source space. This fast option is not available anymore with these averages across runs. * The visualization filters will not be available anymore after we apply a Zscore normalization. If we want to display Zscore source maps that are smoothed in time, we will have to apply explicitly the filters on the file, with the Process1 tab. ==== Lowpass filter ==== * Clear the Process1 list, then drag and drop the new averages in it. * Run process "'''Preprocess > Bandpass filter'''": [0,40] Hz<<BR>><<BR>> {{attachment:average_filter.gifwidth="476",height="338"}} * '''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 [[http://neuroimage.usc.edu/brainstorm/Tutorials/Epoching#Epoch_lengthhere]], we should have been importing longer epochs in order to filter and normalize the averages properly. <<BR>><<BR>> {{attachment:average_filter2.gifwidth="378",height="323"}} ==== Zscore normalization ==== * In Process1, select the two filtered averages. * Run process "'''Standardize > Baseline normalization'''", baseline='''[100,1.7]ms, ''''''Zscore.'''<<BR>><<BR>> {{attachment:average_zscore.gifwidth="539",height="292"}} * Four new files are accessible in the database: two filtered and two filtered+normalized. <<BR>><<BR>> {{attachment:average_zscore_files.gifwidth="341",height="205"}} * Doubleclick on the source averages to display them (deviant=top, standard=bottom).<<BR>><<BR>> {{attachment:average_display_zscore.gifwidth="530",height="295"}} * The Zscore source values at 90ms are higher for the standard condition (~25) than for the deviant condition (~15). We observe this because the two conditions have very different signaltonoise ratios. The standard condition has about 5x more trials, therefore the standard deviation over the baseline is a lot lower, leading to higher Zscore. * '''Delete''' the nonnormalized filtered files, we will not use them in the following tutorials. <<BR>><<BR>> {{attachment:average_zscore_files2.gifwidth="294",height="147"}} == Note for beginners == Everything below is advanced documentation, you can skip it for now. <<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?skip=Tutorials/Scouts")>> <<TAG(Advanced)>> == 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 zscore 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. '''dSPM''' * 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 (minnorm, zscored), we now dropped this option. Thus averaging dSPM maps now results in maps with variance less than 1, and is consistent with handling minnorm, zscored 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. <<BR>><<BR>> {{attachment:dspm_scale.gifwidth="547",height="307"}} '''Zscore''' * The same SNR issues arise while averaging Zscores: the average of the Zscores is lower than the Zscore of the average. * When computing averages at the subject level: Always '''avoid averaging Zscore maps'''. <<BR>>Average the current density maps, then normalize. '''sLORETA''' * 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 minnorm 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". <<TAG(Advanced)>> == Display: Contact sheets and movies == A good way to represent what is happening in time is to generate contact sheets or videos. Rightclick on any figure and go to the menu Snapshot to check out all the possible options. For a nicer result, take some time to adjust the size of the figure, the amplitude threshold and the colormap options (hiding the colorbar can be a good option for contact sheets). A time stamp is added to the captured figure. The size of the text font is fixed, so if you want it to be readable in the contact sheet, you should make you figure very small before starting the capture. The screen captures below where produced with the colormap "hot". * '''Contact sheet:''' Rightclick on any figure > Snapshot > '''Time contact sheet: Figure'''<<BR>><<BR>> {{attachment:contact_options.gifwidth="128",height="209"}} {{attachment:contact_image.gifwidth="405",height="249"}} * '''Movies''': Rightclick on any figure > Snapshot > '''Movie (time): All figures'''<<BR>><<BR>> {{attachment:video_options.gifwidth="126",height="263"}} [[http://neuroimage.usc.edu/wikidocs/average_sources_unconstr.avi{{attachment:video_vlc.gifhttp://neuroimage.usc.edu/wikidocs/average_sources_unconstr.aviwidth="409",height="262"}}]] <<TAG(Advanced)>> == Model evaluation == One way to evaluate the accuracy of the source reconstruction if to simulate recordings using the estimated source maps. This is done simply by multiplying the source time series with the forward model:<<BR>>MEG_simulated [Nmeg x Ntime] = Forward_model [Nmeg x Nsources] * MN_sources [Nsources x Ntime]<<BR>>Then you can compare visually the original MEG recordings with the simulated ones. More formally, you can compute an error measure from the residuals (recordings  simulated). To simulate MEG recordings from a minimum norm source model, rightclick on the source file, then select the menu "'''Model evaluation > Simulate recordings'''". . {{attachment:model_popup.gifwidth="597",height="151"}} Open sidebyside the original and simulated MEG recordings for the same condition: . {{attachment:model_results.gifwidth="332",height="218"}} <<TAG(Advanced)>> == Advanced options: Minimum norm == Rightclick on the '''deviant average''' for '''Run#01''' > '''Compute sources [2018]'''. <<BR>>Click on the button ['''Show details'''] to bring up all the advanced minimum norm options. <<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:minnorm_details.gifwidth="390",height="392"}} <<HTML(</A></div>)>> === Depth weighting === Briefly, the use of various depth weightings was far more debated in the 1990s, before the introduction of MNE normalization via dSPM, sLORETA, and other "zscoring" methods, which mostly cancel the effects of depth weighting (put another way, after normalization min norm results tend to look quite similar whether depth weighting is used or not). By modifying the source covariance model at each point in the source grid, deeper sources are "boosted" to increase their signal strength relative to the shallower dipoles; otherwise, the resulting MNE current density maps are too dominated by the shallower sources. If using dSPM or sLORETA, little difference in using depth weighting should be noted. To understand how to set these parameters, please refer to the [[https://mne.tools/mnecmanual/MNEmanual2.7.3.pdfMNE manual]]. (options depth, weightexp and weightlimit). === 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 can cause the eigenspectrum of the covariance matrix to be illconditioned (i.e. a large eigenvalue spread or matrix condition number). In Brainstorm, we provide several means to "stabilize" or "regularize" the noise covariance matrix, so that source estimation calculations are more robust to small errors. * '''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 interdependencies 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 zeroout the crosscovariances. 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_2018.m" for notes and details. * '''Recommended option''': This author (Mosher) votes for the '''median eigenvalue '''as being generally effective. The other options are useful for comparing with other software packages that generally employ similar regularization methods. '''[TODO]''' === Regularization parameter [TODO] === In minimum norm estimates, as mentioned above in the comparisons among methods, the data covariance matrix is essentially synthesized by adding the noise covariance matrix to a modeled signal covariance matrix. The signal covariance matrix is generated by passing the source prior through the forward model. The source prior is in turn prescribed by the source model orientation and the depth weighting. A final regularization parameter, however, determines how much weight the signal model should be given relative to the noise model, i.e. the "signal to noise ratio" (SNR). In Brainstorm, we follow the definition of SNR as first defined in the original MNE software of Hamalainen. The signal covariance matrix is "whitened" by the noise covariance matrix, such that the whitened eigenspectrum has elements in terms of SNR (power). We find the mean of this spectrum, then take the square root to yield the average SNR (amplitude). The default in MNE and in Brainstorm is "3", i.e. the average SNR (power) is 9. * '''Signaltonoise ratio''': Use SNR of 3 as the classic recommendation, as discussed above. * '''RMS source amplitude''': An alternative definition of SNR, but still under test and may be dropped. '''[TODO]''' 
Line 101:  Line 340: 
* '''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 those 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. * '''Default ?''' * 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 chose 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...). === Source orientation === * '''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 macrocolumns 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? Those 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 combining 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 properties === * '''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 (lambda^2 = 1/SNR^2). 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 offdiagonal elements of the noise covariance matrix (assuming heteroscedastic uncorrelated noise). Corresponds in MNE to the diagnoise option. * '''Regularize noise covariance''': Regularize the noisecovariance 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). === Depth weighting === The minimumnorm estimates have a bias towards superficial currents. This tendency can be alleviated by adjusting those parameters. To understand how to set those 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: 1. '''wMNE, Kernel only, Constrained '''(already computed) * Check the contents of the file (rightclick > File > View file contents): The inverse operator is saved in the field ImagingKernel [nVertices x nSensors] 1. '''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 1. '''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 1. '''dSPM, Kernel only, Unconstrained''' 1. '''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 those values to anything but the exact same type of inverse solution) * Observe around 46ms the respective behavior of those 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 {{attachment:compMNEconstr.gifheight="173px",width="193px"}} {{attachment:compMNEconstr.gifheight="171px",width="190px"}} {{attachment:compMNE.gifheight="170px",width="189px"}} dSPM, sLORETA: {{attachment:compSPM.gifheight="173px",width="192px"}} {{attachment:compLORETA.gifheight="172px",width="191px"}} Now delete all those 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?1535DipolestrengthunitsfordSPMandsLORETAhttp://neuroimage.usc.edu/forums/showthread.php?1535]] * EEG reference: http://neuroimage.usc.edu/forums/showthread.php?1525#post6718 * Sign of the MNE values: http://neuroimage.usc.edu/forums/showthread.php?1649#post7014 = From continuous tutorials = == Inverse model == Rightclick 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 those links mean and the operations performed to display them, please refer to the [[http://neuroimage.usc.edu/brainstorm/Tutorials/TutSourceEstimationtutorial #8 "Source estimation"]]. {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=inverseDb.gifinverseDb.gifclass="attachment"}} == Explore the sources == Rightclick 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). {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=inverse100.gifinverse100.gifclass="attachment"}} 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: rightclick 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. {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=inverse480.gifinverse480.gifclass="attachment"}} 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''' (topleft): First response, primary somatosensory cortex (S1 right) * '''30 ms''' (topright): S1 right * '''60 ms''' (bottomleft): Secondary somatosensory cortex (S2 right) * '''70 ms''' (bottomright): Activity ipsilateral to the stimulus (S2 left + S2 right) {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=sources16.gifsources16.gifclass="attachment"}} {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=sources30.gifsources30.gifclass="attachment"}} {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=sources60.gifsources60.gifclass="attachment"}} {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawAvg?action=AttachFile&do=get&target=sources70.gifsources70.gifclass="attachment"}} 
As mentioned above, these methods create a convenient linear imaging kernel that is "tall" in the number of elemental dipoles (one or three per grid point) and "wide" only in the number of sensors. At subsequent visualization time, we efficiently multiply the kernel with the data matrix to compute the min norm images. For some custom purposes, however, a user may find it convenient to premultiply the data matrix and generate the full source estimation matrix. This would only be recommended in small data sets, since the full results can become quite large. * '''Kernel only''': Saves only the linear inverse operator, a model that converts sensor values into source values. The size of this matrix is: number of sources (15000) x number of MEG sensors (274). The multiplication with the recordings is done on the fly by Brainstorm in a transparent way. For long recordings or numerous epochs, this form of compact storage helps saving a lot of disk space and computation time, and it speeds up significantly the display. Always select this option when possible. * '''Full results''': Saves in one big matrix the values of all the sources (15,000) for all the time samples (361). The size in memory of such a matrix is about 45Mb for 600ms 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 face "Out of memory" errors in Matlab, or fill your hard drive quickly. * 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 [[http://neuroimage.usc.edu/brainstorm/Tutorials/NoiseCovariance#Data_covariancedata covariance matrix]] in order to enable the option "LCMV beamformer" in the interface. Note that prewhitening 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, for the time being an identity matrix can be substituted for the noise covariance (to enable the option "LCMV beamformer" in the interface) until noise prewhitening is implemented. <<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:lcmv_options.gifwidth="390",height="350"}} <<HTML(</A></div>)>> '''Measure''' The only option "Pseudo Neural Activity Index" (PNAI), is named after the definition of the Neural Activity Index (NAI). We have modified Van Veen’s definition to rely strictly on the data covariance, without need for a separate noise covariance matrix, but the basic premise is the same as in dSPM, sLORETA, and other normalizations. Viewing the resulting "map," in an identical manner to that with MNE, dSPM, and sLORETA described above, reveals possibly multiple sources as peaks in the map. The PNAI scores analogously to zscoring. '''Dipole orientations''' We recommend you choose "'''unconstrained'''" and let the later [[Tutorials/TutDipScanDipole scanning]] process, which finds the best fitting dipole at each time point, optimize the orientation with respect to the data. '''Data covariance regularization''' Same definitions as in MNE, only applied to the data covariance matrix, rather than the noise covariance matrix. Our recommendation is to use '''median eigenvalue'''. == Advanced options: Dipole modeling == Dipole modeling fits a single dipole at each potential source location to produce a dipole scanning map. This map can be viewed as a indication of how well, and where, the dipole fits at each time point. However, we recommend using the subsequent bestdipole fitting routine ([[Tutorials/TutDipScandipole scanning]]) to determine the final location and orientation of the dipole (one per time point). Please note that '''this function does not fit multiple simultaneous dipoles'''. Although not widely recognized, dipole modeling and beamforming are more alike than they are different – when comparing the inverse operators required to compute the dipole scanning map (dipole modeling) and the beamformer output map (LCMV), we see that they differ only in that the former uses an inverse noise covariance matrix while the latter replaces this with the inverse of the data covariance. <<HTML(<div style="padding: 0px 0px 0px 10px; float: right;">)>> {{attachment:dipoles_options.gifwidth="390",height="314"}} <<HTML(</A></div>)>> '''Measure''' This field is now missing, but the resulting imaging kernel file is directly analogous to the PNAI result from LCMV beamforming. The user can display this scanning measure just as with the LCMV case, where again the normalization and units are a form of zscoring. '''Dipole orientations''' Use "unconstrained source" modeling and let the process "[[Tutorials/TutDipScandipole scanning]]" optimize the orientation of the dipole for every time instance. '''Noise covariance regularization''' Similarly, use "median eigenvalue". <<HTML(<div style="clear: right;"></div>)>> The tutorial "[[http://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElekta#Dipole_source_estimationMEG current phantom (Elekta)]]" demonstrates dipole modeling of 32 individual dipoles under realistic experimental noise conditions. <<TAG(Advanced)>> == 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 ([[https://www.frontiersin.org/articles/10.3389/fnins.2019.00076/full#B4Baillet 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. <<TAG(Advanced)>> == On the hard drive == ==== Constrained shared kernel ==== Rightclick on a shared inverse file in the database explorer > File > '''View file contents'''. . {{attachment:kernel_contents.gifwidth="669",height="342"}} ==== Structure of the source files: results_*.mat ==== Mandatory fields: * '''ImagingKernel''': [Nsources x Nchannels] Linear inverse operator that must be multiplied by the recordings in order to get the full source time series. If defined, ImageGridAmp must be empty. * '''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. * '''nComponents''': Number of dipoles per grid point: 1=Constrained, 3=Unconstrained, 0=Mixed. In the case of mixed head models, the atlas GridAtlas documents region by region how the list of grid points matches the list of dipoles. * '''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". * '''GridLoc''': [Nvertices x 3], (x,y,z) positions of the grid of source points. In the case of a surface head model, it is empty and you read directly the positions from the surface file. * '''GridOrient''': [Nvertices x 3], direction of the normal to the surface for each vertex point (copy of the 'VertNormals' matrix of the cortex surface). Empty in the case of a volume head model or unconstrained sources. * '''GridAtlas''': Atlas "Source model" used with mixed source models. * '''GoodChannel''': [1 x Nchannels] Array of channel indices used to estimate the sources. * '''DataFile''': Relative path to the recordings file for which the sources where computed. If this field is set, the source file appears as a dependent of the DataFile. * '''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"). Optional fields: * '''Options''': Structure that contains all the options of the inverse calculation. This is saved in the file only for bookkeeping. * '''Whitener''': Noise covariance whitener computed in bst_inverse_linear_2018.m * '''DataWhitener''': Data covariance whitener computed in bst_inverse_linear_2018.m * '''SourceDecompVa''': [3 x Nsources] Concatenated right singular vectors from the SVD decomposition of the whitened leadfield for each source (only for unconstrained sources). * '''SourceDecompSa''': [3 x Nvertices] Vector diagonal of the singular values from the SVD decomposition of the whitened leadfield for each source (only for unconstrained sources). * '''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 ''Intrasubject'', rightclick on one of the normalized averages > File > '''View file contents'''. . {{attachment:zscore_contents.gifwidth="426",height="377"}} 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 Zscore process updated the field Function ('mn' => 'zscore') ==== Source links ==== * The links are not real files on the hard drive, if you select the menu "View file contents" for any of them it would display the structure of the corresponding shared kernel. * They are saved in the database as strings with a specific structure: "'''linkkernel_filedata_file'''". This string associates a shared inverse operator with some recordings. The two files have to be available to load the this file. All the functions in Brainstorm are equipped to reconstruct the full source matrix dynamically. ==== Filename tags ==== * '''_KERNEL_''': Indicates that the file contains only an inverse kernel, it needs to be associated with recordings to be opened. ==== Useful functions ==== * '''in_bst_results'''(ResultsFile, LoadFull, FieldsList): Load a source file and optionally reconstruct the full source time series on the fly (ImagingKernel * recordings). * '''in_bst'''(FileName, TimeBounds, LoadFull): Load any Brainstorm data file with the possibility to load only a specific part of the file. * '''bst_process'''(''''LoadInputFile'''', FileName, Target, TimeWindow, OPTIONS): The most highlevel function for reading data files, can compute scout values on the fly. == Additional documentation == ==== Articles ==== * '''Minimum norm''': Baillet S, Mosher JC, Leahy RM<<BR>>[[http://neuroimage.usc.edu/paperspdf/BailletMosherLeahy_IEEESPMAG_Nov2001.pdfElectromagnetic brain mapping]], IEEE SP MAG 2001. * '''dSPM''': Dale AM, Liu AK, Fischl BR, Buckner RL, Belliveau JW, Lewine JD, Halgren E<<BR>>[[http://www.ncbi.nlm.nih.gov/pubmed/10798392Dynamic statistical parametric mapping: combining fMRI and MEG for highresolution imaging of cortical activity]]. Neuron 2000 Apr, 26(1):5567 * '''sLORETA''': PascualMarqui RD<<BR>>[[http://www.ncbi.nlm.nih.gov/pubmed/12575463Standardized lowresolution brain electromagnetic tomography (sLORETA): technical details]], Methods Find Exp Clin Pharmacol 2002, 24 Suppl D:512 ==== Tutorials ==== * Tutorial: [[Tutorials/TutVolSourceVolume source estimation]] * Tutorial: [[Tutorials/DeepAtlasDeep cerebral structures]] * Tutorial: [[Tutorials/TutDipScanComputing and displaying dipoles]] * Tutorial: [[Tutorials/DipoleFittingDipole fitting with FieldTrip]] * Tutorial: [[Tutorials/TutBEstMaximum Entropy on the Mean (MEM)]] ==== Forum discussions ==== * Forum: Minimum norm units (pA.m): [[http://neuroimage.usc.edu/forums/showthread.php?1246DoubtaboutcurrentdensityunitspA.morpAm2http://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?1535DipolestrengthunitsfordSPMandsLORETAhttp://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#post7014http://neuroimage.usc.edu/forums/showthread.php?1649]] * Forum: Combine MEG+EEG: https://neuroimage.usc.edu/forums/t/combiningeegandmegforsourceanalysis/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/simulatescalprecording/2421/3 * Forum: Prewhitening: 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/donttrustthesourcepowerspectrumresults/21265https://neuroimage.usc.edu/forums/t/21265]] * Forum: Subset of sensors: https://neuroimage.usc.edu/forums/t/26496 <<HTML(<! ENDPAGE >)>> 
Tutorial 22: Source estimation
Authors: Francois Tadel, Elizabeth Bock, Rey R Ramirez, John C Mosher, Richard M Leahy, Sylvain Baillet
You have in your database a forward model that explains how the cortical sources determine the values on the sensors. This is useful for simulations, but what we need next is to solve the inverse problem: how to estimate the sources when we have the recordings. This tutorial introduces the tools available in Brainstorm for solving this inverse problem. (For backward compatibility, see the ?old tutorials)
Contents
 Illposed 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
Illposed problem
Our goal is to estimate the activity of the thousands of dipoles described by our forward model. However we only have a few hundred spatial measurements as input (the number of sensors). This inverse problem is illposed, 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 wellstudied 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: minimumnorm 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.
Method
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 prestim recordings in event related studies) or simply assume a whitenoise 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.
LCMV beamformer
 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.
Recommended option
 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.
Measure [TODO]
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 Am, 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 "depthweighted" linear L2minimum 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: picoampermeter (pAm).
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 zscore statistical map. Units: unitless "z".
sLORETA: Standardized LOw Resolution brain Electromagnetic TomogrAphy (PasqualMarqui, 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 (PasqualMarqui 2002), this theoretical data covariance simplifies sLORETA to an alternative form that results in a "resolution" kernel (eq.(17) of (PasqualMarqui 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 macrocolumns 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.
Sensors
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, crossmodality 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 crosscovariances between different sensor types, since regularization and stability of crossmodalities 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, rightclick on the average response for the deviant stim > Compute sources [2018].
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 crossvalidation 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
Rightclick on the sources for the deviant average > Cortical activations > Display on cortex.
Doubleclick on the recordings for the deviant average to have a time reference.
In the filter tab, add a lowpass 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:
Amplitude: Only the sources that have a value superior to 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 source activity on the cortex surface.
Take a few minutes to understand what the amplitude threshold represents.
The colorbar maximum depends on the way you configured your Sources colormap. If the option "Maximum: Global" is selected, the maximum should be around 150 pA.m. This value is a rough estimate of the maximum amplitude, and this default value is not always adapted to your figure. To edit the maximum value, use the colormap option "Maximum: Custom".
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 orientationconstrained 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 (~510mm 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, orientationconstrained 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
Rightclick on the sources for the deviant average > Cortical activations > Display on MRI (MRI Viewer).
The MRI viewer was introduced in tutorials #2 and #3.
Additionally you can change the current time and amplitude threshold from the Brainstorm window.This figure shows the sources computed on the cortical surface and reinterpolated 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 Anatomy: Checkbox in the MRI Viewer figure. For each slice, display the maximum value over all the slices instead of the original value in the structural MRI ("glass brain" view).
MIP Functional: Same as for MIP Anatomy, but with the layer of functional values.
Smooth level: The sources values can be smoothed after being reinterpolated in the volume. Rightclick 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 topright of the Brainstorm window (or use the time series figure).
Display: MRI 3D
Rightclick 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.
Rightclick 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. Doubleclick 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 timefrequency decompositions and connectivity analysis. For estimating frequency measures on the source maps it is essential that we retain the sign of the time course at each location so that the correct oscillatory frequencies are identified.
Unconstrained orientations
In cases where the orientation constraint imposed on the dipole orientations produces implausible results, it is possible to relax it partially (option "loose constraints") or completely (option "unconstrained"). This produces a vector (3 component) current source at each location which can complicate interpretation, but avoids some of the noisy and discontinuous features in the current map that are often seen in the constrained maps. Unconstrained solutions are particularly appropriate when using the MNI template instead of the subject's anatomy, or when studying deeper or noncortical brain regions for which the normal to the cortical surface obtained with FreeSurfer or BrainSuite is unlikely to match any physiological reality.
In terms of data representation, the option "unconstrained" and "loose constraints" are very similar. Instead of using one dipole at each cortical location, a base of three orthogonal dipoles is used. Here we will only illustrate the fully unconstrained case.
In Run#01, rightclick on the average response for the deviant stim > Compute sources [2018].
Select the options: Minimum norm imaging, Current density map, Unconstrained.Doubleclick on the new source file for the deviant average, open the time series simultaneously. The two brain maps below represent the same file at 91ms, with different colormap options (absolute values on the left, relative values on the right). Explanations below.
We have to be careful with the visual comparisons of constrained and unconstrained source maps displayed on the cortex surface, because they are very different types of data. In unconstrained source maps, we have three dipoles with orthogonal orientations at each cortex location, therefore we cannot represent all the information at once. To display them as an activity map, Brainstorm computes the norm of the vectorial sum of the three orientations at each vertex.
S = sqrt(Sx^{2} + Sy^{2} + Sz^{2})
This explains that we only observe positive values (no blue values when the colormap is set to display positive and negative values): the norm displayed at each vertex is always positive. The underlying values along each orientation (x,y,z) can be positive or negative and oscillate around zero in time, but we cannot get access to this information with these static cortical maps.
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 in the forward model between two adjacent points of the grid for a vector dipole represented in Cartesian coordinates while the normal to the surface for two nearby points can be very different, resulting in rapidly changing forward models for the constrained case.
Delete the unconstrained file, we will not explore this option in the introduction tutorials. You can refer to the tutorial EEG and epilepsy for an example of analysis using unconstrained sources.
Source map normalization
The current density values returned by the minimum norm method have a few problems:
 They depend a lot on the SNR of the signal, which may vary significantly between subjects. 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, prestimulus baseline or resting state recordings) can help with all these issues at the same time. In the case of dSPM and sLORETA, the normalizations are computed as part of the inverse routine and based on noise and data covariances, respectively. While dSPM does produce a Zscore map, we also provide an explicit Zscore normalization that offers the user more flexibility in defining a baseline period over which Brainstorm computes the standard deviation for normalization.
The normalization options do not change the temporal dynamics of your results when considering a single location but they do alter the relative scaling of each point in the min norm map. If you look at the time series associated with one given source, it will be exactly the same for all normalizations, except for a scaling factor. Only the relative weights change between the sources, and these weights do not change over time.
dSPM, sLORETA
In Run#01, rightclick on the average recordings for the deviant stim > Compute sources [2018].
Select successively the two normalization options: dSPM, sLORETA, (constrained).
Doubleclick 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 Zscores, 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.
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. The maps are unitless, but unlike dSPM cannot be interpreted as Zscores so are more difficult to interpret.
Zscore
The Ztransformation converts the current density values to a score that represents the number of standard deviations with respect to a baseline period. We define a baseline period in our file (in this case, the prestimulus 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 location, so the sources with a low variability during baseline will be more salient in the cortical maps poststimulus.
In Process1: Select the constrained current density maps (file MN: MEG(Constr)).
Run process "Standardize > Baseline normalization", [100,1.7]ms, Zscore transformation
Do not select "Use absolute values": We want the sign of the current values.
Doubleclick 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.
A value of 3 in this figure means: at this vertex, the value is 3 times higher than the standard deviation from zero during the baseline. If the values during the baseline follow a normal distribution N(μ,σ^{2}), then the values we computed follow a N(0,1)=Z distribution. We can get a level of significance from this well known distribution, for instance a value Z=1.96 corresponds to a pvalue of 0.05. These questions will be discussed in more details in the statistics tutorial.
The Znormalized source maps are not impacted by the visualization filters. If you open simultaneously the time series and all the files you have now (MN, dSPM, sLORETA, Zscore) and modify the options in the Filter tab, all the figures are updated except for the Zscore one. We can filter easily all the linear models (MN, dSPM, sLORETA), but we would lose the interesting properties of the Zvalues if we were filtering them (the values would not follow a Zdistribution anymore).
If the baseline and the active state are not in the same file, you can use the Process2 tab: place the baseline in the left list (Files A) and the file to normalize in the right list (Files B).
Typical recommendations
Use nonnormalized current density maps for:
 Computing shared kernels applied to single trials.
 Averaging files across MEG runs.
 Computing timefrequency decompositions or connectivity measures on the single trials.
Use normalized maps (dSPM, sLORETA, Zscore) 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:
 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 and inferences. 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 easier to manipulate in Brainstorm. sLORETA maps can be smoother but they are difficult to interpret. dSPMs, as zscore maps, are much easier to understand and interpret.
 Znormalized current density maps are also easy to interpret. They represent explicitly a "deviation from experimental baseline" as defined by the user. In contrast, dSPM indicates the deviation from the data that was used to define the noise covariance used in computing the min norm map.
Delete your experiments
Select all the source files you computed until now and delete them.
Computing sources for single trials
Because the minimum norm model is linear, we can compute an inverse model independently from the recordings and apply it on the recordings when needed. We will now illustrate how to compute a shared inverse model for all the imported epochs.
Rightclick on the head model or the folder for Run#01 > Compute sources [2018].
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 rightclick 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.
The database explorer now shows one source link to this inverse model for each block of recordings available in this folder, single trials and averages. These links are not real files saved on the hard drive, but you can use them exactly like the previous source files we calculated for the deviant average. If you load a link, Brainstorm loads the corresponding MEG recordings, loads the inverse kernel and multiplies the two on the fly before displaying it. This optimized approach saves a lot of computation time and lot of space on the hard drive.
Averaging in source space
Computing the average
First compute the same source model for the the second acquisition run.
In Run#02, rightclick on the head model or the folder > Compute sources [2018].
Select: Minimum norm imaging, Current density map, Constrained: Normal to cortex
Now we have the source maps available for all the recordings, we can average them in source space across runs. This allows us to average MEG recordings that were recorded with different head positions (in this case Run#01 and Run#02 have different channel files so they could potentially have different head positions preventing the direct averaging at the sensor level).
Thanks to the linearity of the minimum norm model: MN(Average(trials)) = Average(MN(trials))
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.
 We will use the second option: using the sources for the sensorlevel averages. 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).
Drag and drop to the Process1 tab the average recordings for Run01 and Run02, then press the [Process sources] button on the left to select the source files instead of the MEG recordings.
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 two averages that are produced (one for each condition) are saved in the folder Intrasubject. This is where all the files computed using information from multiple folders within the same subject are saved. If you prefer to have them somewhere else, you can create new folders and move them there, just like you would do with a regular file explorer.
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.
Rightclick 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.Doubleclick on the source averages to display them (deviant=top, standard=bottom).
Open the sensorlevel averages as a time reference.
Use the predefined view "Left, Right" for looking at the two sides at the same time (shortcut: "7").
Visualization filters
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.
 It was not a problem until now because the source files were saved in the compact form (Kernel*recordings) and the visualization filters were applied on the recordings, then projected to the source space. This fast option is not available anymore with these averages across runs.
 The visualization filters will not be available anymore after we apply a Zscore normalization. If we want to display Zscore source maps that are smoothed in time, we will have to apply explicitly the filters on the file, with the Process1 tab.
Lowpass filter
 Clear the Process1 list, then drag and drop the new averages in it.
Run process "Preprocess > Bandpass 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.
Zscore normalization
 In Process1, select the two filtered averages.
Run process "Standardize > Baseline normalization", baseline=[100,1.7]ms, Zscore.
Four new files are accessible in the database: two filtered and two filtered+normalized.
Doubleclick on the source averages to display them (deviant=top, standard=bottom).
 The Zscore source values at 90ms are higher for the standard condition (~25) than for the deviant condition (~15). We observe this because the two conditions have very different signaltonoise ratios. The standard condition has about 5x more trials, therefore the standard deviation over the baseline is a lot lower, leading to higher Zscore.
Delete the nonnormalized 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 zscore 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.
dSPM
 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 (minnorm, zscored), we now dropped this option. Thus averaging dSPM maps now results in maps with variance less than 1, and is consistent with handling minnorm, zscored 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.
Zscore
 The same SNR issues arise while averaging Zscores: the average of the Zscores is lower than the Zscore of the average.
When computing averages at the subject level: Always avoid averaging Zscore maps.
Average the current density maps, then normalize.
sLORETA
 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 minnorm 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
A good way to represent what is happening in time is to generate contact sheets or videos. Rightclick on any figure and go to the menu Snapshot to check out all the possible options. For a nicer result, take some time to adjust the size of the figure, the amplitude threshold and the colormap options (hiding the colorbar can be a good option for contact sheets).
A time stamp is added to the captured figure. The size of the text font is fixed, so if you want it to be readable in the contact sheet, you should make you figure very small before starting the capture. The screen captures below where produced with the colormap "hot".
Contact sheet: Rightclick on any figure > Snapshot > Time contact sheet: Figure
Movies: Rightclick on any figure > Snapshot > Movie (time): All figures
Model evaluation
One way to evaluate the accuracy of the source reconstruction if to simulate recordings using the estimated source maps. This is done simply by multiplying the source time series with the forward model:
MEG_simulated [Nmeg x Ntime] = Forward_model [Nmeg x Nsources] * MN_sources [Nsources x Ntime]
Then you can compare visually the original MEG recordings with the simulated ones. More formally, you can compute an error measure from the residuals (recordings  simulated).
To simulate MEG recordings from a minimum norm source model, rightclick on the source file, then select the menu "Model evaluation > Simulate recordings".
Open sidebyside the original and simulated MEG recordings for the same condition:
Advanced options: Minimum norm
Rightclick on the deviant average for Run#01 > Compute sources [2018].
Click on the button [Show details] to bring up all the advanced minimum norm options.
Depth weighting
Briefly, the use of various depth weightings was far more debated in the 1990s, before the introduction of MNE normalization via dSPM, sLORETA, and other "zscoring" methods, which mostly cancel the effects of depth weighting (put another way, after normalization min norm results tend to look quite similar whether depth weighting is used or not).
By modifying the source covariance model at each point in the source grid, deeper sources are "boosted" to increase their signal strength relative to the shallower dipoles; otherwise, the resulting MNE current density maps are too dominated by the shallower sources. If using dSPM or sLORETA, little difference in using depth weighting should be noted. To understand how to set these parameters, please refer to the MNE manual. (options depth, weightexp and weightlimit).
Noise covariance regularization [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 can cause the eigenspectrum of the covariance matrix to be illconditioned (i.e. a large eigenvalue spread or matrix condition number). In Brainstorm, we provide several means to "stabilize" or "regularize" the noise covariance matrix, so that source estimation calculations are more robust to small errors.
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 interdependencies 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 zeroout the crosscovariances. 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_2018.m" for notes and details.
Recommended option: This author (Mosher) votes for the median eigenvalue as being generally effective. The other options are useful for comparing with other software packages that generally employ similar regularization methods. [TODO]
Regularization parameter [TODO]
In minimum norm estimates, as mentioned above in the comparisons among methods, the data covariance matrix is essentially synthesized by adding the noise covariance matrix to a modeled signal covariance matrix. The signal covariance matrix is generated by passing the source prior through the forward model. The source prior is in turn prescribed by the source model orientation and the depth weighting.
A final regularization parameter, however, determines how much weight the signal model should be given relative to the noise model, i.e. the "signal to noise ratio" (SNR). In Brainstorm, we follow the definition of SNR as first defined in the original MNE software of Hamalainen. The signal covariance matrix is "whitened" by the noise covariance matrix, such that the whitened eigenspectrum has elements in terms of SNR (power). We find the mean of this spectrum, then take the square root to yield the average SNR (amplitude). The default in MNE and in Brainstorm is "3", i.e. the average SNR (power) is 9.
Signaltonoise ratio: Use SNR of 3 as the classic recommendation, as discussed above.
RMS source amplitude: An alternative definition of SNR, but still under test and may be dropped. [TODO]
Output mode
As mentioned above, these methods create a convenient linear imaging kernel that is "tall" in the number of elemental dipoles (one or three per grid point) and "wide" only in the number of sensors. At subsequent visualization time, we efficiently multiply the kernel with the data matrix to compute the min norm images.
For some custom purposes, however, a user may find it convenient to premultiply the data matrix and generate the full source estimation matrix. This would only be recommended in small data sets, since the full results can become quite large.
Kernel only: Saves only the linear inverse operator, a model that converts sensor values into source values. The size of this matrix is: number of sources (15000) x number of MEG sensors (274). The multiplication with the recordings is done on the fly by Brainstorm in a transparent way. For long recordings or numerous epochs, this form of compact storage helps saving a lot of disk space and computation time, and it speeds up significantly the display. Always select this option when possible.
Full results: Saves in one big matrix the values of all the sources (15,000) for all the time samples (361). The size in memory of such a matrix is about 45Mb for 600ms 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 face "Out of memory" errors in Matlab, or fill your hard drive quickly.
 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 prewhitening 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, for the time being an identity matrix can be substituted for the noise covariance (to enable the option "LCMV beamformer" in the interface) until noise prewhitening is implemented.
Measure
The only option "Pseudo Neural Activity Index" (PNAI), is named after the definition of the Neural Activity Index (NAI). We have modified Van Veen’s definition to rely strictly on the data covariance, without need for a separate noise covariance matrix, but the basic premise is the same as in dSPM, sLORETA, and other normalizations. Viewing the resulting "map," in an identical manner to that with MNE, dSPM, and sLORETA described above, reveals possibly multiple sources as peaks in the map. The PNAI scores analogously to zscoring.
Dipole orientations
We recommend you choose "unconstrained" and let the later Dipole scanning process, which finds the best fitting dipole at each time point, optimize the orientation with respect to the data.
Data covariance regularization
Same definitions as in MNE, only applied to the data covariance matrix, rather than the noise covariance matrix. Our recommendation is to use median eigenvalue.
Advanced options: Dipole modeling
Dipole modeling fits a single dipole at each potential source location to produce a dipole scanning map. This map can be viewed as a indication of how well, and where, the dipole fits at each time point. However, we recommend using the subsequent bestdipole fitting routine (dipole scanning) to determine the final location and orientation of the dipole (one per time point). Please note that this function does not fit multiple simultaneous dipoles.
Although not widely recognized, dipole modeling and beamforming are more alike than they are different – when comparing the inverse operators required to compute the dipole scanning map (dipole modeling) and the beamformer output map (LCMV), we see that they differ only in that the former uses an inverse noise covariance matrix while the latter replaces this with the inverse of the data covariance.
Measure
This field is now missing, but the resulting imaging kernel file is directly analogous to the PNAI result from LCMV beamforming. The user can display this scanning measure just as with the LCMV case, where again the normalization and units are a form of zscoring.
Dipole orientations
Use "unconstrained source" modeling and let the process "dipole scanning" optimize the orientation of the dipole for every time instance.
Noise covariance regularization
Similarly, use "median eigenvalue".
The tutorial "MEG current phantom (Elekta)" demonstrates dipole modeling of 32 individual dipoles under realistic experimental noise conditions.
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
Rightclick on a shared inverse file in the database explorer > File > View file contents.
Structure of the source files: results_*.mat
Mandatory fields:
ImagingKernel: [Nsources x Nchannels] Linear inverse operator that must be multiplied by the recordings in order to get the full source time series. If defined, ImageGridAmp must be empty.
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.
nComponents: Number of dipoles per grid point: 1=Constrained, 3=Unconstrained, 0=Mixed. In the case of mixed head models, the atlas GridAtlas documents region by region how the list of grid points matches the list of dipoles.
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".
GridLoc: [Nvertices x 3], (x,y,z) positions of the grid of source points. In the case of a surface head model, it is empty and you read directly the positions from the surface file.
GridOrient: [Nvertices x 3], direction of the normal to the surface for each vertex point (copy of the 'VertNormals' matrix of the cortex surface). Empty in the case of a volume head model or unconstrained sources.
GridAtlas: Atlas "Source model" used with mixed source models.
GoodChannel: [1 x Nchannels] Array of channel indices used to estimate the sources.
DataFile: Relative path to the recordings file for which the sources where computed. If this field is set, the source file appears as a dependent of the DataFile.
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").
Optional fields:
Options: Structure that contains all the options of the inverse calculation. This is saved in the file only for bookkeeping.
Whitener: Noise covariance whitener computed in bst_inverse_linear_2018.m
DataWhitener: Data covariance whitener computed in bst_inverse_linear_2018.m
SourceDecompVa: [3 x Nsources] Concatenated right singular vectors from the SVD decomposition of the whitened leadfield for each source (only for unconstrained sources).
SourceDecompSa: [3 x Nvertices] Vector diagonal of the singular values from the SVD decomposition of the whitened leadfield for each source (only for unconstrained sources).
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 Intrasubject, rightclick 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 Zscore process updated the field Function ('mn' => 'zscore')
Source links
 The links are not real files on the hard drive, if you select the menu "View file contents" for any of them it would display the structure of the corresponding shared kernel.
They are saved in the database as strings with a specific structure: "linkkernel_filedata_file". This string associates a shared inverse operator with some recordings. The two files have to be available to load the this file. All the functions in Brainstorm are equipped to reconstruct the full source matrix dynamically.
Filename tags
_KERNEL_: Indicates that the file contains only an inverse kernel, it needs to be associated with recordings to be opened.
Useful functions
in_bst_results(ResultsFile, LoadFull, FieldsList): Load a source file and optionally reconstruct the full source time series on the fly (ImagingKernel * recordings).
in_bst(FileName, TimeBounds, LoadFull): Load any Brainstorm data file with the possibility to load only a specific part of the file.
bst_process('LoadInputFile', FileName, Target, TimeWindow, OPTIONS): The most highlevel function for reading data files, can compute scout values on the fly.
Additional documentation
Articles
Minimum norm: Baillet S, Mosher JC, Leahy RM
Electromagnetic brain mapping, IEEE SP MAG 2001.dSPM: Dale AM, Liu AK, Fischl BR, Buckner RL, Belliveau JW, Lewine JD, Halgren E
Dynamic statistical parametric mapping: combining fMRI and MEG for highresolution imaging of cortical activity. Neuron 2000 Apr, 26(1):5567sLORETA: PascualMarqui RD
Standardized lowresolution brain electromagnetic tomography (sLORETA): technical details, Methods Find Exp Clin Pharmacol 2002, 24 Suppl D:512
Tutorials
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 discussions
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/combiningeegandmegforsourceanalysis/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/simulatescalprecording/2421/3
Forum: Prewhitening: 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