非官方中文翻译 本页译自 MNE-Python 官方文档： 查看英文原文 ↗。 文档版权归 MNE Developers 所有（BSD-3-Clause）。

使用时频稀疏先验计算 MxNE

TF-MxNE 求解器是一种分布式逆方法（类似 dSPM 或 sLORETA），可促进聚焦（稀疏）源（类似偶极子拟合技术）[1] [2]。该方法的优点包括：

• 它是时空方法，且不假设平稳性（源属性可随时间变化）

• 一步即可在空间、时间与频率上定位激活。

• 借助基于短时傅里叶变换（STFT）的内置滤波过程，数据无需低通（只需高通使信号均值为零）。

• 求解器求解凸优化问题，因此不会陷入局部极小。

# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
#         Daniel Strohmeier <daniel.strohmeier@tu-ilmenau.de>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import numpy as np

import mne
from mne.datasets import sample
from mne.inverse_sparse import make_stc_from_dipoles, tf_mixed_norm
from mne.minimum_norm import apply_inverse, make_inverse_operator
from mne.viz import (
    plot_dipole_amplitudes,
    plot_dipole_locations,
    plot_sparse_source_estimates,
)

print(__doc__)

data_path = sample.data_path()
subjects_dir = data_path / "subjects"
meg_path = data_path / "MEG" / "sample"
fwd_fname = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
ave_fname = meg_path / "sample_audvis-no-filter-ave.fif"
cov_fname = meg_path / "sample_audvis-shrunk-cov.fif"

# Read noise covariance matrix
cov = mne.read_cov(cov_fname)

# Handling average file
condition = "Left visual"
evoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))
# 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
forward = mne.read_forward_solution(fwd_fname)

    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
Reading /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-no-filter-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 = 64 - 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 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
    Read MEG forward solution (7498 sources, 306 channels, free orientations)
    Desired named matrix (kind = 3523 (FIFF_MNE_FORWARD_SOLUTION_GRAD)) not available
    Read EEG forward solution (7498 sources, 60 channels, free orientations)
    Forward solutions combined: MEG, EEG
    Source spaces transformed to the forward solution coordinate frame

运行求解器

# alpha parameter is between 0 and 100 (100 gives 0 active source)
alpha = 40.0  # 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.0 / 9.0, 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.0,
    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)

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]
info["bads"] and noise_cov["bads"] do not match, excluding bad channels from both
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 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.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
Adjusting 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 (nchan = 364 nzero = 4)
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]
info["bads"] and noise_cov["bads"] do not match, excluding bad channels from both
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.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
Adjusting 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.845725527788 :: d_obj 4408.353440672163

    Iteration 20 :: n_active 3
    dgap 5.67e-01 :: p_obj 4410.85949217969 :: d_obj 4410.292945587396

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.669663490561 :: d_obj 4410.668049322309

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

Convergence reached!

Debiasing converged after 190 iterations max(|D - D0| = 5.549147e-07 < 1.000000e-06)
[done]

	General
	Filename(s)	sample_audvis-no-filter-ave.fif
	MNE object type	Evoked
	Measurement date	2002-12-03 at 19:01:10 UTC
	Participant	Unknown
	Experimenter	Unknown
	Acquisition
	Aggregation	average of 64 epochs
	条件	Left visual
	Time range	-0.050 – 0.300 s
	Baseline	-0.200 – 0.000 s
	Sampling frequency	600.61 Hz
	Time points	211
	Channels
	Magnetometers	102
	Gradiometers	203
	EEG	59
	Head & sensor digitization	146 points
	Filters
	Highpass	0.10 Hz
	Lowpass	172.18 Hz
	Projections	PCA-v1 (on) PCA-v2 (on) PCA-v3 (on) Average EEG reference (on)

绘制偶极子激活

plot_dipole_amplitudes(dipoles)

用 MRI 切片绘制最强偶极子的位置

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')

显示梯度计的诱发响应与残差

ylim = dict(grad=[-120, 120])
evoked.pick(picks="grad", exclude="bads")
evoked.plot(
    titles=dict(grad="Evoked Response: Gradiometers"),
    ylim=ylim,
    proj=True,
    time_unit="s",
)

residual.pick(picks="grad", exclude="bads")
residual.plot(
    titles=dict(grad="Residuals: Gradiometers"), ylim=ylim, proj=True, time_unit="s"
)

• 
• 

由偶极子生成 STC

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

Converting dipoles into a SourceEstimate.
[done]

以 2D 与 3D 查看（类似 3D 的“玻璃脑”图）

plot_sparse_source_estimates(
    forward["src"],
    stc,
    bgcolor=(1, 1, 1),
    opacity=0.1,
    fig_name=f"TF-MxNE (cond {condition})",
    modes=["sphere"],
    scale_factors=[1.0],
)

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",
)
brain.add_label("V1", color="yellow", scalar_thresh=0.5, borders=True)
brain.add_label("V2", color="red", scalar_thresh=0.5, borders=True)

Total number of active sources: [160797 167557]

参考文献

[1]

Alexandre Gramfort, Daniel T. Strohmeier, Jens Haueisen, Matti S. Hämäläinen, and Matthieu Kowalski. Time-frequency mixed-norm estimates: sparse M/EEG imaging with non-stationary source activations. NeuroImage, 70:410–422, 2013. doi:10.1016/j.neuroimage.2012.12.051.

[2]

Alexandre Gramfort, Daniel Strohmeier, Jens Haueisen, Matti S. Hämäläinen, and Matthieu Kowalski. Functional brain imaging with M/EEG using structured sparsity in time-frequency dictionaries. In Gábor Székely and Horst K. Hahn, editors, Information Processing in Medical Imaging, volume 6801, pages 600–611. Springer, Berlin; Heidelberg, 2011. doi:10.1007/978-3-642-22092-0_49.