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

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

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. This computation varies if we are using constrained or unconstrained sources.

Method for constrained sources

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.

width="100%"

Let's compute coherence between the EMG sensor and the sources from the cortical surface/constrained model.

Double-click on the 1xN connectivity file 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.

Method for unconstrained sources

In the case of orientation-unconstrained (either surface or volume) sources, each vertex (source location) is associated with 3 time series, each one corresponding to the X, Y and Z directions. Thus, a dimension reduction across directions is needed. In this tutorial the dimension reduction consisted in flattening the unconstrained maps to a single orientation with the PCA method. This results in a source model similar to that with orientation-constrained sources, thus, coherence is then computed as with orientation-constrained (described above).

PCA Options: If PCA is selected for flattening unconstrained maps (or as scout function), additional options are available through the [Edit] button.

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

Let's compute coherence between the EMG sensor and: the sources from the cortical surface/unconstrained and brain volume/unconstrained models.

Repeat now using the source files form the brain volume/unconstrained model. Do not forget to update the Tag text to (volume)(Unconstr) in the Add tag process.

An alternative to address the issue of the dimension reduction across directions when estimating coherence with unconstrained sources is to flatten the results after computing coherence. Coherence is computed between the EMG signal (one time series) and each of the 3 time series (one for each source orientation: x, y and z) associated with each source location, resulting in 3 coherence spectra per source location. Then, these coherence spectra are collapsed (flattened) into a single coherence spectrum by keeping the maximum coherence value across those in the three directions for each frequency bin. 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 connectivity 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_unconstr.png

Comparison

Close the search tab. If the 3 new connectivity files ) are not featured in the database explorer: refresh the database display by pressing [F5] or by clicking again on the selected button "Functional data".

tre_coh_src_files.gif

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.

If we use the alternative method of reducing dimensions across directions by obtaining the maximum coherence across orientations. Coherence maps look even smoother than with PCA, this is due to the maximum operation. The figure below shows the results of coherence between EMG sensor and: constrained sources (top), unconstrained sources flattened with PCA flattening (middle), and unconstrained with maximum aggregation across coherence orientations (bottom).

res_coh_surf_all.gif

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. To compute scout-wise coherence, it is necessary to specify two parameters to indicate how the data is aggregated in each scout:

1. The '''scout function''' (for instance, the mean), and

2. When the scout aggregation if performed, either before or after the coherence computation).

Method for constrained sources

When the scout function is used with constrained sources, the computation of coherence is carried as follows:

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

https://neuroimage.usc.edu/brainstorm/Tutorials/CorticomuscularCoherence?action=AttachFile&do=get&target=diagram_1xn_coh_sct_aft_c.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 on 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.

Method for unconstrained sources

As described above, when computing coherence with unconstrained sources it is necessary to perform dimension reduction across directions. When this reduction is done by first flattening the unconstrained sources, is performed on unconstrained sources coherence computations is carried out as described in the previous section.

However, if the sources are not flattened before computing coherence is carried out as follows:

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

Coherence Scouts x Scouts

When computing whole-brain connectomes (i.e., NxN connectivity matrices between all brain sources), additional technical considerations need to be taken into account: the number of time series is typically too large to be computationally tractable on most workstations (see introduction above); if using unconstrained source maps, with 3 time series at each brain location, this requires an extra step to derive connectivity estimates.

The last two sections explained the computation of coherence measures between one reference sensor (EMG) and brain source maps, with sensor x sources and sensor x scouts scenarios, both with constrained and unconstrained source orientation models.

This section explains dimension reduction (across number of sources) using scouts to compute scouts x scouts coherence, for constrained and unconstrained source maps. As in the sensor-scout coherence scenario, it is necessary to specify two parameters to indicate how the data is aggregated in each scout:

1. The scout function (for instance, the mean), and

2. When the scout aggregation if performed, either before or after the coherence computation).

The sections below explain the processing steps for the cases of constrained and unconstrained sources, for coherence between two scouts (Scout-1 and Scout-2). In either case, the outcome is one coherence spectrum per pair of scouts.

Method for constrained sources

The computation steps depend of when the scout function is applied:

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

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

Here, we will compute the coherence among three scouts, using mean as scout function with the before option. We will use the Schaefer 100 parcellation atlas on to the orientation-constrained source map. These three scouts are: SomMotA_2R, SomMotA_4R and VisCent_ExStr_2L. The firs two scouts are adjacent on the right motor cortex, while the third is located on the left occipital lobe. We expect a higher coherence between the first around 15 Hz, and low coherence across all frequencies between the third scout and the others.

Open the Pipeline editor:

pro_coh_srfc_nn_bef_sct.png

pro_coh_srfc_nn_bef_sct2.png

Right-click on the new NxN connectivity file (), and display it as image, spectrum ad graph. For the connectivity graph, adjust the intensity threshold to its minimum to display all the connections. Se the frequency slider at 14.65 Hz. Note that MSC is higher between the scouts over the right motor cortex, than between any of these scouts and the scout over the visual cortex.

res_coh_srfc_nn_bef_sct.png

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

Method for unconstrained sources

As described above, when computing coherence with unconstrained sources it is necessary to perform dimension reduction across directions. When this reduction is done by first flattening the unconstrained sources, is performed on unconstrained sources coherence computations is carried out as described in the previous section. However, if the sources are not flattened before computing coherence is carried out as follows:

Before: The scout function (e.g., mean) is applied on each direction on the elementary source time series each scout, which results in 3 time series per elementary direction, per scout. These scouts time series are then used to compute coherence measures, yielding to 9 coherence spectra (1x-2x, 1x-2y, 1x-2z, 1y-2x, ..., 1z-2z) for each pair of scouts. The resulting coherence spectra are finally aggregated across the 9 dimensions with maximum per frequency bin, resulting in one single coherence spectrum for each pair of scouts.

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

After: With this option, coherence is first computed between each pair of source locations (each with 3 time series) in the two scouts, yielding to 9 coherence spectra for each pair of sources in the scouts. Then, these coherence spectra are aggregated across the 9 combinations of dimensions (1x-2x, 1x-2y, 1x-2z, 1y-2x, ..., 1z-2z) to obtain one single coherence spectrum for each pair of sources within the scouts. Finally, the scout function (e.g., mean) is applied twice (once per scout) to the coherence spectra for each source within the scouts to obtain one single coherence spectrum for each pair of scouts. Applying the scout function "after" is a considerably more greedy option that the applying the scout function "before", as it computes coherence measures between all 45000x45000 elementary source signals, and this number is multiplied by the number of frequency bins. This requires more computing resources; see the introduction for an example.

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

Method for mixed source models

Mixed models refer to source models with multiple regions that may include surface and volume regions, and regions with different source orientation constraints. To compute connectivity between surface and volume regions, you must first create volume scouts for the volume regions of the mixed model, and then use the process2 tab with the same files on both sides. On one side you can select surface scouts, and the other, the volume scouts you created. For connectivity between surface regions only, or volume regions only, use the process1 tab.

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-15 19:46:08 by RaymundoCassani)