Compute Phase Slope Index (PSI) in source space for a visual stimulus#

This example demonstrates how the phase slope index (PSI) [1] can be computed in source space based on single trial dSPM source estimates. In addition, the example shows advanced usage of the connectivity estimation routines by first extracting a label time course for each epoch and then combining the label time course with the single trial source estimates to compute the connectivity.

The result clearly shows how the activity in the visual label precedes more widespread activity (as a postivive PSI means the label time course is leading).

# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#
# License: BSD (3-clause)

import mne
import numpy as np
from mne.datasets import sample
from mne.minimum_norm import apply_inverse_epochs, read_inverse_operator

from mne_connectivity import phase_slope_index, seed_target_indices

print(__doc__)

data_path = sample.data_path()
subjects_dir = data_path / "subjects"
fname_inv = data_path / "MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif"
fname_raw = data_path / "MEG/sample/sample_audvis_filt-0-40_raw.fif"
fname_event = data_path / "MEG/sample/sample_audvis_filt-0-40_raw-eve.fif"
fname_label = data_path / "MEG/sample/labels/Vis-lh.label"

event_id, tmin, tmax = 4, -0.2, 0.5
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)

# Load data
inverse_operator = read_inverse_operator(fname_inv)
raw = mne.io.read_raw_fif(fname_raw, preload=True)
events = mne.read_events(fname_event)

# pick MEG channels
picks = mne.pick_types(
    raw.info, meg=True, eeg=False, stim=False, eog=True, exclude="bads"
)

# Read epochs
epochs = mne.Epochs(
    raw,
    events,
    event_id,
    tmin,
    tmax,
    picks=picks,
    baseline=(None, 0),
    reject=dict(mag=4e-12, grad=4000e-13, eog=150e-6),
)

# Compute inverse solution and for each epoch. Note that since we are passing
# the output to both extract_label_time_course and the phase_slope_index
# functions, we have to use "return_generator=False", since it is only possible
# to iterate over generators once.
snr = 1.0  # use lower SNR for single epochs
lambda2 = 1.0 / snr**2
stcs = apply_inverse_epochs(
    epochs, inverse_operator, lambda2, method, pick_ori="normal", return_generator=True
)

# Now, we generate seed time series by averaging the activity in the left
# visual corex
label = mne.read_label(fname_label)
src = inverse_operator["src"]  # the source space used
seed_ts = mne.extract_label_time_course(
    stcs, label, src, mode="mean_flip", verbose="error"
)

# Combine the seed time course with the source estimates. There will be a total
# of 7500 signals:
# index 0: time course extracted from label
# index 1..7499: dSPM source space time courses
stcs = apply_inverse_epochs(
    epochs, inverse_operator, lambda2, method, pick_ori="normal", return_generator=True
)
comb_ts = list(zip(seed_ts, stcs))

# Construct indices to estimate connectivity between the label time course
# and all source space time courses
vertices = [src[i]["vertno"] for i in range(2)]
n_signals_tot = 1 + len(vertices[0]) + len(vertices[1])

indices = seed_target_indices([0], np.arange(1, n_signals_tot))

# Compute the PSI in the frequency range 10Hz-20Hz. We exclude the baseline
# period from the connectivity estimation.
fmin = 10.0
fmax = 20.0
tmin_con = 0.0
sfreq = epochs.info["sfreq"]  # the sampling frequency

psi = phase_slope_index(
    comb_ts,
    mode="multitaper",
    indices=indices,
    sfreq=sfreq,
    fmin=fmin,
    fmax=fmax,
    tmin=tmin_con,
)

# Generate a SourceEstimate with the PSI. This is simple since we used a single
# seed (inspect the indices variable to see how the PSI scores are arranged in
# the output)
psi_stc = mne.SourceEstimate(
    psi.get_data(), vertices=vertices, tmin=0, tstep=1, subject="sample"
)

# Now we can visualize the PSI using the :meth:`~mne.SourceEstimate.plot`
# method. We use a custom colormap to show signed values
v_max = np.max(np.abs(psi.get_data()))
brain = psi_stc.plot(
    surface="inflated",
    hemi="lh",
    time_label="Phase Slope Index (PSI)",
    subjects_dir=subjects_dir,
    clim=dict(kind="percent", pos_lims=(95, 97.5, 100)),
)
brain.show_view("medial")
brain.add_label(str(fname_label), color="green", alpha=0.7)
mne inverse psi visual
Reading inverse operator decomposition from /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif...
    Reading inverse operator info...
    [done]
    Reading inverse operator decomposition...
    [done]
    305 x 305 full covariance (kind = 1) found.
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Noise covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 2) found.
    Source covariance matrix read.
    22494 x 22494 diagonal covariance (kind = 6) found.
    Orientation priors read.
    22494 x 22494 diagonal covariance (kind = 5) found.
    Depth priors read.
    Did not find the desired covariance matrix (kind = 3)
    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
    Read a total of 4 projection items:
        PCA-v1 (1 x 102) active
        PCA-v2 (1 x 102) active
        PCA-v3 (1 x 102) active
        Average EEG reference (1 x 60) active
    Source spaces transformed to the inverse solution coordinate frame
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
Not setting metadata
70 matching events found
Setting baseline interval to [-0.19979521315838786, 0.0] s
Applying baseline correction (mode: mean)
Created an SSP operator (subspace dimension = 3)
3 projection items activated
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 = 3)
    Created the whitener using a noise covariance matrix with rank 302 (3 small eigenvalues omitted)
    Computing noise-normalization factors (dSPM)...
