Day 25: Good news, everyone!

Today I continued my validation of the combined source space using simulated data. This time, I used a source space
that contained the bilateral cortical surface as well as the bilateral amygdala. Then I generated two waveforms and
assigned each to a random voxel within the left or right VLPFC. Finally, I fed this data through the inverse
operator and measured the response at the left and right vlpfc and the left and right amygdala.

The results look very promising. The left and right VLPFC do a good job accurately representing the simulated
waveforms, while the left and right amygdala show very little activation (y-axes are normalized for comparison).

The next step is to simulate a source from within the amygdala and measure the cortical activation. This will require
some legwork, since the simulation functions currently support cortical surfaces. I will need to write functions that
work specifically with volume source spaces.

In [6]:

%matplotlib inline
%load spm_test_subcort.py

In [9]:

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import mne
from mne.datasets import spm_face

# set the random state
r = 2

# set data files
data_path = spm_face.data_path()
labels_path = data_path + '/MEG/spm/labels/%s.label'
fstring = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D'
fwd_fname = fstring + '_oct-6-amyg-vol-fwd.fif'
evo_fname = fstring + '-ave.fif'
cov_fname = fstring + '-cov.fif'

# read evoked data
evo = mne.read_evokeds(evo_fname, condition=0)
sfreq = evo.info['sfreq']
times = evo.times

# read the cortical parcellation for left and right VLPFC
label_names = ['parstriangularis-lh', 'parstriangularis-rh']
labels = [mne.read_label(labels_path % ln) for ln in label_names]

# read the forward model
fwd = mne.read_forward_solution(fwd_fname, force_fixed=True)

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

# read noise covariance
cov = mne.read_cov(cov_fname)

# generate two waveforms
stc_data = np.zeros((2, len(times)))
ws = mne.time_frequency.morlet(sfreq, [3, 10], n_cycles=[1, 1.5])
stc_data[0][:len(ws[0])] = np.real(ws[0])
stc_data[1][:len(ws[1])] = np.real(ws[1])
stc_data[0] = np.roll(stc_data[0], 100)
stc_data[1] = np.roll(stc_data[1], 200)
stc_data *= 100 * 1e-9

# generate simulated sparse source estimates
stc = mne.simulation.generate_sparse_stc(src, labels, stc_data, times.min(), 
                                         1./sfreq, random_state=r)

# set parameters for inverse operator
snr = 10
lambda2 = 1. / snr**2

# generate simulated evoked data
sim = mne.simulation.generate_evoked(fwd, stc, evo, cov, snr, tmin=0, tmax=0.2,
                                     random_state=r)

# make the inverse operator
inv = mne.minimum_norm.make_inverse_operator(evo.info, fwd, cov, loose=None, 
                                             depth=None, fixed=True)

# apply the inverse operator
stc = mne.minimum_norm.apply_inverse(sim, inv, lambda2)

# right hemisphere : 1738 vertices
# left hemisphere : 1840 vertices

# set the plot parameters
t = stc.times*1000
ylim = (-0.1, 0.4)
xlim = (t[0], t[-1])
sz = 'xx-large'

# plot results for bilateral VLPFC and amygdala
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(221)
ix = np.in1d(stc.lh_vertno, labels[0].vertices)
x = stc.lh_data[ix].mean(0)
ax.plot(t, x, 'b')
ax.plot(t, stc_data[0]*1e7, 'g')
ax.set_title('Left VLPFC', size=sz)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_xticks([])
ax.set_ylabel('Source Estimate', size=sz)

ax = fig.add_subplot(222)
ix = np.in1d(stc.rh_vertno, labels[1].vertices)
x = stc.rh_data[ix].mean(0)
l1, = ax.plot(t, x, 'b')
l2, = ax.plot(t, stc_data[1]*1e7, 'g')
ax.set_title('Right VLPFC', size=sz)
ax.legend([l1, l2], ['Evoked', 'Simulated'], fontsize=sz, loc=2)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_xticks([])
ax.set_yticks([])

ax = fig.add_subplot(223)
x = stc.data[9930:].mean(0)
ax.plot(t, x, 'b')
ax.set_title('Left Amygdala', size=sz)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_xlabel('Time (ms)', size=sz)
ax.set_ylabel('Source Estimate', size=sz)

ax = fig.add_subplot(224)
x = stc.data[8192:9930].mean(0)
ax.plot(t, x, 'b')
ax.set_title('Right Amygdala', size=sz)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_yticks([])
ax.set_xlabel('Time (ms)', size=sz)

plt.show()