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





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)