Compute a sparse inverse solution using the Gamma-MAP empirical Bayesian method¶

See 1 for details.

# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#         Daniel Strohmeier <daniel.strohmeier@tu-ilmenau.de>
#
# License: BSD (3-clause)

import numpy as np

import mne
from mne.datasets import sample
from mne.inverse_sparse import gamma_map, make_stc_from_dipoles
from mne.viz import (plot_sparse_source_estimates,
                     plot_dipole_locations, plot_dipole_amplitudes)

print(__doc__)

data_path = sample.data_path()
subjects_dir = data_path + '/subjects'
fwd_fname = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
evoked_fname = data_path + '/MEG/sample/sample_audvis-ave.fif'
cov_fname = data_path + '/MEG/sample/sample_audvis-cov.fif'

# Read the evoked response and crop it
condition = 'Left visual'
evoked = mne.read_evokeds(evoked_fname, condition=condition,
                          baseline=(None, 0))
evoked.crop(tmin=-50e-3, tmax=300e-3)

# Read the forward solution
forward = mne.read_forward_solution(fwd_fname)

# Read noise noise covariance matrix and regularize it
cov = mne.read_cov(cov_fname)
cov = mne.cov.regularize(cov, evoked.info, rank=None)

# Run the Gamma-MAP method with dipole output
alpha = 0.5
dipoles, residual = gamma_map(
    evoked, forward, cov, alpha, xyz_same_gamma=True, return_residual=True,
    return_as_dipoles=True)

Out:

Reading /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-ave.fif ...
    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
    Found the data of interest:
        t =    -199.80 ...     499.49 ms (Left visual)
        0 CTF compensation matrices available
        nave = 67 - aspect type = 100
Projections have already been applied. Setting proj attribute to True.
Applying baseline correction (mode: mean)
Reading forward solution from /home/circleci/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
    366 x 366 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
Computing rank from covariance with rank=None
    Using tolerance 3.3e-13 (2.2e-16 eps * 305 dim * 4.8  max singular value)
    Estimated rank (mag + grad): 302
    MEG: rank 302 computed from 305 data channels with 3 projectors
    Using tolerance 4.7e-14 (2.2e-16 eps * 59 dim * 3.6  max singular value)
    Estimated rank (eeg): 58
    EEG: rank 58 computed from 59 data channels with 1 projector
8 projection items activated
    MAG regularization : 0.1
    Created an SSP operator (subspace dimension = 3)
Computing rank from covariance with rank={'meg': 302, 'eeg': 58}
    Using tolerance 2.5e-14 (2.2e-16 eps * 102 dim * 1.1  max singular value)
    Estimated rank (mag): 99
    MAG: rank 99 computed from 102 data channels with 3 projectors
    Setting small MAG eigenvalues to zero (without PCA)
    GRAD regularization : 0.1
Computing rank from covariance with rank={'meg': 302, 'eeg': 58, 'mag': 99}
    Using tolerance 1.8e-13 (2.2e-16 eps * 203 dim * 3.9  max singular value)
    Estimated rank (grad): 203
    GRAD: rank 203 computed from 203 data channels with 0 projectors
    Setting small GRAD eigenvalues to zero (without PCA)
    EEG regularization : 0.1
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank={'meg': 302, 'eeg': 58, 'mag': 99, 'grad': 203}
    Setting small EEG eigenvalues to zero (without PCA)
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 364 channels.
    364 out of 366 channels remain after picking
Selected 364 channels
Creating the depth weighting matrix...
Applying loose dipole orientations to surface source spaces: 0.2
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 4)
Computing rank from covariance with rank=None
    Using tolerance 3.3e-13 (2.2e-16 eps * 305 dim * 4.8  max singular value)
    Estimated rank (mag + grad): 302
    MEG: rank 302 computed from 305 data channels with 3 projectors
    Using tolerance 4.7e-14 (2.2e-16 eps * 59 dim * 3.6  max singular value)
    Estimated rank (eeg): 58
    EEG: rank 58 computed from 59 data channels with 1 projector
    Setting small MEG eigenvalues to zero (without PCA)
    Setting small EEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Whitening data matrix.
