Corticomuscular coherence (CTF MEG)

[TUTORIAL UNDER DEVELOPMENT: NOT READY FOR PUBLIC USE]

Authors: Raymundo Cassani, Francois Tadel & Sylvain Baillet.

Corticomuscular coherence measures the degree of similarity between electrophysiological signals (MEG, EEG, ECoG sensor traces or source time series, especially over the contralateral motor cortex) and the EMG signal recorded from muscle activity during voluntary movement. This signal similarity is due mainly to the descending communication along corticospinal pathways between primary motor cortex (M1) and muscles. For consistency and reproducibility purposes across major software toolkits, the present tutorial replicates the processing pipeline "Analysis of corticomuscular coherence" by FieldTrip.

Background

Coherence measures the linear relationship between two signals in the frequency domain. Previous studies (Conway et al., 1995, Kilner et al., 2000) have reported corticomuscular coherence effects in the 15–30 Hz range during maintained voluntary contractions.

IMAGE OF EXPERIMENT, SIGNALS and COHERENCE

Dataset description

The dataset comprises MEG recordings (151-channel CTF MEG system) and bipolar EMG recordings (from left and right extensor carpi radialis longus muscles) from one participant who was tasked to lift their hand and exert a constant force against a lever for about 10 seconds. The force was monitored by strain gauges on the lever. The participant performed two blocks of 25 trials using either the left or right wrist. EOG signals were also recorded, which will be useful for detection and attenuation of ocular artifacts. We will analyze the data from the left-wrist trials in the present tutorial. Replicating the pipeline with right-wrist data is a good exercise to do next!

Download and installation

Importing anatomy data

MEG and EMG recordings

MEG-MRI coregistration

Reviewing continuous recordings

Event markers

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

Pre-process

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

Removal of power line artifacts

EMG pre-processing

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

Two new folders SubjectCMC_notch_high and SubjectCMC_notch_high_abs are added to Brainstorm database explorer. We can now safely delete folders that are not needed anymore:

MEG pre-processing

Detection and removal of artifacts with SSP

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

Detection of "bad" data segments:

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

We encourage users to review and validate the segments marked using this procedure. In the present case, the segments detected as bad clearly point at contaminated MEG data segments. We will now label these as "bad".

Importing data epochs

At this point we are finished with the pre-processing of the EMG and MEG recordings. We will now extract and import specific data segments of interest into the Brainstorm database for further derivations. We refer to these segments as epochs or trials. As mentioned previously, we will focus on the Left (wrist) category of events. To follow the same pipeline as the FieldTrip tutorial, only the first 8 1-s epochs for each trial are imported. In addition DC offset is removed only for MEG signals.

pro_import.png

pro_remove_dc.png

A new folder SubjectCMC_notch_high_abs without the 'raw' indication is created for Subject01.

The new folder contains a copy of the channel file from the original raw file, and 192 (24 trials, 8 epochs each) individual trials tagged as Left in a new trial group (). Expand the trial group and note there are trials marked with an exclamation mark in a red circle (). This indicates trials that contain segments of signal occurred in the bad segments we previously identified. All the bad trials will be automatically ignored for further processing, whenever dropped into the Process1 and Process2 tabs.

To have a glance of the signals after the pre-processing, plot the MEG signal for one sensor over the left motor-cortex (MRC21) and the EMG signals, for trial 1 (left image). Note that these traces are similar to the ones obtained in the FieldTrip tutorial (right image).

bst_ft_trial1.png

Coherence 1xN (sensor level)

We will now compute the magnitude square coherence (MSC) between the left EMG signal and each of the MEG sensor data.

Open the Pipeline editor:

coh_meg_emgleft.png

coh_meg_emgleft2.png

We can now average magnitude of the MSC across a frequency band of interest (15-20 Hz):

The resulting file mscohere(0.6Hz,555win): EMGlft | MEG sensors | tfbands has only one MSC value for each sensor (the MSC average in the 15-20 Hz band). You may visualize the topography of this MSC statistics via 3 possible representations: 2D Sensor cap, 2D Sensor cap and 2D Disk, which are all accessible via a right-click over the MSC node. We clicked on sensor MRC21 below; it is shown in red.

We can observe higher MSC values between the EMG signal and MEG sensor signals over the contralateral set of central sensors in the beta band. Unfortunately, sensor-level connectivity presents the disadvantages of not being interpretable, and being subject to spurious results due to volume conduction.

In the next section we will compute coherence in the source level. To do this, we first need to estimate the sources time series from the sensor data.

MRI segmentation

