Compute point-spread functions (PSFs) for MNE/dSPM/sLORETA

PSFs are computed for four labels in the MNE sample data set for linear inverse operators (MNE, dSPM, sLORETA). PSFs describe the spread of activation from one label across the cortical surface.

  • ../../_images/sphx_glr_plot_mne_point_spread_function_001.png
  • ../../_images/sphx_glr_plot_mne_point_spread_function_002.png

Out:

Reading forward solution from /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif...
    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
    Desired named matrix (kind = 3523) not available
    Read MEG forward solution (7498 sources, 306 channels, free orientations)
    Desired named matrix (kind = 3523) not available
    Read EEG forward solution (7498 sources, 60 channels, free orientations)
    MEG and EEG forward solutions combined
    Source spaces transformed to the forward solution coordinate frame
    Converting to surface-based source orientations...
    Average patch normals will be employed in the rotation to the local surface coordinates....
[done]
Reading inverse operator decomposition from /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif...
    Reading inverse operator info...
    [done]
    Reading inverse operator decomposition...
    [done]
    364 x 364 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
Reading inverse operator decomposition from /home/ubuntu/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
About to process 4 labels
    Converting to surface-based source orientations...
    Average patch normals will be employed in the rotation to the local surface coordinates....
[done]
<Label  |  unknown, u'Aud-rh', rh : 883 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 1978.94911214  1218.86857762   528.31045019   157.61122813    75.27112997]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 68.5% variance in label.
<Label  |  unknown, u'Aud-lh', lh : 1097 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 1760.92377927  1439.64644748   743.05280713   156.72701522    80.65306411]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 53.8% variance in label.
<Label  |  unknown, u'Vis-rh', rh : 623 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 979.45596116  579.38306825  218.70036716   45.46351293   14.07251956]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 71.3% variance in label.
<Label  |  unknown, u'Vis-lh', lh : 1253 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 1261.76052917   659.9069819    559.49851826    70.12488601    44.44986223]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 67.7% variance in label.
About to apply inverse operator for method='MNE' and lambda2=0.111111111111
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 full noise covariance matrix (4 small eigenvalues omitted)
Picked 364 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
[done]
About to process 4 labels
    Converting to surface-based source orientations...
    Average patch normals will be employed in the rotation to the local surface coordinates....
[done]
<Label  |  unknown, u'Aud-rh', rh : 883 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 0.00274084  0.00179732  0.00066999  0.00021198  0.00018807]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 66.4% variance in label.
<Label  |  unknown, u'Aud-lh', lh : 1097 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 0.00223286  0.00095637  0.00044671  0.0002232   0.00012375]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 80.6% variance in label.
<Label  |  unknown, u'Vis-rh', rh : 623 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [  1.19709958e-03   8.38498995e-04   2.75738170e-04   1.07129273e-04
   9.00191185e-05]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 64.1% variance in label.
<Label  |  unknown, u'Vis-lh', lh : 1253 vertices>
Computing SVD within labels, using 1 component(s)
First 5 singular values: [ 0.0015879   0.00091905  0.00027931  0.0001392   0.00012593]
(This tells you something about variability of forward solutions in sub-leadfield for label)
Your 1 component(s) explain(s) 72.1% variance in label.
About to apply inverse operator for method='MNE' and lambda2=0.111111111111
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 full noise covariance matrix (3 small eigenvalues omitted)
Picked 305 channels from the data
Computing inverse...
(eigenleads need to be weighted)...
[done]
Updating smoothing matrix, be patient..
Smoothing matrix creation, step 1
Smoothing matrix creation, step 2
Smoothing matrix creation, step 3
Smoothing matrix creation, step 4
Smoothing matrix creation, step 5
Smoothing matrix creation, step 6
Smoothing matrix creation, step 7
Smoothing matrix creation, step 8
Smoothing matrix creation, step 9
Smoothing matrix creation, step 10
colormap: fmin=-3.61e-02 fmid=0.00e+00 fmax=3.61e-02 transparent=0
Updating smoothing matrix, be patient..
Smoothing matrix creation, step 1
Smoothing matrix creation, step 2
Smoothing matrix creation, step 3
Smoothing matrix creation, step 4
Smoothing matrix creation, step 5
Smoothing matrix creation, step 6
Smoothing matrix creation, step 7
Smoothing matrix creation, step 8
Smoothing matrix creation, step 9
Smoothing matrix creation, step 10
colormap: fmin=-2.52e-02 fmid=0.00e+00 fmax=2.52e-02 transparent=0

# Authors: Olaf Hauk <olaf.hauk@mrc-cbu.cam.ac.uk>
#          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#
# License: BSD (3-clause)

from mayavi import mlab

import mne
from mne.datasets import sample
from mne.minimum_norm import read_inverse_operator, point_spread_function

print(__doc__)

data_path = sample.data_path()
subjects_dir = data_path + '/subjects/'
fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
fname_inv_eegmeg = (data_path +
                    '/MEG/sample/sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif')
fname_inv_meg = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_label = [data_path + '/MEG/sample/labels/Aud-rh.label',
               data_path + '/MEG/sample/labels/Aud-lh.label',
               data_path + '/MEG/sample/labels/Vis-rh.label',
               data_path + '/MEG/sample/labels/Vis-lh.label']


# read forward solution
forward = mne.read_forward_solution(fname_fwd, force_fixed=False,
                                    surf_ori=True)

# read inverse operators
inverse_operator_eegmeg = read_inverse_operator(fname_inv_eegmeg)
inverse_operator_meg = read_inverse_operator(fname_inv_meg)

# read label(s)
labels = [mne.read_label(ss) for ss in fname_label]

# regularisation parameter
snr = 3.0
lambda2 = 1.0 / snr ** 2
method = 'MNE'  # can be 'MNE' or 'sLORETA'
mode = 'svd'
n_svd_comp = 1

stc_psf_eegmeg, _ = point_spread_function(
    inverse_operator_eegmeg, forward, method=method, labels=labels,
    lambda2=lambda2, pick_ori='normal', mode=mode, n_svd_comp=n_svd_comp)

stc_psf_meg, _ = point_spread_function(
    inverse_operator_meg, forward, method=method, labels=labels,
    lambda2=lambda2, pick_ori='normal', mode=mode, n_svd_comp=n_svd_comp)

# save for viewing in mne_analyze in order of labels in 'labels'
# last sample is average across PSFs
# stc_psf_eegmeg.save('psf_eegmeg')
# stc_psf_meg.save('psf_meg')

time_label = "EEGMEG %d"
brain_eegmeg = stc_psf_eegmeg.plot(hemi='rh', subjects_dir=subjects_dir,
                                   time_label=time_label,
                                   figure=mlab.figure(size=(500, 500)))

time_label = "MEG %d"
brain_meg = stc_psf_meg.plot(hemi='rh', subjects_dir=subjects_dir,
                             time_label=time_label,
                             figure=mlab.figure(size=(500, 500)))

# The PSF is centred around the right auditory cortex label,
# but clearly extends beyond it.
# It also contains "sidelobes" or "ghost sources"
# in middle/superior temporal lobe.
# For the Aud-RH example, MEG and EEGMEG do not seem to differ a lot,
# but the addition of EEG still decreases point-spread to distant areas
# (e.g. to ATL and IFG).
# The chosen labels are quite far apart from each other, so their PSFs
# do not overlap (check in mne_analyze)

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

Generated by Sphinx-Gallery