Iteration: 0     active set size: 22494  convergence: 9.054e-01
Iteration: 8     active set size: 22236  convergence: 1.864e-01
Iteration: 9     active set size: 20715  convergence: 1.377e-01
Iteration: 10    active set size: 17181  convergence: 1.063e-01
Iteration: 11    active set size: 13002  convergence: 8.293e-02
Iteration: 12    active set size: 10110  convergence: 6.591e-02
Iteration: 13    active set size: 8313   convergence: 5.335e-02
Iteration: 14    active set size: 7209   convergence: 4.383e-02
Iteration: 15    active set size: 6354   convergence: 3.644e-02
Iteration: 16    active set size: 5700   convergence: 3.060e-02
Iteration: 17    active set size: 5145   convergence: 2.591e-02
Iteration: 18    active set size: 4701   convergence: 2.211e-02
Iteration: 19    active set size: 4296   convergence: 1.900e-02
Iteration: 20    active set size: 4005   convergence: 1.644e-02
Iteration: 21    active set size: 3693   convergence: 1.432e-02
Iteration: 22    active set size: 3432   convergence: 1.255e-02
Iteration: 23    active set size: 3153   convergence: 1.107e-02
Iteration: 24    active set size: 2958   convergence: 9.830e-03
Iteration: 25    active set size: 2745   convergence: 8.781e-03
Iteration: 26    active set size: 2562   convergence: 7.891e-03
Iteration: 27    active set size: 2397   convergence: 7.133e-03
Iteration: 28    active set size: 2226   convergence: 6.484e-03
Iteration: 29    active set size: 2070   convergence: 6.061e-03
Iteration: 30    active set size: 1935   convergence: 5.884e-03
Iteration: 31    active set size: 1839   convergence: 5.719e-03
Iteration: 32    active set size: 1725   convergence: 5.563e-03
Iteration: 33    active set size: 1614   convergence: 5.431e-03
Iteration: 34    active set size: 1512   convergence: 5.326e-03
Iteration: 35    active set size: 1437   convergence: 5.220e-03
Iteration: 36    active set size: 1350   convergence: 5.113e-03
Iteration: 37    active set size: 1245   convergence: 5.004e-03
Iteration: 38    active set size: 1200   convergence: 4.891e-03
Iteration: 39    active set size: 1128   convergence: 4.775e-03
Iteration: 40    active set size: 1077   convergence: 4.654e-03
Iteration: 41    active set size: 1011   convergence: 4.530e-03
Iteration: 42    active set size: 975    convergence: 4.402e-03
Iteration: 43    active set size: 927    convergence: 4.271e-03
Iteration: 44    active set size: 876    convergence: 4.138e-03
Iteration: 45    active set size: 840    convergence: 4.003e-03
Iteration: 46    active set size: 798    convergence: 3.867e-03
Iteration: 47    active set size: 744    convergence: 3.731e-03
Iteration: 48    active set size: 705    convergence: 3.595e-03
Iteration: 49    active set size: 687    convergence: 3.459e-03
Iteration: 50    active set size: 657    convergence: 3.326e-03
Iteration: 51    active set size: 627    convergence: 3.194e-03
Iteration: 52    active set size: 600    convergence: 3.065e-03
Iteration: 53    active set size: 588    convergence: 2.940e-03
Iteration: 54    active set size: 570    convergence: 2.817e-03
Iteration: 55    active set size: 552    convergence: 2.698e-03
Iteration: 56    active set size: 528    convergence: 2.583e-03
Iteration: 57    active set size: 501    convergence: 2.472e-03
Iteration: 58    active set size: 486    convergence: 2.365e-03
Iteration: 59    active set size: 462    convergence: 2.263e-03
Iteration: 60    active set size: 444    convergence: 2.164e-03
Iteration: 61    active set size: 417    convergence: 2.070e-03
Iteration: 62    active set size: 408    convergence: 1.980e-03
Iteration: 63    active set size: 402    convergence: 1.894e-03
Iteration: 64    active set size: 387    convergence: 1.813e-03
Iteration: 65    active set size: 375    convergence: 1.735e-03
Iteration: 66    active set size: 345    convergence: 1.660e-03
Iteration: 67    active set size: 339    convergence: 1.590e-03
Iteration: 68    active set size: 333    convergence: 1.523e-03
Iteration: 69    active set size: 330    convergence: 1.459e-03
Iteration: 70    active set size: 327    convergence: 1.398e-03
Iteration: 71    active set size: 312    convergence: 1.341e-03
Iteration: 72    active set size: 306    convergence: 1.286e-03
Iteration: 73    active set size: 297    convergence: 1.235e-03
Iteration: 74    active set size: 288    convergence: 1.185e-03
Iteration: 75    active set size: 273    convergence: 1.139e-03
Iteration: 76    active set size: 261    convergence: 1.094e-03
Iteration: 77    active set size: 252    convergence: 1.052e-03
Iteration: 78    active set size: 240    convergence: 1.012e-03
Iteration: 79    active set size: 234    convergence: 9.744e-04
Iteration: 80    active set size: 222    convergence: 9.383e-04
Iteration: 81    active set size: 219    convergence: 9.040e-04
Iteration: 82    active set size: 213    convergence: 8.714e-04
Iteration: 83    active set size: 210    convergence: 8.404e-04
Iteration: 84    active set size: 207    convergence: 8.108e-04
Iteration: 85    active set size: 198    convergence: 7.827e-04
Iteration: 86    active set size: 192    convergence: 7.559e-04
Iteration: 87    active set size: 186    convergence: 7.303e-04
Iteration: 88    active set size: 180    convergence: 7.060e-04
Iteration: 89    active set size: 174    convergence: 6.827e-04
Iteration: 91    active set size: 168    convergence: 6.393e-04
Iteration: 92    active set size: 159    convergence: 6.191e-04
Iteration: 93    active set size: 153    convergence: 5.997e-04
Iteration: 94    active set size: 150    convergence: 5.812e-04
Iteration: 97    active set size: 147    convergence: 5.302e-04
Iteration: 98    active set size: 141    convergence: 5.146e-04
Iteration: 99    active set size: 132    convergence: 4.996e-04
Iteration: 100   active set size: 129    convergence: 4.853e-04
Iteration: 102   active set size: 126    convergence: 4.582e-04
Iteration: 103   active set size: 123    convergence: 4.455e-04
Iteration: 104   active set size: 120    convergence: 4.332e-04
Iteration: 105   active set size: 117    convergence: 4.215e-04
Iteration: 106   active set size: 111    convergence: 4.101e-04
Iteration: 107   active set size: 108    convergence: 3.992e-04
Iteration: 108   active set size: 105    convergence: 3.887e-04
Iteration: 110   active set size: 102    convergence: 3.688e-04
Iteration: 113   active set size: 99     convergence: 3.414e-04
Iteration: 114   active set size: 93     convergence: 3.329e-04
Iteration: 115   active set size: 90     convergence: 3.247e-04
Iteration: 119   active set size: 87     convergence: 2.945e-04
Iteration: 120   active set size: 84     convergence: 2.875e-04
Iteration: 123   active set size: 75     convergence: 2.679e-04
Iteration: 125   active set size: 69     convergence: 2.558e-04
Iteration: 127   active set size: 66     convergence: 2.444e-04
Iteration: 133   active set size: 63     convergence: 2.139e-04
Iteration: 136   active set size: 60     convergence: 2.005e-04
Iteration: 154   active set size: 57     convergence: 1.384e-04
Iteration: 155   active set size: 54     convergence: 1.357e-04
Iteration: 168   active set size: 51     convergence: 1.054e-04
Iteration: 177   active set size: 48     convergence: 8.890e-05
Iteration: 178   active set size: 45     convergence: 8.725e-05
Iteration: 181   active set size: 42     convergence: 8.248e-05
Iteration: 211   active set size: 39     convergence: 4.750e-05
Iteration: 225   active set size: 36     convergence: 3.685e-05
Iteration: 252   active set size: 33     convergence: 2.265e-05
Iteration: 292   active set size: 30     convergence: 1.105e-05
Iteration: 303   active set size: 27     convergence: 9.077e-06
Iteration: 310   active set size: 24     convergence: 8.008e-06
Iteration: 313   active set size: 21     convergence: 7.589e-06
Iteration: 340   active set size: 18     convergence: 4.680e-06
Iteration: 342   active set size: 15     convergence: 4.515e-06
Iteration: 381   active set size: 12     convergence: 2.247e-06
Iteration: 427   active set size: 12     convergence: 9.865e-07

