Motor cortex microcircuit simulation based on brain activity mapping (Chadderdon et al. 2014)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:146949
"... We developed a computational model based primarily on a unified set of brain activity mapping studies of mouse M1. The simulation consisted of 775 spiking neurons of 10 cell types with detailed population-to-population connectivity. Static analysis of connectivity with graph-theoretic tools revealed that the corticostriatal population showed strong centrality, suggesting that would provide a network hub. ... By demonstrating the effectiveness of combined static and dynamic analysis, our results show how static brain maps can be related to the results of brain activity mapping."
Reference:
1 . Chadderdon GL, Mohan A, Suter BA, Neymotin SA, Kerr CC, Francis JT, Shepherd GM, Lytton WW (2014) Motor cortex microcircuit simulation based on brain activity mapping. Neural Comput 26:1239-62 [PubMed]
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 V1 L6 pyramidal corticothalamic cell; Neocortex M1 L2/6 pyramidal intratelencephalic cell; Neocortex fast spiking (FS) interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Oscillations; Laminar Connectivity;
Implementer(s): Lytton, William [bill.lytton at downstate.edu]; Neymotin, Sam [samn at neurosim.downstate.edu]; Shepherd, Gordon MG [g-shepherd at northwestern.edu]; Chadderdon, George [gchadder3 at gmail.com]; Kerr, Cliff [cliffk at neurosim.downstate.edu];
Search NeuronDB for information about:  Neocortex V1 L6 pyramidal corticothalamic cell; Neocortex M1 L2/6 pyramidal intratelencephalic cell; GabaA; AMPA; NMDA; Gaba; Glutamate;
/
src
README
infot.mod *
intf6.mod *
intfsw.mod *
matrix.mod
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
boxes.hoc *
col.hoc
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
filtutils.hoc *
gcelldata.hoc
gmgs102.nqs
grvec.hoc *
infot.hoc *
init.hoc
intfsw.hoc *
labels.hoc *
load.py
local.hoc *
main.hoc
misc.h *
miscfuncs.py
network.hoc
neuroplot.py *
nload.hoc
nqs.hoc *
nqsnet.hoc
nrnoc.hoc *
params.hoc
run.hoc
samutils.hoc *
saveoutput.hoc
saveweights.hoc
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
wdmaps2.nqs
xgetargs.hoc *
                            
# load.py -- Script for loading and analyzing data for the Gordon-cell
#   M1 model.
#
# Updated: 11/1/12 (georgec)

#
# Global parameters
#

# Use the NEURON GUI?
use_NEURON_GUI = True

# Set up directory for paper figs.
paperfigsdir = '.'

# Index to grvec panel object for latest loaded sim.
curr_grvec_pind = 1

#
# Misc functions
#

def loadsimbyfstem (filestem,grvec_load=False):
   global num_cells
   global cell_types
   global cell_zlocs
   global curr_grvec_pind

   # Get LFPs NQS (nqLFP).
   h('if (nqLFP != nil) nqsdel(nqLFP)')
   h('nqLFP = new NQS()')
   fname = '%s-lfp.nqs' % filestem
   h.nqLFP.rd(fname)

   # Get spikes NQS (snq).
   h('if (snq != nil) nqsdel(snq)')
   h('snq = new NQS()')
   fname = '%s-spk.nqs' % filestem
   h.snq.rd(fname)
   spks_tvec = nqscol2narr('snq','t')
   spks_idvec = nqscol2narr('snq','id')

   # Get cells info NQS (cellsnq).
   h('if (cellsnq != nil) nqsdel(cellsnq)')
   h('cellsnq = new NQS()')
   fname = '%s-loc.nqs' % filestem
   h.cellsnq.rd(fname)

   # Get connections info NQS (connsnq).
   h('if (connsnq != nil) nqsdel(connsnq)')
   h('connsnq = new NQS()')
   fname = '%s-con.nqs' % filestem
   h.connsnq.rd(fname)

   # Set up second connsnq variable for info extraction.
   h('if (connsnq2 != nil) nqsdel(connsnq2)')
   h('connsnq2 = new NQS()')

   # Load information from the cellsnq NQS table.
   num_cells = int(h.cellsnq.size())
   cell_types = nqscol2narr('cellsnq','ty')
   cell_zlocs = nqscol2narr('cellsnq','zloc')

   # If the NEURON GUI is up, load the grvec printlist from the saved data.
   if use_NEURON_GUI and grvec_load:
      ldgrvec('%s-prlist.gvc' % outfilestem)
      curr_grvec_pind += 1

def ctypandind2cind (ctyp,ctind):
   cell_inds = find(cell_types == get_ctyp_num(ctyp))
   return cell_inds[ctind]

def cind2ctypandind (cell_num):
   # Get the cell type number.
   ctypnum = cell_types[cell_num]

   # Get the first cell index for that type.
   ctypfind = find(cell_types == ctypnum)[0]

   # Show the cell type string and the relative cell index.
   return (get_ctyp_str(ctypnum), cell_num - ctypfind)

#
# Data analysis functions
#

def get_ctyp_fire_rate (ctyp,min_t=0.0,max_t=-1.0,allcellrates=False):
   if (max_t == -1.0):
      max_t = sim_duration
   tvec = nqscol2narr('snq','t')
   idvec = nqscol2narr('snq','id')
   tvec,idvec = get_vec_subset(tvec,idvec,mintvec=min_t,maxtvec=max_t)
   freqs = get_ave_fire_freqs(tvec,idvec,num_cells,max_t - min_t)
   cell_inds = find(cell_types == get_ctyp_num(ctyp))
   if (allcellrates):
      return freqs[cell_inds]
   else:
      return freqs[cell_inds].mean()