[done]
Picked 305 channels from the data
Computing inverse...
    Eigenleads need to be weighted ...
Processing epoch : 1 / 70 (at most)
Processing epoch : 2 / 70 (at most)
Processing epoch : 3 / 70 (at most)
Processing epoch : 4 / 70 (at most)
Processing epoch : 5 / 70 (at most)
Processing epoch : 6 / 70 (at most)
Processing epoch : 7 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 8 / 70 (at most)
Processing epoch : 9 / 70 (at most)
Processing epoch : 10 / 70 (at most)
Processing epoch : 11 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 12 / 70 (at most)
Processing epoch : 13 / 70 (at most)
Processing epoch : 14 / 70 (at most)
Processing epoch : 15 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 16 / 70 (at most)
Processing epoch : 17 / 70 (at most)
Processing epoch : 18 / 70 (at most)
Processing epoch : 19 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 20 / 70 (at most)
Processing epoch : 21 / 70 (at most)
Processing epoch : 22 / 70 (at most)
Processing epoch : 23 / 70 (at most)
Processing epoch : 24 / 70 (at most)
Processing epoch : 25 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 26 / 70 (at most)
Processing epoch : 27 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 28 / 70 (at most)
Processing epoch : 29 / 70 (at most)
Processing epoch : 30 / 70 (at most)
Processing epoch : 31 / 70 (at most)
Processing epoch : 32 / 70 (at most)
Processing epoch : 33 / 70 (at most)
Processing epoch : 34 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 35 / 70 (at most)
Processing epoch : 36 / 70 (at most)
Processing epoch : 37 / 70 (at most)
Processing epoch : 38 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 39 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 40 / 70 (at most)
Processing epoch : 41 / 70 (at most)
Processing epoch : 42 / 70 (at most)
Processing epoch : 43 / 70 (at most)
Processing epoch : 44 / 70 (at most)
Processing epoch : 45 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 46 / 70 (at most)
Processing epoch : 47 / 70 (at most)
Processing epoch : 48 / 70 (at most)
Processing epoch : 49 / 70 (at most)
Processing epoch : 50 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 51 / 70 (at most)
Processing epoch : 52 / 70 (at most)
Processing epoch : 53 / 70 (at most)
Processing epoch : 54 / 70 (at most)
Processing epoch : 55 / 70 (at most)
    Rejecting  epoch based on EOG : ['EOG 061']
Processing epoch : 56 / 70 (at most)
Estimating phase slope index (PSI)
Connectivity computation...
    using t=0.000s..0.499s for estimation (76 points)
    frequencies: 11.9Hz..19.8Hz (5 points)
    computing connectivity for 7498 connections
    Using multitaper spectrum estimation with 7 DPSS windows
    the following metrics will be computed: Coherency
    computing cross-spectral density for epoch 1
    computing cross-spectral density for epoch 2
    computing cross-spectral density for epoch 3
    computing cross-spectral density for epoch 4
    computing cross-spectral density for epoch 5
    computing cross-spectral density for epoch 6
    computing cross-spectral density for epoch 7
    computing cross-spectral density for epoch 8
    computing cross-spectral density for epoch 9
    computing cross-spectral density for epoch 10
    computing cross-spectral density for epoch 11
    computing cross-spectral density for epoch 12
    computing cross-spectral density for epoch 13
    computing cross-spectral density for epoch 14
    computing cross-spectral density for epoch 15
    computing cross-spectral density for epoch 16
    computing cross-spectral density for epoch 17
    computing cross-spectral density for epoch 18
    computing cross-spectral density for epoch 19
    computing cross-spectral density for epoch 20
    computing cross-spectral density for epoch 21
    computing cross-spectral density for epoch 22
    computing cross-spectral density for epoch 23
    computing cross-spectral density for epoch 24
    computing cross-spectral density for epoch 25
    computing cross-spectral density for epoch 26
    computing cross-spectral density for epoch 27
    computing cross-spectral density for epoch 28
    computing cross-spectral density for epoch 29
    computing cross-spectral density for epoch 30
    computing cross-spectral density for epoch 31
    computing cross-spectral density for epoch 32
    computing cross-spectral density for epoch 33
    computing cross-spectral density for epoch 34
    computing cross-spectral density for epoch 35
    computing cross-spectral density for epoch 36
    computing cross-spectral density for epoch 37
    computing cross-spectral density for epoch 38
    computing cross-spectral density for epoch 39
    computing cross-spectral density for epoch 40
    computing cross-spectral density for epoch 41
    computing cross-spectral density for epoch 42
    computing cross-spectral density for epoch 43
    computing cross-spectral density for epoch 44
    computing cross-spectral density for epoch 45
    computing cross-spectral density for epoch 46
    computing cross-spectral density for epoch 47
    computing cross-spectral density for epoch 48
    computing cross-spectral density for epoch 49
    computing cross-spectral density for epoch 50
    computing cross-spectral density for epoch 51
    computing cross-spectral density for epoch 52
    computing cross-spectral density for epoch 53
    computing cross-spectral density for epoch 54
    computing cross-spectral density for epoch 55
    computing cross-spectral density for epoch 56
[Connectivity computation done]
Computing PSI from estimated Coherency: <SpectralConnectivity | freq : [11.854243228310033, 19.757072047183385], , nave : 56, nodes, n_estimated : 7499, 7498, ~1.2 MiB>
[PSI Estimation Done]
Using control points [0.03622925 0.05660114 0.13862278]

References#

Total running time of the script: (0 minutes 16.163 seconds)

Gallery generated by Sphinx-Gallery