JDG Lab Notebook

Applying BOLD HRFs to neural time series data

The BOLD monitor in TVB calculates the haemodynamic response (i.e. 'fMRI activity') associated with a given neural activity time series. This can be obtained as an output from a TVB simulation.

But sometimes one wants to calculate BOLD directly from a set of numbers (e.g. non-TVB simulations; from physiological recording data, etc.).

The following notes show how to do this.

Contents

Background

BOLD signal is calculated in TVB using convolution with a haemodynamic response function kernel.

Alternative ways of doing this I have seen by e.g. Gus Deco is to numerically integrate a separate ODE for the BOLD signal.

But here we will be following the convolution approach.

Jarrod Millman's fmri stats lectures here and here are extremely useful and practical background reading on convolution in the context of event-related fMRI analysis. The wikipedia page on convolution is also very good, with some nice visual demonstrations of the concepts.

In terms of the TVB implementation, the following is essentially continuing from the 'exploring the BOLD monitor' demo, which shows how to visualize the HRF kernel shape. So check that out first for a bit more context.

Now we get to an important point, and also the reason why it seems worth writing these notes out in a reasonably coherent form: When I first started looking into this, my initial plan was simply to use directly whatever function TVB uses for calculating BOLD internally.

However, after pouring over the BOLD monitor code, and asking around a bit, it became pretty clear that that probably isn't the most sensible approach.

In a nutshell, the reasons for this are that

  • the BOLD monitor has to do a number of things related to the in-step sampling of neural activity during that would simply be unnecessarily clunky to reproduce directly when one already has a complete time series
  • pandas has some nice moving window average functions that make life a lot easier
  • it seems more straightforward to use numpy.convolve rather than write out the convolution explicitly, as is done in the monitor code

The basic aim of the following notes is to construct a function that takes as input a (simulated) neural time series (temporal average monitor), and returns as output a BOLD signal that matches reasonably closely to the BOLD monitor output from the same simulation.

A secondary aim is to explain a few things here and there along the way.

Notebook Setup

Define some variables

# put some system-specific variables in the namespace ('le' dict)
%run ~/set_localenv_vars.py

nb_name = 'applying_bold_hrfs_to_neural_time_series'

# this is where I look for stuf
outdir = '%s/notebooks/%s' %(le['data_dir'],nb_name)

# this is where stuff actually lives
outdir_loc = '%s/%s' %('/mnt/titania-hecuba/mcintosh_lab/john/Data/notebooks', nb_name)
#outdir_loc = '%s/%s' %('/hecuba/mcintosh_lab/john/Data/notebooks', nb_name)

# symlink the two
import os
if not os.path.isdir(outdir_loc): os.makedirs(outdir_loc)
if not os.path.isdir(outdir): os.system('ln -s %s %s' %(outdir_loc,outdir))
                                        

# root path to tvb-library, tvb-data, etc. libraries. 
tvb_folder = le['code_dir'] + '/libraries_of_mine/github'


# stuff for workdocs-cloudfiles

aws_key = 'drowssaperucesyreva'
aws_secret = '?teytidettopsuoyevah'

Importage

# Generic imports

import os,sys,glob,h5py,itertools,multiprocessing,\
       numpy as np,pandas as pd
from datetime import datetime
from time import time
from scipy import optimize
from scipy.spatial import cKDTree
from IPython.display import display as d

# Visualization stuff

%matplotlib inline
from matplotlib import pyplot as plt
from IPython.display import Image,display as d
import seaborn as sns
sns.set(style='white')


# Neuroimaging stuff

import nibabel as nib
from nilearn.plotting import plot_roi
from dipy.tracking.utils import apply_affine

# TVB stuff

nbso,nbse = sys.stdout,sys.stderr # hack part 1/2 to keep output printing properly
sys.path.append(tvb_folder +'/tvb-library')
sys.path.append(tvb_folder + '/tvb-data')
from tvb.simulator.lab import *
sys.stdout,sys.stderr = nbso,nbse  # ...hack part 2/2
from tvb.simulator.plot.tools import plot_surface_mpl,plot_surface_mpl_mv


from copy import copy # seems that this needs to be done last otherwise it gets 
                      # replaced with numpy copy function
    
    
# aws api and workdocs-cloudfiles stuff

sys.path.append(le['ipynb_workdocs_dir'])
from new_utils import nb_fig,cloudfiles_nb

Initialize aws api and workdocs-cloudfiles folder

bucket_name = nb_name.replace('_', '-')
cnb = cloudfiles_nb('aws', [aws_key,aws_secret])
cnb.initialize_folder(bucket_name)

Load calico document extensions

%%javascript
IPython.load_extensions('calico-spell-check', 
                        'calico-document-tools', 'calico-cell-tools');

Go to output folder

os.chdir(outdir)

Ok, let's get cracking

Generate simulation data

Simple run simulation function