# get_curr_fire_rate() -- get an estimate of a firing rate for a target cell at a 
#    particular time.  By default, the counting window is 100 ms and starts at time t.
#    Requires that snq table is loaded.
def get_curr_fire_rate (targtype='',targind=0,t=0.0,spkwindwid=100.0,spkwindoffset=0.0):
   # Remember the old snq table verbosity.
   oldverbose = h.snq.verbose

   # Turn off the table verbosity.
   h.snq.verbose = 0

   # Default the cell index to targind.
   cell_ind = targind

   # If we have a cell type, get the cell IDs and get the indexed value.
   if (targtype != ''):
      targ_cells = find(cell_types == get_ctyp_num(targtype)) 
      cell_ind = targ_cells[targind]     
   
   # Get number all of the desired cell spikes in the desired time window.
   num_spks = h.snq.select('id',cell_ind,'t','[]',t+spkwindoffset, \
      t+spkwindoffset+spkwindwid)

   # Get the firing rate from the spike count and window width.
   fire_rate = num_spks * (1000.0 / spkwindwid)

   # Reset the snq table.
   h.snq.tog()
   h.snq.verbose = oldverbose

   return fire_rate

def get_avg_ctyp_rates (min_t=0.0,max_t=-1.0):
   ctyps = ['E2','I2','I2L','E5R','E5B','I5','I5L','E6','I6','I6L']
   rates = zeros(len(ctyps))
   for rind in range(len(ctyps)):
      rates[rind] = get_ctyp_fire_rate(ctyps[rind],min_t,max_t)
   return ctyps,rates

def get_ctyp_spk_count (ctyp,min_t=0.0,max_t=-1.0):
   if (max_t == -1.0):
      max_t = sim_duration

   # Get all of the spikes.
   tvec = nqscol2narr('snq','t')
   idvec = nqscol2narr('snq','id')

   # Get all of the cell indices for the desire cell type.
   cell_inds = find(cell_types == get_ctyp_num(ctyp))

   # Get only the spikes of the right cell type and in the time window.
   tvec,idvec = get_vec_subset(tvec,idvec,mintvec=min_t,maxtvec=max_t, \
      minvec=cell_inds[0],maxvec=cell_inds[-1])

   # Return the numbrer of spikes.
   return len(tvec)

def get_bin_integral (F,xs,xstart,xend): 
   # Get the relevant bin indices using xs, xstart, and xend.
   bin_inds = find((xs >= xstart) & (xs <= xend))

   # Start with the sum being the sum all of the bins but the first and last.
   sumval = F[bin_inds[1:-1]].sum()

   # Add the portion of the left bin we have.
   sumval += (F[bin_inds[0]] * (xs[bin_inds[0]] - xstart) / (xs[1] - xs[0]))

   # Add the portion of the right bin we have.
   sumval += (F[bin_inds[-1]] * (xend - xs[bin_inds[-1]]) / (xs[1] - xs[0]))

   return sumval

# get_conv_conns() -- load the (static) convergence connections into a unit into 
#    the NQS table connsnq2.  Requires the connection table connsnq to be loaded.
def get_conv_conns (srctype='E2',targtype='E5R',targind=0):
   # Remember the old connsnq table verbosity, and turn it off.
   oldverbose = h.connsnq.verbose
   h.connsnq.verbose = 0

   # Remember the old connsnq2 table verbosity, and turn it off.
   oldverbose2 = h.connsnq2.verbose
   h.connsnq2.verbose = 0

   # Get source cell IDs and save in v1.
   src_cells = find(cell_types == get_ctyp_num(srctype))
   narr2hv('v1',src_cells)

   # Get target cell IDs.
   targ_cells = find(cell_types == get_ctyp_num(targtype))

   # Select the entries with the desired weights and return the number of connections found.
   nc = h.connsnq.select("id1","EQW",h.v1,"id2",targ_cells[targind])
   
   # If we found connections...
   if (nc > 0):
      # Copy those weights over to connsnq2.
      h.connsnq2.cp(h.connsnq.out)

      # Sort the list.
      h.connsnq2.sort("id1")

   # Reset the connsnq table.
   h.connsnq.tog("db")
   h.connsnq.verbose = oldverbose

   # Reset the connsnq2 table.
   h.connsnq2.verbose = oldverbose2

   # Return the number of connections found.
   return nc

# get_div_conns() -- load the (static) divergence connections out of a unit into 
#    the NQS table connsnq2.  Requires the connection table connsnq to be loaded.
def get_div_conns (srctype='E2',targtype='E5R',srcind=0):
   # Remember the old connsnq table verbosity, and turn it off.
   oldverbose = h.connsnq.verbose
   h.connsnq.verbose = 0

   # Remember the old connsnq2 table verbosity, and turn it off.
   oldverbose2 = h.connsnq2.verbose
   h.connsnq2.verbose = 0

   # Get source cell IDs.
   src_cells = find(cell_types == get_ctyp_num(srctype))

   # Get target cell IDs and save in v1.
   targ_cells = find(cell_types == get_ctyp_num(targtype))
   narr2hv('v1',targ_cells)

   # Select the entries with the desired weights and return the number of connections found.
   nc = h.connsnq.select("id1",src_cells[srcind],"id2","EQW",h.v1)

   # If we found connections...
   if (nc > 0):   
      # Copy those weights over to connsnq2.
      h.connsnq2.cp(h.connsnq.out)

      # Sort the list.
      h.connsnq2.sort("id2")

   # Reset the connsnq table.
   h.connsnq.tog("db")
   h.connsnq.verbose = oldverbose

   # Reset the connsnq2 table.
   h.connsnq2.verbose = oldverbose2

   # Return the number of connections found.
   return nc

