Simulations

Author: Francois Tadel

This tutorial addresses everything that can be called "simulation" in Brainstorm: simulation of simple signals, simulation of full source maps from a few scouts, simulation of MEG/EEG recordings from source maps (real or synthetic). Its goal is to group all the information that was previously spread across multiple tutorials and forum posts, and difficult to access. The examples are based on the introduction tutorials, but can easily be transposed to any study.

Simulate signals

You have no data in your database. The processes in the Simulate menu can generate signals and save them as new files in the database. Because you start from nothing, the Process1 list needs to be empty. Click on Run and try the following processes.

Simulate generic signals

This process allows to define a set of signals with Matlab code and save it in a Brainstorm file. With the text box, you can write a free-form piece of code, which will be evaluated with the function eval(). This code must at least define the variable Data [Nsignals x Ntime], based on in the input time vector t (defined automatically by the process from the number of time points and the sampling frequency). The more rows you add to the Data matrix, the more signals you will have in your file. The example below creates three signals. Another example is available in the tutorial Time-frequency.

simulate_generic.gif

Simulate AR signals [TODO]

For evaluating connectivity metrics, it is useful to create signals that have known connectivity properties. Simulating signals with auto-regressive models is a solution illustrated in the tutorial Connectivity. Two processes are available, using different methods.

Process > Simulate AR signals (ARfit)

simulate_ar.gif

Process > Simulate AR signals (random)

simulate_ar_random.gif

Simulate PAC signals

This process generates synthesized data containing cross-frequency phase-amplitude coupling. It is fully documented in the tutorial Phase-amplitude coupling: Method.

simulate_pac.gif

Simulate MEG/EEG from sources

After estimating the sources from real MEG/EEG recordings with the minimum norm method, you can simulate recordings from the model. This is done by multiplying the source time series with the forward model:

EEG_sim [Nelec x Ntime] = Forward_model [Nelec x Nsources] * MN_sources [Nsources x Ntime]

From full source maps

To simulate MEG/EEG recordings from a minimum norm source model: right-click on the source file, then select the menu Model evaluation > Simulate recordings. This can be useful for evaluating the quality of the model. The process is documented in the tutorial Source estimation.

model_eval_menu.gif

From scouts

If you want to limit this simulation to a specific ROI, one easy solution is to set all the sources outside of the ROI to zero, and calculate the MEG/EEG signal exactly as above. This creates a new file which represents the scalp MEG/EEG recordings if only the selected region of the cortex was activated. Double-click on the source file to open it for display, then in the Scout tab select the scout(s) of interested and use the menu Sources > Simulate recordings.

simulate_scout.gif

simulate_scout2.gif

Generate full source maps from scouts

You have the values corresponding to the activity of one or a few dipoles, and you would like to save this a source file in the Brainstorm database. This source file can later be used to simulate MEG/EEG recordings, or be exported as a .nii volume or .gii surface.

Select the scouts time series in the Process1 tab, then run the process Simulate > Full source maps from scouts. Then select the scouts corresponding to the signals in input. Each input file is processed independently and results in the creation of one source file. You need a channel file and a head model file available in the same folder as the input files.

For constrained dipole orientations, you need one signal in the input file for each selected scout.
For unconstrained dipole orientations, you need three signals for each selected scout.

The order of the signals must match the order of the scouts. In the unconstrained case, the signals must be sorted by scout and then by orientation (scout1_orientX, scout1_orientY, scout1_orientZ, scout2_orientX...). In the process options, you cannot change the order of the scouts, however you can re-arrange them if you generate a script.

This process is the first step executed by the process "Simulate > Simulate recordings from scouts", illustrated in the following section.

simulate_sources.gif

Simulate MEG/EEG from synthetic dipoles

You have the values corresponding to the activity of one or a few dipoles, and you would like to simulate the corresponding MEG/EEG surface recordings. The process Simulate > Simulate recordings from scouts does this as a combination of two mechanisms previously described on this page: 1) it generates a source map with zeros everywhere except for a few scouts where it uses instead the signals provided in the input file (see previous section), 2) it multiplies this source map with the forward model available in the folder to simulate MEG/EEG recordings. Noise can be added at each stage.