We then need to segment the head tissues to obtain the surfaces required to derive a realistic MEG head models (aka "forward models"). Here, we will perform MRI segmentation with CAT12, this process takes between 30 to 60 minutes.

CAT12 computes the scalp surface (), and the cortex surfaces () for white matter, pial envelope and the midpoint between them. The default surface of each type is indicated in green. Moreover, as part of the MRI segmentation with CAT12, the anatomy data is normalized in the MNI space, and several anatomical parcellations () are computed. For further information on the anatomy files see the Display the anatomy tutorial.

MEG source modelling

We will perform source modelling using a distributed model approach for two possible source spaces: the cortex surface and the entire MRI volume. In the surface space, a source grid is made with the source located on the cortical surface obtained from the participant's MRI. For the volume space, the source grid consists of elementary sources uniformly distributed across the entire brain volume. Before estimating the brain sources, we need to derive the sensor noise covariance matrix, and the head model.

Noise covariance

The recommendation for MEG, is to extract basic noise statistics from empty-room recordings. However, when recommended empty-room recordings are not available, as with this tutorial data, resting-state data can be used as proxies for MEG noise covariance. See the noise covariance tutorial for more details.

Head model

The head model, aka forward model, accounts for how neural electrical currents (in a source space) produce magnetic fields captured by sensors outside the head, considering head tissues electromagnetic properties and geometry, independently of actual empirical measurements. As the head model depends on the source space, a distinct head model is required for the surface and volume source spaces. Please refer to the head model tutorial for more in-depth explanations.

Surface

The cortical head model will be derived from each of the 15,000 sources (surface vertices) as defined in the default cortex.

The Overlapping spheres (surface) head model () now appears in the database explorer.

Volume

The Overlapping spheres (volume) head model is now added to the database explorer. The green color indicates this is the default head model for the current folder (this can be changed by simply double clicking over the head model nodes.)

Source estimation

Now that the noise covariance and head model(s) are available, we will perform source estimation, to find the sources that gave origin to the signals registered in the sensors. From the diverse source estimation methods available in Brainstorm, in this tutorial the minimum-norm imaging method is used. The minimum-norm method estimates the linear combination of the current at each point in the source grid that explains the recorded sensor signals favouring minimum energy (L2-norm) solutions. As result, a large matrix called the imaging kernel is obtained. By multiplying the imaging kernel with the sensor data, it is possible to obtain the estimates of brain sources time series. A different imaging kernel is derived for each of the head models we have produced above: surface and volume. See the source estimation tutorial for more details.

Each dipole in the source grid may point arbitrarily in any direction in a 3D space.

Only for surface grids, the dipole orientation can be fixed to be normal to the cortical surface, this approach is based on anatomical observations of the brain cortex. The result is then in a smaller model that is faster to compute and display.
A discussion on constrained vs unconstrained sources is presented here.

Surface

Here we will estimate the sources in the surface space for constrained (normal to the cortex) and unconstrained dipole orientations.

pro_sources_srfc.png

pro_sources_srfu.png

The inversion kernels () MN: MEG (surface)(Constr) 2018 and MN: MEG (surface)(Unconstr) 2018 are now available in the database explorer.

Volume

To compute the imaging kernel for the volume source space:

At this point the imaging kernel () MN: MEG (volume)(Unconstr) 2018 is now also available in the database explorer.


Note that each trial is associated with three source link () nodes, that correspond to each of the imaging kernels obtained above.

The source link nodes are not files containing the sources time series, instead, the links indicate Brainstorm to: load the corresponding MEG recordings, load the respective inverse kernel, and multiply the two on the fly to generate the sources time series. This approach saves space on the hard drive.

Coherence 1xN (source level)

Once we have computed the time series for the sources, it is time to compute coherence between the EMG signal and brain sources obtained with each of the imaging kernels.

Let's start with sources from the MN: MEG (surface)(Constr) kernel:

This will create a new tab in the database explorer. This new tab contains only the files that match the search criteria.

Open the Pipeline editor:

pro_coh_srf.png

pro_coh_srf2.png

Results (Surface)

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

Find the location and frequency with the highest coherence value.

The highest coherence value is located on the right primary motor cortex (precentral gyrus) around 14.65 Hz for the analysis using constrained and unconstrained sources. Set the amplitude threshold to 0% to see the extension of the high coherence values.

res_coh_srfc.png

res_coh_srfu.png

MSC @ 14.65Hz (surface)(Constr)

MSC @ 14.65Hz (surface)(Unconstr)

We observe that results obtained with constrained and unconstrained sources agree in the location and frequency of the peak coherence. The main difference between these results is that coherence values obtained with unconstrained sources appear smoother, this caused by the maximum aggregation performed across directions, explained at detail in the next section.