Convergence reached !

Explained  22.3% variance
[done]

Plot dipole activations

plot_dipole_amplitudes(dipoles)

# Plot dipole location of the strongest dipole with MRI slices
idx = np.argmax([np.max(np.abs(dip.amplitude)) for dip in dipoles])
plot_dipole_locations(dipoles[idx], forward['mri_head_t'], 'sample',
                      subjects_dir=subjects_dir, mode='orthoview',
                      idx='amplitude')

# # Plot dipole locations of all dipoles with MRI slices
# for dip in dipoles:
#     plot_dipole_locations(dip, forward['mri_head_t'], 'sample',
#                           subjects_dir=subjects_dir, mode='orthoview',
#                           idx='amplitude')

Show the evoked response and the residual for gradiometers

ylim = dict(grad=[-120, 120])
evoked.pick_types(meg='grad', exclude='bads')
evoked.plot(titles=dict(grad='Evoked Response Gradiometers'), ylim=ylim,
            proj=True, time_unit='s')

residual.pick_types(meg='grad', exclude='bads')
residual.plot(titles=dict(grad='Residuals Gradiometers'), ylim=ylim,
              proj=True, time_unit='s')

Out:

Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Removing projector <Projection | Average EEG reference, active : True, n_channels : 60>
Removing projector <Projection | PCA-v1, active : True, n_channels : 102>
Removing projector <Projection | PCA-v2, active : True, n_channels : 102>
Removing projector <Projection | PCA-v3, active : True, n_channels : 102>
Removing projector <Projection | Average EEG reference, active : True, n_channels : 60>

Generate stc from dipoles

stc = make_stc_from_dipoles(dipoles, forward['src'])

Out:

Converting dipoles into a SourceEstimate.
[done]

View in 2D and 3D (“glass” brain like 3D plot) Show the sources as spheres scaled by their strength

scale_factors = np.max(np.abs(stc.data), axis=1)
scale_factors = 0.5 * (1 + scale_factors / np.max(scale_factors))

plot_sparse_source_estimates(
    forward['src'], stc, bgcolor=(1, 1, 1),
    modes=['sphere'], opacity=0.1, scale_factors=(scale_factors, None),
    fig_name="Gamma-MAP")

Out:

Total number of active sources: 4

References¶

1: David Wipf and Srikantan Nagarajan. A unified Bayesian framework for MEG/EEG source imaging. NeuroImage, 44(3):947–966, 2009. doi:10.1016/j.neuroimage.2008.02.059.

Total running time of the script: ( 1 minutes 2.475 seconds)

Estimated memory usage: 357 MB