Large-scale neural model of visual short-term memory (Ulloa, Horwitz 2016; Horwitz, et al. 2005,...)

 Download zip file 
Help downloading and running models
Large-scale neural model of visual short term memory embedded into a 998-node connectome. The model simulates electrical activity across neuronal populations of a number of brain regions and converts that activity into fMRI and MEG time-series. The model uses a neural simulator developed at the Brain Imaging and Modeling Section of the National Institutes of Health.
1 . Tagamets MA, Horwitz B (1998) Integrating electrophysiological and anatomical experimental data to create a large-scale model that simulates a delayed match-to-sample human brain imaging study. Cereb Cortex 8:310-20 [PubMed]
2 . Ulloa A, Horwitz B (2016) Embedding Task-Based Neural Models into a Connectome-Based Model of the Cerebral Cortex. Front Neuroinform 10:32 [PubMed]
3 . Horwitz B, Warner B, Fitzer J, Tagamets MA, Husain FT, Long TW (2005) Investigating the neural basis for functional and effective connectivity. Application to fMRI. Philos Trans R Soc Lond B Biol Sci 360:1093-108 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Prefrontal cortex (PFC);
Cell Type(s):
Gap Junctions:
Simulation Environment: Python;
Model Concept(s): Working memory;
Implementer(s): Ulloa, Antonio [antonio.ulloa at];
# ============================================================================
#                            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 ( was created on April 17, 2015.
#   Author: Antonio Ulloa
#   Last updated by Antonio Ulloa on August 11 2015
#   Based on computer code originally developed by Barry Horwitz et al
# **************************************************************************/

# Calculate and plot fMRI BOLD signal using a convolution with hemodynamic response
# function modeled as a poisson time-series, based on data from visual
# delay-match-to-sample simulation. 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:
# V1 ROI (right hemisphere, includes LSNM units and TVB nodes) 
# V4 ROI (right hemisphere, includes LSNM units and TVB nodes)
# IT ROI (right hemisphere, includes LSNM units and TVB nodes)
# FS ROI (right hemisphere, includes LSNM units and TVB nodes)
# D1 ROI (right hemisphere, includes LSNM units and TVB nodes)
# D2 ROI (right hemisphere, includes LSNM units and TVB nodes)
# FR ROI (right hemisphere, includes LSNM units and TVB nodes)
# IT ROI (left hemisphere, contains only  TVB nodes)

import numpy as np
import matplotlib.pyplot as plt

from scipy.stats import poisson

# define the name of the input file where the synaptic activities are stored
SYN_file  = 'synaptic_in_ROI.npy'

# define the name of the output file where the BOLD timeseries will be stored
BOLD_file = 'lsnm_bold_poisson.npy'

# 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 5 ms.
Ti = 0.005 * 10

# define constant needed for hemodynamic function (in milliseconds)
lambda_ = 6

# Total time of scanning experiment in seconds (timesteps X 5)
T = 198

# Time for one complete trial in seconds
Ttrial = 5.5

# 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 = 8

# read the input file that contains the synaptic activities of all ROIs
syn = np.load(SYN_file)

# extract the synaptic activities corresponding to each ROI:
v1_syn = syn[0]
v4_syn = syn[1]
it_syn = syn[2]
fs_syn = syn[3]
d1_syn = syn[4]
d2_syn = syn[5]
fr_syn = syn[6]
lit_syn= syn[7]

# Extract number of timesteps from one of the synaptic activity arrays
synaptic_timesteps = v1_syn.size

# 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, Tr)

# the following calculates a Poisson distribution (that will represent a hemodynamic
# function, given lambda (the Poisson time constant characterizing width and height
# of hemodynamic function), and tau (the time step)
# if you would do it manually you would do the following:
#h = [lambda_ ** tau * m.exp(-lambda_) / m.factorial(tau) for tau in time_in_seconds]
h = poisson.pmf(time_in_seconds, lambda_)

# the following calculates the impulse response (convolution kernel) of the gamma
# function, approximating the BOLD response (Boynton model).

