Compute envelope correlations in source space

Compute envelope correlations of orthogonalized activity [1][2] using pairwise and symmetric orthogonalization [3] in source space using resting state CTF data.

Note that the original procedure for symmetric orthogonalization in [3] is:

1. Extract inverse label data from raw

2. Symmetric orthogonalization

3. Band-pass filter

4. Hilbert transform and absolute value

5. Low-pass (1 Hz)

Here we follow the procedure:

1. Epoch data, then for each

2. Extract inverse label data for each epoch

3. Symmetric orthogonalization for each epoch

4. Band-pass filter each epoch

5. Hilbert transform and absolute value (inside envelope_correlation)

The differences between these two should hopefully be fairly minimal given the pairwise orthogonalization used in [2] used a similar pipeline.

# Authors: Eric Larson <larson.eric.d@gmail.com>
#          Sheraz Khan <sheraz@khansheraz.com>
#          Denis Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)

import os.path as op

import numpy as np
import matplotlib.pyplot as plt

import mne
import mne_connectivity
from mne_connectivity import envelope_correlation
from mne.minimum_norm import make_inverse_operator, apply_inverse_epochs
from mne.preprocessing import compute_proj_ecg, compute_proj_eog

data_path = mne.datasets.brainstorm.bst_resting.data_path()
subjects_dir = op.join(data_path, "subjects")
subject = "bst_resting"
trans = op.join(data_path, "MEG", "bst_resting", "bst_resting-trans.fif")
src = op.join(subjects_dir, subject, "bem", subject + "-oct-6-src.fif")
bem = op.join(subjects_dir, subject, "bem", subject + "-5120-bem-sol.fif")
raw_fname = op.join(
    data_path, "MEG", "bst_resting", "subj002_spontaneous_20111102_01_AUX.ds"
)

Here we do some things in the name of speed, such as crop (which will hurt SNR) and downsample. Then we compute SSP projectors and apply them.

raw = mne.io.read_raw_ctf(raw_fname, verbose="error")
raw.crop(0, 60).pick_types(meg=True, eeg=False).load_data().resample(80)
raw.apply_gradient_compensation(3)
projs_ecg, _ = compute_proj_ecg(raw, n_grad=1, n_mag=2)
projs_eog, _ = compute_proj_eog(raw, n_grad=1, n_mag=2, ch_name="MLT31-4407")
raw.add_proj(projs_ecg + projs_eog)
raw.apply_proj()
raw.filter(0.1, None)  # this helps with symmetric orthogonalization later
cov = mne.compute_raw_covariance(raw)  # compute before band-pass of interest

NOTE: pick_types() is a legacy function. New code should use inst.pick(...).
Reading 0 ... 144000  =      0.000 ...    60.000 secs...
Including 0 SSP projectors from raw file
Running ECG SSP computation
Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 5 - 35 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hann window
- Lower passband edge: 5.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 4.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 800 samples (10.000 s)

Number of ECG events detected : 88 (average pulse 88 / min.)
Computing projector
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 35 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 800 samples (10.000 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
Not setting metadata
88 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 88 events and 49 original time points ...
    Rejecting  epoch based on MAG : ['MLT31-4407', 'MRT31-4407']
    Rejecting  epoch based on MAG : ['MLT31-4407', 'MLT41-4407', 'MRT31-4407', 'MRT41-4407']
    Rejecting  epoch based on MAG : ['MLT31-4407', 'MLT41-4407', 'MRT31-4407', 'MRT41-4407']
4 bad epochs dropped
No channels 'grad' found. Skipping.
Adding projection: axial--0.200-0.400-PCA-01
Adding projection: axial--0.200-0.400-PCA-02
No channels 'eeg' found. Skipping.
Done.
Including 0 SSP projectors from raw file
Running EOG SSP computation
Using EOG channel: MLT31-4407
EOG channel index for this subject is: [137]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Selecting channel MLT31-4407 for blink detection
Setting up band-pass filter from 1 - 10 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hann window
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz)
- Upper passband edge: 10.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 10.25 Hz)
- Filter length: 800 samples (10.000 s)

