Corticomuscular coherence (CTF MEG)

[TUTORIAL UNDER DEVELOPMENT: NOT READY FOR PUBLIC USE]

Authors: Raymundo Cassani, Francois Tadel & Sylvain Baillet.

Corticomuscular coherence measures the degree of similarity between electrophysiological signals (MEG, EEG, ECoG sensor traces or source time series, especially over the contralateral motor cortex) and the EMG signals recorded from muscle activity during voluntary movements. Cortical-muscular signals similarities are conceived as due mainly to the descending communication along corticospinal pathways between primary motor cortex (M1) and the muscles attached to the moving limb(s). For consistency purposes, the present tutorial replicates, with Brainstorm tools, the processing pipeline "Analysis of corticomuscular coherence" of the FieldTrip toolbox.

Dataset description

The dataset is identical to that of the FieldTrip tutorial: Analysis of corticomuscular coherence:

Corticomuscular coherence:

Download and installation

Pre-requisites

Download the dataset

Brainstorm

Importing anatomy

MEG and EMG recordings

MEG-MRI coregistration

Reviewing

Event markers

The colored dots above the data time series indicate event markers (or triggers) saved with this dataset. The trial onset information of the left-wrist and right-wrist trials is saved in an auxiliary channel of the raw data named Stim. To add these markers, these events need to be decoded as follows:

Pre-processing

In this tutorial, we will analyze only the Left trials (left-wrist extensions). In the following sections, we will process only the first 330 s of the recordings, where the left-wrist trials were performed.

Power line artifacts

EMG: Filter and rectify

Two typical pre-processing steps for EMG consist in high-pass filtering and rectifying.

Stereotypical artifacts such eye blinks and heartbeats can be identified from their respective characteristic spatial distributions. Their contamination of MEG signals can then be attenuated specifically using Signal-Space Projections (SSPs). For more details, consult the dedicated tutorials about the detection and removal of artifacts with SSP. The present tutorial dataset features an EOG channel but no ECG. We will perform only the removal of eye blinks.

Detection of "bad" data segments

Here we will use the automatic detection of artifacts to identify data segments contaminated by e.g., large eye and head movements and muscle contractions.

Epoching

We are finished with the pre-processing of the EMG and MEG recordings. We will now extract and import specific data segments of interest into the Brainstorm database for further derivations. As mentioned previously, we will focus on the Left category of events (left wrist movements). To follow the same pipeline as the FieldTrip tutorial: we will consider 8 seconds of recordings after each trigger (out of the 10s of each trial), and split them in epochs of 1 second. In addition DC offset is removed, only for MEG signals.

pro_import.png

pro_remove_dc.png

Comparison with FieldTrip

The figures below represent the EMG and MRC21 channels (sensor over the left motor-cortex) from the epoch #1.1, in Brainstorm (left) and in the FieldTrip tutorial (right).

bst_ft_trial1.png

Coherence: EMG x MEG

Let's compute the magnitude square coherence (MSC) between the left EMG and the MEG channels.

coh_meg_emgleft.png

coh_meg_emgleft2.png

Source estimation

MRI segmentation

In order to estimate the brain sources for these MEG recordings, we first need to reconstruct the cortex surface from the T1 MRI imported at the beginning of this tutorial. For this puropose, we decided to use CAT12 because it is fast (30-60min) and fully integrated with Brainstorm as a plugin.

Head models

We will perform source modeling using a distributed model approach for two different source spaces: the cortex surface and the entire MRI volume. The forward model, labelled head model in Brainstorm, accounts for how neural electrical currents produce magnetic fields captured by sensors outside the head, considering head tissues electromagnetic properties and geometry, independently of actual empirical measurements (more information). As the head model depends on the source space, a distinct head model is required for the surface and volume source spaces: we will compute them both now.

Surface

Volume

Noise covariance

The recommendation for MEG is to extract basic noise statistics from empty-room recordings. When not available, resting-state data can be used as proxies for MEG noise covariance. We will use a segment of the MEG recordings, away from the task or major artifacts: 18s-29s.

Inverse models

We will now compute three inverse models, with different source spaces: cortex surface with constrained (normal to the cortex) dipole orientations, cortex surface with unconstrained orientation, and MRI volume (more information).

Surface

pro_sources_srfc.png

pro_sources_srfu.png

Volume

Coherence: EMG x Sources

We can now compute the coherence between the EMG signal and the brain source time series, for each of the source models. Let's start with the surface/constrained model.

pro_coh_srf.png

pro_coh_srf2.png

Surface

Double-click the 1xN connectivity files for the two (surface) source space to show the results on the cortex. If you are not familiar with the options in the cortex figures, check Display: Cortex surface. Find the location and frequency with the highest coherence value.

Volume

Method

For constrained sources, each vertex in the source grid is associated with ONE time series, as such, when coherence is computed with the EMG signal (also one time series), the result is ONE coherence spectrum per vertex. In other words, for each frequency bin, there is a coherence brain map.

https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_1xn_coh_constr.png

In the case of unconstrained or volume sources, each vertex in the grid is associated with THREE time series, each one corresponding to the X, Y and Z directions. Thus, when coherence is computed with the EMG signal (one time series), there are THREE coherence spectra. To be represented on the cortex, these three values need to be flattened into one, resulting in one coherence spectrum per vertex. For performing this dimension reduction, we decided to take the maximum across three directions, for each frequency bin for each vertex.