# rescale the array containing the poisson to increase its size and match the size of
# the synaptic activity array (using linear interpolation)
scanning_timescale = np.arange(0, synaptic_timesteps, synaptic_timesteps / (T/Tr))
synaptic_timescale = np.arange(0, synaptic_timesteps)

h = np.interp(synaptic_timescale, scanning_timescale, h)

# now, we need to convolve the synaptic activity with a hemodynamic delay
# function and sample the array at Tr regular intervals (back to the scanning
# timescale)
v1_BOLD = np.convolve(v1_syn, h)[scanning_timescale]
v4_BOLD = np.convolve(v4_syn, h)[scanning_timescale]
it_BOLD = np.convolve(it_syn, h)[scanning_timescale]
d1_BOLD = np.convolve(d1_syn, h)[scanning_timescale]
d2_BOLD = np.convolve(d2_syn, h)[scanning_timescale]
fs_BOLD = np.convolve(fs_syn, h)[scanning_timescale]
fr_BOLD = np.convolve(fr_syn, h)[scanning_timescale]
lit_BOLD= np.convolve(lit_syn,h)[scanning_timescale]

# now we are going to remove the first trial
# estimate how many 'synaptic ticks' there are in each trial
synaptic_ticks = Ttrial/Ti
# estimate how many 'MR ticks' there are in each trial
mr_ticks = round(Ttrial/Tr)

# remove first few scans from BOLD signal array (to eliminate edge effects from
# convolution)
v1_BOLD = np.delete(v1_BOLD, np.arange(scans_to_remove))
v4_BOLD = np.delete(v4_BOLD, np.arange(scans_to_remove))
it_BOLD = np.delete(it_BOLD, np.arange(scans_to_remove))
d1_BOLD = np.delete(d1_BOLD, np.arange(scans_to_remove))
d2_BOLD = np.delete(d2_BOLD, np.arange(scans_to_remove))
fs_BOLD = np.delete(fs_BOLD, np.arange(scans_to_remove))
fr_BOLD = np.delete(fr_BOLD, np.arange(scans_to_remove))
lit_BOLD= np.delete(lit_BOLD,np.arange(scans_to_remove))

# ...and normalize the BOLD signal of each module (convert to percentage signal change)
#v1_BOLD = v1_BOLD / np.mean(v1_BOLD) * 100. - 100.
#v4_BOLD = v4_BOLD / np.mean(v4_BOLD) * 100. - 100.
#it_BOLD = it_BOLD / np.mean(it_BOLD) * 100. - 100.
#d1_BOLD = d1_BOLD / np.mean(d1_BOLD) * 100. - 100.
#d2_BOLD = d2_BOLD / np.mean(d2_BOLD) * 100. - 100.
#fs_BOLD = fs_BOLD / np.mean(fs_BOLD) * 100. - 100.
#fr_BOLD = fr_BOLD / np.mean(fr_BOLD) * 100. - 100.

# pack BOLD time-series in preparation for saving to output file
lsnm_BOLD = np.array([v1_BOLD, v4_BOLD, it_BOLD,
                      fs_BOLD, d1_BOLD, d2_BOLD, fr_BOLD,
                      lit_BOLD ])

print lsnm_BOLD.shape

# now, save all BOLD timeseries to a single file 
# Please note that we are saving the original BOLD time-series, before removing the
# edge effects, lsnm_BOLD)

# Set up figure to plot synaptic activity


# Plot synaptic activities

# Set up separate figures to plot fMRI BOLD signal


plt.plot(v1_BOLD, linewidth=3.0, color='yellow')

print v1_BOLD.shape



plt.plot(v4_BOLD, linewidth=3.0, color='green')


plt.plot(it_BOLD, linewidth=3.0, color='blue')


plt.plot(d1_BOLD, linewidth=3.0, color='red')


plt.plot(lit_BOLD, linewidth=3.0, color='pink')

# Show the plots on the screen

Loading data, please wait...