Prosthetic electrostimulation for information flow repair in a neocortical simulation (Kerr 2012)

 Download zip file   Auto-launch 
Help downloading and running models
This model is an extension of a model ( ) recently published in Frontiers in Computational Neuroscience. This model consists of 4700 event-driven, rule-based neurons, wired according to anatomical data, and driven by both white-noise synaptic inputs and a sensory signal recorded from a rat thalamus. Its purpose is to explore the effects of cortical damage, along with the repair of this damage via a neuroprosthesis.
1 . Kerr CC, Neymotin SA, Chadderdon GL, Fietkiewicz CT, Francis JT, Lytton WW (2012) Electrostimulation as a prosthesis for repair of information flow in a computer model of neocortex. IEEE Trans Neural Syst Rehabil Eng 20:153-60 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; Neocortex fast spiking (FS) interneuron; Neocortex spiny stellate cell;
Channel(s): I Chloride; I Sodium; I Potassium;
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA; Gaba;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Deep brain stimulation; Information transfer; Brain Rhythms;
Implementer(s): Lytton, William [bill.lytton at]; Neymotin, Sam [Samuel.Neymotin at]; Kerr, Cliff [cliffk at];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; GabaA; AMPA; NMDA; Gaba; I Chloride; I Sodium; I Potassium; Gaba; Glutamate;
infot.mod *
intf6_.mod *
intfsw.mod *
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
decmat.hoc *
decnqs.hoc *
default.hoc *
drline.hoc *
infot.hoc *
local.hoc *
misc.h *
ratlfp.dat *
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
xgetargs.hoc *
This script generates an output figure for the NEST paper -- comparing the 
rasters in the healthy, damaged, and restored cases. It's mostly copied-and-
pasted out of, and bears an even stronger resemblance to

It also plots LFP time series, so the name is not very informative. Sorry.

Usage: just run!, but choose cortical or thalamic with the switch on line 22.

Version: 2012dec10

# Import useful tools -- not sure how many of these will actually be used
import time; tic=time.clock()
from scipy import loadtxt, size, shape, zeros, mod, floor, array, ceil
from pylab import figure, plot, scatter, xlabel, ylabel, xlim, ylim, hold, figtext, r_, show

whichlfp=4 # Which LFP to use, default 4 = from all layers
calculate=1 # Whether or not to recalculate

if calculate:
    dr='/home/cliffk/bill/ckintf/ck/1109' # Directory
    nfiles=size(inputstems) # Number of files
    killtime=0 # Number of time points to remove from beginning and end
    fs=200 # Sampling rate
    results=list() # Initialize list to store results
    changedindices=1; # WARNING 2011mar21 -- this is needed since Bill/Sam seem to have changed which indices correspond to which cells 
    popdata=list() # For storing spikes per population
    for i in range(nfiles):
        popdata.append([]) # Add an element to popdata to store this sim's entire raster
        print "Loading spike data for file %i..." % i
        allspkdata=loadtxt("%s-spk.txt" % inputstems[i]) # Set spike filename (if used)) # Load the file
        cellscale=ceil(max(allspkdata[:,1])/470) # Calculate how many more cells than normal there are -- 470 is standard, so any more is unusual
        popinfo=[["I2L", "I2", "E2", "I4L", "I4", "E4", "I5L", "I5", "E5R", "E5B", "I6L", "I6", "E6"],array([23, 20, 18, 17, 16, 15, 14, 13, 12, 11, 10, 8, 7])+changedindices,array([13, 25, 150, 14, 20, 30, 13, 25, 65, 17, 13, 25, 60])*cellscale] # List cell population names, numbers, and number of cells in each
        ncells=int(sum(popinfo[2])) # Calculate total number of cells
        npops=size(popinfo,1) # Get the number of cell populations
        maxtime=3 # Maximum time is the time of the latest spike in seconds, rounded up to the nearest integer
        numcols=int(allspkdata[-1,3]+1) # Number of columns is the last number listed in the file, plus one
        nspikes=size(allspkdata,0) # Find out how many spikes across all populations and columns
        tmp2=list() # Re-initialize the spikes-per-population for this file
        for j in range(npops):
            tmp2=allspkdata[allspkdata[:,2]==popinfo[1][j],:] # Pull out the spikes from a particular population
            popdata[i].append(tmp2) # Add one population to popinfo
        # And now for LFPs...
        print "Loading LFP data..."
        lfporig=loadtxt("%s-lfp.txt" % inputstems[i]) # Set spike filename (if used)) # Load the file
        [rows,cols]=shape(lfporig) # Find out how big the array is
        nlayers=int(cols/numcols) # Since layers and columns are interlaced, dividing the total columns by the number of model columns gives layers
        tmp=zeros((rows-2*killtime,nlayers,cols/nlayers)) # Dimensions are: | time | layer | column |
        for j in range(cols):# Reshape the array from 2D to 3D (time, layer, column)
            L=mod(j,nlayers) # Which layer is it? Loop with period nlayers.
            C=int(floor(j/nlayers)) # Which column is it? Increment every nlayers.
            tmp[:,L,C]=lfporig[:,j] # Reshape array, eliminating the first <killtime> points
    # Normalize LFPs
    for i in range(nfiles): lfpmax=max(lfpmax,lfpdata[i][:,whichlfp,0].max())
    for i in range(nfiles): lfpdata[i]/=lfpmax
