Compare simulated and estimated source activity#

This example illustrates how to compare the simulated and estimated source time courses (STC) by computing different metrics. Simulated source is a cortical region or dipole. It is meant to be a brief introduction and only highlights the simplest use case.

# Author: Kostiantyn Maksymenko <kostiantyn.maksymenko@gmail.com>
#
# License: BSD (3-clause)

import numpy as np
import matplotlib.pyplot as plt
from functools import partial

import mne
from mne.datasets import sample
from mne.minimum_norm import make_inverse_operator, apply_inverse
from mne.simulation.metrics import (region_localization_error,
                                    f1_score, precision_score,
                                    recall_score, cosine_score,
                                    peak_position_error,
                                    spatial_deviation_error)

random_state = 42  # set random state to make this example deterministic

# Import sample data
data_path = sample.data_path()
subjects_dir = data_path / 'subjects'
subject = 'sample'
evoked_fname = data_path / 'MEG' / subject / 'sample_audvis-ave.fif'
info = mne.io.read_info(evoked_fname)
tstep = 1. / info['sfreq']

# Import forward operator and source space
fwd_fname = data_path / 'MEG' / subject / 'sample_audvis-meg-eeg-oct-6-fwd.fif'
fwd = mne.read_forward_solution(fwd_fname)
src = fwd['src']

# To select source, we use the caudal middle frontal to grow
# a region of interest.
selected_label = mne.read_labels_from_annot(
    subject, regexp='caudalmiddlefrontal-lh', subjects_dir=subjects_dir)[0]
    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
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)
    Forward solutions combined: MEG, EEG
    Source spaces transformed to the forward solution coordinate frame
Reading labels from parcellation...
   read 1 labels from /home/circleci/mne_data/MNE-sample-data/subjects/sample/label/lh.aparc.annot
   read 0 labels from /home/circleci/mne_data/MNE-sample-data/subjects/sample/label/rh.aparc.annot

In this example we simulate two types of cortical sources: a region and a dipole sources. We will test corresponding performance metrics.

Define main parameters of sources#

First we define both region and dipole sources in terms of Where?, What? and When?.

# WHERE?

# Region
location = 'center'  # Use the center of the label as a seed.
extent = 20.  # Extent in mm of the region.
label_region = mne.label.select_sources(
    subject, selected_label, location=location, extent=extent,
    subjects_dir=subjects_dir, random_state=random_state)

# Dipole
location = 1915  # Use the index of the vertex as a seed
extent = 0.  # One dipole source
label_dipole = mne.label.select_sources(
    subject, selected_label, location=location, extent=extent,
    subjects_dir=subjects_dir, random_state=random_state)

# WHAT?
# Define the time course of the activity
source_time_series = np.sin(2. * np.pi * 18. * np.arange(100) * tstep) * 10e-9

# WHEN?
# Define when the activity occurs using events.
n_events = 50
events = np.zeros((n_events, 3), int)
events[:, 0] = 200 * np.arange(n_events)  # Events sample.
events[:, 2] = 1  # All events have the sample id.
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.1s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.1s finished
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s finished

Create simulated source activity#

Here, SourceSimulator is used.

Simulate raw data#

Project the source time series to sensor space with multivariate Gaussian noise obtained from the noise covariance from the sample data.

Setting up raw simulation: 1 position, "cos2" interpolation
Event information stored on channel:              STI 014
    Interval 0.000-1.665 sec
Setting up forward solutions
Computing gain matrix for transform #1/1
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.498 sec
    10 STC iterations provided
[done]
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>
Adding noise to 59/68 channels (59 channels in cov)
Setting up raw simulation: 1 position, "cos2" interpolation
Event information stored on channel:              STI 014
    Interval 0.000-1.665 sec
Setting up forward solutions
Computing gain matrix for transform #1/1
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.665 sec
    Interval 0.000-1.498 sec
    10 STC iterations provided
[done]
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>
Adding noise to 59/68 channels (59 channels in cov)
Measurement date December 03, 2002 19:01:10 GMT
Experimenter Unknown
Participant Unknown
Digitized points 146 points
Good channels 9 Stimulus, 59 EEG
Bad channels None
EOG channels Not available
ECG channels Not available
Sampling frequency 600.61 Hz
Highpass 0.10 Hz
Lowpass 40.00 Hz
Projections Average EEG reference : on
Duration 00:00:17 (HH:MM:SS)


Compute evoked from raw data#

Averaging epochs corresponding to events.

# Region
events = mne.find_events(raw_region, initial_event=True)
tmax = (len(source_time_series) - 1) * tstep
epochs = mne.Epochs(raw_region, events, 1, tmin=0, tmax=tmax, baseline=None)
evoked_region = epochs.average()

# Dipole
events = mne.find_events(raw_dipole, initial_event=True)
tmax = (len(source_time_series) - 1) * tstep
epochs = mne.Epochs(raw_dipole, events, 1, tmin=0, tmax=tmax, baseline=None)
evoked_dipole = epochs.average()
50 events found
Event IDs: [1]
Not setting metadata
50 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 1)
1 projection items activated
50 events found
Event IDs: [1]
Not setting metadata
50 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 1)
1 projection items activated

Create true stcs corresponding to evoked#

Before we computed stcs corresponding to raw data. To be able to compare it with the reconstruction, based on the evoked, true stc should have the same number of time samples.

# Region
stc_true_region = \
    source_simulator_region.get_stc(start_sample=0,
                                    stop_sample=len(source_time_series))