def run_sim(sim_len=100000,tavg_per=100,bold_per=2000):
    
    tavg_mon = monitors.TemporalAverage(period=tavg_per)
    bold_mon = monitors.Bold(period=bold_per)
    mons = (tavg_mon,bold_mon)
    mod =models.Generic2dOscillator(a=2.)
    conn = connectivity.Connectivity(load_default=True)
    cpl=coupling.Linear()
    solver=integrators.HeunDeterministic(dt=1.)
    
    sim = simulator.Simulator(model=mod,connectivity=conn,
                              integrator=solver,coupling=cpl,
                              monitors=mons)
    sim.configure()

    (tavg_time,tavg_data),(bold_time,bold_data) = sim.run(simulation_length=sim_len)

    df_tavg = pd.DataFrame(np.squeeze(tavg_data),index=tavg_time,columns=conn.region_labels)
    df_bold = pd.DataFrame(np.squeeze(bold_data),index=bold_time,columns=conn.region_labels)
    
    return df_tavg,df_bold

Generate 200 seconds of simulated data with the 2D oscillator model, and also collecting BOLD:

df_tavg,df_bold = run_sim(sim_len=200000)

Looking at a random selection of nodes from this simulation, we can see clearly that the BOLD (red) tracks the low-frequency fluctuations visible in the envelope of the neural activity (black):

(-_-) 
f = outdir + '/bold_vs_tavg_ts.png'


# Make figure

fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,5))

a = ax[0][0]; a.set_title(df_tavg.columns[0])
df_tavg.loc[10000:].iloc[:,0].plot(legend=False,alpha=0.5,ax=a,linewidth=0.5,c='k')
df_bold.loc[10000:].iloc[:,0].plot(secondary_y=True,legend=False,ax=a,c='r')

a = ax[0][1]; a.set_title(df_tavg.columns[5])
df_tavg.loc[10000:].iloc[:,5].plot(legend=False,alpha=0.5,ax=a,linewidth=0.5,c='k')
df_bold.loc[10000:].iloc[:,5].plot(secondary_y=True,legend=False,ax=a,c='r')

a = ax[1][0]; a.set_title(df_tavg.columns[25])
df_tavg.loc[10000:].iloc[:,25].plot(legend=False,alpha=0.5,ax=a,linewidth=0.5,c='k')
df_bold.loc[10000:].iloc[:,25].plot(secondary_y=True,legend=False,ax=a,c='r')

a = ax[1][1]; a.set_title(df_tavg.columns[50])
df_tavg.loc[10000:].iloc[:,50].plot(legend=False,alpha=0.5,ax=a,linewidth=0.5,c='k')
df_bold.loc[10000:].iloc[:,50].plot(secondary_y=True,legend=False,ax=a,c='r')



plt.tight_layout()
plt.savefig(f, bbox_inches='tight')
plt.close()


# Upload to cloud

cap = ''
label  = "Bold vs. temporal average time series examples"

 
im = nb_fig(f,label,cap,cnb,
            size=(1000,800),show_fignum=False)
d(im)
broken link Bold vs. temporal average time series examples.

This is even visible to some extend in the average over nodes:

(-_-) 
f = outdir + '/bold_vs_tavg_ts_mean.png'


# Make figure

fig, ax = plt.subplots(figsize=(12,3))

df_tavg.loc[10000:].mean(axis=1).plot(ax=ax,alpha=0.5,linewidth=0.5,c='k')
df_bold.loc[10000:].mean(axis=1).plot(ax=ax,secondary_y=True,c='r')

plt.tight_layout()

plt.savefig(f, bbox_inches='tight')
plt.close()


# Upload to cloud

cap = ''
label  = "Bold vs. temporal average time series examples - mean over nodes"
im = nb_fig(f,label,cap,cnb,
            size=(1000,800),show_fignum=False)
d(im)
broken link Bold vs. temporal average time series examples - mean over nodes.

The aim of the following is to reproduce the BOLD signal shown above, directly from the neural signal.

Computing BOLD signal post-hoc

First let's take a look at the BOLD monitor and its HRF functions.

The meat of the BOLD monitor is the equations object cfound in the .hrf_kernel attribute:

d(monitors.Bold())

Bold(period=2000.0, hrf_length=20000.0, hrf_kernel=FirstOrderVolterra(bound=False, value=None), variables_of_interest=[], post_expr=, pre_expr=)

The default bold monitor equation is the following first order Volterra kernel:

d(monitors.Bold().hrf_kernel.equation)

'1/3. * exp(-0.5*(var / tau_s)) * (sin(sqrt(1./tau_f - 1./(4.*tau_s**2)) * var)) / (sqrt(1./tau_f - 1./(4.*tau_s**2)))'

Borrowing from the TVB bold monitor demo, we can plot the shape of this function (and ther other available HRF kernels) as follows:

