Day 10: Comparing Approaches

Today I compared the activation of the right amygdala calculated as part of a whole brain volume (Day 9)
with the activation calculated by merging a cortical surface with a subcortical volume (Day 8). The
plots are included below, with grey bars indicating regions of significant difference (fdr corrected).

Maybe its just Monday, but I found it really difficult to identify vertices from the whole brain volume
that belong to the amygdala. The original MRI coordinates are in 256 x 256 x 256 space, while the
vertices select 4608 vertices from a space of 22 x 28 x 23. I finally settled on converting both to MRI
coordinates, then defining a bounding cube for the right amygdala, and averaging together any vertices
that fell inside that cube. The amygdala is not cube shaped, so I probably included more vertices than
actually belong. My addition to the Day 9 code is included below.

Tomorrow I'll consult with the mne-python community and see if the similarities between these two methods
might be due to some source of error. I've also received some Matlab code from a collaborator outlining
another method of combining forward operators, I'd still like to test these methods on simulated data
and I'd like to see if these results can be reproduced using EEG data.

In [3]:

%matplotlib inline
from IPython.display import Image
Image('/Users/Alan/Documents/wholebrain.png')

Out[3]:

In [4]:

Image('/Users/Alan/Documents/merge.png')

Out[4]:

In []:

#########################################
# average the data for the right amygdala

# get positions for the right amygdala
amyg_pos = mne.source_space.get_segment_positions('spm', 'Right-Amygdala',
                                                  subjects_dir=subjects_dir)
rr = amyg_pos['rr']

# get positions for these sources
src_pos = fwd['src'][0]['rr'][fwd['src'][0]['vertno']]
x = src_pos

# find sources that fit inside the right amygdala bounding cube
ix = np.mean((x[:, 0]>rr[:, 0].min(), x[:, 0]<rr[:, 0].max(), 
              x[:, 1]>rr[:, 1].min(), x[:, 1]<rr[:, 1].max(),
              x[:, 2]>rr[:, 2].min(), x[:, 2]<rr[:, 2].max()), 0) == 1

# get data for each event type
x = X[epochs.events[:, 2]==1][:, ix].mean(1)
y = X[epochs.events[:, 2]==2][:, ix].mean(1)

# compare the two conditions
t, p = stats.ttest_ind(x, y)

# correct for multiple comparisons
sig, pcorr = mne.stats.fdr_correction(p)

#########################
# plot the two potentials

dt = epochs.times

ax = plt.axes()

l1, = ax.plot(dt, x.mean(0))
l2, = ax.plot(dt, y.mean(0))

ylim = ax.get_ylim()

ax.fill_between(dt, ylim[1]*np.ones(dt.shape), ylim[0]*np.ones(dt.shape),
        sig, facecolor='k', alpha=0.3)
ax.set_ylim(ylim)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Potential (uV)')
ax.set_title('Right Amygdala Activation')
ax.legend((l1, l2), ('Faces', 'Scrambled'))

plt.show()

# print run time
print time.time()-t0