# get_conv_wts() -- for the convergent weights into a unit (and potentially at a 
#    particular time) return an array of the source cell IDs and the weights (static 
#    or plastic).
def get_conv_wts (srctype='E2',targtype='E5R',targind=0):
   # Get source cell IDs.
   src_cells = find(cell_types == get_ctyp_num(srctype))

   # Get target cell IDs.
   targ_cells = find(cell_types == get_ctyp_num(targtype))

   # Get the appropriate connections into connsnq2 table.
   nc = get_conv_conns(srctype,targtype,targind)

   # If we actually have some connections...
   if (nc > 0):
      # Get IDs and subtract first index for type.
      nqids = nqscol2narr('connsnq2','id1')
      nqids -= src_cells[0]
      nqids = array(nqids,int)

      # Get the primary weights from the table.
      nqwts = nqscol2narr('connsnq2','wt1')
   else:
      nqids = []
      nqwts = []

   return [nqids, nqwts]

def get_avg_conv_wts (srctype='E2',targtype='E5R'):
   # Get target cell IDs.
   targ_cells = find(cell_types == get_ctyp_num(targtype))

   # Set up an array for the weights.
   avg_conv_wts = zeros(len(targ_cells))

   # For all of the target cells...
   for ii in range(len(targ_cells)):
      # Read the convergent weight IDs and weight values.
      nqids, nqwts = get_conv_wts(srctype, targtype, ii)

      # Save the average weight from this.
      if (len(nqwts) > 0):
         avg_conv_wts[ii] = nqwts.mean()
      else:
         avg_conv_wts[ii] = 0

   # Return the average of averages.
   return avg_conv_wts.mean()

#
# Display functions
#

def pravgrates (min_t=0.0,max_t=-1.0):
   print 'Average Firing Rates by Cell Type:'
   print '   E2: %.2f Hz' % get_ctyp_fire_rate('E2',min_t,max_t)
   print '   I2: %.2f Hz' % get_ctyp_fire_rate('I2',min_t,max_t)
   print '  I2L: %.2f Hz' % get_ctyp_fire_rate('I2L',min_t,max_t)
   print '  E5R: %.2f Hz' % get_ctyp_fire_rate('E5R',min_t,max_t)
   print '  E5B: %.2f Hz' % get_ctyp_fire_rate('E5B',min_t,max_t)
   print '   I5: %.2f Hz' % get_ctyp_fire_rate('I5',min_t,max_t)
   print '  I5L: %.2f Hz' % get_ctyp_fire_rate('I5L',min_t,max_t)
   print '   E6: %.2f Hz' % get_ctyp_fire_rate('E6',min_t,max_t)
   print '   I6: %.2f Hz' % get_ctyp_fire_rate('I6',min_t,max_t)
   print '  I6L: %.2f Hz' % get_ctyp_fire_rate('I6L',min_t,max_t)

#
# Plotting functions
#

def plot_raster (startt=0,endt=1000,xmarg=10,ymarg=10,drawpopbounds=True,bwmode=False):
   def plotcells (ctyp,pcol,drawpopbounds=True,drawpoplabels=True):
      # Select only the cells in the right time window and of the right type.
      spk_count = h.snq.select("t","[]",startt,endt,"type",get_ctyp_num(ctyp))

      # If there are spikes...
      if (spk_count):
         # Get the spikes and show them.
         ts = nqscol2narr('snq','t')
         ts /= 1000.0
         cs = nqscol2narr('snq','id')
         marks = scatter(ts, cs, marker='+')
         setp(marks, color=pcol)

      # Reset the snq table.
      h.snq.tog("DB")

      # If we want to draw population boundaries...
      if (drawpopbounds):
         cell_inds = find(cell_types == get_ctyp_num(ctyp))
         axhline(cell_inds[0] - 0.5,color='k')
         axhline(cell_inds[-1] + 0.5,color='k') 

      # If we want to draw population labels...
      if (drawpoplabels):
         cell_inds = find(cell_types == get_ctyp_num(ctyp))
         # If E5B, convert cell type to 'SPI'.
         if (ctyp == 'E5B'):
            ctyp = 'SPI'
         elif (ctyp == 'E5R'):
            ctyp = 'STR'
         text((xmarg + endt) / 1000.0 + ((endt - startt) / 1000.0 / 80.0), (cell_inds[0] + cell_inds[-1]) / 2.0 - 12, ctyp)

   # If not black-and-white mode...
   if (not bwmode):
      # Plot the cells with color-coding.
