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.
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
numpy.convolverather 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.
Define some variables
tvb_lib_folder = '~/Code/libraries_of_others/tvb-library' tvb_dat_folder = '~/Code/libraries_of_others/tvb-data' data_dir = '../data'
import os,sys,glob,numpy as np,pandas as pd %matplotlib inline from matplotlib import pyplot as plt import seaborn as sns sys.path += [tvb_lib_folder,tvb_dat_folder] from tvb.simulator.lab import *
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:
%%time df_tavg,df_bold = run_sim(sim_len=200000)
WARNING File 'hemispheres' not found in ZIP. CPU times: user 2min 1s, sys: 2.44 s, total: 2min 3s Wall time: 2min 11s
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):
fig, ax = plt.subplots(ncols=2,nrows=2, figsize=(12,5)) a = ax; a.set_title(df_tavg.columns) 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; a.set_title(df_tavg.columns) 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; a.set_title(df_tavg.columns) 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; a.set_title(df_tavg.columns) 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()
This is even visible to some extend in the average over nodes:
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()
The aim of the following is to reproduce the BOLD signal shown above, directly from the neural signal.
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||FirstOrderVolterra(bound=False, value=None)||A tvb.datatypes.equation object which describe the haemodynamic response function used to compute the BOLD signal.|
|period||2000.0||For the BOLD monitor, sampling period in milliseconds must be an integral multiple of 500. Typical measurment interval (repetition time TR) is between 1-3 s. If TR is 2s, then Bold period is 2000ms.|
|hrf_length||20000.0||Duration of the hrf kernel|
|variables_of_interest||||Indices of model’s variables of interest (VOI) that this monitor should record. Note that the indices should start at zero, so that if a model offers VOIs V, W and V+W, and W is selected, and this monitor should record W, then the correct index is 0.|