https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_1xn_coh_unconstr.png

An alternative approach in the literature, to address the 3-dimensional nature of the unconstrained sources, consists in flattening the vertex X, Y and Z time series before the coherence computation; resulting in a similar case as the constrained sources. Common methods for this flattening include: PCA (only first component is kept), and (Euclidean) norm. This flattening of the time series can be performed in Brainstorm with the process: Sources > Unconstrained to flat map.

https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_1xn_coh_flattened.png

  • Flattened sources are saved as full rather than recordings+kernel.
  • We have tested this flattening approach with simulations and found detrimental effects on the expected results.

Coherence: EMG x Scouts

So far, we have computed coherence at the source level, thus, a coherence spectrum is computed for each of the 15002 source points. This large dimension hinders later analysis of the results. Therefore, the strategy is to reduce the dimensionality of the source space by using a surface- or volume-parcellation scheme, in Brainstorm jargon this is an atlas that is made of scouts. See the scout tutorial for detail information on atlases and scouts in Brainstorm.

Under this approach, instead of providing one result (coherence spectrum) per source vertex, one result is computed for each scout. When computing coherence (or other connectivity metrics) at the scout level, it is necessary to provide two parameters that define how the data is aggregated per scout: The scout function (mean is often used), and when the within-scout aggregation takes place (before or after the coherence computation).

https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_1xn_coh_sct_bef.png

https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_1xn_coh_sct_aft.png

Let's here compute the coherence using scouts, using mean as scout function alongside with the Before option. We will use the Schaefer 100 parcellation atlas on the results from constrained sources.

Open the Pipeline editor:

pro_coh_srfc_bef_sct.png

pro_coh_srfc_bef_sct2.png

Open the output file by double-clicking on it. This time the coherence spectra are not displayed on the cortex, but they are plotted for each scout. Moreover, 1xN connectivity file can be shown as image.

res_coh_srfc_bef_sct.png

res_coh_srfc_bef_sct2.png

Note that for 14.65 Hz, the highest two peaks correspond to the SomMotA_4 R and SomMotA_2 R scouts, both located over the right primary motor cortex.

The choice of the optimal parcellation scheme for the source space is not easy, this is still an active field of research. The only recommendation we can give at the moment is to use regions of interest that are small and homogenous in size.

Advanced

Additional documentation

Articles

Tutorials

Scripting

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