The input signals can be simulated (as explained in the first sections of this page) or extracted from sources estimated from real MEG/EEG recordings (also explained earlier in this page). You need to provide one signal per scout for constrained sources, or three signals per scout for unconstrained sources. The examples below illustrate three types of source models: cortex surface / constrained orientation, cortex surface / unconstrained orientation, volume grid / unconstrainted orientation.

Surface: Constrained orientation

Simulation of a single dipole located on the cortex surface, oriented along the normal to the cortex surface, pointing outwards.

Surface: Unconstrained orientation

Simulation of three dipoles at the same location on the cortex surface, oriented along three orthogonal axes (X,Y,Z) as defined in the SCS coordinate system.

Volume: Unconstrained orientation

In the first two examples, we showed how to reuse existing recordings already imported and processed. If you don't have any real MEG/EEG recordings to start with, you may not have any channel file a noise covariance to use in your simulation.

In this third example, we will demonstrate how to simulate EEG activity using a template anatomy and a standard 10-10 electrode cap, and from a dipole in a volume grid instead of the cortex surface.

Prepare head model

Create volume scout

Noise covariance

Simulation

Additional documentation

Articles

A simulation article was published by Prof. Debener's neuropsychology lab in Oldenburg. It uses Brainstorm to evaluate and improve the design of ear-EEG devices. All the Matlab/Brainstorm scripts used in this study are available on figshare, together with some videos illustrating the simulations that were done.

Meiser A, Tadel F, Debener S, Bleichner MG
The Sensitivity of Ear-EEG: Evaluating the Source-Sensor Relationship Using Forward Modeling
Brain Topography, Aug 2020

meiser.gif

Forum discussions

Scripting

The following script from the Brainstorm distribution reproduces the third example presented in this tutorial page (Volume: Unconstrained orientation): brainstorm3/toolbox/script/tutorial_simulations.m

