Corticomuscular coherence (CTF MEG)

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 at the top of the figure, above the data time series, indicate event markers (or triggers) saved along with MEG data. The actual onsets of the left- and right-wrist trials are not displayed yet: they are saved in an auxiliary channel of the raw data named Stim. To add these markers to the display, follow this procedure:

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 original 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 specific tutorial sections about the detection and removal of artifacts with SSP. The present tutorial dataset features an EOG channel but no ECG. We will therefore only remove artifacts caused by 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, or muscle contractions.

Epoching

We are now finished with the pre-processing of 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). For consistency with the FieldTrip tutorial, we will analyze 8s of recordings following each movement (from the original 10s around each trial), and split them in 1-s epochs. For each 1-s epoch, the DC offset will be also removed from each MEG channel. And finally, we will uniform the time vector for all 1-s epochs. This last steps is needed as the time vector of each 1-s vector is with respect to the event that was used to import the 8-s trial.

pro_import.png

pro_remove_dc.png

Comparison with FieldTrip

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

bst_ft_trial1.png

Coherence: EMG x MEG

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

coh_meg_emgleft.png

coh_meg_emgleft2.png

Source estimation

MRI segmentation

We first need to extract the cortical surface from the T1 MRI volume we imported at the beginning of this tutorial. CAT12 is a Brainstorm pluing that will perform this task in 30-60min.

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. Forward models are called head models in Brainstorm. They account 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). A distinct head model is required for the cortex surface and head volume source spaces.

Cortical surface

Whole-head volume

Noise covariance

The recommendation for MEG source imaging is to extract basic noise statistics from empty-room recordings. When not available, as here, 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 and major artifacts: 18s-29s.

Inverse models

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

Cortical surface

pro_sources_srfc.png

pro_sources_srfu.png

Whole-head volume

Coherence: EMG x Sources

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

pro_coh_srf.png

pro_coh_srf2.png

Surface

Double-click on the 1xN connectivity files for the two (surface) source spaces (constrained and unconstrained source orientations) to visualize the cortical maps. See main tutorial Display: Cortex surface for all available options. Pick the cortical location and frequency with the highest coherence value.

Volume

Method

For cortical and orientation-constrained sources, each vertex in the source grid is associated with 1 time series. As such, when coherence is computed with the EMG signal (also consisting of one time series), this result is 1 coherence spectrum per vertex. In other words, for each frequency bin, we obtain one coherence brain map.

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

In the case of orientation-unconstrained or volume sources, each vertex in the grid is associated with 3 time series, each one corresponding to the X, Y and Z directions. Thus, when coherence is computed with the EMG signal (one time series), 3 coherence spectra are derived (one for each source orientation: x, y and z). For visualization purposes, these three values need to be flattened into one, yielding only one coherence spectrum per vertex, as in the orientation-constrained case above. At each vertex location, we will select the maximum coherence value across those in the three directions, for each frequency bin.

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

Coherence: EMG x Scouts

We have now obtained coherence spectra at each of the 15002 brain source locations. This is a very large amount of data to shift through. We therefore recommend to further reduce the dimensionality of the source space by using a surface- or volume-parcellation scheme. In Brainstorm vernacular, this can be achieved via an atlas consisting of scout regions. See the Scout Tutorial for detailed information about atlases and scouts in Brainstorm.

* Scout function: With this approach, a coherence spectrum is computed for each scout and the EMG signal, for each elementary source orientation. It is necessary to specify two parameters to indicate how the data is aggregated in each scout: The scout function (for instance, the mean), and the within-scout aggregation procedure (before or after the coherence computation).

* Unconstrained maps / maximum: The graphs below the scout has 3 orientations (x,y,z). Coherence is computed between each elementary sources in all three orientations (x, y,z) and the EMG signal, yielding 3 coherence spectra. From these 3 values at each frequency bin, only the maximum value is kept to summarize the connectivity estimate between the scout and the EMG signal. The connectvity outcome measure is therefore one coherence spectrum for each scout. The choice of the maximum statistic is empirical: it is not rotation invariant (if the specifications of the NAS/LPA/RPA fiducials change, the x,y,z axes also change, which may influence conenctivity measures quantitatively). Nevertheless, our empirical tests showed this option produced smooth and robust (to noise) connectivity estimates.

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

Here, we will compute the coherence from scouts, using mean as scout function with the Before option. We will use the Schaefer 100 parcellation atlas appplied to the orientation-constrained source map.

Open the Pipeline editor:

pro_coh_srfc_bef_sct.png

pro_coh_srfc_bef_sct2.png

Double-click on the new file: the coherence spectra are not displayed on the cortex; they are plotted for each scout. The 1xN connectivity file can also be shown as an image as displayed below:

res_coh_srfc_bef_sct.png

res_coh_srfc_bef_sct2.png

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

The brain parcellation scheme is the user's decision. We recommend the Schaefer100 atlas used here by default.

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 2024-04-12 23:56:58 by RaymundoCassani)