# Dipole
stc_true_dipole = \
    source_simulator_dipole.get_stc(start_sample=0,
                                    stop_sample=len(source_time_series))

Reconstruct simulated sources#

Compute inverse solution using sLORETA.

# Region
snr = 30.0
inv_method = 'sLORETA'
lambda2 = 1.0 / snr ** 2

inverse_operator = make_inverse_operator(evoked_region.info, fwd, cov,
                                         loose='auto', depth=0.8,
                                         fixed=True)

stc_est_region = apply_inverse(evoked_region, inverse_operator, lambda2,
                               inv_method, pick_ori=None)

# Dipole
snr = 3.0
inv_method = 'sLORETA'
lambda2 = 1.0 / snr ** 2

inverse_operator = make_inverse_operator(evoked_dipole.info, fwd, cov,
                                         loose='auto', depth=0.8,
                                         fixed=True)

stc_est_dipole = apply_inverse(evoked_dipole, inverse_operator, lambda2,
                               inv_method, pick_ori=None)
Computing inverse operator with 59 channels.
    59 out of 366 channels remain after picking
Selected 59 channels
Creating the depth weighting matrix...
    59 EEG channels
    limit = 7499/7498 = 2.118742
    scale = 155292 exp = 0.8
    Picked elements from a free-orientation depth-weighting prior into the fixed-orientation one
    Average patch normals will be employed in the rotation to the local surface coordinates....
    Converting to surface-based source orientations...
    [done]
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank=None
    Using tolerance 5.2e-18 (2.2e-16 eps * 59 dim * 0.0004  max singular value)
    Estimated rank (eeg): 58
    EEG: rank 58 computed from 59 data channels with 1 projector
    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 = 4.2931
    scaling factor to adjust the trace = 1.34986e+22 (nchan = 59 nzero = 1)
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 50
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 1)
    Created the whitener using a noise covariance matrix with rank 58 (1 small eigenvalues omitted)
    Computing noise-normalization factors (sLORETA)...
[done]
Applying inverse operator to "1"...
    Picked 59 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained 100.0% variance
    sLORETA...
[done]
Computing inverse operator with 59 channels.
    59 out of 366 channels remain after picking
Selected 59 channels
Creating the depth weighting matrix...
    59 EEG channels
    limit = 7499/7498 = 2.118742
    scale = 155292 exp = 0.8
    Picked elements from a free-orientation depth-weighting prior into the fixed-orientation one
    Average patch normals will be employed in the rotation to the local surface coordinates....
    Converting to surface-based source orientations...
    [done]
Whitening the forward solution.
    Created an SSP operator (subspace dimension = 1)
Computing rank from covariance with rank=None
    Using tolerance 5.2e-18 (2.2e-16 eps * 59 dim * 0.0004  max singular value)
    Estimated rank (eeg): 58
    EEG: rank 58 computed from 59 data channels with 1 projector
    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 = 4.2931
    scaling factor to adjust the trace = 1.34986e+22 (nchan = 59 nzero = 1)
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 50
    Created the regularized inverter
    Created an SSP operator (subspace dimension = 1)
    Created the whitener using a noise covariance matrix with rank 58 (1 small eigenvalues omitted)
    Computing noise-normalization factors (sLORETA)...
[done]
Applying inverse operator to "1"...
    Picked 59 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained  85.9% variance
    sLORETA...
[done]

Compute performance scores for different source amplitude thresholds#

thresholds = [10, 30, 50, 70, 80, 90, 95, 99]

For region#

# create a set of scorers
scorers = {'RLE': partial(region_localization_error, src=src),
           'Precision': precision_score, 'Recall': recall_score,
           'F1 score': f1_score}

# compute results
region_results = {}
for name, scorer in scorers.items():
    region_results[name] = [scorer(stc_true_region, stc_est_region,
                                   threshold=f'{thx}%', per_sample=False)
                            for thx in thresholds]

# Plot the results
f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    2, 2, sharex='col', constrained_layout=True)
for ax, (title, results) in zip([ax1, ax2, ax3, ax4], region_results.items()):
    ax.plot(thresholds, results, '.-')
    ax.set(title=title, ylabel='score', xlabel='Threshold',
           xticks=thresholds)

f.suptitle('Performance scores per threshold')  # Add Super title
ax1.ticklabel_format(axis='y', style='sci', scilimits=(0, 1))  # tweak RLE

# Cosine score with respect to time
f, ax1 = plt.subplots(constrained_layout=True)
ax1.plot(stc_true_region.times, cosine_score(stc_true_region, stc_est_region))
ax1.set(title='Cosine score', xlabel='Time', ylabel='Score')
  • Performance scores per threshold, RLE, Precision, Recall, F1 score
  • Cosine score

For Dipoles#

# create a set of scorers
scorers = {
    'Peak Position Error': peak_position_error,
    'Spatial Deviation Error': spatial_deviation_error,
}


# compute results
dipole_results = {}
for name, scorer in scorers.items():
    dipole_results[name] = [scorer(stc_true_dipole, stc_est_dipole, src=src,
                                   threshold=f'{thx}%', per_sample=False)
                            for thx in thresholds]


# Plot the results
for name, results in dipole_results.items():
    f, ax1 = plt.subplots(constrained_layout=True)
    ax1.plot(thresholds, 100 * np.array(results), '.-')
    ax1.set(title=name, ylabel='Error (cm)', xlabel='Threshold',
            xticks=thresholds)
  • Peak Position Error
  • Spatial Deviation Error

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

Estimated memory usage: 196 MB

Gallery generated by Sphinx-Gallery