= Import and process Neuromag raw recordings =
This tutorial describes how to process raw Neuromag recordings. It is based on median nerve stimulation acquired at MGH in 2005 with a Neuromag Vectorview 306 system. The sample dataset contains the results for one subject for both arms.

This document shows what to do step by step, but does not really explain what is happening, the meaning of the different options or processes, the issues or bugs you can encounter, and does not provide an exhaustive description of the software features. Those topics are introduced in the basic tutorials based on CTF recordings; so make sure that you followed all those initial tutorials before going through this one.

The script file '''tutorial_mind_neuromag.m''' in the brainstorm3/toolbox/script folder performs exactly the same tasks automatically, without any user interaction. Please have a look at this file if you plan to write scripts to process recordings in .fif format.

<<TableOfContents(2,2)>>

== Download and installation ==
 * It is considered that you followed all the basic tutorials and that you already have a fully working copy of Brainstorm installed on your computer.
 * Go to the [[http://neuroimage.usc.edu/brainstorm3_register/download.php|Download]] page of this website and dowload the file "bst_sample_mind_neuromag.zip", unless you already followed the previous tutorial "Review raw recordings and edit markers".
 * Unzip it in a folder that is not in any of the Brainstorm folders (program folder or database folder). This is really important that you always keep your original data files in a separate folder: the program folder can be deleted when updating the software, and the contents of the database folder is supposed to be manipulated only by the program itself.
 * Start Brainstorm (Matlab scripts or stand-alone version)
 * Select the menu File > Create new protocol. Name it "!TutorialNeuromag", and select all the options as illustrated: No default anatomy, Use individual anatomy, Use one channel file per subject (we have only one run, it's possible to share the channels positions and head models between diffrerent conditions).<<BR>><<BR>> {{attachment:createProtocol.gif}}

== Importing anatomy ==
=== Create the subject ===
 * Select the "Anatomy" exploration mode (first button in the toolbar above the database explorer).
 * Right-click on the protocol node and select "New subject".<<BR>><<BR>> {{attachment:createSubjectMenu.gif}}
 * You can leave the default name "Subject01" and the default options. Then click on Save.<<BR>><<BR>> {{attachment:createSubject.gif}}

=== Import the MRI ===
 * Right-click on the subject node, and select "Import MRI". Select the file: sample_mind_neuromag/anatomy/T1.mri<<BR>><<BR>> {{attachment:mriImportMenu.gif}} {{attachment:mriFile.gif}}
 * The orientation of the MRI is already ok, you just have to mark all the six fiducial points, as explained in the following wiki page: [[CoordinateSystems|Coordinate systems]]. Click on save button when you're done.<<BR>><<BR>> {{attachment:fiducials.gif}}
 * If at some point you do not see the files you've just imported in the database explorer, you just need to refresh the display: by pressing F5, or using the menu File > Refresh tree display.

=== Import the surfaces ===
 * Right-click again on the subject node, and select "Import surfaces". Select the "!FreeSurfer surfaces" file type, and select all the surface files at once, holding the SHIFT or CTRL buttons. Click on Open.<<BR>><<BR>> {{attachment:surfaceFiles.gif}}
 * You should see the following nodes in the database.<<BR>><<BR>> {{attachment:surfListInit.gif}}
  * '''outer skin''': Head, used for surfaces/MRI registration
  * '''inner skull''': Internal surface of the skull, used for the computation of the "overlapping spheres" forward model
  * '''lh''': Left hemisphere
  * '''rh''': Right hemisphere
  * '''pial''': External cortical surface = grey matter/CSF interface (used for source estimation)
  * '''white''': External white matter surface = interface between the grey and the white matter (not used in this tutorial)
 * Downsample the two pial and the two white surfaces to 7,500 vertices each (right-click > Less vertices)
 * Merge the two downsampled hemispheres of the pial surface (select both lh.pial_7500 and rh.pial_7500 files, right-click > Merge surfaces). Rename the new surface in "cortex".
 * Do the same with the white matter. Call the result "white".
 * Delete all the intermediate lh and rh surfaces. Rename the head and the inner skull with shorter names. Display all the surfaces and play with the colors and transparency to check that everything what imported correctly. You should obtain something like this:<<BR>><<BR>> {{attachment:surfListFinal.gif}} {{attachment:surfTest.gif}}

== Importing MEG recordings ==
 * Select the "Functional data / sorted by subject" exploration mode (second button in the toolbar on top of the database explorer).
 * Right click on the subject, select Import MEG/EEG. Select the "Neuromag FIFF" file type, and open the only .fif file in the folder sample_mind_neuromag/data.<<BR>><<BR>> {{attachment:megImportMenu.gif}}
 * Select "Event channel" when asked about the events.<<BR>><<BR>> {{attachment:eventsImport.gif}}
 * You have then to select from which technical track you want to read the events. <<BR>><<BR>> {{attachment:eventsChannels.gif}}
 * In this dataset, the tracks of interest are:
  * STI 001: Contains the onsets of the electric stimulation of the '''right '''arm.
  * STI 002: Contains the onsets of the electric stimulation of the '''left '''arm.
  * Those are not standard settings, they depend on your acquisition setup. In most cases, you would rather import "events" defined with the Neuromag acquisition software, which are saved in tracks STI 014, STI 101, STI 201...
 * Select "STI 001" and "STI 002". The full track is read, all the changes of values are detected from those two tracks and a list of events is created.
 * The following figure appears, that asks how to import those recordings in the Brainstorm database.<<BR>><<BR>> {{attachment:importOptions.gif}}
  * '''Time window''': Time range of interest. Now we want to keep the whole time definition, we are interested by all the stimulations.
  * '''Split''': Useful to import continuous recordings without events. We do not need this here.
  * '''Events selection''': Check the "Use events" option, and select both STI001 and STI002. The value between parenthesis represents the number of occurrences of this event.
  * '''Epoch time''': Time instants that will be extracted before an after each event, to create the epochs that will be saved in the database. Set it to [-100, +300] ms
  * '''Use Signal Space Projections''': Save and use the SSP vectors created by !MaxFilter or during other pre-processing steps. Keep this option selected.
  * '''Remove DC Offset''': Check this option, and select: Time range: [-100, 0] ms. <<BR>>For each epoch, this will: compute the average of each channel over the baseline (pre-stimulus interval: -100ms to 0ms), and subtract it from the channel at all the times in [-100,+300]ms.
  * '''Resample recordings''': Keep this unchecked
  * '''Create new conditions for epochs''': If selected, a new condition is created for each event type (here, it will create a condition "Event_5"). If not selected, all the epochs from all the selected events are saved in the same condition, the one that was selected in the database explorer (if no condition selected, create a condition called "New condition").
  * Click on Import and wait.
 * At the end, you're asked if you want to "Refine the registration with the head  points". This operation improves the initial MRI/MEG registration by  applying and ICP algorithm to make fit the head points digitized at the  MEG acquisition and the scalp surface. Answer yes. You'll see that, even  if the result is not perfect, it dramatically improves the positioning  of the head in the MEG helmet.<<BR>><<BR>> {{attachment:refine.gif}}
 * At the end, you get something like this:<<BR>><<BR>> {{attachment:treeAfterImport1.gif}}
 * The files that appear in the folder "''(Common files)''" will be shared for all the conditions of Subject01. This mechanism allows to compute the forward and inverse models only once for all the conditions when there is only one run (= one channel file = one position of sensors recorded).
 * Rename the condition "STI_001_5" and "STI_002_5" into respectively "Right" and "Left" (Right-click > Rename, or press F2, or click twice waiting a while before the two clicks).<<BR>><<BR>> {{attachment:treeAfterImport2.gif}}
 * __Note__: Another way to proceed would be: 1) to create a new empty condition Right, and 2) inside it import the STI001 events only. Same for Left and STI002. It can be a safer way to proceed, depending on the cases.

== Pre-processing ==
=== Select files to process ===
We want to apply the same pre-processing operations to all the trials of both conditions. Drag'n'drop the "Subject01" node from the database explorer to the tab ''Process1''. Select on the "Recordings" button, then click ''Run''.

 . {{attachment:cutStimList.gif}}

=== Create a processing pipeline ===
 * '''Process > Cut stimulation artifact''': The electric stimulation of the median nerve induces a strong artefact  right after 0ms. This process uses a simple trick to remove this  artifact: re-interpolate the values for the artifact time window (linear  interpolation). It doesn't affect much the data but will make all the  displays much better. Set the time window to [0, 3.9]ms. Select "Overwrite initial file".

 * '''Process > Band-pass filter''': Set the frequency range to [2Hz -90Hz]. Select "Overwrite initial file".
 * '''Average > Average files''': Select "By condition".<<BR>><<BR>> {{attachment:pipeline.gif}}

=== Run the pipeline ===
 * With all those process tags added, the comments of the files are getting a little too long. Rename them respectively Right and Left, by renaming directly the list (the node that contains all the epochs in each condition). It takes a while but will improve a lot the readability later.<<BR>><<BR>> {{attachment:treeAfterPreprocess.gif}}

=== Review the epochs ===
It is always very important to keep an eye on the quality of the data at the different steps of the analysis. There is always a few epochs that are too artifacted or noisy to be used, or one bad sensor. Unfortunately, there are no tools yet in Brainstorm for autoamtic artifact detection and rejection. But here is a procedure for reviewing your epochs manually.

 * Double-click on the Right epochs list to expand the list of files contained in it. Then display the time series for the first Right epoch by double-clicking on it.<<BR>><<BR>> {{attachment:reviewEpochs.gif}}
 * Then press F3 to go to the next epoch (or menu "Navigator > Next data file" in the main Brainstorm window). Do that until you reach the last epoch. For each epoch:
  * Check that the amplitude range in consistant with the previous epochs (in this study: between 1000fT and 2500fT).
  * You can tag some channels as bad: select a sensor by clicking on it and press the Delete key (or right-click > Channel > Mark selected as bad)
  * You can completely delete an epoch if you consider it as unusable (delete the files from the tree).
 * Time series can be also displayed in columns (select the "Display mode" toggle button in the toolbar of the main Brainstorm window). To select create your own selection of sensors, or use a pre-defined set: Right-click on the time-series figure > Display setup > ... Example: Right-temporal sensors only.<<BR>><<BR>> {{attachment:selectionsRT.gif}}
 * In the case of this study, there are so many trials for each condition and the recorded signal is so strong that you do not really need to spend too much time on the selection of the bad channels and bad epochs. If you just keep everything as it is, it is just fine.
 * Do the same for the Left condition.

=== Averaging ===
You can use Brainstorm to work on individual trials or on average recordings. But even if you plan to work on single trials, start your exploration of the recordings by computing an average per condition. It would give you a good idea of the quality of your recordings and pre-processing operations. If you do not see anything looking like the effect you are supposed to observe on the average, it is a complete waste of time to go on with source or time-frequency analysis.

 * Drag'n'drop again the node "Subject01" in the Processes tab, and click on Run process.
 * Select "Average > Average by condition (Grand average). Click on Run.
 * Two new files called "GAVE: ...." appeared in the database. They are the averaged response for each condition, locked on the stimulus.<<BR>><<BR>> {{attachment:treeAfterAverage.gif}}
 * Display the time series for each average: Right-click on the file > MEG (All) > Display time series.
 * Display the 2D topographies for each average: Right-click on the file > MEG (All) > 2D Sensor cap (or press CTRL+T from the time series figures)<<BR>><<BR>> {{attachment:reviewAverages.gif}}
 * Everything looks ok, the time series figures show a clean signal with a very high SNR, and the topographies at the first response peak (about 22ms) show an activity on the left side for the Right condition, and on the right side for the Left condition.
 * Other available displays: 3D Sensor cap, 2D Disc, 2D Layout. Note the shared selection of sensors.<<BR>><<BR>>

== Forward model ==
First step of the source estimation process: establishing a model that describes the way the brain electric activities influence the magnetic fields that are recorded by the MEG sensors. This model can be designated in the documentation by the following terms: head model, forward model, lead field matrix.

=== MEG / MRI registration ===
An accurate forward model requires first of all an accurate registration of the anatomy files (MRI+surfaces) and functional recordings (position of the MEG sensors and EEG electrodes). A basic registration is provided by the alignment of the fiducials (Nasion, LPA, RPA), that were both located before the acquisition of the recordings and marked on the MRI in Brainstorm. This registration based on three points only can be very inaccurate, because it is sometimes hard to identify clearly those points, and not everybody identify them the same way. Two methods described in the [[Tutorials/TutImportRecordings|introduction tutorial #3]] may help you to improve this registration.

 * Right-click on the channel file > MRI Registration > Check MEG-MRI registration<<BR>><<BR>> {{attachment:checkRegistration.gif}}  {{attachment:checkRegistrationFig.gif}}
 * The yellow squares represent the Neuromag MEG sensors, and the green dots represent the head points that where saved with the MEG recordings.
 * In this case the registration looks acceptable. If you consider it is not your case, you can try the two menus:
  * Auto-registration: This operation will try to find a better registration of the head points (green dots), on the head surface (grey surface)
  * Manual resgistration: You rotate/translate the sensor array manually respect to the head surface.

=== Compute head model ===
 * Right-click on any node that contains the channel file (including the channel file itself), and select: "Compute head model".<<BR>><<BR>> {{attachment:headmodelStart.gif}}
 * Select "overlapping spheres" and click on Run. Select "Inner skull" when asked what surface to use to estimate the spheres. For more information, read the tutorial "[[Tutorials/TutHeadModel|Computing a head model]]".<<BR>><<BR>> {{attachment:headmodelSurface.gif}}
 * When the computation is done, close the "Check spheres" figure. The lead field matrix is saved in file "Overlapping spheres" in "Common files", in a the "Gain" field.

== Noise covariance matrix ==
To estimate the sources properly, we need an estimation of the noise level for each sensor. A good way to do this is to compute the covariance matrix of the concatenation of the baselines from all the trials in both conditions.

 * Select at the same time the two groups of trials (Right and Left). To do this: hold the Control key and click successively on the Right and the Left trial lists.
 * Right-click on one of them and select: Noise covariance > Compute from recordings. Leave all the options unchanged, and click on Ok.<<BR>><<BR>> {{attachment:noisecovMenu.gif}} {{attachment:noisecovOptions.gif}}
 * This operation computes the noise covariance matrix based on the baseline, [-100, 0]ms, of all the trials (624 files), and keep only the diagonal values (more stable). The result is stored in a new file in the "Common files" folder.

== Source estimation ==
Reconstruction of the cortical sources with a weighted minimum norm estimator (wMNE).

 * Right-click on "Common files", (or on the head model, or on the subject node, this is all equivalent), and select "Compute sources". <<BR>><<BR>> {{attachment:inverseMenu.gif}}
 * Select "Yes" to ignore the warning about noise covariance used at the same time for averages and single trials (if you do not see such a warning, that's fine too). In this case, this warning does not really apply, the results are the same if we compute separately the minimum norm estimator for the single trials and the averages.
 * Select "Normal mode" (not "Expert mode"), wMNE, and both gradiometers and magnetometers, as in the figure below. Click on Run.<<BR>><<BR>> {{attachment:inverseOptions.gif}} {{attachment:treeAfterInverse.gif}}
 * A shared inversion kernel was created in the ''(Common files)'' folder, and a link node is now visible for each recording file, single trials and averages. For more information about what those links mean and the operations performed to display them, please refer to the [[Tutorials/TutSourceEstimation|tutorial #4 "Source estimation"]].
 * For the average of condition Right: display the time series, the 2D sensor topography, the sources on the cortex, the sources on 3D orthogonal slices, and the sources in the MRI Viewer.<<BR>><<BR>> {{attachment:inverseAllDisplay.gif}}
 * Create a scout called LeftSS that covers the most active region of the cortex at 33ms for condition Right. Display its time series for conditions Left and Right at the same time (select both Left and Right source files, right click > Cortical activations > Scouts time series (check "Overlay: conditions").<<BR>><<BR>> {{attachment:inverseScoutLeft.gif}}

== Time-frequency ==