WARNING: This tutorial uses an old version of SimNIBS, which includes the segmentation pipeline headreco. Please see instead the new version of this tutorial, based on SimNIBS 4 and the CHARM pipeline. However, we recommned to use the current version with SimNibs3 headreco. We are working with SimNibs team to keep the Headreco availble for the users [headreco is available on all releases prior to v4.0].


FEM tutorial with MEG/EEG median nerve stimulation (headreco)

Authors: Takfarinas Medani, Juan Garcia-Prieto, Wayne Mead, Michael Funke, Francois Tadel

This tutorial introduces the most advanced FEM modeling options available in the Brainstorm environment, applied to MEG+EEG recordings of a median nerve stimulation. The pipeline presented here includes: FEM mesh reconstruction with SimNIBS3/CAT12, FEM head model including DTI tensors for modeling anisotropic conductivities (using BrainSuite), FEM forward model estimation with DUNEuro. This pipeline requires many third-party programs and very long computation times. For a simpler FEM pipeline, please refer to the tutorial: Realistic head model: FEM with DUNEuro.

Note that the operations used here are not detailed, the goal of this tutorial is not to introduce Brainstorm to new users. For in-depth explanations of the interface and theoretical foundations, please refer to the introduction tutorials.

Dataset description

Experiment

The experiment consists of an unilateral median nerve stimulation conducted in a MEG laboratory with an Elekta Triux (Megin, Finland) scanner.

  • The stimulation signal was a square-wave pulse with 2Hz frequency and duration of 0.2ms.
  • An ISI (inter-stimulus interval) of 500ms with a variation of ±20ms in order to be able to average out time-locked noise to the stimulation, while remaining unnoticeable by the subject.
  • The stimulation was performed on the left hand/wrist, for two minutes, with a Digitimer DS7A stimulator. An electrode was placed on the hand/wrist while current value was tuned to match the motor threshold of the subject on the stimulated hand, with a result of 10~12mA approximately.
  • Recordings were performed with a 1kHz sampling rate. Continuous HPI was disabled during these recordings. High-pass filters were set to DC for MEG channels and 0.03Hz for EEG channels.
  • All MEG files underwent a MaxFilter (version 2.3.13) tsss post-processing.

MRI imaging scanning (Philips Medical Systems):

  • T1w image without contrast: 3T field strength, flip angle 8°, TR 7.9s, dimensions 512 × 512 × 200, and voxel dimensions of 0.469 × 0.469 × 0.939 millimeters
  • T2w spin-echo image: 3T field strength, 78.57% phase FOV, 90° flip angle, SAR 0.327, voxel dimensions 560 × 560 × 55, dimensions of 0.429 × 0.429 × 3 millimeters
  • DWI sequences: voxel dimensions 512 x 512 x 200 and dimensions of 0.469 × 0.469 × 0.939 millimeters, diffusion-sensitizing gradients in 32 non-collinear directions

Files

The dataset we distribute with this tutorial follows the Brain Imaging Data Structure (BIDS) standard for neuroimaging data organization. The files that will be imported in this tutorial are the following:

sample_fem/sub-fem01/: Raw data for subject fem01

  • ses-meg/: Simultaneous recordings of MEG and EEG.

    • sub-fem01_ses-meg_task-mediannerve_run-01_proc-tsss_meg.fif

  • ses-mri/: Imaging exams.

    • anat/sub-fem01_ses-mri_T1w.nii.gz: T1-weighted MRI

    • anat/sub-fem01_ses-mri_T2w.nii.gz: T2-weighted MRI

    • dwi/sub-fem01_ses-mri_dwi.*: Diffusion-Weighted Imaging (DWI)

License

Creative Commons CC0 1.0 Universal. This dataset is distributed in the public domain.

Download and installation

Requirements

Additional programs

You need to install the following programs (listed by order of appearance):

Download the dataset

Import the anatomy

T1 MRI

T2 MRI

Diffusion imaging

The FEM has the ability to incorporate anisotropic conductivity from MRI diffusion imaging, which is particularly interesting for the modeling of the white matter. Brainstorm can load the Diffusion-Weighted Images (DWI), and compute the tensors (DTI) using the BrainSuite Diffusion Pipeline (BDP). This requires BrainSuite to be installed on your computer, with the bdp program available in the system path.

FEM head model

The FEM approach requires a segmentation of the head volume in different tissues, represented as hexahedral or tetrahedral 3D meshes. The methods available within Brainstorm are listed in the tutorial FEM mesh generation. Here we illustrate only the use of SimNIBS.

FEM mesh with SimNIBS

Running SimNIBS:

At the end of this computation, new files are available in the database:

At the end of the process, make sure that the file "cortex_15002V" is selected (downsampled pial surface, which will be used for the source estimation). If it is not, double-click on it to select it as the default cortex surface.

simnibs.gif

Remesh with Iso2mesh

This step is optional, to be considered only when the DUNEuro FEM computation fails. In some cases, the FEM mesh generated with SimNIBS causes issues with the DUNEuro FEM solver, due to the air cavities that are not tesselated. In order to avoid these possible issues in the next steps, we will correct the mesh using the Iso2mesh functions integrated within Brainstorm.

The following figure shows the two model, left is the initial model obtained with SimNibs, right is the second model obtained from the Iso2Mesh remesh: the air cavity in the bottom-right corner is now fully tesselated. Note that you can also use this process to generate FEM models with higher mesh densities.

femMeshSimNibsVSiso2mesh2.jpg