def tvb_bold_hrf(times=None,mon_params={},return_df=True):
    
    """
    times   = array of time steps to evaluate respons (unit = ms)
    mon_params = bold monitor parameters dict
    """
    
    bold_mon = monitors.Bold(**mon_params)

    hrfkeqn = bold_mon.hrf_kernel
    label = hrfkeqn.__class__.__name__
    
    if times is None:
        times = np.linspace(0,bold_mon.hrf_length,100)/1000.
    
    hrfkeqn._set_pattern(times)
    
    res = hrfkeqn.pattern 
    # (note 'pattern' here is a decorator; i.e. this expression
    # is actually evaluating a function
    # (computing hrf equation for the specified times and parameters)
    
    if return_df:
        df = pd.DataFrame(res,index=times,columns=[label])
        return df
    else:
        return res

    # Alternative; gives slightly different result though (?)
    #  hrf_params = bold_mon.hrf_kernel.parameters
    #  eqn_str = bold_mon.hrf_kernel.equation
    #  hrf_params.update({'var': times})
    #  res = evaluate(eqn_str,local_dict=hrf_params)
(-_-) 
f = outdir + '/hrf_kernel_functions.png'


# Make figure

hrf_kernels = [equations.FirstOrderVolterra(),
               equations.DoubleExponential(),
               equations.MixtureOfGammas(),
               equations.Gamma()]

fig, ax = plt.subplots()

for k in hrf_kernels:
    tvb_bold_hrf(mon_params={'hrf_kernel': k}).plot(ax=ax)   
    

plt.tight_layout()
plt.savefig(f, bbox_inches='tight')
plt.close()


# Upload to cloud

cap = ''
label  = "HRF Kernel functions"

im = nb_fig(f,label,cap,cnb,
            size=(1000,800),show_fignum=False)
d(im)
broken link HRF Kernel functions.

We'll now use the above function to construct our convolution function for computing bold from neural activity time series

def compute_bold_from_ts(data,times,hrf_len=500,scale_time=1E-3,n_subsamp=20,winlen=10):
    
    hrf_time = times[:hrf_len]*scale_time
    hrf_at_time = tvb_bold_hrf(times=hrf_time,return_df=False)

    data_conv = convolve(data,hrf_at_time)
    n_to_remove = len(hrf_at_time) - 1
    data_conv = data_conv[:-n_to_remove]
    
    idxs = np.arange(0,times.shape[0],n_subsamp)
    
    df_data_conv = pd.DataFrame(data_conv,index=times)
    
    df_data_conv_rwm = df_data_conv.rolling(window=winlen).mean()
    
    df_data_conv_rwm_ss = df_data_conv_rwm.iloc[idxs,:]
    
    return df_data_conv_rwm,df_data_conv_rwm_ss

Applying this to the simulated data:

times = df_tavg.index.values    # time stamps
stuff = {}
for c in range(df_tavg.columns.shape[0]):
    data = df_tavg.iloc[:,c].values
    stuff[c] = compute_bold_from_ts(data,times)
df_data_conv_rwm_all = pd.concat({k: v[0][0] for k,v in stuff.items()},axis=1)
df_data_conv_rwm_ss_all = pd.concat({k: v[1][0] for k,v in stuff.items()},axis=1)
(-_-) 
f = outdir + '/bold_vs_posthoc_bold_ts.png' 


# Make figure

fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,5))

a = ax[0][0]; a.set_title(df_tavg.columns[0])
df_data_conv_rwm_ss_all.loc[20000:].iloc[:,0].plot(legend=False,alpha=0.5,ax=a,c='k')
df_bold.loc[20000:].iloc[:,0].plot(secondary_y=True,legend=False,ax=a,c='r')

a = ax[0][1]; a.set_title(df_tavg.columns[5])
df_data_conv_rwm_ss_all.loc[20000:].iloc[:,5].plot(legend=False,alpha=0.5,ax=a,c='k')
df_bold.loc[20000:].iloc[:,5].plot(secondary_y=True,legend=False,ax=a,c='r')

a = ax[1][0]; a.set_title(df_tavg.columns[25])
df_data_conv_rwm_ss_all.loc[20000:].iloc[:,25].plot(legend=False,alpha=0.5,ax=a,c='k')
df_bold.loc[20000:].iloc[:,25].plot(secondary_y=True,legend=False,ax=a,c='r')

a = ax[1][1]; a.set_title(df_tavg.columns[50])
df_data_conv_rwm_ss_all.loc[20000:].iloc[:,50].plot(legend=False,alpha=0.5,ax=a,c='k')
df_bold.loc[20000:].iloc[:,50].plot(secondary_y=True,legend=False,ax=a,c='r')


plt.tight_layout()

plt.savefig(f, bbox_inches='tight')
plt.close()


# Upload to cloud

cap = ''
label  = "Post-hoc computed BOLD time series"

im = nb_fig(f,label,cap,cnb,
            size=(1000,800),show_fignum=False)
d(im)
broken link Post-hoc computed BOLD time series.

Excellent. We can see that the new post-hoc computed BOLD activity, and the outputs of the TVB BOLD monitor match up very nicely.

Additional Notes

  • Absolute value of the bold monitor output at the new function here arent currently matching. Hence use of secondary_y axis in above plots. Not quite sure why this is.
  • With this implementation, the post-hoc computed BOLD activity only starts to follow the BOLD monitor activity after the initial transient of a few seconds. This is probably due to either the moving average window length or the convolution window length; neither of which were tuned to any great deal in the above.
  • Effects of moving average window length and convolution window length should perhaps be investigated.