Simulate raw data using subject anatomy

This example illustrates how to generate source estimates and simulate raw data using subject anatomy with the mne.simulation.SourceSimulator class. Once the raw data is simulated, generated source estimates are reconstructed using dynamic statistical parametric mapping (dSPM) inverse operator.

# Author: Ivana Kojcic <ivana.kojcic@gmail.com>
#         Eric Larson <larson.eric.d@gmail.com>
#         Kostiantyn Maksymenko <kostiantyn.maksymenko@gmail.com>
#         Samuel Deslauriers-Gauthier <sam.deslauriers@gmail.com>

# License: BSD (3-clause)

import os.path as op

import numpy as np

import mne
from mne.datasets import sample

print(__doc__)

# To simulate the sample dataset, information of the sample subject needs to be
# loaded. This step will download the data if it not already on your machine.
# Subjects directory is also set so it doesn't need to be given to functions.
data_path = sample.data_path()
subjects_dir = op.join(data_path, 'subjects')
subject = 'sample'
meg_path = op.join(data_path, 'MEG', subject)

# First, we get an info structure from the sample subject.
fname_info = op.join(meg_path, 'sample_audvis_raw.fif')
info = mne.io.read_info(fname_info)
tstep = 1 / info['sfreq']

# To simulate sources, we also need a source space. It can be obtained from the
# forward solution of the sample subject.
fwd_fname = op.join(meg_path, 'sample_audvis-meg-eeg-oct-6-fwd.fif')
fwd = mne.read_forward_solution(fwd_fname)
src = fwd['src']

# To simulate raw data, we need to define when the activity occurs using events
# matrix and specify the IDs of each event.
# Noise covariance matrix also needs to be defined.
# Here, both are loaded from the sample dataset, but they can also be specified
# by the user.

fname_event = op.join(meg_path, 'sample_audvis_raw-eve.fif')
fname_cov = op.join(meg_path, 'sample_audvis-cov.fif')

events = mne.read_events(fname_event)
noise_cov = mne.read_cov(fname_cov)

# Standard sample event IDs. These values will correspond to the third column
# in the events matrix.
event_id = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
            'visual/right': 4, 'smiley': 5, 'button': 32}

Out:

    Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
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

In order to simulate source time courses, labels of desired active regions need to be specified for each of the 4 simulation conditions. Make a dictionary that maps conditions to activation strengths within aparc.a2009s 1 labels. In the aparc.a2009s parcellation:

  • ‘G_temp_sup-G_T_transv’ is the label for primary auditory area

  • ‘S_calcarine’ is the label for primary visual area

In each of the 4 conditions, only the primary area is activated. This means that during the activations of auditory areas, there are no activations in visual areas and vice versa. Moreover, for each condition, contralateral region is more active (here, 2 times more) than the ipsilateral.

activations = {
    'auditory/left':
        [('G_temp_sup-G_T_transv-lh', 100),          # label, activation (nAm)
         ('G_temp_sup-G_T_transv-rh', 200)],
    'auditory/right':
        [('G_temp_sup-G_T_transv-lh', 200),
         ('G_temp_sup-G_T_transv-rh', 100)],
    'visual/left':
        [('S_calcarine-lh', 100),
         ('S_calcarine-rh', 200)],
    'visual/right':
        [('S_calcarine-lh', 200),
         ('S_calcarine-rh', 100)],
}

annot = 'aparc.a2009s'

# Load the 4 necessary label names.
label_names = sorted(set(activation[0]
                         for activation_list in activations.values()
                         for activation in activation_list))
region_names = list(activations.keys())

#  Define the time course of the activity for each region to activate. We use a
#  sine wave and it will be the same for all 4 regions.
source_time_series = np.sin(np.linspace(0, 4 * np.pi, 100)) * 10e-9

Create simulated source activity

Here, SourceSimulator is used, which allows to specify where (label), what (source_time_series), and when (events) event type will occur.

We will add data for 4 areas, each of which contains 2 labels. Since add_data method accepts 1 label per call, it will be called 2 times per area. All activations will contain the same waveform, but the amplitude will be 2 times higher in the contralateral label, as explained before.

When the activity occurs is defined using events. In this case, they are taken from the original raw data. The first column is the sample of the event, the second is not used. The third one is the event id, which is different for each of the 4 areas.

source_simulator = mne.simulation.SourceSimulator(src, tstep=tstep)
for region_id, region_name in enumerate(region_names, 1):
    events_tmp = events[np.where(events[:, 2] == region_id)[0], :]
    for i in range(2):
        label_name = activations[region_name][i][0]
        label_tmp = mne.read_labels_from_annot(subject, annot,
                                               subjects_dir=subjects_dir,
                                               regexp=label_name,
                                               verbose=False)
        label_tmp = label_tmp[0]
        amplitude_tmp = activations[region_name][i][1]
        source_simulator.add_data(label_tmp,
                                  amplitude_tmp * source_time_series,
                                  events_tmp)

# To obtain a SourceEstimate object, we need to use `get_stc()` method of
# SourceSimulator class.
stc_data = source_simulator.get_stc()