FEM tensors

We can now incorporate the diffusion information into the FEM model, and compute anisotropic conductivity tensors for the tetrahedral elements of the white matter (we consider the other tissues to have isotropic conductivities).

Visualization

Brainstorms include the possibilities to display the FEM head models and the tensors, users can also overlay the display with the MRI as well as with the different surfaces.

overlayModalities.png

Right-click on the FEM mesh > Display FEM tensors: The FEM tensors can be displayed on the mesh or MRI, as arrows on the main eigenvector or as ellipsoids on each FEM element.

dispTensorMenu.png

tensorsOnBrain.jpg

tensorsOnMRI.jpg

To configure the display: right-click on the figure > FEM tensors menu. The keyboard shortcuts for changing the size of the displayed tensors are the Up and Down arrows keys. You can also switch the display mode by using the shortcut "Shift + Space".

Advanced

On the hard drive

The DTI-EIG file as the same structure as any MRI file, with 12 volumes stacked along the 4th dimension. From 1 to 9: components of the three eigenvectors; from 10 to 12: the values of their norm to the eigenvalue.

The FEM mesh contains the following fields:

BEM head model

For the purpose of comparison between the FEM and the BEM, we will generate also the BEM surfaces for this subject and we will follow the same step as explained in the BEM tutorial.

Access the recordings

Channel file

Pre-processing

Frequency filters

EEG: Average reference

Epoching

In this experiment, the electric stimulation is sent with a frequency of 2Hz, meaning that the inter-stimulus interval is 500ms. We are going to import epochs of 300ms around the stimulation events (-100 to 200ms).

Averaging

Forward model

We are going to use the realistic FEM model previously generated from the MRI. Go to the "Anatomy" view, and make sure that the FEM head model is highlighted in green color (it should be the case if you have only one model). You may also highlight the cortex to use for the computation (select the cortex_15002V).

You can compute the forward model both for EEG and MEG simultaneously, however, using the high mesh resolution model we recommend to compute separately the head model for each modality (EEG and then MEG). The time required for EEG is around one hour for ~70 channels, for the MEG with 306 sensors it can take up to 4 hours or more (with the integrations points).

The EEG/MEG FEM computation depends on the computation of the FEM transfer matrix, which is related to the resolution of the FEM head mesh (number of vertices) and the number of sensors. In most of the case, the number of EEG sensors is lower than the number of MEG sensors. Furthermore, internally the MEG sensors modeling uses the integrations points, which increase the number of computation points (~ multiplied by 4 for the magnetometers and by 8 for the gradiometers). Therefore, the MEG requires more time than the EEG.

To reduce the MEG computation time, there are some tips:

  1. Use only the inners tissues (wm, gm, and CSF) ==> reduce the number of vertices

  2. Do not use the integration points ==> reduce the number of virtual sensors

  3. These parameters can be tuned from the DUNEuro options panel (see the advanced panel)

EEG with DTI tensors

MEG with DTI tensors

EEG with isotropic conductivity

If the DTIs are not available, it is possible to use the isotropic conductivities instead. In this tutorial, we will compute another forward model in order to compare the sources obtained with anisotropic or isotropic conductivities.

MEG with isotropic conductivity

Source estimation

Noise covariance matrix

Inverse model

It is recommended to study separately the two modalities because Brainstorm does not offer any reliable method for combining MEG and EEG source imaging yet. More explanations in the Source estimation tutorial.

Dipole scanning

Comparison: DTI vs isotropic

The following figures (left EEG, right MEG), show the localization of the dipoles on the MRI at the 21ms. The ISO model is colored in green whereas the DTI model is colored in red.

eegMegDipolesSagittal.jpg eegMegDipolesAxial.jpg eegMegDipolesCoronal.jpg

In this experiment, the anisotropy does not show a significant effect on the source localization, whereas it shows a slight difference in the orientation.

Comparison: MEG vs EEG

When we compare between the EEG and the MEG, there is a difference of 20mm between the localization between the MEG and the EEG as well as a difference in the orientation in the coronal view.

If we compare the difference between the EEG and MEG dipoles, the following figures (left anisotropy, right isotropic), show the localization of the dipoles on the MRI at the 21ms. The EEG dipole is colored in green whereas the MEG dipole is colored in red.

eegMegDipolesSagittal2.jpg eegMegDipolesAxial2.jpg eegMegDipolesCoronal2.jpg

Aniso: Eeg [-7.5 -35.7 89.8] vs Meg [-0.2 -33.6 88.6] ==> distance = 7.7mm

Iso: Eeg [-7.5 -35.7 89.8] vs Meg [-0.2 -33.6 88.6] ==> distance = 7.7mm

In both modalities, the dipoles are located exactly in the same positions, and there is 7.7mm between the two dipoles. Furthermore, we notice a difference on the orientation on the coronal and sagittal views

Comparison: FEM vs BEM/OS

We will now compare qualitatively the FEM results with the default methods implemented in Brainstorm: Overlapping spheres for MEG, BEM for EEG.

In the following figures, left is the dipoles computed from the EEG (BEM in red and FEM in green), right are the dipoles computed from the MEG (OS in red and FEM in green)

eegMegDipolesSagittal3.jpg eegMegDipolesAxial3.jpg eegMegDipolesCoronal3.jpg

EEG: fem [-7.5 -35.7 89.8] ; bem [-10.7 -35.0 89.5] ==> 3.3mm