1 function tutorial_simulations(tutorial_dir, reports_dir) 2 % TUTORIAL_SIMULATIONS: Script that reproduces the results of the online tutorial "Simulations". 3 % 4 % CORRESPONDING ONLINE TUTORIALS: 5 % https://neuroimage.usc.edu/brainstorm/Tutorials/Simulations 6 % 7 % INPUTS: 8 % - tutorial_dir : Directory where the sample_epilepsy.zip file has been unzipped (NOT USED) 9 % - reports_dir : Directory where to save the execution report (instead of displaying it) 10 11 % @============================================================================= 12 % This function is part of the Brainstorm software: 13 % https://neuroimage.usc.edu/brainstorm 14 % 15 % Copyright (c) University of Southern California & McGill University 16 % This software is distributed under the terms of the GNU General Public License 17 % as published by the Free Software Foundation. Further details on the GPLv3 18 % license can be found at http://www.gnu.org/copyleft/gpl.html. 19 % 20 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 21 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 22 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 23 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 24 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 25 % 26 % For more information type "brainstorm license" at command prompt. 27 % =============================================================================@ 28 % 29 % Author: Francois Tadel, 2020 30 31 32 % ===== FILES TO IMPORT ===== 33 % Output folder for reports 34 if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir) 35 reports_dir = []; 36 end 37 38 % ===== CREATE PROTOCOL ===== 39 % The protocol name has to be a valid folder name (no spaces, no weird characters...) 40 ProtocolName = 'TutorialSimulation'; 41 % Start brainstorm without the GUI 42 if ~brainstorm('status') 43 brainstorm nogui 44 end 45 % Delete existing protocol 46 gui_brainstorm('DeleteProtocol', ProtocolName); 47 % Create new protocol 48 gui_brainstorm('CreateProtocol', ProtocolName, 1, 0); 49 % Start a new report 50 bst_report('Start'); 51 % Create new subject 52 SubjectName = 'Simulation'; 53 [sSubject, iSubject] = db_add_subject(SubjectName, [], 1, 0); 54 % Create new folder 55 FolderName = 'VolumeEEG'; 56 iStudy = db_add_condition(SubjectName, FolderName); 57 58 % ===== EXAMPLE 3: VOLUME/UNCONSTRAINED ===== 59 60 % === 3. FORWARD MODEL === 61 % Get EEG template ICBM152/10-10 (65 channels) 62 eegTemplate = bst_get('EegDefaults', 'icbm152', '10-10 65'); 63 % Set channel file for the simulation folder 64 ChannelFile = db_set_channel(iStudy, eegTemplate.contents.fullpath, 1, 0); 65 % Process: Compute head model 66 bst_process('CallProcess', 'process_headmodel', [], [], ... 67 'Comment', '', ... 68 'sourcespace', 2, ... % MRI volume 69 'volumegrid', struct(... 70 'Method', 'isotropic', ... 71 'nLayers', 17, ... 72 'Reduction', 3, ... 73 'nVerticesInit', 4000, ... 74 'Resolution', 0.005, ... 75 'FileName', []), ... 76 'meg', 1, ... % 77 'eeg', 3, ... % OpenMEEG BEM 78 'ecog', 1, ... % 79 'seeg', 1, ... % 80 'openmeeg', struct(... 81 'BemFiles', {{}}, ... 82 'BemNames', {{'Scalp', 'Skull', 'Brain'}}, ... 83 'BemCond', [1, 0.0125, 1], ... 84 'BemSelect', [1, 1, 1], ... 85 'isAdjoint', 0, ... 86 'isAdaptative', 1, ... 87 'isSplit', 0, ... 88 'SplitLength', 4000), ... 89 'channelfile', ChannelFile); 90 % Set identity noise covariance 91 import_noisecov(iStudy, 'Identity'); 92 93 % === 3. CREATE SCOUT === 94 % Load MRI file 95 sSubject = bst_get('Subject', iSubject); 96 sMri = in_mri_bst(sSubject.Anatomy(1).FileName); 97 % Load head model file 98 sStudy = bst_get('Study', iStudy); 99 HeadModelFile = sStudy.HeadModel(1).FileName; 100 HeadModelMat = in_bst_headmodel(HeadModelFile, 0, 'GridLoc'); 101 102 % Convert grid locations to MNI coordinates 103 GridLocMni = cs_convert(sMri, 'scs', 'mni', HeadModelMat.GridLoc); 104 % Find closest grid point to a target MNI point 105 MniTarget = [-48, -1, -5]; 106 [dist, iVertex] = min(sqrt(sum(bst_bsxfun(@minus, GridLocMni, MniTarget./1000).^2,2))); 107 108 % Create new scout structure 109 DipoleName = 'dip1'; 110 sScout = db_template('scout'); 111 sScout.Vertices = iVertex; 112 sScout.Seed = iVertex; 113 sScout.Color = [0 1 0]; 114 sScout.Label = DipoleName; 115 sScout.Region = 'LT'; 116 % Load cortex surface 117 TessFile = file_fullpath(sSubject.Surface(sSubject.iCortex).FileName); 118 TessMat = in_tess_bst(TessFile); 119 % Get or create a volume atlas and add scout to it 120 AtlasName = sprintf('Volume %d', size(HeadModelMat.GridLoc,1)); 121 iAtlas = find(strcmpi({TessMat.Atlas.Name}, AtlasName)); 122 if isempty(iAtlas) 123 iAtlas = length(TessMat.Atlas) + 1; 124 end 125 TessMat.Atlas(iAtlas).Name = AtlasName; 126 TessMat.Atlas(iAtlas).Scouts = sScout; 127 TessMat.iAtlas = iAtlas; 128 % Update cortex surface 129 bst_save(TessFile, TessMat, 'v7'); 130 131 % === 3. SIMULATE SIGNALS === 132 % Dipole X: (1,0,0) 133 % Process: Simulate generic signals 134 sFilesMatrixX = bst_process('CallProcess', 'process_simulate_matrix', [], [], ... 135 'subjectname', SubjectName, ... 136 'condition', FolderName, ... 137 'samples', 600, ... 138 'srate', 600, ... 139 'matlab', ['Data(1,:) = sin(2*pi*t);' 10 'Data(2,:) = 0 * t;' 10 'Data(3,:) = 0 * t;']); 140 % Process: Set name: Simulated X 141 sFilesMatrixX = bst_process('CallProcess', 'process_set_comment', sFilesMatrixX, [], ... 142 'tag', 'Simulated X', ... 143 'isindex', 1); 144 145 % Dipole Y: (0,1,0) 146 % Process: Simulate generic signals 147 sFilesMatrixY = bst_process('CallProcess', 'process_simulate_matrix', [], [], ... 148 'subjectname', SubjectName, ... 149 'condition', FolderName, ... 150 'samples', 600, ... 151 'srate', 600, ... 152 'matlab', ['Data(1,:) = 0 * t;' 10 'Data(2,:) = sin(2*pi*t);' 10 'Data(3,:) = 0 * t;']); 153 % Process: Set name: Simulated X 154 sFilesMatrixY = bst_process('CallProcess', 'process_set_comment', sFilesMatrixY, [], ... 155 'tag', 'Simulated Y', ... 156 'isindex', 1); 157 158 % Dipole Z: (0,0,1) 159 % Process: Simulate generic signals 160 sFilesMatrixZ = bst_process('CallProcess', 'process_simulate_matrix', [], [], ... 161 'subjectname', SubjectName, ... 162 'condition', FolderName, ... 163 'samples', 600, ... 164 'srate', 600, ... 165 'matlab', ['Data(1,:) = 0 * t;' 10 'Data(2,:) = 0 * t;' 10 'Data(3,:) = sin(2*pi*t);']); 166 % Process: Set name: Simulated X 167 sFilesMatrixZ = bst_process('CallProcess', 'process_set_comment', sFilesMatrixZ, [], ... 168 'tag', 'Simulated Z', ... 169 'isindex', 1); 170 171 % === 3. SIMULATE EEG === 172 % Process: Simulate recordings from scouts 173 sFilesRec = bst_process('CallProcess', 'process_simulate_recordings', [sFilesMatrixX, sFilesMatrixY, sFilesMatrixZ], [], ... 174 'scouts', {AtlasName, {DipoleName}}, ... 175 'savedata', 1, ... 176 'savesources', 1, ... 177 'isnoise', 1, ... 178 'noise1', 0.01, ... 179 'noise2', 0.01); 180 181 % === 3. SCREEN CAPTURES === 182 % Process: Snapshot: Sensors/MRI registration 183 bst_process('CallProcess', 'process_snapshot', sFilesRec(1), [], ... 184 'target', 1, ... % Sensors/MRI registration 185 'modality', 4, ... % EEG 186 'orient', 1, ... % left 187 'Comment', '10-10 electrode cap, 65 electrodes'); 188 % Process: Snapshot: Recordings time series 189 bst_process('CallProcess', 'process_snapshot', sFilesRec, [], ... 190 'target', 5, ... % Recordings time series 191 'modality', 4, ... % EEG 192 'time', 0.223, ... 193 'Comment', ''); 194 % Process: Snapshot: Recordings topography (one time) 195 bst_process('CallProcess', 'process_snapshot', sFilesRec, [], ... 196 'target', 6, ... % Recordings topography (one time) 197 'modality', 4, ... % EEG 198 'time', 0.223, ... 199 'Comment', ''); 200 201 202 % ===== EXAMPLE 4: SINGLE DIPOLES ===== 203 % Create new folder 204 FolderName = 'Dipoles'; 205 iStudy = db_add_condition(SubjectName, FolderName); 206 % Get EEG template ICBM152/10-10 (65 channels) 207 eegTemplate = bst_get('EegDefaults', 'icbm152', '10-10 65'); 208 % Set channel file for the simulation folder 209 db_set_channel(iStudy, eegTemplate.contents.fullpath, 1, 0); 210 % Set identity noise covariance 211 import_noisecov(iStudy, 'Identity'); 212 213 % Process: Simulate generic signals 214 sFileSim = bst_process('CallProcess', 'process_simulate_matrix', [], [], ... 215 'subjectname', SubjectName, ... 216 'condition', FolderName, ... 217 'samples', 600, ... 218 'srate', 600, ... 219 'matlab', ['Data(1,:) = sin(2*pi*t);' 10 'Data(2,:) = sin(2*pi*t + pi/4);']); 220 % Process: Set name: Simulated bilateral 221 sFileSim = bst_process('CallProcess', 'process_set_comment', sFileSim, [], ... 222 'tag', 'Simulated bilateral', ... 223 'isindex', 0); 224 % Process: Simulate recordings from dipoles 225 sFileEeg = bst_process('CallProcess', 'process_simulate_dipoles', sFileSim, [], ... 226 'dipoles', ['-48, -2, -4, 0.2, 1, 0' 10 '48, -2, -4, -0.2, 1, 0'], ... 227 'cs', 'mni', ... % MNI 228 'meg', {''}, ... % 229 'eeg', {'openmeeg'}, ... % OpenMEEG BEM 230 'ecog', {''}, ... % 231 'seeg', {''}, ... % 232 'openmeeg', struct(... 233 'BemFiles', {{}}, ... 234 'BemNames', {{'Scalp', 'Skull', 'Brain'}}, ... 235 'BemCond', [1, 0.0125, 1], ... 236 'BemSelect', [1, 1, 1], ... 237 'isAdjoint', 0, ... 238 'isAdaptative', 1, ... 239 'isSplit', 0, ... 240 'SplitLength', 4000), ... 241 'isnoise', 1, ... 242 'noise1', 0.1, ... 243 'noise2', 0.5, ... 244 'savedip', 1, ... 245 'savedata', 1); 246 247 % Process: Snapshot: Recordings time series 248 bst_process('CallProcess', 'process_snapshot', sFileSim, [], ... 249 'target', 5, ... % Recordings time series 250 'time', 0); 251 % Process: Snapshot: Recordings time series 252 bst_process('CallProcess', 'process_snapshot', sFileEeg, [], ... 253 'target', 5, ... % Recordings time series 254 'modality', 4, ... % EEG 255 'time', 0); 256 % Process: Snapshot: Recordings topography (contact sheet) 257 bst_process('CallProcess', 'process_snapshot', sFileEeg, [], ... 258 'target', 7, ... % Recordings topography (contact sheet) 259 'modality', 4, ... % EEG 260 'contact_time', [0, 0.998], ... 261 'contact_nimage', 12); 262 263 % Process: Select files using search query 264 sFileDip = bst_process('CallProcess', 'process_select_search', [], [], ... 265 'search', '([type EQUALS "Dipoles"])'); 266 % Display dipoles 267 hFig = view_dipoles(sFileDip.FileName, 'mri3d'); 268 figure_3d('SetStandardView', hFig, 'top'); 269 bst_report('Snapshot', hFig, sFileDip.FileName, 'Two dipoles: MNI(-48,-2,-4) and MNI(+48,-2,-4)', [200, 200, 640, 400]); 270 figure_3d('SetStandardView', hFig, 'left'); 271 bst_report('Snapshot', hFig, sFileDip.FileName, 'Two dipoles: MNI(-48,-2,-4) and MNI(+48,-2,-4)', [200, 200, 640, 400]); 272 close(hFig); 273 274 275 % ===== SAVE REPORT ===== 276 % Save and display report 277 ReportFile = bst_report('Save', []); 278 if ~isempty(reports_dir) && ~isempty(ReportFile) 279 bst_report('Export', ReportFile, reports_dir); 280 else 281 bst_report('Open', ReportFile); 282 end 283 284 disp([10 'BST> tutorial_simulations: Done.' 10]); 285





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


Tutorials/Simulations (last edited 2020-09-29 17:59:35 by FrancoisTadel)