Compute source power using DICS beamformer

Compute a Dynamic Imaging of Coherent Sources (DICS) [1] filter from single-trial activity to estimate source power across a frequency band. This example demonstrates how to source localize the event-related synchronization (ERS) of beta band activity in the somato dataset.

# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
#         Roman Goj <roman.goj@gmail.com>
#         Denis Engemann <denis.engemann@gmail.com>
#         Stefan Appelhoff <stefan.appelhoff@mailbox.org>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import numpy as np

import mne
from mne.beamformer import apply_dics_csd, make_dics
from mne.datasets import somato
from mne.time_frequency import csd_morlet

print(__doc__)

Reading the raw data and creating epochs:

data_path = somato.data_path()
subject = "01"
task = "somato"
raw_fname = data_path / f"sub-{subject}" / "meg" / f"sub-{subject}_task-{task}_meg.fif"

# Use a shorter segment of raw just for speed here
raw = mne.io.read_raw_fif(raw_fname)
raw.crop(0, 120)  # one minute for speed (looks similar to using all ~800 s)

# Read epochs
events = mne.find_events(raw)

epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, preload=True)
del raw

# Paths to forward operator and FreeSurfer subject directory
fname_fwd = (
    data_path / "derivatives" / f"sub-{subject}" / f"sub-{subject}_task-{task}-fwd.fif"
)

subjects_dir = data_path / "derivatives" / "freesurfer" / "subjects"

Opening raw data file /home/circleci/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
    Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
Ready.
15 events found on stim channel STI 014
Event IDs: [1]
Not setting metadata
15 matching events found
Setting baseline interval to [-1.498464098687909, 0.0] s
Applying baseline correction (mode: mean)
0 projection items activated
Loading data for 15 events and 1052 original time points ...
0 bad epochs dropped

We are interested in the beta band. Define a range of frequencies, using a log scale, from 12 to 30 Hz.

freqs = np.logspace(np.log10(12), np.log10(30), 9)

Computing the cross-spectral density matrix for the beta frequency band, for different time intervals. We use a decim value of 20 to speed up the computation in this example at the loss of accuracy.

csd = csd_morlet(epochs, freqs, tmin=-1, tmax=1.5, decim=20)
csd_baseline = csd_morlet(epochs, freqs, tmin=-1, tmax=0, decim=20)
# ERS activity starts at 0.5 seconds after stimulus onset
csd_ers = csd_morlet(epochs, freqs, tmin=0.5, tmax=1.5, decim=20)
info = epochs.info
del epochs

Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/15 [00:00<?,       ?it/s]
  7%|▋         | CSD epoch blocks : 1/15 [00:00<00:03,    4.30it/s]
 13%|█▎        | CSD epoch blocks : 2/15 [00:00<00:03,    4.32it/s]
 20%|██        | CSD epoch blocks : 3/15 [00:00<00:02,    4.31it/s]
 27%|██▋       | CSD epoch blocks : 4/15 [00:00<00:02,    4.31it/s]
 33%|███▎      | CSD epoch blocks : 5/15 [00:01<00:02,    4.30it/s]
 40%|████      | CSD epoch blocks : 6/15 [00:01<00:02,    4.31it/s]
 47%|████▋     | CSD epoch blocks : 7/15 [00:01<00:01,    4.31it/s]
 53%|█████▎    | CSD epoch blocks : 8/15 [00:01<00:01,    4.31it/s]
 60%|██████    | CSD epoch blocks : 9/15 [00:02<00:01,    4.31it/s]
 67%|██████▋   | CSD epoch blocks : 10/15 [00:02<00:01,    4.31it/s]
 73%|███████▎  | CSD epoch blocks : 11/15 [00:02<00:00,    4.31it/s]
 80%|████████  | CSD epoch blocks : 12/15 [00:02<00:00,    4.31it/s]
 87%|████████▋ | CSD epoch blocks : 13/15 [00:03<00:00,    4.31it/s]
 93%|█████████▎| CSD epoch blocks : 14/15 [00:03<00:00,    4.31it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:03<00:00,    4.31it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:03<00:00,    4.31it/s]
[done]
Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/15 [00:00<?,       ?it/s]
  7%|▋         | CSD epoch blocks : 1/15 [00:00<00:02,    5.88it/s]
 13%|█▎        | CSD epoch blocks : 2/15 [00:00<00:02,    5.89it/s]
 20%|██        | CSD epoch blocks : 3/15 [00:00<00:02,    5.87it/s]
 27%|██▋       | CSD epoch blocks : 4/15 [00:00<00:01,    5.88it/s]
 33%|███▎      | CSD epoch blocks : 5/15 [00:00<00:01,    5.88it/s]
 40%|████      | CSD epoch blocks : 6/15 [00:01<00:01,    5.88it/s]
 47%|████▋     | CSD epoch blocks : 7/15 [00:01<00:01,    5.88it/s]
 53%|█████▎    | CSD epoch blocks : 8/15 [00:01<00:01,    5.88it/s]
 60%|██████    | CSD epoch blocks : 9/15 [00:01<00:01,    5.87it/s]
 67%|██████▋   | CSD epoch blocks : 10/15 [00:01<00:00,    5.88it/s]
 73%|███████▎  | CSD epoch blocks : 11/15 [00:01<00:00,    5.88it/s]
 80%|████████  | CSD epoch blocks : 12/15 [00:02<00:00,    5.88it/s]
 87%|████████▋ | CSD epoch blocks : 13/15 [00:02<00:00,    5.88it/s]
 93%|█████████▎| CSD epoch blocks : 14/15 [00:02<00:00,    5.88it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:02<00:00,    5.87it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:02<00:00,    5.87it/s]