#      plotcells('E2','b',drawpopbounds)
#      plotcells('I2','r',drawpopbounds)  
#      plotcells('I2L','orange',drawpopbounds)
##      plotcells('E4','b',drawpopbounds) 
##      plotcells('I4','r',drawpopbounds)  
##      plotcells('I4L','orange',drawpopbounds)   
#      plotcells('E5R','b',drawpopbounds) 
#      plotcells('E5B','g',drawpopbounds)    
#      plotcells('I5','r',drawpopbounds)  
#      plotcells('I5L','orange',drawpopbounds)
#      plotcells('E6','b',drawpopbounds)  
#      plotcells('I6','r',drawpopbounds) 
#      plotcells('I6L','orange',drawpopbounds)

      plotcells('E2','b',drawpopbounds)
      plotcells('I2','orange',drawpopbounds)  
      plotcells('I2L','m',drawpopbounds)
#      plotcells('E4','b',drawpopbounds) 
#      plotcells('I4','orange',drawpopbounds)  
#      plotcells('I4L','m',drawpopbounds)   
      plotcells('E5R','g',drawpopbounds) 
      plotcells('E5B','r',drawpopbounds)    
      plotcells('I5','orange',drawpopbounds)  
      plotcells('I5L','m',drawpopbounds)
      plotcells('E6','c',drawpopbounds)  
      plotcells('I6','orange',drawpopbounds) 
      plotcells('I6L','m',drawpopbounds)
   else:
      # Plot the cells with black.
      plotcells('E2','k',drawpopbounds)
      plotcells('I2','k',drawpopbounds)  
      plotcells('I2L','k',drawpopbounds)
#      plotcells('E4','k',drawpopbounds) 
#      plotcells('I4','k',drawpopbounds)  
#      plotcells('I4L','k',drawpopbounds)   
      plotcells('E5R','k',drawpopbounds) 
      plotcells('E5B','k',drawpopbounds)    
      plotcells('I5','k',drawpopbounds)  
      plotcells('I5L','k',drawpopbounds)
      plotcells('E6','k',drawpopbounds)  
      plotcells('I6','k',drawpopbounds) 
      plotcells('I6L','k',drawpopbounds)

   # Set the x limits according to the min and max time.
   xlim((startt - xmarg) / 1000.0, (endt + xmarg) / 1000.0)

   # Set the y limits according to which cells exist.
   ylim(-ymarg, num_cells + ymarg - 1)

   # Label the axes.
   xlabel('Time (s)')
   ylabel('Cell #')

def plot_cell_densities ():
   plot(lyrbins,M1dens)
   title('Cell Density vs. Depth')
   xlabel('Distance from pia (microns)')
   ylabel('Cell Density (E and I)')
   axvline(143.0,color='k') # L1 / L2/3 boundary
   axvline(451.8,color='k') # L2/3 / L5A boundary
   axvline(663.6,color='k') # L5A / L5B boundary
   axvline(1059.0,color='k') # L5B / L6 boundary
   text(70,0.11,'L1')
   text(300,0.11,'L2/3')
   text(525,0.11,'L5A')
   text(850,0.11,'L5B')
   text(1250,0.11,'L6')

def plot_vtrace (tvec,vec,ctyp='E2'):
   plot(tvec,vec)
   ylim(-80,60)
   if (ctyp in ('E2','E5R','E6')):
      axhline(-25,color='r')
      axhline(-40,color='k')
      axhline(-65,color='k',linestyle='--')
   elif (ctyp == 'E5B'):
      axhline(-25,color='r')
      axhline(-37.5,color='k')
      axhline(-65,color='k',linestyle='--')
   elif (ctyp in ('I2','I5','I6')):
      axhline(-40,color='k')
      axhline(-63,color='k',linestyle='--')
   elif (ctyp in ('I2L','I5L','I6L')):
      axhline(-47,color='k')
      axhline(-65,color='k',linestyle='--')
   title('%s Voltage Trace' % ctyp)
   xlabel('Time (ms)')
   ylabel('Memb. Pot. (mV)')

def plot_inter_ctyp_matrix (ctyp_names, vals):
   # Make the reverse list of names.
   rnames=[]
   rnames.extend(ctyp_names)
   rnames.reverse() 

   npops = size(ctyp_names)
   figh = figure(figsize=(9,8)) # Open figure and store its handle
   axh = subplot(111)
   figh.subplots_adjust(left=0.1) # Less space on right
   figh.subplots_adjust(right=1.1) # Less space on right
   ex = array(range(npops+1))
   why = npops - ex
   tmp = axh.pcolor(ex,why,vals.transpose(),linestyle='solid',linewidth=0.5,color='k')
   xlim([0,npops])
   ylim([0,npops])
   axh.xaxis.set_ticks(array(range(npops))+0.5)
   axh.xaxis.set_ticklabels(ctyp_names)
   axh.xaxis.set_ticks_position('top')
   axh.yaxis.set_ticks(array(range(npops))+0.5)
   axh.yaxis.set_ticklabels(rnames)
   tmp.cmap = cm.bone_r # Blues_r # summer_r#cm.gist_earth # The _r reverses the order -- beautiful solution!
   colorbar(tmp, fraction=0.2)
   xlabel('Pre-Synaptic Cell Type')
   ylabel('Post-Synaptic Cell Type')