1 function tutorial_coherence(tutorial_dir, reports_dir) 2 % TUTORIAL_COHERENCE: Script that runs the Brainstorm corticomuscular coherence tutorial 3 % https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence 4 % 5 % INPUTS: 6 % - tutorial_dir : Directory where the SubjectCMC.zip file has been unzipped 7 % - reports_dir : Directory where to save the execution report (instead of displaying it) 8 9 % @============================================================================= 10 % This function is part of the Brainstorm software: 11 % https://neuroimage.usc.edu/brainstorm 12 % 13 % Copyright (c) University of Southern California & McGill University 14 % This software is distributed under the terms of the GNU General Public License 15 % as published by the Free Software Foundation. Further details on the GPLv3 16 % license can be found at http://www.gnu.org/copyleft/gpl.html. 17 % 18 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 19 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 20 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 21 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 22 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 23 % 24 % For more information type "brainstorm license" at command prompt. 25 % =============================================================================@ 26 % 27 % Author: Raymundo Cassani & Francois Tadel, 2022 28 29 % Output folder for reports 30 if (nargin < 2) || isempty(reports_dir) || ~isfolder(reports_dir) 31 reports_dir = []; 32 end 33 % You have to specify the folder in which the tutorial dataset is unzipped 34 if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir) 35 error('The first argument must be the full path to the dataset folder.'); 36 end 37 38 % Protocol name 39 ProtocolName = 'TutorialCMC'; 40 % Subject name 41 SubjectName = 'Subject01'; 42 % Channel selection 43 emg_channel = 'EMGlft'; % Name of EMG channel 44 meg_sensor = 'MRC21'; % MEG sensor over the left motor-cortex (MRC21) 45 % Coherence parameters 46 cohmeasure = 'mscohere'; % Magnitude-squared Coherence |C|^2 = |Cxy|^2/(Cxx*Cyy) 47 win_length = 0.5; % 500ms 48 overlap = 50; % 50% 49 maxfreq = 80; % 80Hz 50 51 % Build the path of the files to import 52 MriFilePath = fullfile(tutorial_dir, 'SubjectCMC', 'SubjectCMC.mri'); 53 MegFilePath = fullfile(tutorial_dir, 'SubjectCMC', 'SubjectCMC.ds'); 54 % Check if the folder contains the required files 55 if ~file_exist(MriFilePath) || ~file_exist(MegFilePath) 56 error(['The folder ' tutorial_dir ' does not contain the folder from the file SubjectCMC.zip.']); 57 end 58 59 60 %% ===== CREATE PROTOCOL ===== 61 % Start brainstorm without the GUI 62 if ~brainstorm('status') 63 brainstorm nogui 64 end 65 % Delete existing protocol 66 gui_brainstorm('DeleteProtocol', ProtocolName); 67 % Create new protocol 68 gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); 69 % Start a new report 70 bst_report('Start'); 71 % Reset colormaps 72 bst_colormaps('RestoreDefaults', 'meg'); 73 % Hide scouts 74 panel_scout('SetScoutShowSelection', 'none'); 75 76 77 %% ===== IMPORT ANATOMY ===== 78 % Process: Import MRI 79 bst_process('CallProcess', 'process_import_mri', [], [], ... 80 'subjectname', SubjectName, ... 81 'mrifile', {MriFilePath, 'ALL'}, ... 82 'nas', [0, 0, 0], ... 83 'lpa', [0, 0, 0], ... 84 'rpa', [0, 0, 0], ... 85 'ac', [0, 0, 0], ... 86 'pc', [0, 0, 0], ... 87 'ih', [0, 0, 0]); 88 % Process: Generate head surface 89 bst_process('CallProcess', 'process_generate_head', [], [], ... 90 'subjectname', SubjectName, ... 91 'nvertices', 10000, ... 92 'erodefactor', 0, ... 93 'fillfactor', 2); 94 95 96 %% ===== LINK TO RAW FILE AND DISPLAY REGISTRATION ===== 97 % Process: Create link to raw files 98 sFileRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 99 'subjectname', SubjectName, ... 100 'datafile', {MegFilePath, 'CTF'}, ... 101 'channelalign', 1); 102 % Process: Snapshot: Sensors/MRI registration 103 bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ... 104 'target', 1, ... % Sensors/MRI registration 105 'modality', 1, ... % MEG (All) 106 'orient', 1, ... % left 107 'Comment', 'MEG/MRI Registration'); 108 % Process: Convert to continuous (CTF): Continuous 109 bst_process('CallProcess', 'process_ctf_convert', sFileRaw, [], ... 110 'rectype', 2); % Continuous 111 112 113 %% ===== EVENT MARKERS ===== 114 % Process: Read from channel 115 sFileRaw = bst_process('CallProcess', 'process_evt_read', sFileRaw, [], ... 116 'stimchan', 'Stim', ... 117 'trackmode', 1, ... % Value: detect the changes of channel value 118 'zero', 0, ... 119 'min_duration', 12); 120 % Load all Event group labels 121 DataMat = in_bst_data(sFileRaw.FileName, 'F'); 122 eventList = {DataMat.F.events.label}; 123 % Labels for Event groups to keep 124 eventKeep = cellfun(@(c)sprintf('U%d',c), num2cell([1:6, 8:25]), 'UniformOutput', 0); 125 % Find useless Events 126 eventDelete = setdiff(eventList, eventKeep); 127 % Process: Delete events 128 sFileRaw = bst_process('CallProcess', 'process_evt_delete', sFileRaw, [], ... 129 'eventname', strjoin(eventDelete, ', ')); 130 % Process: Merge events 131 sFileRaw = bst_process('CallProcess', 'process_evt_merge', sFileRaw, [], ... 132 'evtnames', strjoin(eventKeep, ', '), ... 133 'newname', 'Left'); 134 135 136 %% ===== REMOVAL OF POWER LINE ARTIFACTS ===== 137 % Process: Notch filter: 50Hz 100Hz 150Hz 138 sFileRawNotch = bst_process('CallProcess', 'process_notch', sFileRaw, [], ... 139 'freqlist', [50, 100, 150], ... 140 'sensortypes', 'MEG, EMG', ... 141 'read_all', 1); 142 % Process: Power spectrum density (Welch) 143 sFilesPsd = bst_process('CallProcess', 'process_psd', [sFileRaw, sFileRawNotch], [], ... 144 'timewindow', [0 330], ... 145 'win_length', 10, ... 146 'win_overlap', 50, ... 147 'clusters', {}, ... 148 'sensortypes', 'MEG, EMG', ... 149 'edit', struct(... 150 'Comment', 'Power', ... 151 'TimeBands', [], ... 152 'Freqs', [], ... 153 'ClusterFuncTime', 'none', ... 154 'Measure', 'power', ... 155 'Output', 'all', ... 156 'SaveKernel', 0)); 157 % Process: Snapshot: Frequency spectrum 158 bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ... 159 'target', 10, ... % Frequency spectrum 160 'Comment', 'Power spectrum density'); 161 162 163 %% ===== EMG PRE-PROCESSING ===== 164 % Process: High-pass:10Hz 165 sFileRawNotchHigh = bst_process('CallProcess', 'process_bandpass', sFileRawNotch, [], ... 166 'sensortypes', 'EMG', ... 167 'highpass', 10, ... 168 'lowpass', 0, ... 169 'tranband', 0, ... 170 'attenuation', 'strict', ... % 60dB 171 'ver', '2019', ... % 2019 172 'mirror', 0, ... 173 'read_all', 0); 174 % Process: Absolute values 175 sFileRawNotchHighAbs = bst_process('CallProcess', 'process_absolute', sFileRawNotchHigh, [], ... 176 'sensortypes', 'EMG'); 177 % Process: Delete folders 178 bst_process('CallProcess', 'process_delete', [sFileRawNotch, sFileRawNotchHigh], [], ... 179 'target', 2); % Delete folders 180 181 182 %% ===== MEG PRE-PROCESSING ===== 183 % Process: Detect eye blinks 184 bst_process('CallProcess', 'process_evt_detect_eog', sFileRawNotchHighAbs, [], ... 185 'channelname', 'EOG', ... 186 'timewindow', [0 330], ... 187 'eventname', 'blink'); 188 % Process: SSP EOG: blink 189 bst_process('CallProcess', 'process_ssp_eog', sFileRawNotchHighAbs, [], ... 190 'eventname', 'blink', ... 191 'sensortypes', 'MEG', ... 192 'usessp', 1, ... 193 'select', 1); 194 % Process: Snapshot: SSP projectors 195 bst_process('CallProcess', 'process_snapshot', sFileRawNotchHighAbs, [], ... 196 'target', 2, ... % SSP projectors 197 'Comment', 'SSP projectors'); 198 199 % Process: Detect other artifacts 200 bst_process('CallProcess', 'process_evt_detect_badsegment', sFileRawNotchHighAbs, [], ... 201 'timewindow', [0 330], ... 202 'sensortypes', 'MEG', ... 203 'threshold', 3, ... % 3 204 'isLowFreq', 1, ... 205 'isHighFreq', 1); 206 % Process: Rename event (1-7Hz > bad_1-7Hz) 207 bst_process('CallProcess', 'process_evt_rename', sFileRawNotchHighAbs, [], ... 208 'src', '1-7Hz', ... 209 'dest', 'bad_1-7Hz'); 210 % Process: Rename event (40-240Hz > bad_40-240Hz) 211 bst_process('CallProcess', 'process_evt_rename', sFileRawNotchHighAbs, [], ... 212 'src', '40-240Hz', ... 213 'dest', 'bad_40-240Hz'); 214 215 216 %% ===== IMPORTING DATA EPOCHS ===== 217 % Process: Import MEG/EEG: Events 218 sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFileRawNotchHighAbs, [], ... 219 'subjectname', SubjectName, ... 220 'condition', '', ... 221 'eventname', 'Left', ... 222 'timewindow', [0, 330], ... 223 'epochtime', [0, 7.9992], ... 224 'split', 1, ... 225 'createcond', 0, ... 226 'ignoreshort', 1, ... 227 'usectfcomp', 1, ... 228 'usessp', 1, ... 229 'freq', [], ... 230 'baseline', 'all', ... 231 'blsensortypes', 'MEG'); 232 % Process: Uniform epoch time 233 sFilesEpochs = bst_process('CallProcess', 'process_stdtime', sFilesEpochs, [], ... 234 'method', 'spline', ... % spline 235 'overwrite', 1); 236 % Process: Snapshot: Recordings time series 237 bst_process('CallProcess', 'process_snapshot', sFilesEpochs(1), [], ... 238 'type', 'data', ... % Recordings time series 239 'modality', 1, ... % MEG (All) 240 'rowname', meg_sensor, ... 241 'Comment', ['Epoch #1: ' meg_sensor]); 242 % Process: Snapshot: Recordings time series 243 bst_process('CallProcess', 'process_snapshot', sFilesEpochs(1), [], ... 244 'type', 'data', ... % Recordings time series 245 'modality', 10, ... % EMG 246 'rowname', '', ... 247 'Comment', 'Epoch #1: EMG'); 248 249 250 %% ===== COHERENCE: EMG x MEG ===== 251 % Process: Coherence NxN [2023] 252 sFileCoh1N = bst_process('CallProcess', 'process_cohere1', {sFilesEpochs.FileName}, [], ... 253 'timewindow', [], ... 254 'src_channel', emg_channel, ... 255 'dest_sensors', 'MEG', ... 256 'includebad', 0, ... 257 'removeevoked', 0, ... 258 'cohmeasure', cohmeasure, ... % Magnitude-squared coherence 259 'tfmeasure', 'stft', ... % Fourier transform 260 'tfedit', struct(... 261 'Comment', 'Complex', ... 262 'TimeBands', [], ... 263 'Freqs', [], ... 264 'StftWinLen', win_length, ... 265 'StftWinOvr', overlap, ... 266 'StftFrqMax', maxfreq, ... 267 'ClusterFuncTime', 'none', ... 268 'Measure', 'none', ... 269 'Output', 'all', ... 270 'SaveKernel', 0), ... 271 'timeres', 'none', ... % None 272 'avgwinlength', 1, ... 273 'avgwinoverlap', 0, ... 274 'outputmode', 'avg'); % across combined files/epochs 275 % Process: Add tag 276 sFileCoh1N = bst_process('CallProcess', 'process_add_tag', sFileCoh1N, [], ... 277 'tag', 'MEG sensors', ... 278 'output', 1); % Add to file name 279 % Set selected sensors 280 bst_figures('SetSelectedRows', {meg_sensor}); 281 % Process: Snapshot: Frequency spectrum 282 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 283 'target', 10, ... % Frequency spectrum 284 'freq', 17.58, ... 285 'Comment', ['MSC ' emg_channel ' x MEG']); 286 % Process: Snapshot: Frequency spectrum 287 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 288 'target', 10, ... % Frequency spectrum 289 'freq', 17.58, ... 290 'rowname', meg_sensor, ... 291 'Comment', ['MSC ' emg_channel ' x ', meg_sensor]); 292 % Process: Snapshot: Recordings topography (one time) 293 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 294 'type', 'topo', ... % Recordings topography (one time) 295 'modality', 1, ... % MEG (All) 296 'freq', 17.58, ... 297 'Comment', ['2D sensor cap MSC ' emg_channel ' x MEG']); 298 299 300 %% ===== COHERENCE: EMG x MEG (FREQUENCY BAND) ===== 301 % Process: Group in time or frequency bands 302 sFileCoh1NBand = bst_process('CallProcess', 'process_tf_bands', sFileCoh1N, [], ... 303 'isfreqbands', 1, ... 304 'freqbands', {'cmc_band', '15, 20', 'mean'}, ... 305 'istimebands', 0, ... 306 'timebands', '', ... 307 'overwrite', 0); 308 % Process: Snapshot: Recordings topography (one time) 309 bst_process('CallProcess', 'process_snapshot', sFileCoh1NBand, [], ... 310 'type', 'topo', ... % Recordings topography (one time) 311 'modality', 1, ... % MEG (All) 312 'Comment', ['2D topography MSC 15-20Hz ' emg_channel ' x MEG']); 313 314 315 %% ===== SOURCE ESTIMATION ===== 316 % Process: Segment MRI with CAT12 317 bst_process('CallProcess', 'process_segment_cat12', [], [], ... 318 'subjectname', SubjectName, ... 319 'nvertices', 15000, ... 320 'tpmnii', {'', 'Nifti1'}, ... 321 'sphreg', 1, ... 322 'vol', 1, ... 323 'extramaps', 0, ... 324 'cerebellum', 0); 325 % Process: Compute covariance (noise or data) 326 bst_process('CallProcess', 'process_noisecov', sFileRawNotchHighAbs, [], ... 327 'baseline', [18, 29], ... 328 'datatimewindow', [], ... 329 'sensortypes', 'MEG', ... 330 'target', 1, ... % Noise covariance (covariance over baseline time window) 331 'dcoffset', 1, ... % Block by block, to avoid effects of slow shifts in data 332 'identity', 0, ... 333 'copycond', 1, ... 334 'copysubj', 0, ... 335 'copymatch', 0, ... 336 'replacefile', 1); % Replace 337 338 % Process: Compute head model (surface) 339 bst_process('CallProcess', 'process_headmodel', sFilesEpochs(1).FileName, [], ... 340 'Comment', 'Overlapping spheres (surface)', ... 341 'sourcespace', 1, ... % Cortex 342 'meg', 3); % Overlapping spheres 343 % Process: Compute sources [2018] 344 bst_process('CallProcess', 'process_inverse_2018', sFilesEpochs(1).FileName, [], ... 345 'output', 1, ... % Kernel only: shared 346 'inverse', struct(... 347 'Comment', 'MN: MEG (surface)', ... 348 'InverseMethod', 'minnorm', ... 349 'InverseMeasure', 'amplitude', ... 350 'SourceOrient', {{'fixed'}}, ... 351 'Loose', 0.2, ... 352 'UseDepth', 1, ... 353 'WeightExp', 0.5, ... 354 'WeightLimit', 10, ... 355 'NoiseMethod', 'reg', ... 356 'NoiseReg', 0.1, ... 357 'SnrMethod', 'fixed', ... 358 'SnrRms', 1e-06, ... 359 'SnrFixed', 3, ... 360 'ComputeKernel', 1, ... 361 'DataTypes', {{'MEG'}})); 362 % Process: Compute sources [2018] 363 bst_process('CallProcess', 'process_inverse_2018', sFilesEpochs(1).FileName, [], ... 364 'output', 1, ... % Kernel only: shared 365 'inverse', struct(... 366 'Comment', 'MN: MEG (surface)', ... 367 'InverseMethod', 'minnorm', ... 368 'InverseMeasure', 'amplitude', ... 369 'SourceOrient', {{'free'}}, ... 370 'Loose', 0.2, ... 371 'UseDepth', 1, ... 372 'WeightExp', 0.5, ... 373 'WeightLimit', 10, ... 374 'NoiseMethod', 'reg', ... 375 'NoiseReg', 0.1, ... 376 'SnrMethod', 'fixed', ... 377 'SnrRms', 1e-06, ... 378 'SnrFixed', 3, ... 379 'ComputeKernel', 1, ... 380 'DataTypes', {{'MEG'}})); 381 382 % Process: Compute head model (volume) 383 bst_process('CallProcess', 'process_headmodel', sFilesEpochs(1).FileName, [], ... 384 'Comment', 'Overlapping spheres (volume)', ... 385 'sourcespace', 2, ... % MRI volume 386 'volumegrid', struct(... 387 'Method', 'isotropic', ... 388 'nLayers', 17, ... 389 'Reduction', 3, ... 390 'nVerticesInit', 4000, ... 391 'Resolution', 0.005, ... 392 'FileName', []), ... 393 'meg', 3 ); % Overlapping spheres 394 % Process: Compute sources [2018] 395 bst_process('CallProcess', 'process_inverse_2018', sFilesEpochs(1).FileName, [], ... 396 'output', 1, ... % Kernel only: shared 397 'inverse', struct(... 398 'Comment', 'MN: MEG (volume)', ... 399 'InverseMethod', 'minnorm', ... 400 'InverseMeasure', 'amplitude', ... 401 'SourceOrient', {{'free'}}, ... 402 'Loose', 0.2, ... 403 'UseDepth', 1, ... 404 'WeightExp', 0.5, ... 405 'WeightLimit', 10, ... 406 'NoiseMethod', 'reg', ... 407 'NoiseReg', 0.1, ... 408 'SnrMethod', 'fixed', ... 409 'SnrRms', 1e-06, ... 410 'SnrFixed', 3, ... 411 'ComputeKernel', 1, ... 412 'DataTypes', {{'MEG'}})); 413 414 415 %% ===== COHERENCE: EMG x SOURCES ===== 416 % Process: Select data files 417 sFilesRecEmg = bst_process('CallProcess', 'process_select_files_data', [], [], ... 418 'subjectname', SubjectName, ... 419 'condition', '', ... 420 'tag', 'Left', ... 421 'includebad', 0, ... 422 'includeintra', 0, ... 423 'includecommon', 0); 424 % Coherence between EMG signal and sources (for different source types) 425 sourceTypes = {'(surface)(Constr)', '(surface)(Unconstr)', '(volume)(Unconstr)'}; 426 for ix = 1 : length(sourceTypes) 427 sourceType = sourceTypes{ix}; 428 % Process: Select results files 429 sFilesResMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ... 430 'subjectname', SubjectName, ... 431 'condition', '', ... 432 'tag', sourceType, ... 433 'includebad', 0, ... 434 'includeintra', 0, ... 435 'includecommon', 0); 436 % Process: Coherence NxN [2023] 437 sFileCoh1N = bst_process('CallProcess', 'process_cohere2', sFilesRecEmg, sFilesResMeg, ... 438 'timewindow', [], ... 439 'src_channel', emg_channel, ... 440 'dest_scouts', {}, ... 441 'flatten', 1, ... 442 'scouttime', 'after', ... % after connectivity metric 443 'scoutfunc', 'mean', ... % Mean 444 'scoutfuncaft', 'mean', ... % Mean 445 'pcaedit', struct(... 446 'Method', 'pcaa', ... 447 'Baseline', [], ... 448 'DataTimeWindow', [], ... 449 'RemoveDcOffset', 'file'), ... 450 'removeevoked', 0, ... 451 'cohmeasure', cohmeasure, ... % Magnitude-squared coherence 452 'tfmeasure', 'stft', ... % Fourier transform 453 'tfedit', struct(... 454 'Comment', 'Complex', ... 455 'TimeBands', [], ... 456 'Freqs', [], ... 457 'StftWinLen', win_length, ... 458 'StftWinOvr', overlap, ... 459 'StftFrqMax', maxfreq, ... 460 'ClusterFuncTime', 'none', ... 461 'Measure', 'none', ... 462 'Output', 'all', ... 463 'SaveKernel', 0), ... 464 'timeres', 'none', ... % None 465 'avgwinlength', 1, ... 466 'avgwinoverlap', 50, ... 467 'outputmode', 'avg'); % across combined files/epochs 468 % Process: Add tag 469 sFileCoh1N = bst_process('CallProcess', 'process_add_tag', sFileCoh1N, [], ... 470 'tag', sourceType, ... 471 'output', 1); % Add to file name 472 % View surface 473 if ~isempty(strfind(sourceType, 'surface')) 474 % Process: Snapshot: Sources (one time) 475 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 476 'type', 'sources', ... % Sources (one time) 477 'orient', 3, ... % top 478 'threshold', 0, ... 479 'surfsmooth', 30, ... 480 'freq', 14.65, ... 481 'Comment', ['MSC 14.65Hz ', sourceType]); 482 % View volume 483 elseif ~isempty(strfind(sourceType, 'volume')) 484 % Process: Snapshot: Sources (MRI Viewer) 485 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 486 'type', 'mriviewer', ... % MRI viewer 487 'threshold', 0, ... 488 'freq', 14.65, ... 489 'mni', [26, -11, 73], ... 490 'Comment', ['MSC 14.65Hz ', sourceType]); 491 end 492 end 493 494 495 %% ===== COHERENCE: EMG x SCOUTS (BEFORE) ===== 496 % Process: Select data files 497 sFilesRecEmg = bst_process('CallProcess', 'process_select_files_data', [], [], ... 498 'subjectname', SubjectName, ... 499 'condition', '', ... 500 'tag', 'Left', ... 501 'includebad', 0, ... 502 'includeintra', 0, ... 503 'includecommon', 0); 504 % Only performed for (surface)(Constrained) 505 sourceType = '(surface)(Constr)'; 506 sFilesResSrfCon = bst_process('CallProcess', 'process_select_files_results', [], [], ... 507 'subjectname', SubjectName, ... 508 'condition', '', ... 509 'tag', sourceType, ... 510 'includebad', 0, ... 511 'includeintra', 0, ... 512 'includecommon', 0); 513 % Process: Coherence NxN [2023] 514 sFileCoh1N = bst_process('CallProcess', 'process_cohere2', sFilesRecEmg, sFilesResSrfCon, ... 515 'timewindow', [], ... 516 'src_channel', emg_channel, ... 517 'dest_scouts', {'Schaefer_100_17net', {'Background+FreeSurfer_Defined_Medial_Wall L', 'Background+FreeSurfer_Defined_Medial_Wall R', 'ContA_IPS_1 L', 'ContA_IPS_1 R', 'ContA_PFCl_1 L', 'ContA_PFCl_1 R', 'ContA_PFCl_2 L', 'ContA_PFCl_2 R', 'ContB_IPL_1 R', 'ContB_PFCld_1 R', 'ContB_PFClv_1 L', 'ContB_PFClv_1 R', 'ContB_Temp_1 R', 'ContC_Cingp_1 L', 'ContC_Cingp_1 R', 'ContC_pCun_1 L', 'ContC_pCun_1 R', 'ContC_pCun_2 L', 'DefaultA_IPL_1 R', 'DefaultA_PFCd_1 L', 'DefaultA_PFCd_1 R', 'DefaultA_PFCm_1 L', 'DefaultA_PFCm_1 R', 'DefaultA_pCunPCC_1 L', 'DefaultA_pCunPCC_1 R', 'DefaultB_IPL_1 L', 'DefaultB_PFCd_1 L', 'DefaultB_PFCd_1 R', 'DefaultB_PFCl_1 L', 'DefaultB_PFCv_1 L', 'DefaultB_PFCv_1 R', 'DefaultB_PFCv_2 L', 'DefaultB_PFCv_2 R', 'DefaultB_Temp_1 L', 'DefaultB_Temp_2 L', 'DefaultC_PHC_1 L', 'DefaultC_PHC_1 R', 'DefaultC_Rsp_1 L', 'DefaultC_Rsp_1 R', 'DorsAttnA_ParOcc_1 L', 'DorsAttnA_ParOcc_1 R', 'DorsAttnA_SPL_1 L', 'DorsAttnA_SPL_1 R', 'DorsAttnA_TempOcc_1 L', 'DorsAttnA_TempOcc_1 R', 'DorsAttnB_FEF_1 L', 'DorsAttnB_FEF_1 R', 'DorsAttnB_PostC_1 L', 'DorsAttnB_PostC_1 R', 'DorsAttnB_PostC_2 L', 'DorsAttnB_PostC_2 R', 'DorsAttnB_PostC_3 L', 'LimbicA_TempPole_1 L', 'LimbicA_TempPole_1 R', 'LimbicA_TempPole_2 L', 'LimbicB_OFC_1 L', 'LimbicB_OFC_1 R', 'SalVentAttnA_FrMed_1 L', 'SalVentAttnA_FrMed_1 R', 'SalVentAttnA_Ins_1 L', 'SalVentAttnA_Ins_1 R', 'SalVentAttnA_Ins_2 L', 'SalVentAttnA_ParMed_1 L', 'SalVentAttnA_ParMed_1 R', 'SalVentAttnA_ParOper_1 L', 'SalVentAttnA_ParOper_1 R', 'SalVentAttnB_IPL_1 R', 'SalVentAttnB_PFCl_1 L', 'SalVentAttnB_PFCl_1 R', 'SalVentAttnB_PFCmp_1 L', 'SalVentAttnB_PFCmp_1 R', 'SomMotA_1 L', 'SomMotA_1 R', 'SomMotA_2 L', 'SomMotA_2 R', 'SomMotA_3 R', 'SomMotA_4 R', 'SomMotB_Aud_1 L', 'SomMotB_Aud_1 R', 'SomMotB_Cent_1 L', 'SomMotB_Cent_1 R', 'SomMotB_S2_1 L', 'SomMotB_S2_1 R', 'SomMotB_S2_2 L', 'SomMotB_S2_2 R', 'TempPar_1 L', 'TempPar_1 R', 'TempPar_2 R', 'TempPar_3 R', 'VisCent_ExStr_1 L', 'VisCent_ExStr_1 R', 'VisCent_ExStr_2 L', 'VisCent_ExStr_2 R', 'VisCent_ExStr_3 L', 'VisCent_ExStr_3 R', 'VisCent_Striate_1 L', 'VisPeri_ExStrInf_1 L', 'VisPeri_ExStrInf_1 R', 'VisPeri_ExStrSup_1 L', 'VisPeri_ExStrSup_1 R', 'VisPeri_StriCal_1 L', 'VisPeri_StriCal_1 R'}}, ... 518 'flatten', 0, ... 519 'scouttime', 'before', ... % before connectivity metric 520 'scoutfunc', 'mean', ... % Mean 521 'scoutfuncaft', 'mean', ... % Mean 522 'pcaedit', struct(... 523 'Method', 'pcaa', ... 524 'Baseline', [], ... 525 'DataTimeWindow', [], ... 526 'RemoveDcOffset', 'file'), ... 527 'removeevoked', 0, ... 528 'cohmeasure', cohmeasure, ... % Magnitude-squared coherence 529 'tfmeasure', 'stft', ... % Fourier transform 530 'tfedit', struct(... 531 'Comment', 'Complex', ... 532 'TimeBands', [], ... 533 'Freqs', [], ... 534 'StftWinLen', win_length, ... 535 'StftWinOvr', overlap, ... 536 'StftFrqMax', maxfreq, ... 537 'ClusterFuncTime', 'none', ... 538 'Measure', 'none', ... 539 'Output', 'all', ... 540 'SaveKernel', 0), ... 541 'timeres', 'none', ... % None 542 'avgwinlength', 1, ... 543 'avgwinoverlap', 50, ... 544 'outputmode', 'avg'); % across combined files/epochs 545 % Process: Add tag 546 sFileCoh1N = bst_process('CallProcess', 'process_add_tag', sFileCoh1N, [], ... 547 'tag', sourceType, ... 548 'output', 1); % Add to file name 549 % Highlight scout of interest 550 bst_figures('SetSelectedRows', {'SomMotA_2 R', 'SomMotA_4 R'}); 551 % Process: Snapshot: Frequency spectrum 552 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 553 'target', 10, ... % Frequency spectrum 554 'freq', 14.65, ... 555 'Comment', ['MSC ,' sourceType, ' Before']); 556 % Process: Snapshot: Connectivity matrix 557 bst_process('CallProcess', 'process_snapshot', sFileCoh1N, [], ... 558 'type', 'connectimage', ... % Connectivity matrix 559 'freq', 14.65, ... 560 'Comment', ['MSC 14.65Hz,' sourceType, ' Before']); 561 562 563 %% ===== COHERENCE: SCOUTS x SCOUTS (BEFORE) ===== 564 % Only performed for (surface)(Constrained) 565 sourceType = '(surface)(Constr)'; 566 sFilesResSrfCon = bst_process('CallProcess', 'process_select_files_results', [], [], ... 567 'subjectname', SubjectName, ... 568 'condition', '', ... 569 'tag', sourceType, ... 570 'includebad', 0, ... 571 'includeintra', 0, ... 572 'includecommon', 0); 573 % Process: Coherence NxN [2023] 574 sFileCohNN = bst_process('CallProcess', 'process_cohere1n', sFilesResSrfCon, [], ... 575 'timewindow', [], ... 576 'scouts', {'Schaefer_100_17net', {'SomMotA_2 R', 'SomMotA_4 R', 'VisCent_ExStr_2 L'}}, ... 577 'flatten', 0, ... 578 'scouttime', 'before', ... % before 579 'scoutfunc', 'mean', ... % Mean 580 'scoutfuncaft', 'mean', ... % Mean 581 'pcaedit', struct(... 582 'Method', 'pcaa', ... 583 'Baseline', [], ... 584 'DataTimeWindow', [], ... 585 'RemoveDcOffset', 'file'), ... 586 'removeevoked', 0, ... 587 'cohmeasure', 'mscohere', ... % Magnitude-squared coherence 588 'tfmeasure', 'stft', ... % Fourier transform 589 'tfedit', struct(... 590 'Comment', 'Complex', ... 591 'TimeBands', [], ... 592 'Freqs', [], ... 593 'StftWinLen', 0.5, ... 594 'StftWinOvr', 50, ... 595 'StftFrqMax', 80, ... 596 'ClusterFuncTime', 'none', ... 597 'Measure', 'none', ... 598 'Output', 'all', ... 599 'SaveKernel', 0), ... 600 'timeres', 'none', ... % None 601 'avgwinlength', 1, ... 602 'avgwinoverlap', 50, ... 603 'outputmode', 'avg'); % across combined files/epochs 604 % Process: Add tag 605 sFileCohNN = bst_process('CallProcess', 'process_add_tag', sFileCohNN, [], ... 606 'tag', sourceType, ... 607 'output', 1); % Add to file name 608 % Highlight scout of interest 609 bst_figures('SetSelectedRows', {'SomMotA_2 R', 'SomMotA_4 R'}); 610 % Process: Snapshot: Frequency spectrum 611 bst_process('CallProcess', 'process_snapshot', sFileCohNN, [], ... 612 'target', 10, ... % Frequency spectrum 613 'freq', 14.65, ... 614 'Comment', ['MSC ,' sourceType, ' Before']); 615 % Process: Snapshot: Connectivity matrix 616 bst_process('CallProcess', 'process_snapshot', sFileCohNN, [], ... 617 'type', 'connectimage', ... % Connectivity matrix 618 'freq', 14.65, ... 619 'Comment', ['MSC 14.65Hz,' sourceType, ' Before']); 620 % Process: Snapshot: Connectivity graph 621 bst_process('CallProcess', 'process_snapshot', sFileCohNN, [], ... 622 'type', 'connectgraph', ... % Connectivity graph 623 'freq', 14.65, ... 624 'threshold', 0, ... 625 'Comment', ['MSC 14.65Hz,' sourceType, ' Before']); 626 627 628 %% ===== SAVE REPORT ===== 629 % Save and display report 630 ReportFile = bst_report('Save', []); 631 if ~isempty(reports_dir) && ~isempty(ReportFile) 632 bst_report('Export', ReportFile, reports_dir); 633 else 634 bst_report('Open', ReportFile); 635 end 636 637 disp([10 'DEMO> Corticomuscular coherence tutorial completed' 10]); 638





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


Tutorials/CorticomuscularCoherence (last edited 2022-05-25 15:28:15 by SylvainBaillet)