MEG: fem [-0.2 -33.6 88.6] ; os [-10.7 -39.3 83.2] ==> 13.1 mm

As expected, the slight difference on the localization can be explained by the difference on the head shape, the conductivity values as well as the resolution method.

Advanced

Additional documentation

Articles

Scripting

The following script from the Brainstorm distribution reproduces the analysis presented in this tutorial page: brainstorm3/toolbox/script/tutorial_fem_headreco.m

1 function tutorial_fem_headreco(tutorial_dir) 2 % TUTORIAL_FEM_HEADRECO: Script that reproduces the online tutorial "FEM tutorial: MEG/EEG Median nerve stimulation (headreco)" 3 % 4 % REFERENCE: 5 % https://neuroimage.usc.edu/brainstorm/Tutorials/FemMedianNerve 6 % 7 % INPUTS: 8 % tutorial_dir: Directory where the sample_fem.zip file has been unzipped 9 10 % @============================================================================= 11 % This function is part of the Brainstorm software: 12 % https://neuroimage.usc.edu/brainstorm 13 % 14 % Copyright (c) University of Southern California & McGill University 15 % This software is distributed under the terms of the GNU General Public License 16 % as published by the Free Software Foundation. Further details on the GPLv3 17 % license can be found at http://www.gnu.org/copyleft/gpl.html. 18 % 19 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 20 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 21 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 22 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 23 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 24 % 25 % For more information type "brainstorm license" at command prompt. 26 % =============================================================================@ 27 % 28 % Author: Francois Tadel, 2021 29 30 31 % ===== FILES TO IMPORT ===== 32 % You have to specify the folder in which the tutorial dataset is unzipped 33 if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir) 34 error('The first argument must be the full path to the tutorial dataset folder.'); 35 end 36 % Build the path of the files to import 37 T1Nii = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'anat', 'sub-fem01_ses-mri_T1w.nii.gz'); 38 T2Nii = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'anat', 'sub-fem01_ses-mri_T2w.nii.gz'); 39 DwiNii = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'dwi', 'sub-fem01_ses-mri_dwi.nii.gz'); 40 DwiBval = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'dwi', 'sub-fem01_ses-mri_dwi.bval'); 41 DwiBvec = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-mri', 'dwi', 'sub-fem01_ses-mri_dwi.bvec'); 42 FifFile = fullfile(tutorial_dir, 'sample_fem', 'sub-fem01', 'ses-meg', 'meg', 'sub-fem01_ses-meg_task-mediannerve_run-01_proc-tsss_meg.fif'); 43 % Check if the folder contains the required files 44 if ~file_exist(T1Nii) || ~file_exist(FifFile) 45 error(['The folder ' tutorial_dir ' does not contain the folder from the file sample_fem.zip.']); 46 end 47 % Subject name 48 SubjectName = 'Subject01'; 49 % Latency of interest 50 Latency = 0.022; 51 52 53 % ===== CHECK SOFTWARE DEPENDENCIES ===== 54 % Start brainstorm without the GUI 55 if ~brainstorm('status') 56 brainstorm nogui 57 end 58 % SimNIBS 3 / headreaco 59 status = system('headreco --version'); 60 if (status ~= 0) 61 error('SimNIBS is not installed or not added to the system path: the command "headreco" could not be found.'); 62 end 63 % BrainSuite 64 if ~file_exist(bst_fullfile(bst_get('BrainSuiteDir'), 'bin')) 65 error('BrainSuite is not configured in the Brainstorm preferences.'); 66 end 67 % SPM12 68 [isInstalled, errMsg] = bst_plugin('Install', 'spm12', 0); 69 if ~isInstalled 70 error(['Could not install plugin: spm12' 10 errMsg]); 71 end 72 % CAT12 73 [isInstalled, errMsg] = bst_plugin('Install', 'cat12', 0); 74 if ~isInstalled 75 error(['Could not install plugin: cat12' 10 errMsg]); 76 end 77 % Iso2mesh 78 [isInstalled, errMsg] = bst_plugin('Install', 'iso2mesh', 0); 79 if ~isInstalled 80 error(['Could not install plugin: iso2mesh' 10 errMsg]); 81 end 82 83 84 % ===== CREATE PROTOCOL ===== 85 % The protocol name has to be a valid folder name (no spaces, no weird characters...) 86 ProtocolName = 'TutorialFem'; 87 % Delete existing protocol 88 gui_brainstorm('DeleteProtocol', ProtocolName); 89 % Create new protocol 90 gui_brainstorm('CreateProtocol', ProtocolName, 0, 1); 91 % Start a new report 92 bst_report('Start'); 93 94 95 %% ===== IMPORT ANATOMY ===== 96 % ===== IMPORT MRI VOLUMES ===== 97 % Create subject 98 [sSubject, iSubject] = db_add_subject(SubjectName, [], 0, 0); 99 % Import T1 MRI 100 T1File = import_mri(iSubject, T1Nii, 'ALL', 0, 0); 101 % Compute the MNI normalization 102 bst_normalize_mni(T1File); 103 % Import T2 MRI 104 T2File = import_mri(iSubject, T2Nii, 'ALL', 0, 0); 105 % Volumes are not registered: Register with SPM 106 mri_coregister(T2File, T1File, 'spm', 1); 107 % Delete the non-registered T2 108 file_delete(file_fullpath(T2File), 1); 109 db_reload_subjects(iSubject); 110 111 % ===== IMPORT DTI ===== 112 % Process: Convert DWI to DTI (BrainSuite) 113 bst_process('CallProcess', 'process_dwi2dti', [], [], ... 114 'subjectname', SubjectName, ... 115 'dwifile', {DwiNii, 'DWI-NII'}, ... 116 'bvalfile', {DwiBval, 'DWI-BVAL'}, ... 117 'bvecfile', {DwiBvec, 'DWI-BVEC'}); 118 119 % ===== RUN SIMNIBS ===== 120 % Process: Generate FEM mesh 121 bst_process('CallProcess', 'process_fem_mesh', [], [], ... 122 'subjectname', SubjectName, ... 123 'method', 'simnibs3', ... % SimNIBS:Call SimNIBS to segment and mesh the T1 (and T2) MRI. 124 'vertexdensity', 0.5, ... 125 'nvertices', 15000, ... 126 'zneck', 0); 127 % Select default cortex: Low-resolution 128 [sSubject, iSubject, iCortex] = bst_get('SurfaceFile', bst_fullfile(SubjectName, 'tess_cortex_pial_low.mat')); 129 db_surface_default(iSubject, 'Cortex', iCortex); 130 131 % ===== REMESH WITH ISO2MESH ===== 132 % Part skipped because: 133 % 1) Not needed: DUNEuro works with cavities in this case 134 % 2) No process call for: Extract surfaces (import_femlayers) 135 % 3) Remeshing with iso2mesh doesn't work... 136 137 % ===== COMPUTE FEM TENSORS ===== 138 % Process: Compute FEM tensors 139 bst_process('CallProcess', 'process_fem_tensors', [], [], ... 140 'subjectname', SubjectName, ... 141 'femcond', struct(... 142 'FemCond', [0.14, 0.33, 1.79, 0.008, 0.43], ... 143 'isIsotropic', [0, 1, 1, 1, 1], ... 144 'AnisoMethod', 'ema+vc', ... 145 'SimRatio', 10, ... 146 'SimConstrMethod', 'wolters')); 147 148 % ===== COMPUTE BEM SURFACES ===== 149 % Process: Generate BEM surfaces 150 bst_process('CallProcess', 'process_generate_bem', [], [], ... 151 'subjectname', SubjectName, ... 152 'nscalp', 1922, ... 153 'nouter', 1922, ... 154 'ninner', 1922, ... 155 'thickness', 4, ... 156 'method', 'brainstorm'); 157 158 159 %% ===== ACCESS THE RECORDINGS ===== 160 % Process: Create link to raw file 161 sFilesRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 162 'subjectname', SubjectName, ... 163 'datafile', {FifFile, 'FIF'}, ... 164 'channelreplace', 0, ... 165 'channelalign', 1, ... % Automatic registration with digitized head points 166 'evtmode', 'value'); 167 168 % Process: Events: Read from channel 169 sFilesRaw = bst_process('CallProcess', 'process_evt_read', sFilesRaw, [], ... 170 'stimchan', 'STI101', ... 171 'trackmode', 1, ... % Value: detect the changes of channel value 172 'zero', 0); 173 174 % Process: Project electrodes on scalp 175 bst_process('CallProcess', 'process_channel_project', sFilesRaw, []); 176 177 % Process: Snapshot: Sensors/MRI registration 178 bst_process('CallProcess', 'process_snapshot', sFilesRaw, [], ... 179 'target', 1, ... % Sensors/MRI registration 180 'modality', 1, ... % MEG (All) 181 'orient', 1, ... % left 182 'Comment', 'MEG/MRI Registration'); 183 % Process: Snapshot: Sensors/MRI registration 184 bst_process('CallProcess', 'process_snapshot', sFilesRaw, [], ... 185 'target', 1, ... % Sensors/MRI registration 186 'modality', 4, ... % EEG 187 'orient', 1, ... % left 188 'Comment', 'EEG/MRI Registration'); 189 190 191 %% ===== PRE-PROCESSING ===== 192 % Process: Power spectrum density (Welch) 193 sFilesPsd = bst_process('CallProcess', 'process_psd', sFilesRaw, [], ... 194 'timewindow', [18, 148.999], ... 195 'win_length', 5, ... 196 'win_overlap', 50, ... 197 'units', 'physical', ... % Physical: U2/Hz 198 'sensortypes', 'MEG, EEG', ... 199 'win_std', 0, ... 200 'edit', struct(... 201 'Comment', 'Power', ... 202 'TimeBands', [], ... 203 'Freqs', [], ... 204 'ClusterFuncTime', 'none', ... 205 'Measure', 'power', ... 206 'Output', 'all', ... 207 'SaveKernel', 0)); 208 209 % Process: Snapshot: Frequency spectrum 210 bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ... 211 'target', 10, ... % Frequency spectrum 212 'Comment', 'Power spectrum density'); 213 214 % Process: Band-pass:20Hz-250Hz 215 sFilesBand = bst_process('CallProcess', 'process_bandpass', sFilesRaw, [], ... 216 'sensortypes', 'MEG, EEG', ... 217 'highpass', 20, ... 218 'lowpass', 250, ... 219 'tranband', 0, ... 220 'attenuation', 'strict', ... % 60dB 221 'ver', '2019', ... % 2019 222 'mirror', 0, ... 223 'read_all', 0); 224 225 % Process: Notch filter: 60Hz 120Hz 180Hz 226 sFilesClean = bst_process('CallProcess', 'process_notch', sFilesBand, [], ... 227 'sensortypes', 'MEG, EEG', ... 228 'freqlist', [60, 120, 180], ... 229 'cutoffW', 2, ... 230 'useold', 0, ... 231 'read_all', 0); 232 233 % Process: Delete selected files 234 bst_process('CallProcess', 'process_delete', [sFilesRaw, sFilesBand], [], ... 235 'target', 2); % Delete selected folders 236 237 % Process: Re-reference EEG 238 bst_process('CallProcess', 'process_eegref', sFilesClean, [], ... 239 'eegref', 'AVERAGE', ... 240 'sensortypes', 'EEG'); 241 242 243 %% ===== IMPORT RECORDINGS ===== 244 % Process: Import MEG/EEG: Events 245 sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFilesClean, [], ... 246 'subjectname', SubjectName, ... 247 'condition', '', ... 248 'eventname', '2', ... 249 'timewindow', [18, 148.999], ... 250 'epochtime', [-0.1, 0.2], ... 251 'createcond', 0, ... 252 'ignoreshort', 1, ... 253 'usectfcomp', 1, ... 254 'usessp', 1, ... 255 'freq', [], ... 256 'baseline', []); 257 258 % Process: Average: By trial group (folder average) 259 sFilesAvg = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ... 260 'avgtype', 5, ... % By trial group (folder average) 261 'avg_func', 1, ... % Arithmetic average: mean(x) 262 'weighted', 0, ... 263 'keepevents', 0); 264 265 % Process: Snapshot: Recordings time series 266 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ... 267 'target', 5, ... % Recordings time series 268 'modality', 4, ... % EEG 269 'Comment', 'EEG ERP'); 270 % Process: Snapshot: Recordings topography (contact sheet) 271 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ... 272 'target', 6, ... % Recordings topography (ont time) 273 'modality', 4, ... % EEG 274 'time', Latency, ... 275 'Comment', 'EEG ERP'); 276 % Process: Snapshot: Recordings time series 277 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ... 278 'target', 5, ... % Recordings time series 279 'modality', 1, ... % MEG (All) 280 'Comment', 'MEG ERF'); 281 % Process: Snapshot: Recordings topography (contact sheet) 282 bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ... 283 'target', 6, ... % Recordings topography (ont time) 284 'modality', 1, ... % MEG (All) 285 'time', Latency, ... 286 'Comment', 'MEG ERF'); 287 288 289 %% ===== NOISE COVARIANCE ===== 290 % Process: Compute covariance (noise or data) 291 bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ... 292 'baseline', [-0.1, -0.01], ... 293 'sensortypes', 'MEG, EEG', ... 294 'target', 1, ... % Noise covariance (covariance over baseline time window) 295 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 296 'replacefile', 1); % Replace 297 298 299 %% ===== FORWARD: FEM EEG DTI ===== 300 % Process: Compute head model 301 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 302 'Comment', 'DUNEuro FEM EEG DTI', ... 303 'sourcespace', 1, ... % Cortex surface 304 'meg', 1, ... % 305 'eeg', 4, ... % DUNEuro FEM 306 'ecog', 1, ... % 307 'seeg', 1, ... % 308 'duneuro', struct(... 309 'FemCond', [], ... 310 'FemSelect', [1, 1, 1, 1, 1], ... 311 'UseTensor', 1, ... 312 'Isotropic', 1, ... 313 'SrcShrink', 0, ... 314 'SrcForceInGM', 0, ... 315 'FemType', 'fitted', ... 316 'SolverType', 'cg', ... 317 'GeometryAdapted', 0, ... 318 'Tolerance', 1e-08, ... 319 'ElecType', 'normal', ... 320 'MegIntorderadd', 0, ... 321 'MegType', 'physical', ... 322 'SolvSolverType', 'cg', ... 323 'SolvPrecond', 'amg', ... 324 'SolvSmootherType', 'ssor', ... 325 'SolvIntorderadd', 0, ... 326 'DgSmootherType', 'ssor', ... 327 'DgScheme', 'sipg', ... 328 'DgPenalty', 20, ... 329 'DgEdgeNormType', 'houston', ... 330 'DgWeights', 1, ... 331 'DgReduction', 1, ... 332 'SolPostProcess', 1, ... 333 'SolSubstractMean', 0, ... 334 'SolSolverReduction', 1e-10, ... 335 'SrcModel', 'venant', ... 336 'SrcIntorderadd', 0, ... 337 'SrcIntorderadd_lb', 2, ... 338 'SrcNbMoments', 3, ... 339 'SrcRefLen', 20, ... 340 'SrcWeightExp', 1, ... 341 'SrcRelaxFactor', 6, ... 342 'SrcMixedMoments', 1, ... 343 'SrcRestrict', 1, ... 344 'SrcInit', 'closest_vertex', ... 345 'BstSaveTransfer', 0, ... 346 'BstEegTransferFile', 'eeg_transfer.dat', ... 347 'BstMegTransferFile', 'meg_transfer.dat', ... 348 'BstEegLfFile', 'eeg_lf.dat', ... 349 'BstMegLfFile', 'meg_lf.dat', ... 350 'UseIntegrationPoint', 1, ... 351 'EnableCacheMemory', 0, ... 352 'MegPerBlockOfSensor', 0), ... 353 'channelfile', ''); 354 % Process: Compute sources [2018] 355 sSrcEegFemDti = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 356 'output', 1, ... % Kernel only: shared 357 'inverse', struct(... 358 'Comment', 'Dipoles: EEG FEM DTI', ... 359 'InverseMethod', 'gls', ... 360 'InverseMeasure', 'performance', ... 361 'SourceOrient', {{'free'}}, ... 362 'Loose', 0.2, ... 363 'UseDepth', 1, ... 364 'WeightExp', 0.5, ... 365 'WeightLimit', 10, ... 366 'NoiseMethod', 'median', ... 367 'NoiseReg', 0.1, ... 368 'SnrMethod', 'rms', ... 369 'SnrRms', 1e-06, ... 370 'SnrFixed', 3, ... 371 'ComputeKernel', 1, ... 372 'DataTypes', {{'EEG'}})); 373 % Process: Dipole scanning 374 sDipEegFemDti = bst_process('CallProcess', 'process_dipole_scanning', sSrcEegFemDti, [], ... 375 'timewindow', [Latency, Latency], ... 376 'scouts', {}); 377 378 379 %% ===== FORWARD: FEM MEG DTI ===== 380 % Process: Compute head model 381 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 382 'Comment', 'DUNEuro FEM MEG DTI', ... 383 'sourcespace', 1, ... % Cortex surface 384 'meg', 5, ... % DUNEuro FEM 385 'eeg', 1, ... % 386 'ecog', 1, ... % 387 'seeg', 1, ... % 388 'duneuro', struct(... 389 'FemCond', [], ... 390 'FemSelect', [1, 1, 1, 0, 0], ... 391 'UseTensor', 1, ... 392 'Isotropic', 1, ... 393 'SrcShrink', 0, ... 394 'SrcForceInGM', 0, ... 395 'FemType', 'fitted', ... 396 'SolverType', 'cg', ... 397 'GeometryAdapted', 0, ... 398 'Tolerance', 1e-08, ... 399 'ElecType', 'normal', ... 400 'MegIntorderadd', 0, ... 401 'MegType', 'physical', ... 402 'SolvSolverType', 'cg', ... 403 'SolvPrecond', 'amg', ... 404 'SolvSmootherType', 'ssor', ... 405 'SolvIntorderadd', 0, ... 406 'DgSmootherType', 'ssor', ... 407 'DgScheme', 'sipg', ... 408 'DgPenalty', 20, ... 409 'DgEdgeNormType', 'houston', ... 410 'DgWeights', 1, ... 411 'DgReduction', 1, ... 412 'SolPostProcess', 1, ... 413 'SolSubstractMean', 0, ... 414 'SolSolverReduction', 1e-10, ... 415 'SrcModel', 'venant', ... 416 'SrcIntorderadd', 0, ... 417 'SrcIntorderadd_lb', 2, ... 418 'SrcNbMoments', 3, ... 419 'SrcRefLen', 20, ... 420 'SrcWeightExp', 1, ... 421 'SrcRelaxFactor', 6, ... 422 'SrcMixedMoments', 1, ... 423 'SrcRestrict', 1, ... 424 'SrcInit', 'closest_vertex', ... 425 'BstSaveTransfer', 0, ... 426 'BstEegTransferFile', 'eeg_transfer.dat', ... 427 'BstMegTransferFile', 'meg_transfer.dat', ... 428 'BstEegLfFile', 'eeg_lf.dat', ... 429 'BstMegLfFile', 'meg_lf.dat', ... 430 'UseIntegrationPoint', 1, ... 431 'EnableCacheMemory', 0, ... 432 'MegPerBlockOfSensor', 0), ... 433 'channelfile', ''); 434 % Process: Compute sources [2018] 435 sSrcMegFemDti = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 436 'output', 1, ... % Kernel only: shared 437 'inverse', struct(... 438 'Comment', 'Dipoles: MEG FEM DTI', ... 439 'InverseMethod', 'gls', ... 440 'InverseMeasure', 'performance', ... 441 'SourceOrient', {{'free'}}, ... 442 'Loose', 0.2, ... 443 'UseDepth', 1, ... 444 'WeightExp', 0.5, ... 445 'WeightLimit', 10, ... 446 'NoiseMethod', 'median', ... 447 'NoiseReg', 0.1, ... 448 'SnrMethod', 'rms', ... 449 'SnrRms', 1e-06, ... 450 'SnrFixed', 3, ... 451 'ComputeKernel', 1, ... 452 'DataTypes', {{'MEG GRAD', 'MEG MAG'}})); 453 % Process: Dipole scanning 454 sDipMegFemDti = bst_process('CallProcess', 'process_dipole_scanning', sSrcMegFemDti, [], ... 455 'timewindow', [Latency, Latency], ... 456 'scouts', {}); 457 458 459 %% ===== FORWARD: FEM EEG ISO ===== 460 % Get updated subject structure 461 sSubject = bst_get('Subject', iSubject); 462 % Get FEM file 463 FemFile = file_fullpath(sSubject.Surface(sSubject.iFEM).FileName); 464 % Remove tensors 465 process_fem_tensors('ClearTensors', FemFile); 466 % Process: Compute head model 467 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 468 'Comment', 'DUNEuro FEM EEG ISO', ... 469 'sourcespace', 1, ... % Cortex surface 470 'meg', 1, ... % 471 'eeg', 4, ... % DUNEuro FEM 472 'ecog', 1, ... % 473 'seeg', 1, ... % 474 'duneuro', struct(... 475 'FemCond', [0.14, 0.33, 1.79, 0.008, 0.43], ... 476 'FemSelect', [1, 1, 1, 1, 1], ... 477 'UseTensor', 0, ... 478 'Isotropic', 1, ... 479 'SrcShrink', 0, ... 480 'SrcForceInGM', 0, ... 481 'FemType', 'fitted', ... 482 'SolverType', 'cg', ... 483 'GeometryAdapted', 0, ... 484 'Tolerance', 1e-08, ... 485 'ElecType', 'normal', ... 486 'MegIntorderadd', 0, ... 487 'MegType', 'physical', ... 488 'SolvSolverType', 'cg', ... 489 'SolvPrecond', 'amg', ... 490 'SolvSmootherType', 'ssor', ... 491 'SolvIntorderadd', 0, ... 492 'DgSmootherType', 'ssor', ... 493 'DgScheme', 'sipg', ... 494 'DgPenalty', 20, ... 495 'DgEdgeNormType', 'houston', ... 496 'DgWeights', 1, ... 497 'DgReduction', 1, ... 498 'SolPostProcess', 1, ... 499 'SolSubstractMean', 0, ... 500 'SolSolverReduction', 1e-10, ... 501 'SrcModel', 'venant', ... 502 'SrcIntorderadd', 0, ... 503 'SrcIntorderadd_lb', 2, ... 504 'SrcNbMoments', 3, ... 505 'SrcRefLen', 20, ... 506 'SrcWeightExp', 1, ... 507 'SrcRelaxFactor', 6, ... 508 'SrcMixedMoments', 1, ... 509 'SrcRestrict', 1, ... 510 'SrcInit', 'closest_vertex', ... 511 'BstSaveTransfer', 0, ... 512 'BstEegTransferFile', 'eeg_transfer.dat', ... 513 'BstMegTransferFile', 'meg_transfer.dat', ... 514 'BstEegLfFile', 'eeg_lf.dat', ... 515 'BstMegLfFile', 'meg_lf.dat', ... 516 'UseIntegrationPoint', 1, ... 517 'EnableCacheMemory', 0, ... 518 'MegPerBlockOfSensor', 0), ... 519 'channelfile', ''); 520 % Process: Compute sources [2018] 521 sSrcEegFemIso = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 522 'output', 1, ... % Kernel only: shared 523 'inverse', struct(... 524 'Comment', 'Dipoles: EEG FEM ISO', ... 525 'InverseMethod', 'gls', ... 526 'InverseMeasure', 'performance', ... 527 'SourceOrient', {{'free'}}, ... 528 'Loose', 0.2, ... 529 'UseDepth', 1, ... 530 'WeightExp', 0.5, ... 531 'WeightLimit', 10, ... 532 'NoiseMethod', 'median', ... 533 'NoiseReg', 0.1, ... 534 'SnrMethod', 'rms', ... 535 'SnrRms', 1e-06, ... 536 'SnrFixed', 3, ... 537 'ComputeKernel', 1, ... 538 'DataTypes', {{'EEG'}})); 539 % Process: Dipole scanning 540 sDipEegFemIso = bst_process('CallProcess', 'process_dipole_scanning', sSrcEegFemIso, [], ... 541 'timewindow', [Latency, Latency], ... 542 'scouts', {}); 543 544 545 %% ===== FORWARD: FEM MEG ISO ===== 546 % Process: Compute head model 547 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 548 'Comment', 'DUNEuro FEM MEG ISO', ... 549 'sourcespace', 1, ... % Cortex surface 550 'meg', 5, ... % DUNEuro FEM 551 'eeg', 1, ... % 552 'ecog', 1, ... % 553 'seeg', 1, ... % 554 'duneuro', struct(... 555 'FemCond', [0.14, 0.33, 1.79, 0.008, 0.43], ... 556 'FemSelect', [1, 1, 1, 0, 0], ... 557 'UseTensor', 0, ... 558 'Isotropic', 1, ... 559 'SrcShrink', 0, ... 560 'SrcForceInGM', 0, ... 561 'FemType', 'fitted', ... 562 'SolverType', 'cg', ... 563 'GeometryAdapted', 0, ... 564 'Tolerance', 1e-08, ... 565 'ElecType', 'normal', ... 566 'MegIntorderadd', 0, ... 567 'MegType', 'physical', ... 568 'SolvSolverType', 'cg', ... 569 'SolvPrecond', 'amg', ... 570 'SolvSmootherType', 'ssor', ... 571 'SolvIntorderadd', 0, ... 572 'DgSmootherType', 'ssor', ... 573 'DgScheme', 'sipg', ... 574 'DgPenalty', 20, ... 575 'DgEdgeNormType', 'houston', ... 576 'DgWeights', 1, ... 577 'DgReduction', 1, ... 578 'SolPostProcess', 1, ... 579 'SolSubstractMean', 0, ... 580 'SolSolverReduction', 1e-10, ... 581 'SrcModel', 'venant', ... 582 'SrcIntorderadd', 0, ... 583 'SrcIntorderadd_lb', 2, ... 584 'SrcNbMoments', 3, ... 585 'SrcRefLen', 20, ... 586 'SrcWeightExp', 1, ... 587 'SrcRelaxFactor', 6, ... 588 'SrcMixedMoments', 1, ... 589 'SrcRestrict', 1, ... 590 'SrcInit', 'closest_vertex', ... 591 'BstSaveTransfer', 0, ... 592 'BstEegTransferFile', 'eeg_transfer.dat', ... 593 'BstMegTransferFile', 'meg_transfer.dat', ... 594 'BstEegLfFile', 'eeg_lf.dat', ... 595 'BstMegLfFile', 'meg_lf.dat', ... 596 'UseIntegrationPoint', 1, ... 597 'EnableCacheMemory', 0, ... 598 'MegPerBlockOfSensor', 0), ... 599 'channelfile', ''); 600 % Process: Compute sources [2018] 601 sSrcMegFemIso = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 602 'output', 1, ... % Kernel only: shared 603 'inverse', struct(... 604 'Comment', 'Dipoles: MEG FEM ISO', ... 605 'InverseMethod', 'gls', ... 606 'InverseMeasure', 'performance', ... 607 'SourceOrient', {{'free'}}, ... 608 'Loose', 0.2, ... 609 'UseDepth', 1, ... 610 'WeightExp', 0.5, ... 611 'WeightLimit', 10, ... 612 'NoiseMethod', 'median', ... 613 'NoiseReg', 0.1, ... 614 'SnrMethod', 'rms', ... 615 'SnrRms', 1e-06, ... 616 'SnrFixed', 3, ... 617 'ComputeKernel', 1, ... 618 'DataTypes', {{'MEG GRAD', 'MEG MAG'}})); 619 % Process: Dipole scanning 620 sDipMegFemIso = bst_process('CallProcess', 'process_dipole_scanning', sSrcMegFemIso, [], ... 621 'timewindow', [Latency, Latency], ... 622 'scouts', {}); 623 624 625 %% ===== FORWARD: BEM EEG ===== 626 % Process: Compute head model 627 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 628 'Comment', 'OpenMEEG BEM EEG', ... 629 'sourcespace', 1, ... % Cortex surface 630 'meg', 1, ... % 631 'eeg', 3, ... % OpenMEEG BEM 632 'ecog', 1, ... % 633 'seeg', 1, ... % 634 'openmeeg', struct(... 635 'BemSelect', [1, 1, 1], ... 636 'BemCond', [1, 0.0125, 1], ... 637 'BemNames', {{'Scalp', 'Skull', 'Brain'}}, ... 638 'BemFiles', {{}}, ... 639 'isAdjoint', 0, ... 640 'isAdaptative', 1, ... 641 'isSplit', 0, ... 642 'SplitLength', 4000), ... 643 'channelfile', ''); 644 % Process: Compute sources [2018] 645 sSrcEegBem = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 646 'output', 1, ... % Kernel only: shared 647 'inverse', struct(... 648 'Comment', 'Dipoles: EEG BEM', ... 649 'InverseMethod', 'gls', ... 650 'InverseMeasure', 'performance', ... 651 'SourceOrient', {{'free'}}, ... 652 'Loose', 0.2, ... 653 'UseDepth', 1, ... 654 'WeightExp', 0.5, ... 655 'WeightLimit', 10, ... 656 'NoiseMethod', 'median', ... 657 'NoiseReg', 0.1, ... 658 'SnrMethod', 'rms', ... 659 'SnrRms', 1e-06, ... 660 'SnrFixed', 3, ... 661 'ComputeKernel', 1, ... 662 'DataTypes', {{'EEG'}})); 663 % Process: Dipole scanning 664 sDipEegBem = bst_process('CallProcess', 'process_dipole_scanning', sSrcEegBem, [], ... 665 'timewindow', [Latency, Latency], ... 666 'scouts', {}); 667 668 669 %% ===== FORWARD: OS MEG ===== 670 % Process: Compute head model 671 bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ... 672 'Comment', 'OS MEG', ... 673 'sourcespace', 1, ... % Cortex surface 674 'meg', 3, ... % Overlapping spheres 675 'eeg', 1, ... % 676 'ecog', 1, ... % 677 'seeg', 1, ... % 678 'channelfile', ''); 679 % Process: Compute sources [2018] 680 sSrcMegOs = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ... 681 'output', 1, ... % Kernel only: shared 682 'inverse', struct(... 683 'Comment', 'Dipoles: MEG OS', ... 684 'InverseMethod', 'gls', ... 685 'InverseMeasure', 'performance', ... 686 'SourceOrient', {{'free'}}, ... 687 'Loose', 0.2, ... 688 'UseDepth', 1, ... 689 'WeightExp', 0.5, ... 690 'WeightLimit', 10, ... 691 'NoiseMethod', 'median', ... 692 'NoiseReg', 0.1, ... 693 'SnrMethod', 'rms', ... 694 'SnrRms', 1e-06, ... 695 'SnrFixed', 3, ... 696 'ComputeKernel', 1, ... 697 'DataTypes', {{'MEG GRAD', 'MEG MAG'}})); 698 % Process: Dipole scanning 699 sDipMegOs = bst_process('CallProcess', 'process_dipole_scanning', sSrcMegOs, [], ... 700 'timewindow', [Latency, Latency], ... 701 'scouts', {}); 702 703 704 %% ===== MERGE DIPOLES ===== 705 % sDipEegFemDti, sDipMegFemDti, sDipEegFemIso, sDipMegFemIso, sDipEegBem, sDipMegOs 706 % FEM DTI: EEG vs MEG 707 % sDipFemDti = dipoles_merge({sDipEegFemDti.FileName, sDipMegFemDti.FileName}); 708 709 710 % Save and display report 711 ReportFile = bst_report('Save', []); 712 bst_report('Open', ReportFile); 713 714 715





Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):


Tutorials/FemMedianNerve (last edited 2024-08-10 06:16:08 by TakfarinasMedani)