Simulate raw data

Project the source time series to sensor space. Three types of noise will be added to the simulated raw data:

  • multivariate Gaussian noise obtained from the noise covariance from the sample data

  • blink (EOG) noise

  • ECG noise

The SourceSimulator can be given directly to the simulate_raw() function.

raw_sim = mne.simulation.simulate_raw(info, source_simulator, forward=fwd,
                                      cov=None)
raw_sim.set_eeg_reference(projection=True).crop(0, 60)  # for speed

mne.simulation.add_noise(raw_sim, cov=noise_cov, random_state=0)
mne.simulation.add_eog(raw_sim, random_state=0)
mne.simulation.add_ecg(raw_sim, random_state=0)

# Plot original and simulated raw data.
raw_sim.plot(title='Simulated raw data')
../../_images/sphx_glr_plot_simulated_raw_data_using_subject_anatomy_001.png

Out:

Setting up raw simulation: 1 position, "cos2" interpolation
Event information stored on channel:              STI 014
Setting up forward solutions
Computing gain matrix for transform #1/1
   Simulating data for forward operator 1/0
   Interval 0.000-1.665 sec
   Interval 1.665-3.330 sec
   Interval 3.330-4.995 sec
   Interval 4.995-6.660 sec
   Interval 6.660-8.325 sec
   Interval 8.325-9.990 sec
   Interval 9.990-11.655 sec
   Interval 11.655-13.320 sec
   Interval 13.320-14.985 sec
   Interval 14.985-16.650 sec
   Interval 16.650-18.315 sec
   Interval 18.315-19.980 sec
   Interval 19.980-21.644 sec
   Interval 21.644-23.309 sec
   Interval 23.309-24.974 sec
   Interval 24.974-26.639 sec
   Interval 26.639-28.304 sec
   Interval 28.304-29.969 sec
   Interval 29.969-31.634 sec
   Interval 31.634-33.299 sec
   Interval 33.299-34.964 sec
   Interval 34.964-36.629 sec
   Interval 36.629-38.294 sec
   Interval 38.294-39.959 sec
   Interval 39.959-41.624 sec
   Interval 41.624-43.289 sec
   Interval 43.289-44.954 sec
   Interval 44.954-46.619 sec
   Interval 46.619-48.284 sec
   Interval 48.284-49.949 sec
   Interval 49.949-51.614 sec
   Interval 51.614-53.279 sec
   Interval 53.279-54.944 sec
   Interval 54.944-56.609 sec
   Interval 56.609-58.274 sec
   Interval 58.274-59.939 sec
   Interval 59.939-61.604 sec
   Interval 61.604-63.268 sec
   Interval 63.268-64.933 sec
   Interval 64.933-66.598 sec
   Interval 66.598-68.263 sec
   Interval 68.263-69.928 sec
   Interval 69.928-71.593 sec
   Interval 71.593-73.258 sec
   Interval 73.258-74.923 sec
   Interval 74.923-76.588 sec
   Interval 76.588-78.253 sec
   Interval 78.253-79.918 sec
   Interval 79.918-81.583 sec
   Interval 81.583-83.248 sec
   Interval 83.248-84.913 sec
   Interval 84.913-86.578 sec
   Interval 86.578-88.243 sec
   Interval 88.243-89.908 sec
   Interval 89.908-91.573 sec
   Interval 91.573-93.238 sec
   Interval 93.238-94.903 sec
   Interval 94.903-96.568 sec
   Interval 96.568-98.233 sec
   Interval 98.233-99.898 sec
   Interval 99.898-101.563 sec
   Interval 101.563-103.228 sec
   Interval 103.228-104.892 sec
   Interval 104.892-106.557 sec
   Interval 106.557-108.222 sec
   Interval 108.222-109.887 sec
   Interval 109.887-111.552 sec
   Interval 111.552-113.217 sec
   Interval 113.217-114.882 sec
   Interval 114.882-116.547 sec
   Interval 116.547-118.212 sec
   Interval 118.212-119.877 sec
   Interval 119.877-121.542 sec
   Interval 121.542-123.207 sec
   Interval 123.207-124.872 sec
   Interval 124.872-126.537 sec
   Interval 126.537-128.202 sec
   Interval 128.202-129.867 sec
   Interval 129.867-131.532 sec
   Interval 131.532-133.197 sec
   Interval 133.197-134.862 sec
   Interval 134.862-136.527 sec
   Interval 136.527-138.192 sec
   Interval 138.192-139.857 sec
   Interval 139.857-141.522 sec
   Interval 141.522-143.187 sec
   Interval 143.187-144.852 sec
   Interval 144.852-146.516 sec
   Interval 146.516-148.181 sec
   Interval 148.181-149.846 sec
   Interval 149.846-151.511 sec
   Interval 151.511-153.176 sec
   Interval 153.176-154.841 sec
   Interval 154.841-156.506 sec
   Interval 156.506-158.171 sec
   Interval 158.171-159.836 sec
   Interval 159.836-161.501 sec
   Interval 161.501-163.166 sec
   Interval 163.166-164.831 sec
   Interval 164.831-166.496 sec
   Interval 166.496-168.161 sec
   Interval 168.161-169.826 sec
   Interval 169.826-171.491 sec
   Interval 171.491-173.156 sec
   Interval 173.156-174.821 sec
   Interval 174.821-176.486 sec
   Interval 176.486-178.151 sec
   Interval 178.151-179.816 sec
   Interval 179.816-181.481 sec
   Interval 181.481-183.146 sec
   Interval 183.146-184.811 sec
   Interval 184.811-186.476 sec
   Interval 186.476-188.140 sec
   Interval 188.140-189.805 sec
   Interval 189.805-191.470 sec
   Interval 191.470-193.135 sec
   Interval 193.135-194.800 sec
   Interval 194.800-196.465 sec
   Interval 196.465-198.130 sec
   Interval 198.130-199.795 sec
   Interval 199.795-201.460 sec
   Interval 201.460-203.125 sec
   Interval 203.125-204.790 sec
   Interval 204.790-206.455 sec
   Interval 206.455-208.120 sec
   Interval 208.120-209.785 sec
   Interval 209.785-211.450 sec
   Interval 211.450-213.115 sec
   Interval 213.115-214.780 sec
   Interval 214.780-216.445 sec
   Interval 216.445-218.110 sec
   Interval 218.110-219.775 sec
   Interval 219.775-221.440 sec
   Interval 221.440-223.105 sec
   Interval 223.105-224.770 sec
   Interval 224.770-226.435 sec
   Interval 226.435-228.100 sec
   Interval 228.100-229.764 sec
   Interval 229.764-231.429 sec
   Interval 231.429-233.094 sec
   Interval 233.094-234.759 sec
   Interval 234.759-236.424 sec
   Interval 236.424-238.089 sec
   Interval 238.089-239.754 sec
   Interval 239.754-241.419 sec
   Interval 241.419-243.084 sec
   Interval 243.084-244.749 sec
   Interval 244.749-246.414 sec
   Interval 246.414-248.079 sec
   Interval 248.079-249.744 sec
   Interval 249.744-251.409 sec
   Interval 251.409-253.074 sec
   Interval 253.074-254.739 sec
   Interval 254.739-256.404 sec
   Interval 256.404-257.378 sec
    155 STC iterations provided
