Visualising spatial contributions to multivariate connectivity

This example demonstrates how the spatial filters and patterns of connectivity obtained from the decomposition tools in the decoding module can be visualised and interpreted.

# Author: Thomas S. Binns <t.s.binns@outlook.com>
# License: BSD (3-clause)

# sphinx_gallery_multi_image = "single"

import mne
from mne import make_fixed_length_epochs
from mne.datasets.fieldtrip_cmc import data_path

from mne_connectivity import CoherencyDecomposition

Background

Multivariate forms of signal analysis allow you to simultaneously consider the activity of multiple signals. In the case of connectivity, the interaction between multiple sensors can be analysed at once and the strongest components of this interaction captured in a lower-dimensional set of connectivity spectra. This approach brings not only practical benefits (e.g. easier interpretability of results from the dimensionality reduction), but can also offer methodological improvements (e.g. enhanced signal-to-noise ratio and reduced bias).

Coherency-based methods are popular approaches for analysing connectivity, capturing correlations between signals in the frequency domain. Various coherency-based multivariate methods exist, including: canonical coherency (CaCoh; multivariate measure of coherency/coherence) [1]; and maximised imaginary coherency (MIC; multivariate measure of the imaginary part of coherency) [2].

These methods are described in detail in the following examples:

• comparison of coherency-based methods - Comparison of coherency-based methods

• CaCoh - Compute multivariate coherency/coherence

• MIC - Compute multivariate measures of the imaginary part of coherency

The CaCoh and MIC methods work by finding spatial filters that decompose the data into components of connectivity, and applying them to the data. Connectivity can then be computed on this transformed data (see Multivariate decomposition for efficient connectivity analysis for more information).

However, in addition to the connectivity scores, useful insights about the data can be gained by visualising the topographies of the spatial filters and their corresponding spatial patterns. These provide important information about the spatial distributions of connectivity information, and represent two complementary aspects:

• The filters represent how the connectivity sources are extracted from the channel data, akin to an inverse model.

• The patterns represent how the channel data is formed by the connectivity sources, akin to a forward model.

This distinction is discussed further in Haufe et al. (2014) [3], but in short: the patterns should be used to interpret the contribution of distinct brain regions/sensors to a given component of connectivity. Accordingly, keep in mind that the filters and patterns are not a replacement for source reconstruction, as without this the patterns will still only tell you about the spatial contributions of sensors, not underlying brain regions, to connectivity.

Generating the filters and patterns

We will first load some example MEG data collected during a hand movement task, which we will generate the spatial filters and patterns for (see here for more information on the data). We divide the data into continuous epochs.

# Load example MEG data
raw = mne.io.read_raw_ctf(data_path() / "SubjectCMC.ds")
raw.pick("mag")
raw.crop(50.0, 110.0).load_data()
raw.notch_filter(50)
raw.resample(100)

# Create epochs
epochs = make_fixed_length_epochs(raw, duration=2.0).load_data()

ds directory : /home/circleci/mne_data/MNE-fieldtrip_cmc-data/SubjectCMC.ds
    res4 data read.
    hc data read.
    Separate EEG position data file not present.
    Quaternion matching (desired vs. transformed):
       0.33   78.32    0.00 mm <->    0.33   78.32    0.00 mm (orig :  -71.62   40.46 -256.48 mm) diff =    0.000 mm
      -0.33  -78.32   -0.00 mm <->   -0.33  -78.32   -0.00 mm (orig :   39.27  -70.16 -258.60 mm) diff =    0.000 mm
     114.65    0.00   -0.00 mm <->  114.65    0.00   -0.00 mm (orig :   64.35   66.64 -262.01 mm) diff =    0.000 mm
    Coordinate transformations established.
    Polhemus data for 3 HPI coils added
    Device coordinate locations for 3 HPI coils added
Picked positions of 4 EEG channels from channel info
    4 EEG locations added to Polhemus data.
    Measurement info composed.
Finding samples for /home/circleci/mne_data/MNE-fieldtrip_cmc-data/SubjectCMC.ds/SubjectCMC.meg4:
    System clock channel is available, checking which samples are valid.
    75 x 12000 = 911610 samples from 191 chs
    390 samples omitted at the end
