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

An alternative approach for addressing the issue of the 3-dimensional nature of coherence with unconstrained sources, consists in flattening (reducing the dimension) the vertex timeseries in the X, Y and Z directions before estimating coherence. This would result in a source model similar to that with orientation-constrained sources (one timeseries per vertex). Dimension reduction techniques such as PCA (keeping only the first PCA component, see (Euclidean) norm) can be applied. In Brainstorm, use the process: Sources > Unconstrained to flat map. However, our empirical tests with this approach showed a reduction of sensitivity, therefore this is not the method we recommend.

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

  • Note also that the timeseries of flattened sources are saved as a full array of values, not as a kernel; file storage volume on disk may be large.

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 % 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 2024-04-12 22:29:00 by RaymundoCassani)