Compute MxNE with time-frequency sparse prior¶

The TF-MxNE solver is a distributed inverse method (like dSPM or sLORETA) that promotes focal (sparse) sources (such as dipole fitting techniques) 1 2. The benefit of this approach is that:

• it is spatio-temporal without assuming stationarity (sources properties can vary over time)

• activations are localized in space, time and frequency in one step.

• with a built-in filtering process based on a short time Fourier transform (STFT), data does not need to be low passed (just high pass to make the signals zero mean).

• the solver solves a convex optimization problem, hence cannot be trapped in local minima.

References¶

1

A. Gramfort, D. Strohmeier, J. Haueisen, M. Hamalainen, M. Kowalski “Time-Frequency Mixed-Norm Estimates: Sparse M/EEG imaging with non-stationary source activations”, Neuroimage, Volume 70, pp. 410-422, 15 April 2013. DOI: 10.1016/j.neuroimage.2012.12.051

2

A. Gramfort, D. Strohmeier, J. Haueisen, M. Hamalainen, M. Kowalski “Functional Brain Imaging with M/EEG Using Structured Sparsity in Time-Frequency Dictionaries”, Proceedings Information Processing in Medical Imaging Lecture Notes in Computer Science, Volume 6801/2011, pp. 600-611, 2011. DOI: 10.1007/978-3-642-22092-0_49

# Author: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#         Daniel Strohmeier <daniel.strohmeier@tu-ilmenau.de>
#

import numpy as np

import mne
from mne.datasets import sample
from mne.minimum_norm import make_inverse_operator, apply_inverse
from mne.inverse_sparse import tf_mixed_norm, 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'
ave_fname = data_path + '/MEG/sample/sample_audvis-no-filter-ave.fif'
cov_fname = data_path + '/MEG/sample/sample_audvis-shrunk-cov.fif'

# Handling average file
condition = 'Left visual'
evoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))
evoked = mne.pick_channels_evoked(evoked)
# We make the window slightly larger than what you'll eventually be interested
# in ([-0.05, 0.3]) to avoid edge effects.
evoked.crop(tmin=-0.1, tmax=0.4)

# Handling forward solution


Out:

    365 x 365 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 59) active
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 = 64 - aspect type = 100
Projections have already been applied. Setting proj attribute to True.
Applying baseline correction (mode: mean)
Computing patch statistics...
[done]
Computing patch statistics...
[done]
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


Run solver

# alpha parameter is between 0 and 100 (100 gives 0 active source)
alpha = 40.  # general regularization parameter
# l1_ratio parameter between 0 and 1 promotes temporal smoothness
# (0 means no temporal regularization)
l1_ratio = 0.03  # temporal regularization parameter

loose, depth = 0.2, 0.9  # loose orientation & depth weighting

# Compute dSPM solution to be used as weights in MxNE
inverse_operator = make_inverse_operator(evoked.info, forward, cov,
loose=loose, depth=depth)
stc_dspm = apply_inverse(evoked, inverse_operator, lambda2=1. / 9.,
method='dSPM')

# Compute TF-MxNE inverse solution with dipole output
dipoles, residual = tf_mixed_norm(
evoked, forward, cov, alpha=alpha, l1_ratio=l1_ratio, loose=loose,
depth=depth, maxit=200, tol=1e-6, weights=stc_dspm, weights_min=8.,
debias=True, wsize=16, tstep=4, window=0.05, return_as_dipoles=True,
return_residual=True)

# Crop to remove edges
for dip in dipoles:
dip.crop(tmin=-0.05, tmax=0.3)
evoked.crop(tmin=-0.05, tmax=0.3)
residual.crop(tmin=-0.05, tmax=0.3)


Out:

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...
203 planar channels
limit = 7262/7498 = 10.020865
scale = 2.58122e-08 exp = 0.9
Applying loose dipole orientations. Loose value of 0.2.
Whitening the forward solution.
Created an SSP operator (subspace dimension = 4)
Computing data rank from covariance with rank=None
Using tolerance 3.5e-13 (2.2e-16 eps * 305 dim * 5.2  max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Using tolerance 1.1e-13 (2.2e-16 eps * 59 dim * 8.7  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
Computing SVD of whitened and weighted lead field matrix.
largest singular value = 5.96729
scaling factor to adjust the trace = 9.38524e+18
Preparing the inverse operator for use...
Scaled noise and source covariance from nave = 1 to nave = 64
Created the regularized inverter
Created an SSP operator (subspace dimension = 4)
Created the whitener using a noise covariance matrix with rank 360 (4 small eigenvalues omitted)
Computing noise-normalization factors (dSPM)...
[done]
Applying inverse operator to "Left visual"...
Picked 364 channels from the data
Computing inverse...
Eigenleads need to be weighted ...
Computing residual...
Explained  60.0% variance
Combining the current components...
dSPM...
[done]
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. Loose value of 0.2.
Whitening the forward solution.
Created an SSP operator (subspace dimension = 4)
Computing data rank from covariance with rank=None
Using tolerance 3.5e-13 (2.2e-16 eps * 305 dim * 5.2  max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Using tolerance 1.1e-13 (2.2e-16 eps * 59 dim * 8.7  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
Reducing source space to 985 sources
Whitening data matrix.
Using block coordinate descent with active set approach

Iteration 10 :: n_active 3
dgap 3.49e+00 :: p_obj 4411.845726 :: d_obj 4408.353441

Iteration 20 :: n_active 3
dgap 5.67e-01 :: p_obj 4410.859492 :: d_obj 4410.292946

dgap 1.51e-01 :: p_obj 4410.670058 :: d_obj 4410.519426 :: n_active 2

Iteration 10 :: n_active 2
dgap 1.61e-03 :: p_obj 4410.669663 :: d_obj 4410.668049

Iteration 20 :: n_active 2
dgap 3.02e-05 :: p_obj 4410.669663 :: d_obj 4410.669633

Iteration 30 :: n_active 2
dgap 5.66e-07 :: p_obj 4410.669663 :: d_obj 4410.669663

dgap 5.66e-07 :: p_obj 4410.669663 :: d_obj 4410.669663 :: n_active 2

Convergence reached!

Debiasing converged after 190 iterations max(|D - D0| = 4.753993e-07 < 1.000000e-06)
4 projection items deactivated
Created an SSP operator (subspace dimension = 4)
4 projection items activated
SSP projectors applied...
0 projection items deactivated
[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])
subjects_dir=subjects_dir, mode='orthoview',
idx='amplitude')

# # Plot dipole locations of all dipoles with MRI slices
# for dip in dipoles:
#                           subjects_dir=subjects_dir, mode='orthoview',
#                           idx='amplitude')


Show the evoked response and the residual for gradiometers

ylim = dict(grad=[-120, 120])
proj=True, time_unit='s')

proj=True, time_unit='s')


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)

plot_sparse_source_estimates(forward['src'], stc, bgcolor=(1, 1, 1),
opacity=0.1, fig_name="TF-MxNE (cond %s)"
% condition, modes=['sphere'], scale_factors=[1.])

time_label = 'TF-MxNE time=%0.2f ms'
clim = dict(kind='value', lims=[10e-9, 15e-9, 20e-9])
brain = stc.plot('sample', 'inflated', 'rh', views='medial',
clim=clim, time_label=time_label, smoothing_steps=5,
subjects_dir=subjects_dir, initial_time=150, time_unit='ms')

Total number of active sources: 2