= Tutorial 22: Source estimation = '''[UNDER CONSTRUCTION]''' ''Authors: Francois Tadel, Elizabeth Bock, Rey R Ramirez, John C Mosher, Richard Leahy, Sylvain Baillet'' You have in your database a forward model matrix that explains how the cortical sources determine the values on the sensors. This is useful for simulations, but what we need is to build the inverse information: how to estimate the sources when we have the recordings. This tutorials introduces the tools available in Brainstorm for solving this inverse problem. <> == Ill-posed problem == Our goal is to estimate the activity of the 45,000 dipoles described by our forward model. However we only have a few hundreds of variables (the number of sensors). This inverse problem is ill-posed, there is an '''infinity of combinations''' of source activity that can generate exactly the same sensor topography. Inverting the forward problem directly is impossible, unless we add some strong priors in our model. Wikipedia says: "Inverse problems are some of the most important and well-studied mathematical problems in science and mathematics because they tell us about parameters that we cannot directly observe. They have wide application in optics, radar, acoustics, communication theory, signal processing, medical imaging, computer vision, geophysics, oceanography, astronomy, remote sensing, natural language processing, machine learning, nondestructive testing, and many other fields." . {{http://neuroimage.usc.edu/brainstorm/Tutorials/HeadModel?action=AttachFile&do=get&target=forward_inverse.gif|forward_inverse.gif|class="attachment"}} Many solutions have been proposed in the literature, based on different assumptions on the way the brain works and depending on the amount of information we already have on the effects we are studying. Among the hundreds of methods available, two classes of inverse models have been widely used in MEG/EEG source imaging in the past years: '''minimum-norm solutions''' and '''beamformers'''. Both approaches have the advantage of being '''linear''': the activity of the sources is a linear recombination of the MEG/EEG recordings. It is possible to solve the inverse problem independently of the recordings, making the data manipulation a lot easier and faster. Both are available in Brainstorm, so you can use the one the most adapted to your recordings or to your own personal expertise. Only the minimum norm estimates will be described in this tutorial, but the other solutions work exactly in the same way. == Source estimation options [TODO] == Before we start estimating the sources for the recordings available in our database, let's start with an overview of the options available. The screen capture below represents the basic options for the minimum norm estimates. The options for the other methods will be described in advanced tutorials. {{attachment:minnorm_options.gif}} === Method === * '''Minimum norm''': __Priors, justification, application case?__<
> Require an estimation of the noise at the level of the sensors (noise covariance matrix). * '''Dipole modeling''': __?__ * '''LCMV beamformer''': __?__<
>Require both a noise covariance matrix and a data covariance matrix (representation of the effect we are trying to localize in the brain, covariance of the latencies of interest). * __'''Recommended option'''__: Provided that we know at what latencies to look at, we can compute a correct data covariance matrix and may obtain a better spatial accuracy with a beamformer. However, in many cases we don't exactly know what we are looking at, the risks of misinterpretation of badly designed beamforming results are high. Brainstorm tends to favor minimum norm solutions, which have the advantage of needing less manual tuning for getting acceptable results. === Measure === The minimum norm estimates gives a measure of the current density flowing at the surface of the cortex. To visualize these results and compare them between subjects, we can normalize the MNE values to get a standardized level of activation with respect to the noise or baseline level (dSPM, sLORETA, MNp). * '''Current density map''': Whitened and depth-weigthed linear L2-minimum norm estimates algorithm inspired from Matti Hamalainen's MNE software. For a full description of this method, please refer to the [[http://www.nmr.mgh.harvard.edu/meg/manuals/MNE-manual-2.7.pdf|MNE manual]], section 6, "The current estimates". * '''dSPM''': Noise-normalized estimate (dynamical Statistical Parametric Mapping [Dale, 2000]). Its computation is based on the MNE solution. * '''sLORETA''': Noise-normalized estimate using the sLORETA approach (standardized LOw Resolution brain Electromagnetic TomogrAphy [Pasqual-Marqui, 2002]). sLORETA solutions have in general a smaller location bias than either the expected current (MNE) or the dSPM. * '''MNp''': ? * __'''Recommended option'''__: Discussed in the section "Source map normalization" === Source orientation === * '''Constrained: Normal to cortex''': Only one dipole at each vertex of the cortex surface, oriented normally to the surface. This is based on the anatomical observation that in the cortex, the neurons are mainly organized in macro-columns that are perpendicular to the cortex surface.<
>Size of the inverse operator: [nVertices x nChannel]. * '''Constrained: Optimal orientation''': Only one dipole at each vertex of the cortex surface, oriented normally to the surface. This is based on the anatomical observation that in the cortex, the neurons are mainly organized in macro-columns that are perpendicular to the cortex surface.<
>Size of the inverse operator: [nVertices x nChannel]. * '''Unconstrained''': At each vertex of the cortex surface, we define a base of three dipoles with orthogonal directions, then we estimate the sources for the three orientations independently. <
>Size of the inverse operator: [3*nVertices x nChannel]. * '''Loose''': A version of the "unconstrained" option with a weak orientation constraint that emphasizes the importance of the sources with orientations that are 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). <
>Size of the inverse operator: [3*nVertices x nChannel]. * __'''Recommended option'''__: ?<
>The constrained options use one dipole per orientation instead of three, therefore the source maps are smaller, faster to compute and display, and much more intuitive to process because we don't have to think about recombining the three values in one. However the normal orientation constraint is most of the time too strong and not realistic. <
>Unconstrained sources look smoother and nicer but are not necessarily more accurate. == Computing sources for a single data file == * In Run#01, right-click on the average response for the '''deviant''' stim > '''Compute sources [2015]'''.<
>Select the options: '''Minimum norm''' imaging, '''Current density''' map, '''Constrained''': Normal to cortex. <
><
> {{attachment:minnorm_single.gif||height="462",width="492"}} * The other menu "Compute sources" brings the interface that was used previously in Brainstorm. We are going to keep maintaining the two implementations in parallel for a while for compatibility and cross-validation purposes. * The result of the computation is displayed as a dependent of the deviant average because it is related only to this file. In the file comment, "MN" stands minimum norm and "Constr" stands for "Constrained: normal orientation". <
><
> {{attachment:minnorm_single_tree.gif}} == Display: Cortex surface == * Right-click on the sources for the deviant average > Cortical activations > '''Display on cortex'''.<
><
> {{attachment:minnorm_single_popup.gif||height="167",width="380"}} * Double-click on the '''recordings '''for the deviant average to have time reference. <
>In the filter tab, add a '''low-pass filter at 100Hz'''.<
><
> {{attachment:display_cortex.gif||height="163",width="482"}} * 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 many display properties from the Surface tab: * '''Amplitude''': Only the sources that have a value superior than a given percentage of the colorbar maximum are displayed. * '''Min. size''': Hide all the small activated regions, ie. the connected color patches that contain a number of vertices smaller than this "min size" value. * '''Transparency''': Change the transparency of the sources on the cortex. * Take a few minutes to understand what 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 30 pA.m. This value is a rough estimate of the maximum amplitude, sometimes you may have to redefine it manually. * On the screen capture below, the threshold value is set to 90%. It means that only the sources that have a value over 0.90*30 = 27 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 somatosensory cortex, which is what we are expecting for an auditory stimulation. <
><
> {{attachment:display_sliders.gif||height="225",width="478"}} == Display: MRI 3D == * Right-click on the source file > Cortical activations > Display on MRI (3D). * 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.<
><
> {{attachment:sources3D.gif}} {{attachment:popupFigMri.gif}} * 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 re-interpolated in the volume. These menus define the size of the smoothing kernel (fig2: smooth=2; fig3: smooth=5).<
><
> {{attachment:mriMipAnat.gif||height="143px",width="176px"}} {{attachment:mriMipFunc.gif||height="142px",width="175px"}} {{attachment:mriSmooth.gif||height="141px",width="173px"}} == Display: MRI Viewer == 1. Right-click 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.<
><
> {{attachment:sourcesMriViwer.gif||height="331px",width="359px"}} == Display: Contact sheets and movies == * '''Standard:''' (Right-click on the 3D figures > Snapshot > Time contact sheet) . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=average_source_standard_left.gif|average_source_standard_left.gif|height="263",width="486",class="attachment"}} * '''Deviant:''' . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=average_source_deviant_left.gif|average_source_deviant_left.gif|height="263",width="486",class="attachment"}} * Contact sheets: in time or in space, for each orientation. You can try all the menus. Example: Right-click on the figure > Snapshot > Volume contact sheet: axial: <
><
> {{attachment:popupSnapshot.gif}} {{attachment:contactAxial.gif||height="288px",width="322px"}} * '''Movies''': Right-click on any figure > Snapshot > '''Movie (time): All figures''' (click to download video) [[http://neuroimage.usc.edu/wikidocs/average_sources.avi|{{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=average_source_video.gif|http://neuroimage.usc.edu/brainstorm/wikidocs/average_sources.avi|height="258",width="484",class="attachment"}}]] == Minimum norm values are not only positive == You should pay attention to a property of the current amplitudes that are given by the wMNE method: they can be positive of negative, and they oscillate around zero. It's not easy to figure out what is the exact meaning of a negative value respect with a positive value, and most of the time we are only interested in knowing what is activated at what time, and therefore we look only at the absolute values of the sources. In some other cases, mainly when doing frequency analysis, we need to pay attention to the sign of these values. Because we cannot do a frequency decomposition of the absolute values of the sources, we need to keep the sign all along our processes. Display again the sources for Right / ERF on the cortex surface (double-click on the source file), and uncheck the Absolute option for the colormap "Sources" (right-click on the figure > Colormap Sources > Absolute values). Decrease the threshold to observe the pattern of alternance between positive and negative values on the surface. Then double click on the colorbar to reset it to its default. {{attachment:relValues.gif}} == Source map normalization == * MNp, dSPM, sLORETA * Z-score : baseline correction Typical recommendations: * Use non-normalized current density maps for shared kernels applied to single trials. Example cases: * Averaging across MEG runs, * Computing time-frequency decompositions or connectivity measures on the single trials. * Use normalized maps (dSPM, sLORETA, MNp, Z-score) for calculating sources for average responses: * Exploring visually the ERP/ERF at the source level, * Normalizing the subjects condition averages before a group analysis. * '''wMNE, Kernel only, Constrained '''(already computed) * Check the contents of the file (right-click > File > View file contents): The inverse operator is saved in the field ImagingKernel [nVertices x nSensors] * '''wMNE, Full results, Constrained''': * Check the contents of the file: The full results matrix is saved in the field ImageGridAmp [nVertices x nTime] * Open solutions #1 and #2 at the same time and check visually at different time points that the results are exactly the same * '''wMNE, Kernel only, Unconstrained''': * Check the contents of the file: The inverse operator contains now the information for three dipoles (three orientations) per vertex, [3*nVertices, nSensors] * Open solutions #1 and #3 at the same time, and observe that the unconstrained solution is much smoother * '''dSPM, Kernel only, Unconstrained''' * '''sLORETA, Kernel only, Unconstrained''' * Open solutions #3 (wMNE), #4 (dSPM) and #5 (sLORETA), all unconstrained. * Notice that the units are different: wMNE values are in pAm, dSPM and sLORETA are in arbitrary units (never try to compare these values to anything but the exact same type of inverse solution) * Observe around 46ms the respective behavior of these three solutions: * 3) wMNE tends to highlight the top of the gyri and the superficial sources, * 4) dSPM tends to correct that behavior and may give higher values in deeper areas, * 5) sLORETA produces a very smooth solution where all the potentially activated area of the brain (given to the low spatial resolution of the source localization with MEG/EEG) is shown as connected, regardless of the depth of the sources. wMNE: constrained (kernel and full), and unconstrained {{attachment:compMNEconstr.gif||height="173px",width="193px"}} {{attachment:compMNEconstr.gif||height="171px",width="190px"}} {{attachment:compMNE.gif||height="170px",width="189px"}} dSPM, sLORETA: {{attachment:compSPM.gif||height="173px",width="192px"}} {{attachment:compLORETA.gif||height="172px",width="191px"}} Now delete all these files when you're done, and keep only the initial solution: wMNE, Constrained. == Computing sources for multiple data files == '''UNCONSTRAINED MN''' 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:<
><