Example of audiovisual stimulus in TVB with MEG evoked analysis in MNE

avevoked.py

import matplotlib as mpl
mpl.use('Agg')
import numpy as np
from scipy.io import loadmat
import mne
from mne.io.array import RawArray
try:
    from mne.io.meas_info import create_info
except ImportError:
    from mne.io.array import create_info
from tvb.simulator.lab import *

# see github.com/the-virtual-brain/tvb-data
tvb_data_path = '/home/duke/proj/vibes/tvb-data/tvb_data'

# read lead field & map region labels (because lead field was computed
# in Brainstorm, ch order doesn't match)
import json
with open('map_4d_to_mne.json', 'r') as fd:
    map_4d_to_mne = json.load(fd)
labels = loadmat(tvb_data_path + '/../tvb-meg-channel-names.mat')
mnenames = [map_4d_to_mne.get(nm[0], '?') for nm in labels['meg_names'][0]]
megnames = [c for c in mnenames if c.startswith('MEG')]
G = loadmat(tvb_data_path + '/../tvb-lead-fields.mat')['meg']
G = G[np.isfinite(G).all(axis=1)]
assert G.shape[0] == 248

# read connectivity
conn = connectivity.Connectivity.from_file(tvb_data_path + '/connectivity/connectivity_74.zip')
conn.speed = 1.0

# setup region lead field (as w/ freesurfer labels)
rmap = np.loadtxt(tvb_data_path + '/regionMapping/original_region_mapping.txt').astype(np.int32)
Gr = []
for ri in r_[:conn.weights.shape[0]]:
    Gr.append(G[:, rmap==ri].sum(axis=1))
Gr = np.array(Gr)
Gr.shape

# stimulus
eqn_t = equations.PulseTrain()
eqn_t.parameters.update({
    'T': 1e3/2.0, # 2 Hz
    'tau': 5.0,
    'onset': 500.0
    })

stimpatt = np.zeros((conn.weights.shape[0], ))
stimpatt[[35, 72]] = 0.05 # r & l V1
stimpatt[[0, 37]] = 0.5 # r & l A1

# oscillator parameters
regime = {'a': -2.5, 'b': -10.0, 'c': 0.0, 'd': 0.02,# * 8,
          'I': 0.0}

# create simulator
sim = simulator.Simulator(
    model            = models.Generic2dOscillator(**regime),
    connectivity     = conn,
    coupling         = coupling.Linear(a=0.1),
    integrator       = integrators.HeunStochastic(
        dt           = 0.5,
        noise        = noise.Additive(
            nsig     = np.array([2**-14.0,]))),
    monitors         = (monitors.TemporalAverage(period=1e3/2034.5),),
    stimulus         = patterns.StimuliRegion(
        temporal     = eqn_t,
        connectivity = conn, 
        weight       = stimpatt)
)
sim.configure()

# perform the simulation
tf = 6e4
_=next(sim(simulation_length=tf));
ys = []
ts = []
for (t,tavg), in sim(simulation_length=tf):
    if t > 0.0:
        ts.append(t)
        ys.append(tavg)
        if len(ts)%1000==0:
            print t
ts = np.array(ts)
ys = np.array(ys)

# convert stimulus to mne event array
pt = sim.stimulus(spatial_indices=0)[0]
_ts = r_[:len(pt)]
onsets = _ts[(pt[:-1]==0.)*(pt[1:]>0.)]
_ = np.ones((onsets.shape[0],))
events = c_[onsets, _*0, _]
event_id = {'stim': 1}
events.shape

# setup mne raw
meg = ys[:, 0, :, 0].dot(Gr).T
nchan = meg.shape[0]
raw = RawArray(meg/1e9,
        create_info(
            megnames,
            1e3/(ts[1] - ts[0]),
            ['mag' for _ in range(nchan)])
        )
raw.filter(5, 100)

# epoch on stimulus
epochs = mne.Epochs(raw, events, event_id, -0.1, 0.5)

# evoked
evoked = epochs.average()
fig = evoked.plot(show=True)
fig.savefig('avevoked-time-course.png')

# topomap
lay = mne.layouts.read_layout('magnesWH3600')
for i, ix in enumerate(np.r_[:0.1:20j].reshape((4, 5))):
    fig = evoked.plot_topomap(times=ix, layout=lay)
    fig.savefig('avevoked-topo-%d.png' % (i,))