In [1]:
#import needed functions
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from nipy.modalities.fmri.hemodynamic_models import spm_hrf, compute_regressor
from statsmodels.tsa.arima_process import arma_generate_sample
from nipy.modalities.fmri.design_matrix import make_dmtx
from nipy.modalities.fmri.experimental_paradigm import (EventRelatedParadigm,
BlockParadigm)
In [2]:
# frame times
TR= 2
nscans_trial= 30
trial_duration=TR*nscans_trial # 60
ntrials_run= 12
nscans_run=nscans_trial*ntrials_run # 12*30=360
run_duration= TR*nscans_run # 360*2=720
frametimes = np.linspace(0, (nscans_run - 1) * TR, nscans_run)
frametimes
Out[2]:
In [3]:
TR= 2
nscans_trial= 30
trial_duration=TR*nscans_trial # 60
ntrials_run= 12
nscans_run=nscans_trial*ntrials_run # 12*30=360
run_duration= TR*nscans_run # 360*2=720
In [4]:
# experimental paradigm
conditions = ['a', 'b', 'c','c', 'a', 'b', 'c', 'a', 'b', 'c','a','b']
#conditions = ['c', 'c', 'c', 'c','c', 'c', 'c', 'c', 'c', 'c', 'c','c']
#onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
onsets=np.arange(0,run_duration,trial_duration)
motion = np.cumsum(np.random.randn(nscans_run, 6), 0)
add_reg_names = ['tx', 'ty', 'tz', 'rx', 'ry', 'rz']
paradigm = EventRelatedParadigm(conditions, onsets)
#hrf_model = 'fir'
#design_matrix = make_dmtx(frametimes,
# paradigm,
# hrf_model=hrf_model,
# drift_model='polynomial',
# drift_order=3,
# fir_delays=np.arange(1, 6))
hrf_model = 'spm'
design_matrix = make_dmtx(frametimes,
paradigm,
hrf_model=hrf_model,
drift_model='polynomial',
drift_order=3,
add_regs=motion,
add_reg_names=add_reg_names)
# plot the matrix
fig = plt.figure(figsize=(14, 16))
ax = plt.subplot(1, 1, 1)
design_matrix.show(ax=ax)
ax.set_title('Event-related design matrix', fontsize=12)
Out[4]:
In [6]:
# generate the errors (AR1)
ar = [1, 0.5]
ma = [1, 0.0]
epsilon_ar1= arma_generate_sample(ar,ma,nscans_run)
# generate the data
noise_scale = 0.3
activation_param = [90.0,60.0,30.0]
motion_param = [.00,.00,.00,.00,.00,.00] # ignore
drift_param = [0.0,0.0,0.0] # ignore
constant_param= [0.0] # zero mean centered
#concatenate cofficients
Beta =np.array([activation_param +motion_param+ drift_param+ constant_param])
# generate BOLD response timecourse for a single voxel
y=np.dot(Beta,design_matrix.matrix.T) + epsilon_ar1*noise_scale
In [8]:
#plot timecourse for one voxel
plt.figure(figsize=(16,4))
plt.axis([0,nscans_run,-1.,3.])
plt.plot(y.T,'--')
plt.xlabel('Scan',fontsize=12)
plt.ylabel('BOLD response',fontsize=12)
for i in onsets/2:
plt.plot(np.arange(i,i+nscans_trial),y.T[i:i+nscans_trial])
In [10]:
# conditions a
time_onsets =onsets[[i for i, j in enumerate(conditions) if j == 'a']]/2
y_cond_a=y.T[time_onsets[0]:time_onsets[0]+nscans_trial]
for i in time_onsets[1:]:
y_cond_a=np.hstack((y_cond_a,y.T[i:i+nscans_trial]))
plt.figure(figsize=(12,9))
plt.axis([0,nscans_run,-1.,3.])
plt.subplot(2,1,1)
plt.plot(y_cond_a)
plt.xlabel('Scans',fontsize=14)
plt.ylabel('BOLD response',fontsize=12)
plt.subplot(2,1,2)
plt.plot(np.mean(y_cond_a,1))
plt.xlabel('Scans ',fontsize=12)
plt.ylabel('BOLD ',fontsize=12)
plt.tight_layout()
In [13]:
# conditions b
time_onsets =onsets[[i for i, j in enumerate(conditions) if j == 'b']]/2
y_cond_b=y.T[time_onsets[0]:time_onsets[0]+nscans_trial]
for i in time_onsets[1:]:
y_cond_b=np.hstack((y_cond_b,y.T[i:i+nscans_trial]))
plt.figure(figsize=(12,9))
plt.axis([0,nscans_run,-1.,3.])
plt.subplot(2,1,1)
plt.plot(y_cond_b)
plt.xlabel('Scans',fontsize=14)
plt.ylabel('BOLD response',fontsize=12)
plt.subplot(2,1,2)
plt.plot(np.mean(y_cond_b,1))
plt.xlabel('Scans ',fontsize=12)
plt.ylabel('BOLD ',fontsize=12)
plt.tight_layout()
In [12]:
# conditions c
time_onsets =onsets[[i for i, j in enumerate(conditions) if j == 'c']]/2
y_cond_c=y.T[time_onsets[0]:time_onsets[0]+nscans_trial]
for i in time_onsets[1:]:
y_cond_c=np.hstack((y_cond_c,y.T[i:i+nscans_trial]))
plt.figure(figsize=(12,9))
plt.axis([0,nscans_run,-1.,3.])
plt.subplot(2,1,1)
plt.plot(y_cond_c)
plt.xlabel('Scans',fontsize=14)
plt.ylabel('BOLD response',fontsize=12)
plt.subplot(2,1,2)
plt.plot(np.mean(y_cond_c,1))
plt.xlabel('Scans ',fontsize=12)
plt.ylabel('BOLD ',fontsize=12)
plt.tight_layout()