Note
Click here to download the full example code
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
import os.path as op
import numpy as np
import mne
from mne.datasets import somato
from mne.time_frequency import csd_morlet
from mne.beamformer import make_dics, apply_dics_csd
print(__doc__)
Reading the raw data and creating epochs:
data_path = somato.data_path()
subject = '01'
task = 'somato'
raw_fname = op.join(data_path, 'sub-{}'.format(subject), 'meg',
'sub-{}_task-{}_meg.fif'.format(subject, task))
# 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 sec)
# 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 = op.join(data_path, 'derivatives', 'sub-{}'.format(subject),
'sub-{}_task-{}-fwd.fif'.format(subject, task))
subjects_dir = op.join(data_path, 'derivatives', 'freesurfer', 'subjects')
Out:
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
Event IDs: [1]
Not setting metadata
Not setting metadata
15 matching events found
Setting baseline interval to [-1.498464098687909, 0.0] sec
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
Out:
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/15 [00:00<?, ?it/s]
7%|6 | CSD epoch blocks : 1/15 [00:00<00:03, 4.45it/s]
13%|#3 | CSD epoch blocks : 2/15 [00:00<00:02, 5.09it/s]
20%|## | CSD epoch blocks : 3/15 [00:00<00:02, 5.36it/s]
27%|##6 | CSD epoch blocks : 4/15 [00:00<00:02, 5.49it/s]
33%|###3 | CSD epoch blocks : 5/15 [00:00<00:01, 5.58it/s]
40%|#### | CSD epoch blocks : 6/15 [00:01<00:01, 5.58it/s]
47%|####6 | CSD epoch blocks : 7/15 [00:01<00:01, 5.61it/s]
53%|#####3 | CSD epoch blocks : 8/15 [00:01<00:01, 5.63it/s]
60%|###### | CSD epoch blocks : 9/15 [00:01<00:01, 5.60it/s]
67%|######6 | CSD epoch blocks : 10/15 [00:01<00:00, 5.63it/s]
73%|#######3 | CSD epoch blocks : 11/15 [00:01<00:00, 5.66it/s]
80%|######## | CSD epoch blocks : 12/15 [00:02<00:00, 5.67it/s]
87%|########6 | CSD epoch blocks : 13/15 [00:02<00:00, 5.67it/s]
93%|#########3| CSD epoch blocks : 14/15 [00:02<00:00, 5.64it/s]
100%|##########| CSD epoch blocks : 15/15 [00:02<00:00, 5.60it/s]
100%|##########| CSD epoch blocks : 15/15 [00:02<00:00, 5.60it/s]
[done]
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/15 [00:00<?, ?it/s]
7%|6 | CSD epoch blocks : 1/15 [00:00<00:01, 8.00it/s]
13%|#3 | CSD epoch blocks : 2/15 [00:00<00:01, 8.29it/s]
20%|## | CSD epoch blocks : 3/15 [00:00<00:01, 8.39it/s]
27%|##6 | CSD epoch blocks : 4/15 [00:00<00:01, 8.38it/s]
33%|###3 | CSD epoch blocks : 5/15 [00:00<00:01, 8.40it/s]
40%|#### | CSD epoch blocks : 6/15 [00:00<00:01, 8.39it/s]
47%|####6 | CSD epoch blocks : 7/15 [00:00<00:00, 8.40it/s]
53%|#####3 | CSD epoch blocks : 8/15 [00:00<00:00, 8.41it/s]
60%|###### | CSD epoch blocks : 9/15 [00:01<00:00, 8.42it/s]
67%|######6 | CSD epoch blocks : 10/15 [00:01<00:00, 8.44it/s]
73%|#######3 | CSD epoch blocks : 11/15 [00:01<00:00, 8.44it/s]
80%|######## | CSD epoch blocks : 12/15 [00:01<00:00, 8.45it/s]
87%|########6 | CSD epoch blocks : 13/15 [00:01<00:00, 8.43it/s]
93%|#########3| CSD epoch blocks : 14/15 [00:01<00:00, 8.29it/s]
100%|##########| CSD epoch blocks : 15/15 [00:01<00:00, 8.29it/s]
100%|##########| CSD epoch blocks : 15/15 [00:01<00:00, 8.32it/s]
[done]
Computing cross-spectral density from epochs...
0%| | CSD epoch blocks : 0/15 [00:00<?, ?it/s]
7%|6 | CSD epoch blocks : 1/15 [00:00<00:01, 8.00it/s]
13%|#3 | CSD epoch blocks : 2/15 [00:00<00:01, 7.97it/s]
20%|## | CSD epoch blocks : 3/15 [00:00<00:01, 7.82it/s]
27%|##6 | CSD epoch blocks : 4/15 [00:00<00:01, 7.94it/s]
33%|###3 | CSD epoch blocks : 5/15 [00:00<00:01, 8.01it/s]
40%|#### | CSD epoch blocks : 6/15 [00:00<00:01, 8.01it/s]
47%|####6 | CSD epoch blocks : 7/15 [00:00<00:00, 8.03it/s]
53%|#####3 | CSD epoch blocks : 8/15 [00:00<00:00, 8.08it/s]
60%|###### | CSD epoch blocks : 9/15 [00:01<00:00, 8.08it/s]
67%|######6 | CSD epoch blocks : 10/15 [00:01<00:00, 8.07it/s]
73%|#######3 | CSD epoch blocks : 11/15 [00:01<00:00, 8.08it/s]
80%|######## | CSD epoch blocks : 12/15 [00:01<00:00, 8.08it/s]
87%|########6 | CSD epoch blocks : 13/15 [00:01<00:00, 8.08it/s]
93%|#########3| CSD epoch blocks : 14/15 [00:01<00:00, 8.11it/s]
100%|##########| CSD epoch blocks : 15/15 [00:01<00:00, 8.13it/s]
100%|##########| CSD epoch blocks : 15/15 [00:01<00:00, 8.11it/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
Out:
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) 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 014', 'STI 004', 'STI 006', 'EOG 061', 'STI 003', 'STI 015', 'STI 016', 'STI 002', 'STI 005', 'STI 001']
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.
Out:
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)
Out:
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.
Total running time of the script: ( 0 minutes 20.937 seconds)
Estimated memory usage: 209 MB