Current compensation grade : 0
Removing 5 compensators from info because not all compensation channels were picked.
Reading 0 ... 72000  =      0.000 ...    60.000 secs...
Filtering raw data in 1 contiguous segment
Setting up band-stop filter from 49 - 51 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 49.38
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 49.12 Hz)
- Upper passband edge: 50.62 Hz
- Upper transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 50.88 Hz)
- Filter length: 7921 samples (6.601 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.3s
Not setting metadata
30 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 30 events and 200 original time points ...
0 bad epochs dropped

We designate the left hemisphere sensors as the seeds and the right hemisphere sensors as the targets. Since this is sensor-space data, we will use the MIC method to analyse connectivity, given its resilience to zero time-lag interactions (see Comparison of coherency-based methods for more information).

# Left hemisphere sensors
seeds = [idx for idx, ch_info in enumerate(epochs.info["chs"]) if ch_info["loc"][0] < 0]

# Right hemisphere sensors
targets = [
    idx for idx, ch_info in enumerate(epochs.info["chs"]) if ch_info["loc"][0] > 0
]

# Define indices
indices = (seeds, targets)

To fit the filters (and in turn compute the corresponding patterns), we instantiate the CoherencyDecomposition object and call the fit() method. We also define our connectivity frequency band of interest to be 20-30 Hz. See Multivariate decomposition for efficient connectivity analysis for more information.

# Instantiate decomposition object
mic = CoherencyDecomposition(
    info=epochs.info,
    method="mic",
    indices=indices,
    mode="multitaper",
    fmin=20,
    fmax=30,
    rank=(3, 3),
)

# Fit filters & generate patterns
mic.fit(epochs.get_data())

    Using multitaper spectrum estimation with 7 DPSS windows
Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/30 [00:00<?,       ?it/s]
  3%|▎         | CSD epoch blocks : 1/30 [00:00<00:02,   11.67it/s]
  7%|▋         | CSD epoch blocks : 2/30 [00:00<00:02,   11.21it/s]
 10%|█         | CSD epoch blocks : 3/30 [00:00<00:02,   11.27it/s]
 13%|█▎        | CSD epoch blocks : 4/30 [00:00<00:02,   10.91it/s]
 17%|█▋        | CSD epoch blocks : 5/30 [00:00<00:02,   10.69it/s]
 20%|██        | CSD epoch blocks : 6/30 [00:00<00:02,   10.52it/s]
 23%|██▎       | CSD epoch blocks : 7/30 [00:00<00:02,   10.66it/s]
 27%|██▋       | CSD epoch blocks : 8/30 [00:00<00:02,   10.76it/s]
 30%|███       | CSD epoch blocks : 9/30 [00:00<00:01,   10.87it/s]
 33%|███▎      | CSD epoch blocks : 10/30 [00:00<00:01,   10.96it/s]
 37%|███▋      | CSD epoch blocks : 11/30 [00:01<00:01,   11.03it/s]
 40%|████      | CSD epoch blocks : 12/30 [00:01<00:01,   11.08it/s]
 43%|████▎     | CSD epoch blocks : 13/30 [00:01<00:01,   11.09it/s]
 47%|████▋     | CSD epoch blocks : 14/30 [00:01<00:01,   10.95it/s]
 50%|█████     | CSD epoch blocks : 15/30 [00:01<00:01,   11.01it/s]
 53%|█████▎    | CSD epoch blocks : 16/30 [00:01<00:01,   11.07it/s]
 57%|█████▋    | CSD epoch blocks : 17/30 [00:01<00:01,   11.11it/s]
 60%|██████    | CSD epoch blocks : 18/30 [00:01<00:01,   11.13it/s]
 63%|██████▎   | CSD epoch blocks : 19/30 [00:01<00:00,   11.13it/s]
 67%|██████▋   | CSD epoch blocks : 20/30 [00:01<00:00,   11.17it/s]
 70%|███████   | CSD epoch blocks : 21/30 [00:01<00:00,   11.20it/s]
 73%|███████▎  | CSD epoch blocks : 22/30 [00:01<00:00,   11.23it/s]
 77%|███████▋  | CSD epoch blocks : 23/30 [00:02<00:00,   11.26it/s]
 80%|████████  | CSD epoch blocks : 24/30 [00:02<00:00,   11.28it/s]
 83%|████████▎ | CSD epoch blocks : 25/30 [00:02<00:00,   11.30it/s]
 87%|████████▋ | CSD epoch blocks : 26/30 [00:02<00:00,   11.33it/s]
 90%|█████████ | CSD epoch blocks : 27/30 [00:02<00:00,   11.34it/s]
 93%|█████████▎| CSD epoch blocks : 28/30 [00:02<00:00,   11.35it/s]
 97%|█████████▋| CSD epoch blocks : 29/30 [00:02<00:00,   11.32it/s]
100%|██████████| CSD epoch blocks : 30/30 [00:02<00:00,   11.29it/s]
100%|██████████| CSD epoch blocks : 30/30 [00:02<00:00,   11.21it/s]
[done]
Computing MIC for connection 1 of 1

CoherencyDecomposition(fmax=30, fmin=20,
                       indices=(array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
        13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
        26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
        39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
        52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
        65,  66,  67,  68,  69,  70, 147, 148, 149, 150], dtype=int32),
                                array([ 71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
        84,  85,  86,  87,  88,  8...
 dev_ctf_t: MEG device -> CTF/4D/KIT head transform
 dev_head_t: MEG device -> head transform
 dig: 7 items (3 Cardinal, 4 EEG)
 highpass: 0.0 Hz
 hpi_results: 1 item (list)
 lowpass: 50.0 Hz
 meas_date: 2005-01-31 14:43:00 UTC
 meas_id: 4 items (dict)
 nchan: 151
 projs: []
 sfreq: 100.0 Hz
 subject_info: <subject_info | his_id: Anonymized25263_1299153182_0>
>,
                       method='mic', n_components=np.int64(3), rank=(3, 3))

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Visualising the patterns

Visualising the patterns as topomaps can be done using the plot_patterns() method.

When interpreting patterns, note that the absolute value reflects the strength of the contribution to connectivity, and that the sign differences can be used to visualise the orientation of the underlying dipole sources. The spatial patterns are not bound between \([-1, 1]\).

Plotting the patterns for 20-30 Hz connectivity below, we find the strongest connectivity (‘MIC0’, i.e. 1st component) between the left and right hemispheres comes from centromedial left and frontolateral right sensors, based on the areas with the largest absolute values. Patterns for the weaker connectivity components (‘MIC1’ & ‘MIC2’ are also shown). As these patterns come from decomposition on sensor-space data, we make no assumptions about the underlying brain regions involved in this connectivity.

# Plot patterns
mic.plot_patterns(info=epochs.info, sensors="m.", size=2)

• 
• 

[<MNEFigure size 900x355 with 4 Axes>, <MNEFigure size 900x355 with 4 Axes>]

Visualising the filters

We can also visualise the filters as topomaps using the plot_filters() method.

Here we see that the filters show a similar topography to the patterns. However, this is not always the case, and you should never confuse the information represented by the filters (i.e. an inverse model) and patterns (i.e. a forward model), which can lead to very incorrect interpretations of the data [3].

# Plot filters
mic.plot_filters(info=epochs.info, sensors="m.", size=2)

• 
• 

[<MNEFigure size 900x355 with 4 Axes>, <MNEFigure size 900x355 with 4 Axes>]

References

[1]

Carmen Vidaurre, Guido Nolte, Ingmar E.J. de Vries, M. Gómez, Tjeerd W. Boonstra, K.-R. Müller, Arno Villringer, and Vadim V. Nikulin. Canonical maximization of coherence: a novel tool for investigation of neuronal interactions between two datasets. NeuroImage, 201:116009, 2019. doi:10.1016/j.neuroimage.2019.116009.

[2]

Arne Ewald, Laura Marzetti, Filippo Zappasodi, Frank C. Meinecke, and Guido Nolte. Estimating true brain connectivity from EEG/MEG data invariant to linear and static transformations in sensor space. NeuroImage, 60(1):476–488, 2012. doi:10.1016/j.neuroimage.2011.11.084.

[3] (1,2)

Stefan Haufe, Frank Meinecke, Kai Görgen, Sven Dähne, John-Dylan Haynes, Benjamin Blankertz, and Felix Bießmann. On the interpretation of weight vectors of linear models in multivariate neuroimaging. NeuroImage, 87:96–110, 2014. doi:10.1016/j.neuroimage.2013.10.067.