Done
Adding average EEG reference projection.
1 projection items deactivated
Average reference projection was added, but has not been applied yet. Use the apply_proj method to apply it.
Adding noise to 366/376 channels (366 channels in cov)
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 90.0 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
2 sources
blink simulated and trace stored on channel:      EOG 061
Setting up forward solutions
Computing gain matrix for transform #1/1
Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 90.0 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
1 sources
ecg simulated and trace not stored
Setting up forward solutions
Computing gain matrix for transform #1/1

Reconstruct simulated source time courses using dSPM inverse operator

Here, source time courses for auditory and visual areas are reconstructed separately and their difference is shown. This was done merely for better visual representation of source reconstruction. As expected, when high activations appear in primary auditory areas, primary visual areas will have low activations and vice versa.

method, lambda2 = 'dSPM', 1. / 9.
epochs = mne.Epochs(raw_sim, events, event_id)
inv = mne.minimum_norm.make_inverse_operator(epochs.info, fwd, noise_cov)
stc_aud = mne.minimum_norm.apply_inverse(
    epochs['auditory/left'].average(), inv, lambda2, method)
stc_vis = mne.minimum_norm.apply_inverse(
    epochs['visual/right'].average(), inv, lambda2, method)
stc_diff = stc_aud - stc_vis

brain = stc_diff.plot(subjects_dir=subjects_dir, initial_time=0.1,
                      hemi='split', views=['lat', 'med'])
  • ../../_images/sphx_glr_plot_simulated_raw_data_using_subject_anatomy_002.png
  • ../../_images/sphx_glr_plot_simulated_raw_data_using_subject_anatomy_003.png
  • ../../_images/sphx_glr_plot_simulated_raw_data_using_subject_anatomy_004.png
  • ../../_images/sphx_glr_plot_simulated_raw_data_using_subject_anatomy_005.png

Out:

320 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
Created an SSP operator (subspace dimension = 4)
4 projection items activated
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.8
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.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.
Computing SVD of whitened and weighted lead field matrix.
    largest singular value = 5.49264
    scaling factor to adjust the trace = 1.64e+19
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 5
    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 "auditory/left"...
    Picked 364 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained  99.4% variance
    Combining the current components...
    dSPM...
[done]
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 3
    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 "visual/right"...
    Picked 364 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained  99.5% variance
    Combining the current components...
    dSPM...
[done]
Using control points [ 204.9476386   323.62258808 1588.3393872 ]

References

1

Destrieux C, Fischl B, Dale A, Halgren E (2010). Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature, vol. 53(1), 1-15, NeuroImage.

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

Estimated memory usage: 1126 MB

Gallery generated by Sphinx-Gallery