Now detecting blinks and generating corresponding events
Found 12 significant peaks
Number of EOG events detected: 12
Computing projector
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 35 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 800 samples (10.000 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
Not setting metadata
12 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 12 events and 33 original time points ...
    Rejecting  epoch based on MAG : ['MRT41-4407']
1 bad epochs dropped
No channels 'grad' found. Skipping.
Adding projection: axial--0.200-0.200-PCA-01
Adding projection: axial--0.200-0.200-PCA-02
No channels 'eeg' found. Skipping.
Done.
4 projection items deactivated
Created an SSP operator (subspace dimension = 4)
4 projection items activated
SSP projectors applied...
Filtering raw data in 1 contiguous segment
Setting up high-pass filter at 0.1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.10
- Lower transition bandwidth: 0.10 Hz (-6 dB cutoff frequency: 0.05 Hz)
- Filter length: 2641 samples (33.013 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.2s
Using up to 300 segments
Number of samples used : 4800
[done]

Compute the forward and inverse

src = mne.read_source_spaces(src)
fwd = mne.make_forward_solution(raw.info, trans, src, bem, verbose=True)
del src
inv = make_inverse_operator(raw.info, fwd, cov)
del fwd

    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    Reading a source space...
    Computing patch statistics...
    Patch information added...
    Distance information added...
    [done]
    2 source spaces read
Source space          : <SourceSpaces: [<surface (lh), n_vertices=156869, n_used=4098>, <surface (rh), n_vertices=155654, n_used=4098>] MRI (surface RAS) coords, subject 'bst_resting', ~27.5 MB>
MRI -> head transform : /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/MEG/bst_resting/bst_resting-trans.fif
Measurement data      : instance of Info
Conductor model   : /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/bem/bst_resting-5120-bem-sol.fif
Accurate field computations
Do computations in head coordinates
Free source orientations

Read 2 source spaces a total of 8196 active source locations

Coordinate transformation: MRI (surface RAS) -> head
     0.999797 -0.005775 -0.019288       2.71 mm
     0.011390  0.952195  0.305279      16.66 mm
     0.016602 -0.305437  0.952068      28.47 mm
     0.000000  0.000000  0.000000       1.00

Read 298 MEG channels from info
Read 26 MEG compensation channels from info
5 compensation data sets in info
Setting up compensation data...
    Desired compensation data (3) found.
    All compensation channels found.
    Preselector created.
    Compensation data matrix created.
    Postselector created.
105 coil definitions read
Coordinate transformation: MEG device -> head
     0.998490 -0.050225 -0.022235       1.90 mm
     0.052235  0.993447  0.101656      13.13 mm
     0.016984 -0.102664  0.994571      66.69 mm
     0.000000  0.000000  0.000000       1.00
MEG coil definitions created in head coordinates.
Removing 5 compensators from info because not all compensation channels were picked.
Source spaces are now in head coordinates.

Setting up the BEM model using /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/bem/bst_resting-5120-bem-sol.fif...

Loading surfaces...

Loading the solution matrix...

Homogeneous model surface loaded.
Loaded linear collocation BEM solution from /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/bem/bst_resting-5120-bem-sol.fif
Employing the head->MRI coordinate transform with the BEM model.
BEM model bst_resting-5120-bem-sol.fif is now set up

Source spaces are in head coordinates.
Checking that the sources are inside the surface (will take a few...)
Checking surface interior status for 4098 points...
    Found 1239/4098 points inside  an interior sphere of radius   49.6 mm
    Found    0/4098 points outside an exterior sphere of radius  105.3 mm
    Found    0/2859 points outside using surface Qhull
    Found    0/2859 points outside using solid angles
    Total 4098/4098 points inside the surface
Interior check completed in 3039.2 ms
Checking surface interior status for 4098 points...
    Found 1235/4098 points inside  an interior sphere of radius   49.6 mm
    Found    0/4098 points outside an exterior sphere of radius  105.3 mm
    Found    0/2863 points outside using surface Qhull
    Found    0/2863 points outside using solid angles
    Total 4098/4098 points inside the surface
Interior check completed in 3127.1 ms

Checking surface interior status for 298 points...
    Found   0/298 points inside  an interior sphere of radius   49.6 mm
    Found 283/298 points outside an exterior sphere of radius  105.3 mm
    Found  15/ 15 points outside using surface Qhull
    Found   0/  0 points outside using solid angles
    Total 0/298 points inside the surface
Interior check completed in 216.9 ms

Composing the field computation matrix...
Computing MEG at 8196 source locations (free orientations)...

Finished.
Converting forward solution to surface orientation
    Average patch normals will be employed in the rotation to the local surface coordinates....
    Converting to surface-based source orientations...
    [done]
Computing inverse operator with 272 channels.
    272 out of 272 channels remain after picking
Removing 5 compensators from info because not all compensation channels were picked.
Selected 272 channels
Creating the depth weighting matrix...
    272 magnetometer or axial gradiometer channels
    limit = 8074/8196 = 10.006553
    scale = 5.92189e-11 exp = 0.8
Applying loose dipole orientations to surface source spaces: 0.2
Whitening the forward solution.
Removing 5 compensators from info because not all compensation channels were picked.
    Created an SSP operator (subspace dimension = 4)
Computing rank from covariance with rank=None
    Using tolerance 4e-14 (2.2e-16 eps * 272 dim * 0.66  max singular value)
    Estimated rank (mag): 268
    MAG: rank 268 computed from 272 data channels with 4 projectors
    Setting small MAG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing SVD of whitened and weighted lead field matrix.
    largest singular value = 3.41214
    scaling factor to adjust the trace = 1.25505e+19 (nchan = 272 nzero = 4)

Now we create epochs and prepare to band-pass filter them.

duration = 10.0
events = mne.make_fixed_length_events(raw, duration=duration)
tmax = duration - 1.0 / raw.info["sfreq"]
epochs = mne.Epochs(
    raw, events=events, tmin=0, tmax=tmax, baseline=None, reject=dict(mag=20e-13)
)
sfreq = epochs.info["sfreq"]
del raw, projs_ecg, projs_eog

Not setting metadata
6 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 4)
4 projection items activated

Do pairwise-orthogonalized envelope correlation

# sphinx_gallery_thumbnail_number = 2

labels = mne.read_labels_from_annot(subject, "aparc_sub", subjects_dir=subjects_dir)
stcs = apply_inverse_epochs(
    epochs, inv, lambda2=1.0 / 9.0, pick_ori="normal", return_generator=True
)
label_ts = mne.extract_label_time_course(
    stcs, labels, inv["src"], return_generator=False
)
del stcs


def bp_gen(label_ts):
    """Make a generator that band-passes on the fly."""
    for ts in label_ts:
        yield mne.filter.filter_data(ts, sfreq, 14, 30)


corr_obj = envelope_correlation(bp_gen(label_ts), orthogonalize="pairwise")
corr = corr_obj.combine()
corr = corr.get_data(output="dense")[:, :, 0]


def plot_corr(corr, title):
    fig, ax = plt.subplots(figsize=(4, 4), constrained_layout=True)
    ax.imshow(corr, cmap="viridis", clim=np.percentile(corr, [5, 95]))
    fig.suptitle(title)


plot_corr(corr, "Pairwise")


def plot_degree(corr, title):
    threshold_prop = 0.15  # percentage of strongest edges to keep in the graph
    degree = mne_connectivity.degree(corr, threshold_prop=threshold_prop)
    stc = mne.labels_to_stc(labels, degree)
    stc = stc.in_label(
        mne.Label(inv["src"][0]["vertno"], hemi="lh")
        + mne.Label(inv["src"][1]["vertno"], hemi="rh")
    )
    return stc.plot(
        clim=dict(kind="percent", lims=[75, 85, 95]),
        colormap="gnuplot",
        subjects_dir=subjects_dir,
        views="dorsal",
        hemi="both",
        smoothing_steps=25,
        time_label=title,
    )


brain = plot_degree(corr, "Beta (pairwise, aparc_sub)")

Reading labels from parcellation...
   read 226 labels from /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/label/lh.aparc_sub.annot
   read 224 labels from /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/label/rh.aparc_sub.annot
Removing 5 compensators from info because not all compensation channels were picked.
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 1
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 4)
    Created the whitener using a noise covariance matrix with rank 268 (4 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 272 channels from the data
Computing inverse...
    Eigenleads need to be weighted ...
Processing epoch : 1 / 6 (at most)
Extracting time courses for 450 labels (mode: mean_flip)
Processing epoch : 2 / 6 (at most)
Extracting time courses for 450 labels (mode: mean_flip)
Processing epoch : 3 / 6 (at most)
Extracting time courses for 450 labels (mode: mean_flip)
Processing epoch : 4 / 6 (at most)
Extracting time courses for 450 labels (mode: mean_flip)
Processing epoch : 5 / 6 (at most)
Extracting time courses for 450 labels (mode: mean_flip)
Processing epoch : 6 / 6 (at most)
Extracting time courses for 450 labels (mode: mean_flip)
[done]
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 449 tasks      | elapsed:    0.1s
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 449 tasks      | elapsed:    0.1s
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 449 tasks      | elapsed:    0.1s
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 449 tasks      | elapsed:    0.1s
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 449 tasks      | elapsed:    0.1s
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 161 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done 287 tasks      | elapsed:    0.1s
[Parallel(n_jobs=1)]: Done 449 tasks      | elapsed:    0.1s
Using control points [73. 80. 89.]

Do symmetric-orthogonalized envelope correlation

Here we need the number of labels to be less than the rank of the data (here around 200), because all label time courses are orthogonalized relative to one another. 'aparc_sub' has over 400 labels, so here we use 'aparc.a2009s', which has fewer than 200.

labels = mne.read_labels_from_annot(subject, "aparc.a2009s", subjects_dir=subjects_dir)
stcs = apply_inverse_epochs(
    epochs, inv, lambda2=1.0 / 9.0, pick_ori="normal", return_generator=True
)
label_ts = mne.extract_label_time_course(
    stcs, labels, inv["src"], return_generator=True
)
del stcs, epochs

label_ts_orth = mne_connectivity.envelope.symmetric_orth(label_ts)
corr_obj = envelope_correlation(  # already orthogonalized earlier
    bp_gen(label_ts_orth), orthogonalize=False
)

# average over epochs, take absolute value, and plot
corr = corr_obj.combine()
corr = corr.get_data(output="dense")[:, :, 0]
corr.flat[:: corr.shape[0] + 1] = 0  # zero out the diagonal
corr = np.abs(corr)

plot_corr(corr, "Symmetric")
plot_degree(corr, "Beta (symmetric, aparc.a2009s)")

Reading labels from parcellation...
   read 75 labels from /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/label/lh.aparc.a2009s.annot
   read 75 labels from /home/circleci/mne_data/MNE-brainstorm-data/bst_resting/subjects/bst_resting/label/rh.aparc.a2009s.annot
Removing 5 compensators from info because not all compensation channels were picked.
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 1
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 4)
    Created the whitener using a noise covariance matrix with rank 268 (4 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 272 channels from the data
Computing inverse...
    Eigenleads need to be weighted ...
Processing epoch : 1 / 6 (at most)
Extracting time courses for 150 labels (mode: mean_flip)
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
Processing epoch : 2 / 6 (at most)
Extracting time courses for 150 labels (mode: mean_flip)
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
Processing epoch : 3 / 6 (at most)
Extracting time courses for 150 labels (mode: mean_flip)
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
Processing epoch : 4 / 6 (at most)
Extracting time courses for 150 labels (mode: mean_flip)
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
Processing epoch : 5 / 6 (at most)
Extracting time courses for 150 labels (mode: mean_flip)
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
Processing epoch : 6 / 6 (at most)
Extracting time courses for 150 labels (mode: mean_flip)
Setting up band-pass filter from 14 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 14.00
- Lower transition bandwidth: 3.50 Hz (-6 dB cutoff frequency: 12.25 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 77 samples (0.963 s)

[Parallel(n_jobs=1)]: Done  17 tasks      | elapsed:    0.0s
[Parallel(n_jobs=1)]: Done  71 tasks      | elapsed:    0.0s
[done]
Using control points [26. 29. 32.]

<mne.viz._brain._brain.Brain object at 0x7f70f3537a30>

References

[1]

Joerg F Hipp, David J Hawellek, Maurizio Corbetta, Markus Siegel, and Andreas K Engel. Large-scale cortical correlation structure of spontaneous oscillatory activity. Nature Neuroscience, 15(6):884–890, 2012. doi:10.1038/nn.3101.

[2] (1,2)

Sheraz Khan, Javeria A. Hashmi, Fahimeh Mamashli, Konstantinos Michmizos, Manfred G. Kitzbichler, Hari Bharadwaj, Yousra Bekhti, Santosh Ganesan, Keri-Lee A. Garel, Susan Whitfield-Gabrieli, Randy L. Gollub, Jian Kong, Lucia M. Vaina, Kunjan D. Rana, Steven M. Stufflebeam, Matti S. Hämäläinen, and Tal Kenet. Maturation trajectories of cortical resting-state networks depend on the mediating frequency band. NeuroImage, 174:57–68, 2018. doi:10.1016/j.neuroimage.2018.02.018.

[3] (1,2)

G. L. Colclough, M. J. Brookes, S. M. Smith, and M. W. Woolrich. A symmetric multivariate leakage correction for MEG connectomes. NeuroImage, 117:439–448, August 2015. doi:10.1016/j.neuroimage.2015.03.071.