print "Plotting..."
figh=figure(figsize=(8,11)) # Open figure and store its handle
figh.subplots_adjust(left=0.16) # Less space on left
figh.subplots_adjust(right=0.98) # Less space on right
figh.subplots_adjust(top=0.95) # Less space on top
figh.subplots_adjust(bottom=0.06) # Less space on bottom
figh.subplots_adjust(wspace=0.25) # More space between
figh.subplots_adjust(hspace=0.2) # More space between

# Plot rasters
colors=array([[0.42,0.67,0.84],[0.50,0.80,1.00],[0.90,0.32,0.00],[0.34,0.67,0.67],[0.42,0.82,0.83],[0.90,0.59,0.00],[0.33,0.67,0.47],[0.42,0.83,0.59],[0.90,0.76,0.00],[1.00,0.85,0.00],[0.71,0.82,0.41],[0.57,0.67,0.33],[1.00,0.38,0.60]]) # Colors for each cell population
spc=0 # Only one column, so pick it
for i in range(nfiles):
    for j in range(npops):
        if size(popdata[i][j][:,3]==spc):
            thisdata=popdata[i][j][popdata[i][j][:,3]==spc,:] # Pull out spikes from column spc
            # Reorder spikes so the E2 is at the top of the plot and I6 is at the bottom
            xlim(0,maxtime) # To convert to seconds
            ylim(0,ncells+10) # Just larger than the number of cells in the model
            ylabel("Neuron number")
    thisaxis.set_xticks(r_[0:maxtime+1]) # Don't show half seconds

# Plot LFPs
lfpcolors=[[0,0,1],[1,0,0],[0,0.5,0]] # Use colors from comparecausality
for i in range(nfiles):
    offset=20 # Don't show the first few points
    timeaxis=r_[0:maxtime*fs]/float(fs) # Set time axis
    tmp=plot(timeaxis,lfpdata[i][0+offset:npts+offset,whichlfp,0],linewidth=1,c=lfpcolors[i]) # The middle index -- default 4 -- indicates which LFP to use
    ylim(0.5,1.0) # Cortex: big voltage change
    if i==nfiles-1: xlabel('Time (s)')
    ylabel('Normalized voltage')
    thisaxis.set_xticks(r_[0:maxtime+1]) # Don't show half seconds

print 'Done; elapsed time was %0.1f seconds.' % (toc-tic)