Added MEG, NAI, PSC, BOLD computation for auditory simulations, and json output data files

Author: Antonio Ulloa

analysis/compute_BOLD_balloon_model_auditory.py

@@ -0,0 +1,342 @@
+# ============================================================================
+#
+#                            PUBLIC DOMAIN NOTICE
+#
+#       National Institute on Deafness and Other Communication Disorders
+#
+# This software/database is a "United States Government Work" under the 
+# terms of the United States Copyright Act. It was written as part of 
+# the author's official duties as a United States Government employee and 
+# thus cannot be copyrighted. This software/database is freely available 
+# to the public for use. The NIDCD and the U.S. Government have not placed 
+# any restriction on its use or reproduction. 
+#
+# Although all reasonable efforts have been taken to ensure the accuracy 
+# and reliability of the software and data, the NIDCD and the U.S. Government 
+# do not and cannot warrant the performance or results that may be obtained 
+# by using this software or data. The NIDCD and the U.S. Government disclaim 
+# all warranties, express or implied, including warranties of performance, 
+# merchantability or fitness for any particular purpose.
+#
+# Please cite the author in any work or product based on this material.
+# 
+# ==========================================================================
+
+
+
+# ***************************************************************************
+#
+#   Large-Scale Neural Modeling software (LSNM)
+#
+#   Section on Brain Imaging and Modeling
+#   Voice, Speech and Language Branch
+#   National Institute on Deafness and Other Communication Disorders
+#   National Institutes of Health
+#
+#   This file (compute_BOLD_balloon_model_auditory.py) was created on April 1, 2016.
+#
+#
+#   Author: Antonio Ulloa
+#
+#   Last updated by Antonio Ulloa on April 1 2016
+#
+# **************************************************************************/
+
+# compute_BOLD_balloon_model_auditory.py
+#
+# Calculate and plot fMRI BOLD signal, using the
+# Balloon model, as described by Stephan et al (2007)
+# and Friston et al (2000), and
+# the BOLD signal model, as described by Stephan et al (2007) and
+# Friston et al (2000).
+# Parameters for both were taken from Friston et al (2000) and they were
+# estimated using a 2T scanner, with a TR of 2 seconds, 
+# 
+# ... using data from auditory delay-match-to-sample simulation stored in 4
+# different synaptic activity files
+# It also saves the BOLD timeseries for each and all modules in a python data file
+# (*.npy)
+#
+# The input data (synaptic activities) and the output (BOLD time-series) are numpy arrays
+# with columns in the following order:
+#
+# A1 ROI (right hemisphere, includes LSNM units and TVB nodes) 
+# A2 ROI (right hemisphere, includes LSNM units and TVB nodes)
+# ST ROI (right hemisphere, includes LSNM units and TVB nodes)
+# PF ROI (right hemisphere, includes LSNM units and TVB nodes)
+#
+# Note: only the BOLD activity of one of the synaptic activity files is computed and stored.
+
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+from scipy.integrate import odeint
+
+# define the name of the input files where the synaptic activities are stored
+# we have four files: TC-PSL, Tones-PSL, TC-DMS, Tones-DMS
+SYN_file_2 = 'output.TC_PSL/synaptic_in_ROI.npy'
+SYN_file_3 = 'output.Tones_PSL/synaptic_in_ROI.npy'
+SYN_file_4 = 'output.TC_DMS/synaptic_in_ROI.npy'
+SYN_file_1 = 'output.Tones_DMS/synaptic_in_ROI.npy'
+
+# define the name of the output file where the BOLD timeseries will be stored
+BOLD_file = 'output.Tones_DMS/lsnm_bold_balloon.npy'
+
+# define balloon model parameters...
+tau_s = 1.5           # rate constant of vasodilatory signal decay in seconds
+                      # from Friston et al, 2000
+
+tau_f = 4.5           # Time of flow-dependent elimination or feedback regulation
+                      # in seconds, from Friston et al, 2000
+
+alpha = 0.2           # Grubb's vessel stiffness exponent, from Friston et al, 2000
+
+tau_0 = 1.0           # Hemodynamic transit time in seconds
+                      # from Friston et al, 2000
+
+epsilon = 0.1         # efficacy of synaptic activity to induce the signal,
+                      # from Friston et al, 2000
+                      
+# define BOLD model parameters...
+r_0 = 25.0            # Slope of intravascular relaxation rate (Hz)
+                      # (Obata et al, 2004)
+
+nu_0 = 40.3           # Frequency offset at the outer surface of magnetized
+                      # vessels (Hz) (Obata et al, 2004)
+
+e = 1.43        # Ratio of intra- and extravascular BOLD signal at rest
+                       # (Obata et al, 2004)
+
+V_0 = 0.02            # Resting blood volume fraction, from Friston et al, 2000
+
+E_0 = 0.8             # Resting oxygen extraction fraction (Friston et al, 2000)
+
+TE = 0.040            # Echo time for a 1.5T scanner (Friston et al, 2000)
+
+# calculate BOLD model coefficients, from Friston et al (2000)
+k1 = 4.3 * nu_0 * E_0 * TE
+k2 = e * r_0 * E_0 * TE
+k3 = 1.0 - epsilon
+
+def balloon_function(y, t, syn):
+    ''' 
+    Balloon model of hemodynamic change
+    
+    '''
+    
+    # unpack initial values
+    s = y[0]
+    f = y[1]
+    v = y[2]
+    q = y[3]
+
+    x = syn[np.floor(t * synaptic_timesteps / T)]
+
+    # the balloon model equations
+    ds = epsilon * x - (1. / tau_s) * s - (1. / tau_f) * (f - 1)
+    df = s
+    dv = (1. / tau_0) * (f - v ** (1. / alpha))
+    dq = (1. / tau_0) * ((f * (1. - (1. - E_0) ** (1. / f)) / E_0) -
+                          (v ** (1. / alpha)) * (q / v))
+
+    return [ds, df, dv, dq]
+
+# define neural synaptic time interval in seconds. The simulation data is collected
+# one data point at synaptic intervals (10 simulation timesteps). Every simulation
+# timestep is equivalent to 3.5 ms.
+Ti = 0.0035 * 10
+
+# Total time of scanning experiment in seconds (timesteps X 5)
+T = 34.3
+
+# Time for one complete trial in milliseconds
+Ttrial = 2.8
+
+# the scanning happened every Tr interval below (in milliseconds). This
+# is the time needed to sample hemodynamic activity to produce
+# each fMRI image.
+Tr = 2
+
+# how many scans do you want to remove from beginning of BOLD timeseries?
+scans_to_remove = 0
+
+# read the input file that contains the synaptic activities of all ROIs
+syn1 = np.load(SYN_file_1)
+syn2 = np.load(SYN_file_2)
+syn3 = np.load(SYN_file_3)
+syn4 = np.load(SYN_file_4)
+
+print 'LENGTH OF SYNAPTIC TIME-SERIES: ', syn1[0].size
+
+# Throw away first value of each synaptic array (it is always zero)
+a1_syn1 = np.delete(syn1[0], 0)
+a2_syn1 = np.delete(syn1[1], 0)
+st_syn1 = np.delete(syn1[2], 0)
+pf_syn1 = np.delete(syn1[3], 0)
+a1_syn2 = np.delete(syn2[0], 0)
+a2_syn2 = np.delete(syn2[1], 0)
+st_syn2 = np.delete(syn2[2], 0)
+pf_syn2 = np.delete(syn2[3], 0)
+a1_syn3 = np.delete(syn3[0], 0)
+a2_syn3 = np.delete(syn3[1], 0)
+st_syn3 = np.delete(syn3[2], 0)
+pf_syn3 = np.delete(syn3[3], 0)
+a1_syn4 = np.delete(syn4[0], 0)
+a2_syn4 = np.delete(syn4[1], 0)
+st_syn4 = np.delete(syn4[2], 0)
+pf_syn4 = np.delete(syn4[3], 0)
+
+# put together all the synaptic activity array that form part of a study, only
+# for normalization purposes
+# put together synaptic arrays that correspond to the same brain region
+a1 = np.append(a1_syn1, [a1_syn2, a1_syn3, a1_syn4])
+a2 = np.append(a2_syn1, [a2_syn2, a2_syn3, a2_syn4])
+st = np.append(st_syn1, [st_syn2, st_syn3, st_syn4])
+pf = np.append(pf_syn1, [pf_syn2, pf_syn3, pf_syn4])
+
+# extract the synaptic activities corresponding to each ROI, and normalize to (0,1):
+a1_syn = (a1_syn1-a1.min()) / (a1.max() - a1.min())
+a2_syn = (a2_syn1-a2.min()) / (a2.max() - a2.min())
+st_syn = (st_syn1-st.min()) / (st.max() - st.min())
+pf_syn = (pf_syn1-pf.min()) / (pf.max() - pf.min())
+
+# Extract number of timesteps from one of the synaptic activity arrays
+synaptic_timesteps = a1_syn.size
+print 'Size of synaptic arrays: ', synaptic_timesteps
+
+# Given neural synaptic time interval and total time of scanning experiment,
+# construct a numpy array of time points (data points provided in data files)
+time_in_seconds = np.arange(0, T, Ti)
+
+# Hard coded initial conditions
+s = 0   # s, blood flow
+f = 1.  # f, blood inflow
+v = 1.  # v, venous blood volume
+q = 1.  # q, deoxyhemoglobin content
+
+# initial conditions vectors
+y_0_a1 = [s, f, v, q]      
+y_0_a2 = [s, f, v, q]      
+y_0_st = [s, f, v, q]      
+y_0_pf = [s, f, v, q]      
+
+# generate synaptic time array
+t_syn = np.arange(0, synaptic_timesteps)
+
+# generate time array for solution
+t = time_in_seconds
+
+# Use the following lines to test the balloon BOLD model
+#v1_syn[0:300] = .01          # 15 seconds do nothing
+#v1_syn[300:320] = .1       # One-second stimulus
+#v1_syn[320:920] = .01        # 30-second do nothing
+#v1_syn[920:940] = .1       # two-second stimulus
+#v1_syn[940:1540] = .01       # 30-second do nothing
+#v1_syn[1540:1560]= .1      # three-second stimulus
+#v1_syn[1560:2160]= .01       # 30-second do nothing
+#v1_syn[2160:2180]= .1      # four-second stimulus
+#v1_syn[2180:2780]= .01       # 30-second do nothing
+#v1_syn[2780:2800]= .1      # five-second stimulus
+#v1_syn[2800:3400]= .01       # 30-second do nothing
+#v1_syn[3400:3420]= .1      # one-second stimulus
+#v1_syn[3420:3959]= .01       # 30-second do nothing
+
+# solve the ODEs for given initial conditions, parameters, and timesteps
+state_a1 = odeint(balloon_function, y_0_a1, t, args=(a1_syn,) )
+state_a2 = odeint(balloon_function, y_0_a2, t, args=(a2_syn,) )
+state_st = odeint(balloon_function, y_0_st, t, args=(st_syn,) )
+state_pf = odeint(balloon_function, y_0_pf, t, args=(pf_syn,) )
+
+# Unpack the state variables used in the BOLD model
+s_a1 = state_a1[:, 0]
+f_a1 = state_a1[:, 1]
+v_a1 = state_a1[:, 2]
+q_a1 = state_a1[:, 3]        
+
+s_a2 = state_a2[:, 0]
+f_a2 = state_a2[:, 1]
+v_a2 = state_a2[:, 2]
+q_a2 = state_a2[:, 3]        
+
+s_st = state_st[:, 0]
+f_st = state_st[:, 1]
+v_st = state_st[:, 2]
+q_st = state_st[:, 3]        
+
+s_pf = state_pf[:, 0]
+f_pf = state_pf[:, 1]
+v_pf = state_pf[:, 2]
+q_pf = state_pf[:, 3]        
+
+# now, we need to calculate BOLD signal at each timestep, based on v and q obtained from solving
+# balloon model ODE above.
+a1_BOLD = np.array(V_0 * (k1 * (1. - q_a1) + k2 * (1. - q_a1 / v_a1) + k3 * (1. - v_a1)) )
+a2_BOLD = np.array(V_0 * (k1 * (1. - q_a2) + k2 * (1. - q_a2 / v_a2) + k3 * (1. - v_a2)) )
+st_BOLD = np.array(V_0 * (k1 * (1. - q_st) + k2 * (1. - q_st / v_st) + k3 * (1. - v_st)) )
+pf_BOLD = np.array(V_0 * (k1 * (1. - q_pf) + k2 * (1. - q_pf / v_pf) + k3 * (1. - v_pf)) )
+
+# The following is for display purposes only:
+# Construct a numpy array of timesteps (data points provided in data file)
+# to convert from timesteps to time in seconds we do the following:
+# Each simulation time-step equals 5 milliseconds
+# However, we are recording only once every 10 time-steps
+# Therefore, each data point in the output files represents 50 milliseconds.
+# Thus, we need to multiply the datapoint times 50 ms...
+# ... and divide by 1000 to convert to seconds
+#t = np.linspace(0, 659*50./1000., num=660)
+t = np.linspace(0, synaptic_timesteps * 35.0 / 1000., num=synaptic_timesteps)
+
+# downsample the BOLD signal arrays to produce scan rate of 2 per second
+scanning_timescale = np.arange(0, synaptic_timesteps, synaptic_timesteps / (T/Tr))
+scanning_timescale = scanning_timescale.astype(int)      # don't forget to convert indices to integer
+mr_time = t[scanning_timescale]
+a1_BOLD = a1_BOLD[scanning_timescale]
+a2_BOLD = a2_BOLD[scanning_timescale]
+st_BOLD = st_BOLD[scanning_timescale]
+pf_BOLD = pf_BOLD[scanning_timescale]
+
+print 'Size of BOLD arrays: ', a1_BOLD.size
+
+# round of mr time for display purposes
+mr_time = np.round(mr_time, decimals=0)
+
+# create a numpy array of timeseries
+lsnm_BOLD = np.array([a1_BOLD, a2_BOLD, st_BOLD, pf_BOLD])
+
+# now, save all BOLD timeseries to a single file 
+np.save(BOLD_file, lsnm_BOLD)
+
+# increase font size for display purposes
+#plt.rcParams.update({'font.size': 30})
+
+# Set up figure to plot synaptic activity
+plt.figure()
+
+plt.suptitle('SIMULATED SYNAPTIC ACTIVITY')
+
+# Plot synaptic activities
+plt.plot(t, a1_syn, label='A1')
+plt.plot(t, a2_syn, label='A2')
+plt.plot(t, st_syn, label='ST')
+plt.plot(t, pf_syn, label='PFC')
+
+plt.legend()
+
+# Set up separate figures to plot fMRI BOLD signal
+plt.figure()
+
+plt.suptitle('SIMULATED fMRI BOLD SIGNAL')
+
+plt.plot(mr_time, a1_BOLD, label='A1')
+
+plt.plot(mr_time, a2_BOLD, label='A2')
+
+plt.plot(mr_time, st_BOLD, label='ST')
+
+plt.plot(mr_time, pf_BOLD, label='PFC')
+
+plt.legend()
+
+# Show the plots on the screen
+plt.show()