[done]
Computing cross-spectral density from epochs...

  0%|          | CSD epoch blocks : 0/15 [00:00<?,       ?it/s]
  7%|▋         | CSD epoch blocks : 1/15 [00:00<00:02,    5.85it/s]
 13%|█▎        | CSD epoch blocks : 2/15 [00:00<00:02,    5.87it/s]
 20%|██        | CSD epoch blocks : 3/15 [00:00<00:02,    5.88it/s]
 27%|██▋       | CSD epoch blocks : 4/15 [00:00<00:01,    5.89it/s]
 33%|███▎      | CSD epoch blocks : 5/15 [00:00<00:01,    5.88it/s]
 40%|████      | CSD epoch blocks : 6/15 [00:01<00:01,    5.89it/s]
 47%|████▋     | CSD epoch blocks : 7/15 [00:01<00:01,    5.89it/s]
 53%|█████▎    | CSD epoch blocks : 8/15 [00:01<00:01,    5.89it/s]
 60%|██████    | CSD epoch blocks : 9/15 [00:01<00:01,    5.89it/s]
 67%|██████▋   | CSD epoch blocks : 10/15 [00:01<00:00,    5.90it/s]
 73%|███████▎  | CSD epoch blocks : 11/15 [00:01<00:00,    5.86it/s]
 80%|████████  | CSD epoch blocks : 12/15 [00:02<00:00,    5.78it/s]
 87%|████████▋ | CSD epoch blocks : 13/15 [00:02<00:00,    5.77it/s]
 93%|█████████▎| CSD epoch blocks : 14/15 [00:02<00:00,    5.79it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:02<00:00,    5.79it/s]
100%|██████████| CSD epoch blocks : 15/15 [00:02<00:00,    5.81it/s]
[done]

To compute the source power for a frequency band, rather than each frequency separately, we average the CSD objects across frequencies.

csd = csd.mean()
csd_baseline = csd_baseline.mean()
csd_ers = csd_ers.mean()

Computing DICS spatial filters using the CSD that was computed on the entire timecourse.

fwd = mne.read_forward_solution(fname_fwd)
filters = make_dics(
    info,
    fwd,
    csd,
    noise_csd=csd_baseline,
    pick_ori="max-power",
    reduce_rank=True,
    real_filter=True,
)
del fwd

Reading forward solution from /home/circleci/mne_data/MNE-somato-data/derivatives/sub-01/sub-01_task-somato-fwd.fif...
    Reading a source space...
    [done]
    Reading a source space...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
    Read MEG forward solution (8155 sources, 306 channels, free orientations)
    Source spaces transformed to the forward solution coordinate frame
Identifying common channels ...
Dropped the following channels:
['STI 015', 'STI 006', 'STI 001', 'STI 016', 'EOG 061', 'STI 004', 'STI 005', 'STI 014', 'STI 002', 'STI 003']
Identifying common channels ...
Computing inverse operator with 306 channels.
    306 out of 306 channels remain after picking
Selected 306 channels
Creating the depth weighting matrix...
Whitening the forward solution.
Computing rank from covariance with rank=None
    Using tolerance 3.1e-14 (2.2e-16 eps * 306 dim * 0.46  max singular value)
    Estimated rank (mag + grad): 64
    MEG: rank 64 computed from 306 data channels with 0 projectors
    Setting small MEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing rank from covariance with rank=None
    Using tolerance 3.1e-14 (2.2e-16 eps * 306 dim * 0.46  max singular value)
    Estimated rank (mag + grad): 64
    MEG: rank 64 computed from 306 data channels with 0 projectors
Computing rank from covariance with rank=None
    Using tolerance 3.2e-14 (2.2e-16 eps * 306 dim * 0.47  max singular value)
    Estimated rank (mag + grad): 64
    MEG: rank 64 computed from 306 data channels with 0 projectors
Computing DICS spatial filters...
Computing beamformer filters for 8155 sources
Filter computation complete

Applying DICS spatial filters separately to the CSD computed using the baseline and the CSD computed during the ERS activity.

baseline_source_power, freqs = apply_dics_csd(csd_baseline, filters)
beta_source_power, freqs = apply_dics_csd(csd_ers, filters)

Computing DICS source power...
[done]
Computing DICS source power...
[done]

Visualizing source power during ERS activity relative to the baseline power.

stc = beta_source_power / baseline_source_power
message = "DICS source power in the 12-30 Hz frequency band"
brain = stc.plot(
    hemi="both",
    views="axial",
    subjects_dir=subjects_dir,
    subject=subject,
    time_label=message,
)

Using control points [1.40188307 1.5102619  2.34636424]

References

[1]

Joachim Groß, Jan Kujala, Matti S. Hämäläinen, Lars Timmermann, Alfons Schnitzler, and Riitta Salmelin. Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proceedings of the National Academy of Sciences, 98(2):694–699, 2001. doi:10.1073/pnas.98.2.694.