Day 24: Cortical spread for a point of reference

Today as a sanity check, I wanted to see how well simulated sources are localized on the cortex alone, in the absence
of subcortical volumes. The code I used is below. In each case, simulated waveforms were assigned to a random vertex
in both the left (blue waveform) and right (green waveform) parstriangularis (subregions of the
ventrolateral prefrontal cortex, which I'm doing my thesis work on). The resulting source maps are included below.
I couldn't figure out how to integrate mayavi windows into ipython, so I commented out the code and included screen
shots.

In each case, it looks like there's some spreading outside the region of interest. The magnitude and extent of this
spread will serve as a good reference for comparing activations in subcortical regions.

In [3]:

%matplotlib inline
%load spm_test_cortical.py

In [4]:

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

mne.set_log_level(False)

data_path = spm_face.data_path()
fwd_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D-oct-6-fwd.fif'
evo_fname = fwd_fname.replace('-oct-6-fwd.fif', '-ave.fif')
cov_fname = evo_fname.replace('-ave.fif', '-cov.fif')

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

# read the cortical parcellation
label_names = ['parstriangularis-lh', 'parstriangularis-rh']
labels = [mne.read_label(data_path + '/MEG/spm/labels/%s.label' % 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

stc = mne.simulation.generate_sparse_stc(src, labels, stc_data, times.min(), 
                                         1./sfreq)

plt.plot(np.vstack((times, times)).T, stc_data.T)

snr = 10
lambda2 = 1. / snr**2

sim = mne.simulation.generate_evoked(fwd, stc, evo, cov, snr, tmin=0, tmax=0.2)

inv = mne.minimum_norm.make_inverse_operator(evo.info, fwd, cov, loose=None, 
                                             depth=None, fixed=True)

stc = mne.minimum_norm.apply_inverse(sim, inv, lambda2)

'''
rh = stc.plot(hemi='rh', time_viewer=True)
rh.add_label(data_path + '/MEG/spm/labels/parstriangularis-rh.label',
             alpha=0.3, color='green')

lh = stc.plot(hemi='lh', time_viewer=True)
lh.add_label(data_path + '/MEG/spm/labels/parstriangularis-lh.label',
             alpha=0.3, color='green')
'''
plt.show()

/home/alan/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/mne/utils.py:493: DeprecationWarning: Function read_evoked is deprecated; 'read_evoked' will be removed in v0.9. Use 'read_evokeds.'
  warnings.warn(msg, category=DeprecationWarning)

In [6]:

from IPython.display import Image
Image('lh.png')

Out[6]:

In [7]:

Image('rh.png')

Out[7]: