Day 12: Simulating amygdala activity

Today I generated simulated evoked data for right amygdala activation. There are several tests I need to perform
with simulated data.

1) Make sure that simulated data from a subcortical space can be accurately located.  
2) Make sure that a subcortical activation is not detected as a cortical activation.  
3) Make sure that a cortical activation is not detected as a subcortical activation.  

The data I generated today will help me with the first two tests.

How does this work? Source localization breaks the brain into a series of vertices in 3D space, generally a
subsample of the actual number of voxels available from the MRI. It then calculates a gain matrix which determines
the relationship between each vertex and each EEG or MEG sensor on the scalp. To simulate data, you can choose
one or more vertices, assign a time series to it, and then multiply it by the gain matrix to see what that source
would look like if it were measured by an array of sensors. The code to do this in MNE is below, as well as a
plot of the evoked data (one black line per sensor). For comparison, I include a plot of the simulated time series.

In [1]:

%matplotlib inline
%load '/Users/Alan/PythonEEG/gen_sim_subcortical.py'

In [3]:

import numpy as np
import matplotlib.pyplot as plt
import nibabel as nib
import mne
from mne import read_forward_solution
from mne.datasets import spm_face
from mne.time_frequency import morlet
from mne.viz import plot_evoked

mne.set_log_level(False)

# Use spm data
data_path = spm_face.data_path()

# Specify relevant files
epo_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D_epo.fif'
fwd_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D_vol-7-fwd.fif'
seg_fname = data_path + '/subjects/spm/mri/aparc+aseg.mgz'

# Read epoch data
epo = mne.read_epochs(epo_fname)

# Compute the noise covariance
cov = mne.compute_covariance(epo.crop(None, 0, copy=True))

# Read the forward solution
fwd = read_forward_solution(fwd_fname)

# Extract the source space
src = fwd['src']

# Set up a template for simulated evoked data
evoked_template = epo['faces'].average()

###############################################################################
# Read the segmentation data
seg = nib.load(seg_fname)
seg_data = seg.get_data()
seg_hdr = seg.get_header()

# Get the indices to the right amygdala
mri_ix = seg_data.flatten() == 54

# Use the interpolation matrix to transform into source space
interp = src[0]['interpolator']
src_val = interp.transpose().dot(mri_ix)

# Use only vertices that are in use
src_val *= src[0]['inuse']
src_ix = np.where(src_val > 0)[0]

# Store vertices as a Label
amg_lab = mne.Label(src_ix, src[0]['rr'][src_ix], None, 'rh')

###############################################################################
# Generate source time courses and the correspond evoked data
snr = 6  # dB
times = evoked_template.times
tmin = times.min()
sfreq = evoked_template.info['sfreq']
tstep = 1. / sfreq

# Generate times series from 2 Morlet wavelets
stc_data = np.zeros((1, len(times)))
Ws = morlet(sfreq, [10], n_cycles=2)
stc_data[0, :len(Ws[0])] = np.real(Ws[0])
stc_data *= 100 * 1e-9  # use nAm as unit

# Use a random amygdala vertex as a source
vertno = np.random.randint(0, len(amg_lab.vertices))
vertno = np.array([amg_lab.vertices[vertno]])

# Set up the source estimate
stc = mne.SourceEstimate(stc_data, vertno, tmin, tstep, 'spm')

###############################################################################
# Generate noisy evoked data

# Compute the data in sensor space
gain = np.array([fwd['sol']['data'][:, vertno[0]]])
data = np.dot(gain.T, stc.data)

# Set up the evoked data
evoked = evoked_template.copy()
evoked.pick_channels(fwd['sol']['row_names'])

# Add filtered noise to the evoked data
noise = mne.simulation.generate_noise_evoked(evoked, cov)
noise.data = mne.filter.low_pass_filter(noise.data, sfreq, 40)
evoked = mne.simulation.add_noise_evoked(evoked, noise, snr)

# Save the evoked data
sim_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D_' + \
                        'sim-rh-amyg_evo.fif'
evoked.save(sim_fname)

# Plot the evoked data
plot_evoked(evoked)

plt.show()

In [11]:

plt.plot(stc.data[0])

Out[11]:

[<matplotlib.lines.Line2D at 0x1104bd910>]