Finally, right-click on any of the cortex figures and select Get coordinates. Then click on the right motor cortex with the crosshair cursor that appears. The SCS coordinates will be useful in the next section.
SCS coordinates X:39.2, Y:-22.3 and Z: 113.0

res_get_coordinates.png

Results (Volume)

Double-click the 1xN connectivity file for the (volume) source space. Note that this time the results are shown in the MRI viewer rather than the cortical surface.

res_coh_vol.png

MSC @ 14.65Hz (volume)(Unconstr)

We note that all the results obtained with constrained (surface) and unconstrained (surface and volume) sources agree with each other, and in the location and frequency of the peak coherence. Moreover, they agree with our hypothesis, previous results in the literature [REFS], and the results presented in the FieldTrip tutorial.

Coherence with constrained and unconstrained sources

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

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

In the case of unconstrained sources, each vertex in the grid is associated with THREE time series, each one corresponding to the X, Y and Z directions. Thus, when coherence is computed with the EMG signal (one time series), there are THREE coherence spectra. This complicates its representation in the brain, thus; these THREE coherence spectra need to be flattened into one, resulting in one coherence spectrum per vertex, the maximum across directions is found for each frequency bin for each vertex.

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

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

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

  • Flattened sources are saved as full rather than recordings+kernel.
  • We have tested this flattening approach with simulations and the real data (from this tutorial) and we have found decrimental effects on the expected results.

Coherence 1xN (scout level)

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

Under this approach, instead of providing one result (coherence spectrum) per source vertex, one result is computed for each scout. When computing coherence (or other connectivity metrics) in the scout level, it is necessary to provide two parameters that define how the data is aggregated per scout:

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

As it can be seen, the After option takes longer and used more resources as it computes the coherence spectrum for each vertex in the scouts, and then, the coherence spectra are aggregated.


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

Open the Pipeline editor:

pro_coh_srfc_bef_sct.png

pro_coh_srfc_bef_sct2.png

Open the file: mscohere(0.6Hz,555win): EMGlft x 102 scouts, mean before | (surface)(Constr) by double-clicking on it. This time the coherence spectra are not displayed on the cortex, but they are plotted for each scout. Moreover, 1xN connectivity file can be shown as image.

res_coh_srfc_bef_sct.png

res_coh_srfc_bef_sct2.png

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

The choice of the optimal parcellation scheme for the source space is not easy.
The optimal choice is to choose a parcellation based on anatomy, for example the Brodmann parcellation.
In Brainstorm these atlases are imported in Brainstorm as scouts (cortical regions of interest), and saved directly in the surface files as explained in this tutorial here.

BRAINSTORM TEAM
Due to the current implementation of the bst_connectivity, the full source map for each trial (185) are loaded in memory, thus replicate the After, only for the (surface)(Constr) option ~30GB of RAM are needed! (Unconstrained take 3 times that).

Coherence NxN (scout level)

In previous sections, we have computed coherence between a reference signal (EMG) and sources (or scouts). However, depending on the experimental setup and hypotheses, we may want to study the brain as network, this is to say, to compute the NxN connectivity. Due to the large number of sources (often several thousands) while in theory it is possible to compute NxN connectivity, it is not practical, as such a common approach is to use of scouts to reflect the functional connections among its cell assemblies, resulting in a connectome.

Diverse studies have shown some overlap between connectomes derived from electrophysiological signals (MEG/EEG) and the ones derived from fMRI, which is reasonably expected as both are the result of the undergoing biological system. However, due to its nature, the electrophysiological connectomes provide unique insights on how functional communication is implemented in the brain (Sadaghiani et al., 2022).

Similar to the coherence 1xN with scouts, we need to define the scout function and when it is going to be applied. Note that computing coherence between two unconstrained leads to 9 (3x3) coherence spectra.

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

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

As in the coherence 1xN with scouts case, the After option takes longer and uses way more resources as it computes the coherence spectrum for each vertex in all the scouts, leading to arrays as big as [sources, sources, frequency bins]!


The result is a NxN connectivity file (), It contains (Scouts x Scouts) coherence spectra.

diagram_nxn_coh_sct_end.png

Such a visualization is not practical, thus the connectivity graph or the adjacent matrix are displayed for each frequency bin in the coherence spectra. For more details, see the connectivity graph tutorial.

A coherence NxN example can be found in the main connectivituy tutorial.

Advanced

Scripting

Additional documentation

Articles

Tutorials

Forum discussions

[TODO] Find relevant Forum posts.





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


Tutorials/CorticomuscularCoherence (last edited 2022-03-11 13:34:56 by FrancoisTadel)