def plot_inter_ctyp_pmat ():
   def cellpops ():
      I2L=0; I2=1;  E2=2; I5L=3; I5=4; E5R=5; E5B=6; I6L=7 ;I6=8; E6=9;
      return I2L, I2, E2, I5L, I5, E5R, E5B, I6L, I6, E6

   # Define connectivity matrix (copied out of network.py)
   def definepmat (npops):
      I2L, I2, E2, I5L, I5, E5R, E5B, I6L, I6, E6 = cellpops()
      pmat = zeros((npops,npops))

      pmat[E2][E2]=0.2     # weak by wiring matrix in (Weiler et al., 2008)
      pmat[E2][E5B]=0.8    # strong by wiring matrix in (Weiler et al., 2008)
      pmat[E2][E5R]=0.8    # strong by wiring matrix in (Weiler et al., 2008)
      pmat[E2][I5L]=0.51   # L2/3 E -> L5 LTS I (justified by (Apicella et al., 2012)
      pmat[E2][E6]=0       # none by wiring matrix in (Weiler et al., 2008)
      pmat[E2][I2L]=0.51
      pmat[E2][I2]=0.43

      pmat[E5B][E2]=0          # none by wiring matrix in (Weiler et al., 2008)
      pmat[E5B][E5B]=0.04 * 4  # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 
      pmat[E5B][E5R]=0         # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 
      pmat[E5B][E6]=0          # none by suggestion of Ben and Gordon over phone
      pmat[E5B][I5L]=0         # ruled out by (Apicella et al., 2012) Fig. 7
      pmat[E5B][I5]=0.43 

      pmat[E5R][E2]=0.2        # weak by wiring matrix in (Weiler et al., 2008)
      pmat[E5R][E5B]=0.21 * 4  # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4
      pmat[E5R][E5R]=0.11 * 4  # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 
      pmat[E5R][E6]=0.2        # weak by wiring matrix in (Weiler et al., 2008)
      pmat[E5R][I5L]=0         # ruled out by (Apicella et al., 2012) Fig. 7
      pmat[E5R][I5]=0.43

      pmat[E6][E2]=0           # none by wiring matrix in (Weiler et al., 2008)
      pmat[E6][E5B]=0.2        # weak by wiring matrix in (Weiler et al., 2008)
      pmat[E6][E5R]=0.2        # weak by wiring matrix in (Weiler et al., 2008)
      pmat[E6][E6]=0.2         # weak by wiring matrix in (Weiler et al., 2008)
      pmat[E6][I6L]=0.51
      pmat[E6][I6]=0.43

      pmat[I2L][E2]=0.35
      pmat[I2L][I2L]=0.09
      pmat[I2L][I2]=0.53
      pmat[I2][E2]=0.44
      pmat[I2][I2L]=0.34
      pmat[I2][I2]=0.62
      pmat[I5L][E5B]=0.35
      pmat[I5L][E5R]=0.35
      pmat[I5L][I5L]=0.09
      pmat[I5L][I5]=0.53
      pmat[I5][E5B]=0.44
      pmat[I5][E5R]=0.44
      pmat[I5][I5L]=0.34
      pmat[I5][I5]=0.62
      pmat[I6L][E6]=0.35
      pmat[I6L][I6L]=0.09
      pmat[I6L][I6]=0.53
      pmat[I6][E6]=0.44
      pmat[I6][I6L]=0.34
      pmat[I6][I6]=0.62

      return pmat

   ctyp_strs = ["I2L", "I2", "E2", "I5L", "I5", "STR", "SPI", "I6L", "I6", "E6"]
   pmat = definepmat(len(ctyp_strs))
   plot_inter_ctyp_matrix(ctyp_strs, pmat)

def plot_inter_ctyp_avg_wts ():
   def cellpops ():
      I2L=0; I2=1;  E2=2; I5L=3; I5=4; E5R=5; E5B=6; I6L=7 ;I6=8; E6=9;
      return I2L, I2, E2, I5L, I5, E5R, E5B, I6L, I6, E6

   I2L, I2, E2, I5L, I5, E5R, E5B, I6L, I6, E6 = cellpops()
   ctyp_labels = ["I2L", "I2", "E2", "I5L", "I5", "STR", "SPI", "I6L", "I6", "E6"]
   ctyp_strs = ["I2L", "I2", "E2", "I5L", "I5", "E5R", "E5B", "I6L", "I6", "E6"]
   wmat = zeros((10,10))

   # still under construction...
   for ii in range(10):
      for jj in range(10):
         wmat[ii][jj] = get_avg_conv_wts(ctyp_strs[ii], ctyp_strs[jj])
   
   plot_inter_ctyp_matrix(ctyp_labels, wmat)
   
#
# Batch functions
#

def batchA (simdatadir='data/12oct11_sims1',savebatch=True):
   # Set up ctyp range to use (useful to make smaller sometimes due to memory problems).
   ctyp_list = ('E2','I2','I2L','E5R','E5B','I5','I5L','E6','I6','I6L')
#   ctyp_list = ('E2','I2','I2L')
#   ctyp_list = ('E5R','E5B','I5','I5L')
#   ctyp_list = ('E6','I6','I6L')
#   ctyp_list = ('E5R','E5B','I5','I5L','E6','I6','I6L')

   # Set up upper bound for the firing response rate calculations.
   short_resp_dur = 10   # ms
   long_resp_dur = 1000  # ms
   short_resp_end = 1002 + short_resp_dur  # start is always at 1002 ms so excluding shock
   long_resp_end = 1002 + long_resp_dur
   if (short_resp_end > 2000):
      short_resp_end = 2000
   if (long_resp_end > 2000):
      long_resp_end = 2000

   for ctyp in ctyp_list:
      exec('global shock_%s_spk_counts' % ctyp)
      exec('global shock_%s_resp_rates' % ctyp)
      exec('global shock_%s_resp_rates2' % ctyp)

   # For each of the shock percentages...
   for ii in range(len(shock_pcts)):
      # Load L2/3 shock data.
      outfilestem = '%s/L2shock%dp' % (simdatadir, shock_pcts[ii])
#      outfilestem = '%s/L2isolshock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,0] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,0] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,short_resp_end))
         exec("shock_%s_resp_rates2[ii,0] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,long_resp_end))

      # Load L5A shock data.
      outfilestem = '%s/L5Ashock%dp' % (simdatadir, shock_pcts[ii])
#      outfilestem = '%s/L5Aisolshock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,1] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,1] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,short_resp_end))
         exec("shock_%s_resp_rates2[ii,1] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,long_resp_end))

      # Load L5B shock data.
      outfilestem = '%s/L5Bshock%dp' % (simdatadir, shock_pcts[ii])
#      outfilestem = '%s/L5Bisolshock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,2] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,2] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,short_resp_end))
         exec("shock_%s_resp_rates2[ii,2] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,long_resp_end))

      # Load L6 shock data.
      outfilestem = '%s/L6shock%dp' % (simdatadir, shock_pcts[ii])
#      outfilestem = '%s/L6isolshock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,3] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,3] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,short_resp_end))
         exec("shock_%s_resp_rates2[ii,3] = get_ctyp_fire_rate('%s',1002,%d)" % \
            (ctyp,ctyp,long_resp_end))

      # Save the data if we want to.
      if (savebatch):
         dirname = '%s/batch_results' % simdatadir

         # If the directory doesn't exist, make it.
         if (not os.path.exists(dirname)):
            os.system('mkdir %s' % dirname)

         # Save all of the numpy arrays.
         for ctyp in ctyp_list:
            exec("fname = '%s/%s-spk_counts.txt'" % (dirname, ctyp))
            exec("savetxt('%s', shock_%s_spk_counts)" % (fname, ctyp))
            exec("fname = '%s/%s-resp_rates.txt'" % (dirname, ctyp))
            exec("savetxt('%s', shock_%s_resp_rates)" % (fname, ctyp))
            exec("fname = '%s/%s-resp_rates2.txt'" % (dirname, ctyp))
            exec("savetxt('%s', shock_%s_resp_rates2)" % (fname, ctyp))            

def loadbatchA (simdatadir='data/12oct11_sims1'):
   # This block doesn't seem to work for this function.
#   for ctyp in ('E2','I2','I2L','E5R','E5B','I5','I5L','E6','I6','I6L'):
#      exec('global shock_%s_spk_counts' % ctyp)
#      exec('global shock_%s_resp_rates' % ctyp)
#      exec('global shock_%s_resp_rates2' % ctyp)
   global shock_E2_spk_counts
   global shock_E2_resp_rates
   global shock_E2_resp_rates2
   global shock_I2_spk_counts
   global shock_I2_resp_rates
   global shock_I2_resp_rates2
   global shock_I2L_spk_counts
   global shock_I2L_resp_rates
   global shock_I2L_resp_rates2
   global shock_E5R_spk_counts
   global shock_E5R_resp_rates
   global shock_E5R_resp_rates2
   global shock_E5B_spk_counts
   global shock_E5B_resp_rates
   global shock_E5B_resp_rates2
   global shock_I5_spk_counts
   global shock_I5_resp_rates
   global shock_I5_resp_rates2
   global shock_I5L_spk_counts
   global shock_I5L_resp_rates
   global shock_I5L_resp_rates2
   global shock_E6_spk_counts
   global shock_E6_resp_rates
   global shock_E6_resp_rates2
   global shock_I6_spk_counts
   global shock_I6_resp_rates
   global shock_I6_resp_rates2
   global shock_I6L_spk_counts
   global shock_I6L_resp_rates
   global shock_I6L_resp_rates2

   dirname = '%s/batch_results' % simdatadir

   # Load all of the numpy arrays.
   # This is not working in this function for some reason.  Seems like the for loop
   # is its own scope and the information loaded in it is lost.
#   for ctyp in ('E2','I2','I2L','E5R','E5B','I5','I5L','E6','I6','I6L'):
#      exec("fname = '%s/%s-spk_counts.txt'" % (dirname, ctyp))
#      exec("shock_%s_spk_counts = loadtxt('%s')" % (ctyp, fname))
#      exec("fname = '%s/%s-resp_rates.txt'" % (dirname, ctyp))
#      exec("shock_%s_resp_rates = loadtxt('%s')" % (ctyp, fname))
#      exec("fname = '%s/%s-resp_rates2.txt'" % (dirname, ctyp))
#      exec("shock_%s_resp_rates2 = loadtxt('%s')" % (ctyp, fname))  
   fname = '%s/E2-spk_counts.txt' % dirname
   shock_E2_spk_counts = loadtxt(fname)
   fname = '%s/E2-resp_rates.txt' % dirname
   shock_E2_resp_rates = loadtxt(fname)
   fname = '%s/E2-resp_rates2.txt' % dirname
   shock_E2_resp_rates2 = loadtxt(fname)

   fname = '%s/I2-spk_counts.txt' % dirname
   shock_I2_spk_counts = loadtxt(fname)
   fname = '%s/I2-resp_rates.txt' % dirname
   shock_I2_resp_rates = loadtxt(fname)
   fname = '%s/I2-resp_rates2.txt' % dirname
   shock_I2_resp_rates2 = loadtxt(fname)

   fname = '%s/I2L-spk_counts.txt' % dirname
   shock_I2L_spk_counts = loadtxt(fname)
   fname = '%s/I2L-resp_rates.txt' % dirname
   shock_I2L_resp_rates = loadtxt(fname)
   fname = '%s/I2L-resp_rates2.txt' % dirname
   shock_I2L_resp_rates2 = loadtxt(fname)

   fname = '%s/E5R-spk_counts.txt' % dirname
   shock_E5R_spk_counts = loadtxt(fname)
   fname = '%s/E5R-resp_rates.txt' % dirname
   shock_E5R_resp_rates = loadtxt(fname)
   fname = '%s/E5R-resp_rates2.txt' % dirname
   shock_E5R_resp_rates2 = loadtxt(fname)

   fname = '%s/E5B-spk_counts.txt' % dirname
   shock_E5B_spk_counts = loadtxt(fname)
   fname = '%s/E5B-resp_rates.txt' % dirname
   shock_E5B_resp_rates = loadtxt(fname)
   fname = '%s/E5B-resp_rates2.txt' % dirname
   shock_E5B_resp_rates2 = loadtxt(fname)

   fname = '%s/I5-spk_counts.txt' % dirname
   shock_I5_spk_counts = loadtxt(fname)
   fname = '%s/I5-resp_rates.txt' % dirname
   shock_I5_resp_rates = loadtxt(fname)
   fname = '%s/I5-resp_rates2.txt' % dirname
   shock_I5_resp_rates2 = loadtxt(fname)

   fname = '%s/I5L-spk_counts.txt' % dirname
   shock_I5L_spk_counts = loadtxt(fname)
   fname = '%s/I5L-resp_rates.txt' % dirname
   shock_I5L_resp_rates = loadtxt(fname)
   fname = '%s/I5L-resp_rates2.txt' % dirname
   shock_I5L_resp_rates2 = loadtxt(fname)

   fname = '%s/E6-spk_counts.txt' % dirname
   shock_E6_spk_counts = loadtxt(fname)
   fname = '%s/E6-resp_rates.txt' % dirname
   shock_E6_resp_rates = loadtxt(fname)
   fname = '%s/E6-resp_rates2.txt' % dirname
   shock_E6_resp_rates2 = loadtxt(fname)

   fname = '%s/I6-spk_counts.txt' % dirname
   shock_I6_spk_counts = loadtxt(fname)
   fname = '%s/I6-resp_rates.txt' % dirname
   shock_I6_resp_rates = loadtxt(fname)
   fname = '%s/I6-resp_rates2.txt' % dirname
   shock_I6_resp_rates2 = loadtxt(fname)

   fname = '%s/I6L-spk_counts.txt' % dirname
   shock_I6L_spk_counts = loadtxt(fname)
   fname = '%s/I6L-resp_rates.txt' % dirname
   shock_I6L_resp_rates = loadtxt(fname)
   fname = '%s/I6L-resp_rates2.txt' % dirname
   shock_I6L_resp_rates2 = loadtxt(fname)

def batchB (simdatadir='data/12oct10_sims2',savebatch=True):
   # Set up ctyp range to use (useful to make smaller sometimes due to memory problems).
   ctyp_list = ('E2','I2','I2L','E5R','E5B','I5','I5L','E6','I6','I6L')

   for ctyp in ctyp_list:
      exec('global shock_%s_spk_counts' % ctyp)
      exec('global shock_%s_resp_rates' % ctyp)
      exec('global shock_%s_resp_rates2' % ctyp)
      exec('global shock_%s_resp_rates_mean' % ctyp)
      exec('global shock_%s_resp_rates_std' % ctyp)
      exec('global shock_%s_resp_rates2_mean' % ctyp)
      exec('global shock_%s_resp_rates2_std' % ctyp)

   # For each of the shock percentages...
   for ii in range(len(shock_pcts)):
      # Load L2/3 shock data.
      outfilestem = '%s/L2shock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,0] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,0] = get_ctyp_fire_rate('%s',1002,1012)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates2[ii,0] = get_ctyp_fire_rate('%s',1002,2000)" % \
            (ctyp,ctyp))

      # Load L5A shock data.
      outfilestem = '%s/L5Ashock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,1] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,1] = get_ctyp_fire_rate('%s',1002,1012)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates2[ii,1] = get_ctyp_fire_rate('%s',1002,2000)" % \
            (ctyp,ctyp))

      # Load L5B shock data.
      outfilestem = '%s/L5Bshock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,2] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,2] = get_ctyp_fire_rate('%s',1002,1012)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates2[ii,2] = get_ctyp_fire_rate('%s',1002,2000)" % \
            (ctyp,ctyp))

      # Load L6 shock data.
      outfilestem = '%s/L6shock%dp' % (simdatadir, shock_pcts[ii])
      loadsimbyfstem(outfilestem)
      for ctyp in ctyp_list:
         exec("shock_%s_spk_counts[ii,3] = get_ctyp_spk_count('%s',1000,1001)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates[ii,3] = get_ctyp_fire_rate('%s',1002,1012)" % \
            (ctyp,ctyp))
         exec("shock_%s_resp_rates2[ii,3] = get_ctyp_fire_rate('%s',1002,2000)" % \
            (ctyp,ctyp))

def rasterdemo ():
   plot_raster(980,1060,xmarg=0,ymarg=0)
   xticks((0.96,0.98,1.0,1.02,1.04,1.06,1.08),('','-20','0','20','40','60',''))
   xlim(0.98,1.06)
   xlabel('Peri-stimulus Time (ms)')
   ylabel('')
   yticks([])

#
# Script code (run always)
#

# Load the sys and os packages.
import sys, os

# Set up the file stem pointing to the simulation output files.
# NOTE: When "ipython -pylab load.py" is called, only 1 arg is loaded, the last.
if (len(sys.argv) > 1):
   outfilestem = sys.argv[1]
else:
   outfilestem = './runsim'

# Load the interpreter object.
from neuron import h

# Load the NEURON GUI.
if use_NEURON_GUI:
   from neuron import gui

# Load the numpy namespace.
from numpy import *

# Load the neuroplot namespace.
from neuroplot import *

# Set up neurosim run-neuron code.
h.xopen("setup.hoc")
h.xopen("nrnoc.hoc")
h.load_file("init.hoc")

# Load functions for expediting NQS processing in Python interpreter.
execfile("miscfuncs.py")

# If the NEURON GUI is up, load the E5R and E5B cell data from the saved grvec file.
#if use_NEURON_GUI:
#   ldgrvec('celldata.gvc')

# Set up hoc objects for NQS tables for loadsimbyfstem()
h('objref nqLFP')
h('objref snq')
h('objref cellsnq')
h('objref connsnq')
h('objref connsnq2')

# Load the first simulation.
loadsimbyfstem(outfilestem,grvec_load=False)

#
# Set layer information.
#

# Set the M1 cell density array.
M1dens = [0.0044,0.0297,0.0253,0.0121,0.0141,0.0106,0.0195,0.0372,0.0361,0.0591,0.0748,0.0835,0.0947,0.0793,0.0923,0.0796,0.0702,0.0809,0.0749,0.0827,0.0792,0.0893,0.0475,0.0811,0.0755,0.0754,0.083,0.0688,0.0892,0.0765,0.1058,0.0791,0.091,0.1018,0.0931,0.0847,0.0919,0.0661,0.0896,0.0686,0.073,0.0774,0.0588,0.0661,0.0707,0.0496,0.0714,0.065,0.0463,0.0766,0.0546,0.0619,0.0584,0.0541,0.0612,0.047,0.069,0.0611,0.0519,0.0601,0.0646,0.0459,0.0634,0.0348,0.0705,0.0546,0.0485,0.0434,0.0659,0.0585,0.0655,0.0638,0.0922,0.0753,0.0833,0.0798,0.0751,0.0831,0.0893,0.1004,0.0951,0.0734,0.0845,0.0918,0.093,0.0867,0.0828,0.0971,0.088,0.0822,0.0925,0.0884,0.0944,0.0892,0.0659,0.0875,0.1046,0.123,0.1187,0.0656]
M1dens = array(M1dens)

# Set up the layer left bin boundary array.
lyrbins = arange(0,1,0.01)
lyrbins *= 1412.0

# Remember the simulation duration (in ms).
sim_duration = 2000.0

#
# Batch trend collection
#

# data for usual (old batch)
#shock_pcts = [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,30,35,\
#   40,45,50,51,52,53,54,55]
# data for new batch
shock_pcts = [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,30,35,\
   40,45,50,51,52,53,54,55,60,65,70,75,80,85,90,95,100]
# data for isolated layer models batch
#shock_pcts = [5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100]
inseeds = [1,2,3,4,5]
wseeds = [1,2,3,4,5]

# Set up the counts for shock spikes and responses.
for ctyp in ('E2','I2','I2L','E5R','E5B','I5','I5L','E6','I6','I6L'):
   # Batch A (single-seed version)
   exec('shock_%s_spk_counts = zeros((%d,4))' % (ctyp, len(shock_pcts)))
   exec('shock_%s_resp_rates = zeros((%d,4))' % (ctyp, len(shock_pcts)))
   exec('shock_%s_resp_rates2 = zeros((%d,4))' % (ctyp, len(shock_pcts)))

   # Batch B (multi-seed version)
#   exec('shock_%s_spk_counts = zeros((%d,4,%d,%d))' % \
#      (ctyp, len(shock_pcts), len(wseeds), len(inseeds)))
#   exec('shock_%s_resp_rates = zeros((%d,4,%d,%d))' % \
#      (ctyp, len(shock_pcts), len(wseeds), len(inseeds)))
#   exec('shock_%s_resp_rates2 = zeros((%d,4,%d,%d))' % \
#      (ctyp, len(shock_pcts), len(wseeds), len(inseeds)))
#   exec('shock_%s_resp_rates_mean = zeros((%d,4))' % (ctyp, len(shock_pcts)))
#   exec('shock_%s_resp_rates_std = zeros((%d,4))' % (ctyp, len(shock_pcts)))
#   exec('shock_%s_resp_rates2_mean = zeros((%d,4))' % (ctyp, len(shock_pcts)))
#   exec('shock_%s_resp_rates2_std = zeros((%d,4))' % (ctyp, len(shock_pcts)))

from numpy import *
from matplotlib.pyplot import *
from matplotlib.mlab import *
rasterdemo()

Loading data, please wait...