Ca+/HCN channel-dependent persistent activity in multiscale model of neocortex (Neymotin et al 2016)

 Download zip file 
Help downloading and running models
Accession:185858
"Neuronal persistent activity has been primarily assessed in terms of electrical mechanisms, without attention to the complex array of molecular events that also control cell excitability. We developed a multiscale neocortical model proceeding from the molecular to the network level to assess the contributions of calcium regulation of hyperpolarization-activated cyclic nucleotide-gated (HCN) channels in providing additional and complementary support of continuing activation in the network. ..."
Reference:
1 . Neymotin SA, McDougal RA, Bulanova AS, Zeki M, Lakatos P, Terman D, Hines ML, Lytton WW (2016) Calcium regulation of HCN channels supports persistent activity in a multiscale model of neocortex. Neuroscience 316:344-66 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell; Synapse; Channel/Receptor; Molecular 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 spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron; Neocortex layer 2-3 interneuron; Neocortex layer 5 interneuron; Neocortex layer 6a interneuron;
Channel(s): I Na,t; I L high threshold; I T low threshold; I A; I K; I M; I h; I K,Ca; I CAN; I Calcium; I_AHP; I_KD; Ca pump;
Gap Junctions:
Receptor(s): mGluR1; GabaA; GabaB; AMPA; NMDA; mGluR; Glutamate; Gaba; IP3;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Activity Patterns; Ion Channel Kinetics; Oscillations; Spatio-temporal Activity Patterns; Signaling pathways; Working memory; Attractor Neural Network; Calcium dynamics; Laminar Connectivity; G-protein coupled; Rebound firing; Brain Rhythms; Dendritic Bistability; Reaction-diffusion; Beta oscillations; Persistent activity; Multiscale;
Implementer(s): Neymotin, Sam [Samuel.Neymotin at nki.rfmh.org]; McDougal, Robert [robert.mcdougal at yale.edu];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; mGluR1; GabaA; GabaB; AMPA; NMDA; mGluR; Glutamate; Gaba; IP3; I Na,t; I L high threshold; I T low threshold; I A; I K; I M; I h; I K,Ca; I CAN; I Calcium; I_AHP; I_KD; Ca pump; Gaba; Glutamate;
/
CaHDemo
readme.html
cagk.mod
cal.mod *
calts.mod *
can.mod *
cat.mod *
gabab.mod
IC.mod *
icalts.mod *
Ih.mod
ihlts.mod *
IKM.mod *
kap.mod
kcalts.mod *
kdmc.mod
kdr.mod
kdrbwb.mod
km.mod *
mglur.mod *
misc.mod
MyExp2SynBB.mod *
MyExp2SynNMDABB.mod
nafbwb.mod
nax.mod
vecst.mod *
aux_fun.inc *
conf.py
declist.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
geom.py
ghk.inc *
grvec.hoc
init.hoc
labels.hoc
labels.py *
local.hoc *
misc.h
mpisim.py
netcfg.cfg
nqs.hoc
nqs.py
nrnoc.hoc *
onepyr.cfg
onepyr.py
pyinit.py *
python.hoc *
pywrap.hoc *
screenshot.png
screenshot1.png
simctrl.hoc *
simdat.py
syncode.hoc *
xgetargs.hoc *
                            
from neuron import *
import sys
import os
import string
h("strdef simname, allfiles, simfiles, output_file, datestr, uname, comment")
h.simname=simname = "onepyr"
h.allfiles=allfiles = "pyinit.py geom.py onepyr.py"
h.simfiles=simfiles = "pyinit.py geom.py onepyr.py"
h("runnum=1")
h.datestr=datestr = "10dec15"
h.output_file=output_file = "data/15dec31.2"
h.uname=uname = "x86_64"
h("templates_loaded=0")
h("xwindows=1.0")
h.xopen("nrnoc.hoc")
h.xopen("init.hoc")
h("proc setMemb () { }") # so e_pas will not get modified
from neuron import h
from pyinit import *
from labels import *
from geom import *
from nqs import *
import random
from pylab import *
from datetime import datetime
from neuron import rxd
import shutil
Vector = h.Vector
# nsubseg == number of subsegments per segment, subsegum == subseg size in microns
# nsubseg > 0 takes precedence over subsegum
nsubseg,subsegum = dconf['nsubseg'],dconf['subsegum'] 

for k in ['simstr', 'saveout', 'recdt', 'recvdt', 'verbose', 'indir', 'outdir']:
  if type(dconf[k]) == str: exec(k + ' = \'' + str(dconf[k]) + '\'')
  else: exec(k + ' = ' + str(dconf[k]))  

# make dir, catch exceptions
def safemkdir (dn):
  try:
    os.mkdir(dn)
    return True
  except OSError:
    if not os.path.exists(dn):
      print 'could not create', dn
      return False
    else:
      return True

safemkdir('data')
safemkdir('meminput')

# backup the config file
def backupcfg (simstr):
  safemkdir('backupcfg')
  fout = 'backupcfg/' + simstr + '.cfg'
  if os.path.exists(fout): os.remove(fout)  
  try:
    shutil.copy2(fcfg, fout) # fcfg created in geom.py via conf.py
  except:
    pass

backupcfg(simstr) # backup the config file

# return True if s exists in file fname
def grepstr (fname, s):
  found = False
  try:
    fp = open(fname,'r')
    lns = fp.readlines()
    for ln in lns:
      if ln.count(s) > 0:
        found = True
        break
    fp.close()
  except:
    pass
  return found

h.tstop = tstop = dconf['tstop']
tstart = 0 # only use previous end time if loading state
h.dt = dconf['dt']
h.steps_per_ms = 1/h.dt
h.celsius = 37
h.fracca_MyExp2SynNMDABB = dconf['nmfracca'] # fraction of NMDA current that is from calcium
wmatex = numpy.zeros( ( STYPi,) ) # external weights
ratex = numpy.zeros( (STYPi,) )  # external rates
mGLURR = dconf['mGLURR'] # ratio of mGLUR weights to AM2 weights
mGLURRWM = dconf['mGLURRWM']; 
NMAMRWM = dconf['NMAMRWM']; AMRWM = dconf['AMRWM'] # WM related weight ratios
GARWM = dconf['GARWM']
EXGain = dconf['EXGain'] # gain for noise inputs , default 1.0
sgrhzEE = dconf['sgrhzEE'] # external E inputs to E cells; 1000 is default
sgrhzEI = dconf['sgrhzEI'] # external E inputs to I cells
sgrhzIE = dconf['sgrhzIE'] # external I inputs to E cells
sgrhzII = dconf['sgrhzII'] # external I inputs to I cells
sgrhzNME = dconf['sgrhzNME'] # external NM inputs to E cells; 10 is default
sgrhzMGLURE = dconf['sgrhzMGLURE'] # external mGLUR inputs to E cells
sgrhzGB2 = dconf['sgrhzGB2'] # external inputs onto E cell GB2 synapses

ratex[AM2]=sgrhzEE
ratex[NM2]=sgrhzNME
ratex[GA]=sgrhzIE
ratex[GA2]=sgrhzIE
ratex[GB2]=sgrhzGB2
wmatex[AM2] = 0.02e-3 
wmatex[NM2] = 0.02e-3 
wmatex[GA] =  0.2e-3 
wmatex[GA2] = 0.2e-3 
wmatex[GB2] = 5e-3
for sy in xrange(STYPi): wmatex[sy] *= EXGain # apply gain control

cell = PyrAdr(0,0,0,0,E5R)
ce = [cell]
rxdsec = [s for s in cell.all_sec] # Section list for use with rxd

#
def setnsubseg ():
  global mGLURRWM
  if nsubseg > 0:
    print 'setting rxd nsubseg to ', dconf['nsubseg']
    rxd.set_solve_type(nsubseg=dconf['nsubseg'])
    mGLURRWM /= float(nsubseg)
    print 'adjusted mGLURRWM to ' , mGLURRWM
  elif subsegum > 0.0:
    for sec in cell.all_sec:
      ns = int(sec.L / subsegum)
      if ns > 1:
        if ns % 2 == 0: ns += 1
        print 'L = ', sec.L, ' subseg = ', ns
        rxd.set_solve_type(sec, nsubseg=ns)
        if sec == cell.Adend3:
          mGLURRWM /= float(ns)
          print 'adjusted mGLURRWM to ' , mGLURRWM
      else:
        print 'L = ', sec.L, ' not using subseg'
  rxd.options.subseg_interpolation = 1
  rxd.options.subseg_averaging = 1

if nsubseg > 0 or subsegum > 0.0: setnsubseg()

# get a recording vector for the synaptic current
def recveccurr (sy):
  try:
    return saferecord(sy.syn._ref_i, recvdt)
  except:
    return saferecord(sy.syn._ref_iNMDA, recvdt)

cell.dsy = {}; dsstr = ['AMPA', 'NMDA', 'GABAS', 'mGLUR', 'GABAB']
cell.vsy = {}
#
def addsyns ():
  ds = {}; ds[cell.Adend3]='Adend3'; ds[cell.Adend2]='Adend2'; ds[cell.Adend1]='Adend1'; ds[cell.Bdend]='Bdend'
  for sec in [cell.Adend3, cell.Adend2, cell.Adend1, cell.Bdend]:
    llsy = [];
    nloc = sec.nseg
    llvsy = []; # for recording currents
    for i,seg in enumerate(sec):
      if seg.x == 0.0 or seg.x == 1.0: continue # skip endpoints
      lsy = []; loc = seg.x; lvsy = [] #AMPA, NMDA, GABAA_slow, GABAB
      lsy.append(Synapse(sect=sec,loc=loc,tau1=0.05,tau2=5.3,e=0)); lvsy.append(recveccurr(lsy[-1]))#AMPA
      lsy.append(SynapseNMDA(sect=sec,loc=loc,tau1NMDA=15,tau2NMDA=150,r=1,e=0)) # NMDA
      lvsy.append(recveccurr(lsy[-1]))
      lsy.append(Synapse(sect=sec,loc=loc,tau1=0.2,tau2=20,e=-80)) # GABAA_slow
      lvsy.append(recveccurr(lsy[-1]))
      lsy.append(SynapsemGLUR(sect=sec,loc=loc)) # mGLUR
      for node in ip3.nodes(seg): node.include_flux(lsy[-1].syn._ref_rip3 ) # all the sub-segments get flux
      lsy.append(SynapseGABAB(sect=sec,loc=loc)) # GABAB
      lvsy.append(recveccurr(lsy[-1]))
      llsy.append(lsy); llvsy.append(lvsy)
    cell.dsy[sec] = llsy; cell.vsy[sec] = llvsy
  sec = cell.soma; llsy = []; nloc = sec.nseg; llvsy = []
  for i,seg in enumerate(sec):
    if seg.x == 0.0 or seg.x == 1.0: continue # skip endpoints
    lsy = []; loc = seg.x; lvsy = []
    lsy.append(Synapse(sect=sec,loc=loc,tau1=0.07,tau2=9.1,e=-80)) # GABAA_fast
    lvsy.append(recveccurr(lsy[-1]))
    lsy.append(Synapse(sect=sec,loc=loc,tau1=0.05,tau2=5.3,e=0) ) # AMPA
    lvsy.append(recveccurr(lsy[-1]))
    lsy.append(SynapseNMDA(sect=sec,loc=loc,tau1NMDA=15,tau2NMDA=150,r=1,e=0)) # NMDA
    lvsy.append(recveccurr(lsy[-1]))
    llsy.append(lsy); llvsy.append(lvsy);
  cell.dsy[sec] = llsy; cell.vsy[sec] = llvsy;

[cyt,er,cyt_er_membrane,ca,caextrude,serca,leak,CB,caCB,buffering]=[None for i in xrange(10)]

vt=Vector(); vt.record(h._ref_t);
vdt=Vector(); vdt.record(h._ref_dt);

lncrec,ltimevec,lidvec=[],[],[] # spike recorders
lX,lY,lZ,lctyID=[0],[0],[0],[E2] # E2 pyramidal cell - nothing particular about the E2
timevec,idvec = h.Vector(),h.Vector()
ncrec = h.NetCon(ce[-1].soma(0.5)._ref_v, None, sec=ce[-1].soma)
ncrec.record(timevec,idvec,0)
ncrec.threshold = 0 # 10 mV is default, lower it
ltimevec.append(timevec); lidvec.append(idvec); lncrec.append(ncrec) # record the spikes at soma

rxd.options.use_reaction_contribution_to_jacobian = False # faster (checked a few days before 10/16/13)
fc, fe = 0.83, 0.17 # cytoplasmic, er volume fractions
cyt = rxd.Region(rxdsec, nrn_region='i', geometry=rxd.FractionalVolume(fc, surface_fraction=1))
er  = rxd.Region(rxdsec, geometry=rxd.FractionalVolume(fe))
cyt_er_membrane = rxd.Region(rxdsec, geometry=rxd.ScalableBorder(1))
caDiff = 0.233
ca = rxd.Species([cyt, er], d=caDiff, name='ca', charge=2, initial=dconf['cacytinit'])
caexinit = dconf['caexinit']
caextrude = rxd.Rate(ca, (caexinit-ca[cyt])/taurcada, regions=cyt, membrane_flux=False)

ip3 = rxd.Species(cyt, d=0.283, name='ip3', initial=0.0)
# action of IP3 receptor
Kip3=0.13; Kact=0.4
minf = ip3[cyt] * 1000. * ca[cyt] / (ip3[cyt] + Kip3) / (1000. * ca[cyt] + Kact)
ip3r_gate_state = rxd.State(cyt_er_membrane, initial=0.8)
h_gate = ip3r_gate_state[cyt_er_membrane]
k = dconf['gip3'] * (minf * h_gate) ** 3 
ip3r = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], k, k, membrane=cyt_er_membrane)    
# IP3 receptor gating
ip3rg = rxd.Rate(h_gate, (1. / (1 + 1000. * ca[cyt] / (0.4)) - h_gate) / 400.0)
# IP3 degradation - moves towards baseline level (ip3_init)
ip3degTau = 1000 # 1000 ms
ip3deg = rxd.Rate(ip3, (0.0-ip3[cyt])/ip3degTau, regions=cyt, membrane_flux=False)

def setmGLURflux (): # mGLUR synapses generate ip3 that is fed into rxd ip3
  for c in ce:
    if ice(c.ty): continue
    for syn,seg in zip([c.Adend3mGLUR.syn,c.Adend2mGLUR.syn,c.Adend1mGLUR.syn],[c.Adend3(0.5), c.Adend2(0.5), c.Adend1(0.5)]):
      for node in ip3.nodes(seg): 
        node.include_flux(syn._ref_rip3)

def setrecip3 ():
  for c in ce:
    if ice(c.ty): continue
    c.soma_ip3cyt = Vector(tstop/h.dt)
    c.soma_ip3cyt.record( ip3[cyt].nodes(c.soma)(0.5)[0]._ref_concentration)
    c.Adend3_ip3cyt = Vector(tstop/h.dt)
    c.Adend3_ip3cyt.record( ip3[cyt].nodes(c.Adend3)(0.5)[0]._ref_concentration)

# SERCA pump: pumps ca from cyt -> ER
Kserca = 0.1 # Michaelis constant for SERCA pump
gserca = dconf['gserca']
serca = rxd.MultiCompartmentReaction(ca[cyt]>ca[er],gserca*(1e3*ca[cyt])**2/(Kserca**2+(1e3*ca[cyt])**2),membrane=cyt_er_membrane,custom_dynamics=True)
gleak = dconf['gleak']   # leak channel: bidirectional ca flow btwn cyt <> ER
leak = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], gleak, gleak, membrane=cyt_er_membrane)

def setreccaer (): # setup recording of ca[er] for each pyramidal cell in Adend3,soma center
  for c in ce:
    if ice(c.ty): continue
    c.soma_caer = Vector(tstop/h.dt)
    c.soma_caer.record( ca[er].nodes(c.soma)(0.5)[0]._ref_concentration )
    c.Adend3_caer = Vector(tstop/h.dt)
    c.Adend3_caer.record( ca[er].nodes(c.Adend3)(0.5)[0]._ref_concentration)

CB_init = dconf["CB_init"]
CB_frate = dconf["CB_frate"]
CB_brate = dconf["CB_brate"]
CBDiff = 0.043   # um^2 / msec
CB = rxd.Species(cyt,d=CBDiff,name='CB',charge=0,initial=CB_init) # CalBindin (Anwar)
caCB = rxd.Species(cyt,d=CBDiff,name='caCB',charge=0,initial=0.0) # Calcium-CB complex
kCB = [CB_frate, CB_brate] # forward,backward buffering rates
buffering = rxd.Reaction(ca+CB <> caCB, kCB[0], kCB[1], regions=cyt)

def setreccacb (): # setup recording of caCB for each pyramidal cell in Adend3,soma center
  for c in ce:
    if ice(c.ty): continue
    c.soma_caCB = Vector(tstop/h.dt)
    c.soma_caCB.record( caCB.nodes(c.soma)(0.5)[0]._ref_concentration)
    c.Adend3_caCB = Vector(tstop/h.dt)
    c.Adend3_caCB.record( caCB.nodes(c.Adend3)(0.5)[0]._ref_concentration)

setreccaer() # NB: only record from RXD variables after ALL species setup!
setreccacb() # otherwise, the pointers get messed up.
setrecip3()
setmGLURflux()

addsyns() # add syns here since mGLUR depends on ip3 creation

# stimuli

#creates NetStims (and associated NetCon,Random) - provide 'noise' inputs
#returns next useable value of sead
def makeNoiseNetStim (cel,nsl,ncl,nrl,nrlsead,syn,w,ISI,time_limit,sead):
  rd2 = h.Random(); rd2.ACG(sead); rd2.uniform(0,1e3)
  ns = h.NetStim()
  ns.interval = ISI
  ns.noise = 1			
  ns.number = 2 * time_limit / ISI  # create enough spikes for extra time, in case goes over limit
  ns.start = 0 # tstart
  if type(syn) == str: nc = h.NetCon(ns,cel.__dict__[syn].syn)
  else: nc = h.NetCon(ns,syn)
  nc.delay = h.dt * 2 # 0
  nc.weight[0] = w
  rds = h.Random()
  rds.negexp(1)            # set random # generator using negexp(1) - avg interval in NetStim
  rds.MCellRan4(sead,sead) # seeds are in order, shouldn't matter			
  ns.noiseFromRandom(rds)  # use random # generator for this NetStim                
  ns.start = 0. # rd2.repick() # start inputs random time btwn 0-1e3 ms to avoid artificial sync
  nsl.append(ns)
  ncl.append(nc)
  nrl.append(rds)
  nrlsead.append(sead)
  cel.infod[syn] = [ns,nc] #store pointers to NetStim,NetCon for easier manipulation

#
def makeNoiseNetStims (simdur,rdmseed):
  nsl = [] #NetStim List
  ncl = [] #NetCon List
  nrl = [] #Random List for NetStims
  nrlsead = [] #List of seeds for NetStim randoms
  for cell in ce: # go through cell types, check weights,rates of inputs
    ct = cell.ty # get cell type code  
    for sec in [cell.Adend3, cell.Adend2, cell.Adend1]:
      llsy = cell.dsy[sec]
      for lsy in llsy:
        for i,sy in enumerate([AM2,NM2,GA2,AM,GB2]):
          if ratex[sy] > 0. and wmatex[sy] > 0. and sy != AM: # AM now for soma AMPA
            makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,lsy[i].syn,wmatex[sy],(1e3/ratex[sy]),simdur,rdmseed); rdmseed+=1
    sec = cell.Bdend; llsy = cell.dsy[sec];
    for lsy in llsy:
      for i,sy in enumerate([AM2,NM2,GA2]):
        if ratex[sy] > 0. and wmatex[sy] > 0. and sy != AM: # AM now for soma AMPA
          makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,lsy[i].syn,wmatex[sy],(1e3/ratex[sy]),simdur,rdmseed); rdmseed+=1
    sec = cell.soma; llsy = cell.dsy[sec];
    for i,sy in enumerate([GA,AM,NM]):
      if ratex[sy] > 0. and wmatex[sy] > 0.:
        for lsy in llsy:
          makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,lsy[i].syn,wmatex[sy],(1e3/ratex[sy]),simdur,rdmseed); rdmseed+=1
  return nsl,ncl,nrl,nrlsead

nsl,ncl,nrl,nrlsead = makeNoiseNetStims(tstart+tstop,dconf['iseed'])

# setup wm-related stimuli
useWM = dconf['useWM']
rnnsl = [] # based on row number -- should always have the same number of entries as rows in the nqm 
rnncl = [] # based on row number times number of cells in the row
nqm=None; rnrds = []; rnseed = []; fihSIGns=None; 
def setupWMInputs ():
  global nqm,rnnsl,rnncl,rnrds,rnseed,fihSIGns
  if os.path.exists(dconf['nqm']): # read nqm if it exists
    nqm = h.NQS(dconf['nqm'])
  else: # setup nqm here if nqm file path non-existant
    nqm = h.NQS('vid','startt','endt','rate','w'); nqm.odec('vid')
    lvid = []; nMem = dconf['nMem']; 
    memstartt = dconf['memstartt']; memintert = dconf['memintert']; memdurt = dconf['memdurt']
    memstartt += tstart # this only has an effect when loadstate != 0
    startt = memstartt; endt = startt+memdurt
    if verbose: print 'startt:',startt,',memstartt:',memstartt,',tstart:',tstart,',loadstate:',loadstate
    for i in xrange(nMem): # number of stimuli
      lvid.append(h.Vector()); vtmp=Vector()
      vtmp.append(0) # only one cell
      lvid[-1].append(vtmp)
      nqm.append(lvid[-1],startt,endt,dconf['memrate'],dconf['memW'])
      startt += (memdurt+memintert);
      if i == nMem-2 and dconf['lastmemdurt'] != memdurt:
        endt = startt + dconf['lastmemdurt']
      else:
        endt = startt + memdurt;
  # backup the nqm file for later retrieval
  fnqwm = 'meminput/' + simstr + '_nqm.nqs'
  if os.path.exists(fnqwm):
    os.remove(fnqwm)
  nqm.sv(fnqwm)
  if verbose:
    print 'this is nqm:'
    nqm.pr()

  def createNetStims (vid,rate,w,startt,endt,seed=1234,syn='Adend3AMPA'):
    global rnnsl,rnncl,rnrds,rnseed # based on row number in nqm
    rnnsl.append( [] ) # a list for each row of nqm
    rnncl.append( [] )
    rnrds.append( [] )
    rnseed.append( [] )
    sidx = -1
    slist=[cell.Adend3, cell.Adend2, cell.Adend1]
    if syn.count('AMPA'): sidx=0
    elif syn.count('NMDA'): sidx=1
    elif syn.count('GABAss'): sidx=2
    elif syn.count('mGLUR'): sidx=3
    if syn.count('0'): slist=[cell.Adend3, cell.Adend2, cell.Adend1]
    elif syn.count('1'): slist=[cell.Adend1]
    elif syn.count('2'): slist=[cell.Adend2]
    elif syn.count('3'): slist=[cell.Adend3]
    if sidx == -1:
      ns = h.NetStim()
      ns.start = startt
      ns.number = (endt-startt) * rate / 1e3
      if verbose: print 'createNetStims:',startt,endt,ns.number,rate,w
      ns.interval = 1e3 / rate
      ns.noise = 1
      rnnsl[-1].append(ns)
      nc = h.NetCon(ns,ce[0].__dict__[syn].syn)
      nc.delay = h.dt * 2
      nc.weight[0] = w
      rnncl[-1].append(nc)
      rds = h.Random()
      rds.negexp(1)
      rds.MCellRan4(seed,seed)
      ns.noiseFromRandom(rds)
      rnrds[-1].append(rds)
      rnseed[-1].append(seed)
    else:
      tmpseed = seed     
      for sec in slist:
        llsy = cell.dsy[sec]
        for lsy in llsy:
          ns = h.NetStim()
          ns.start = startt
          ns.number = (endt-startt) * rate / 1e3
          if verbose: print 'createNetStims:',startt,endt,ns.number,rate,w
          ns.interval = 1e3 / rate
          ns.noise = 1
          rnnsl[-1].append(ns)
          nc = h.NetCon(ns,lsy[sidx].syn)
          nc.delay = h.dt * 2
          nc.weight[0] = w
          rnncl[-1].append(nc)
          rds = h.Random()
          rds.negexp(1)
          rds.MCellRan4(tmpseed,tmpseed)
          ns.noiseFromRandom(rds)
          rnrds[-1].append(rds)
          rnseed[-1].append(tmpseed)
          tmpseed += 1

  def createNetStimsFromNQ (nqm,row,seed=1234,syn='Adend3AMPA',wfctr=1.0):
    nqm.tog("DB")
    vid = h.Vector()
    rate = nqm.getcol("rate").x[row]
    w = nqm.getcol("w").x[row] * wfctr # wfctr is weight scaling factor
    vid.copy(nqm.get("vid",row).o[0])
    startt = float(nqm.getcol("startt").x[row]) # + tstart
    endt = float(nqm.getcol("endt").x[row]) # + tstart
    createNetStims(vid,rate,w,startt,endt,seed,syn)

  def setStims ():
    global nqm
    lapicIDX = []
    try:
      lapicIDX = [int(dconf['apicIDX'])]
    except:
      for i in dconf['apicIDX'].split(','): lapicIDX.append(int(i))     
    nqm.tog("DB"); sz = int(nqm.v[0].size()); lastmGLURON = dconf['lastmGLURON']; nMemInhib = dconf['nMemInhib']
    if lastmGLURON:
      for i in xrange(sz):
        if i == sz-1:
          for j in lapicIDX: createNetStimsFromNQ(nqm,i,seed=(i+1)*12345,syn='Adend'+str(j)+'mGLUR',wfctr=mGLURRWM)
        else:
          for j in lapicIDX:
            createNetStimsFromNQ(nqm,i,seed=(i+1)*12345,syn='Adend'+str(j)+'AMPA',wfctr=AMRWM)
            createNetStimsFromNQ(nqm,i,seed=(i+1)*12345,syn='Adend'+str(j)+'NMDA',wfctr=NMAMRWM)
    elif nMemInhib > 0:
      for i in xrange(0,sz-nMemInhib,1):
        for j in lapicIDX:
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'AMPA',wfctr=AMRWM)
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'NMDA',wfctr=NMAMRWM)
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'mGLUR',wfctr=mGLURRWM)
        for i in xrange(sz-nMemInhib,sz,1):
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'GABAss',wfctr=GARWM)
    else:
      for i in xrange(sz):
        for j in lapicIDX:
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'AMPA',wfctr=AMRWM)
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'NMDA',wfctr=NMAMRWM)
          createNetStimsFromNQ(nqm,i,seed=j*(i+1)*12345,syn='Adend'+str(j)+'mGLUR',wfctr=mGLURRWM)

  def initSigNetStims ():
    for i in xrange(len(rnnsl)):
      for j in xrange(len(rnnsl[i])):
        rds = rnrds[i][j]
        sead = rnseed[i][j]
        rds.MCellRan4(sead,sead)
        rds.negexp(1)			

  fihSIGns = h.FInitializeHandler(0, initSigNetStims)
  setStims() # create the inputs based on contents of nqm

#this should be called @ beginning of each sim - done in an FInitializeHandler
def init_NetStims ():
  for i in xrange(len(nrl)):
    rds = nrl[i]
    sead = nrlsead[i]
    rds.MCellRan4(sead,sead)
    rds.negexp(1)			

fihns = h.FInitializeHandler(0, init_NetStims)

#
def printt ():
  sys.stdout.write('\rt: {0} ...'.format(h.t))
  sys.stdout.flush()

# handler for printing out time during simulation run
def fi():
  for i in xrange(int(tstart),int(tstart+tstop),100): h.cvode.event(i, printt)

fih = h.FInitializeHandler(1, fi)

##### RECORDINGS

# get x locations in a line, ordered from Bdend to Adend3
def getlx ():
  lx = []; startx = 0
  for sdx in [1,0,2,3,4]:
    sec = cell.all_sec[sdx]
    LX = [node.x for node in ca[cyt].nodes(sec)]
    if sdx == 1:
      LX.reverse()
      for x in LX: lx.append(-x*sec.L + startx)
    else:
      for x in LX: lx.append(x*sec.L + startx)
    if sdx != 1: startx += sec.L
  return lx

# cell.snames # ['soma', 'Bdend', 'Adend1', 'Adend2', 'Adend3']
#
lcai = []; lcaer = []; lv = []; lip3 = []
for sdx in [1,0,2,3,4]:
  sec = cell.all_sec[sdx]
  lx = [node.x for node in ca[cyt].nodes(sec)]
  if sdx == 1: lx.reverse() # since Bdend connected at other end
  for x in lx:
    if len(ca[er].nodes(sec)(x)) > 0:
      vc = Vector(); vc.record(ca[cyt].nodes(sec)(x)[0]._ref_concentration)
      lcai.append(vc)
      vc = Vector(); vc.record(ca[er].nodes(sec)(x)[0]._ref_concentration)
      lcaer.append(vc)
      vc = Vector(); vc.record(sec(x)._ref_v)
      lv.append(vc)
      vc = Vector(); vc.record(ip3.nodes(sec)(x)[0]._ref_concentration)
      lip3.append(vc)

# plot currents
def plotcurr (dstr='Adend3',dsec=cell.Adend3):
  figure(); subplot(1,3,1)
  tv = arange(0,tstop,recvdt)
  llvsy = cell.vsy[dsec]
  idx = int(len(llvsy) / 2)
  clr = ['r','g','b','k']
  for i in xrange(4): plot(tv,llvsy[idx][i],clr[i])
  legend(['AM','NM','GA','GB'],loc='best')
  xlabel('Time (ms)'); ylabel('i')  
  subplot(1,3,2)
  li = [cell.__dict__[dstr + s] for s in ['_ik','_ina','_ica','_ih']]
  for vec,clr in zip(li,['b','g','r', 'c']): plot(tv,vec.as_numpy(),clr)
  legend(['ik','ina','ica','ih'],loc='best'); xlabel('Time (ms)'); ylabel(dstr + ' i')
  subplot(1,3,3)
  for vec,clr in zip([cell.soma_ik,cell.soma_ina,cell.soma_ica,cell.soma_ih],['b','g','r', 'c']): plot(tv,vec.as_numpy(),clr)
  legend(['ik','ina','ica','ih'],loc='best'); xlabel('Time (ms)'); ylabel('soma i')

# get sag - for hyperpolarizing pulses - difference between minimum and steadystate
def getsag (dat,sampr,stimon=0.5,verbose=False):
  vpeak = amin(dat[int(stimon*sampr):int((stimon+0.1)*sampr)])
  vsteady = mean(dat[int((stimon+0.89)*sampr):int((stimon+0.99)*sampr)])
  sag = (vpeak-vsteady)/vsteady
  if verbose: print 'sag:',sag,'peak:',vpeak, 'steady:',vsteady
  return sag

#
def getrin (dat,dt,ic):
  sidx,eidx = int(ic.delay/dt),int((ic.delay+ic.dur)/dt)
  if ic.amp < 0: dv = dat[sidx] - amin(dat[sidx+1:eidx])
  else: dv = amax(dat[sidx+1:eidx]) - dat[sidx]
  di = abs(ic.amp)
  print 'dv:',dv,'di:',di,'rin:',dv/di
  return dv/di

#
def plotout (dstr='Adend3',xl=((0,tstop/1e3))):
  lx = getlx(); cellL = [lx[0],lx[-1]]
  tv = arange(0,tstop,h.dt)
  figure(); npa=numpy.array(lcai);
  subplot(2,2,1); imshow(npa,origin='lower',aspect='auto',extent=(0,tstop/1e3,lx[0],lx[-1]))
  title(r'$Ca^{2+}_{cyt}$'); colorbar(); xlim(xl)
  subplot(2,2,2); npb=numpy.array(lcaer);
  imshow(npb,origin='lower',aspect='auto',extent=(0,tstop/1e3,lx[0],lx[-1]))
  title(r'$Ca^{2+}_{ER}$'); colorbar(); xlim(xl)
  subplot(2,2,3); npv=numpy.array(lv);
  imshow(npv,origin='lower',aspect='auto',extent=(0,tstop/1e3,lx[0],lx[-1]))
  title('v (mV)'); colorbar(); xlim(xl)
  subplot(2,2,4); npi=numpy.array(lip3);
  imshow(npi,origin='lower',aspect='auto',extent=(0,tstop/1e3,lx[0],lx[-1]))
  title(r'$IP_3 (mM)$'); colorbar(); xlim(xl)
  tvec = vt.as_numpy()
  if len(tv) == cell.soma_volt.size(): tvec = tv
  figure(); 
  subplot(1,4,1); plot(tvec/1e3,cell.soma_cai,'b'); plot(tvec/1e3,cell.__dict__[dstr+'_cai'],'r');title(r'$Ca^{2+}_{cyt}$');
  legend(['soma','dend'],loc='best'); xlim(xl)
  subplot(1,4,2);plot(tvec/1e3,cell.soma_Ihp1,'b');plot(tvec/1e3,cell.__dict__[dstr+'_Ihp1'],'r');title(r'$cAMP$');legend(['soma','dend'],loc='best'); xlim(xl)
  subplot(1,4,3); plot(tvec/1e3,cell.soma_Ihm,'b'); plot(tvec/1e3,cell.__dict__[dstr+'_Ihm'],'r');title(r'$g_h$');legend(['soma','dend'],loc='best'); xlim(xl)
  subplot(1,4,4);
  plot(tvec/1e3,cell.soma_volt,'b'); plot(tvec/1e3,cell.__dict__[dstr+'_volt'],'r'); title('v'); legend(['soma','dend'],loc='best');  xlim(xl)

def plotfixeddt (vec,clr,lw): plot(linspace(0,tstop/1e3,vec.size()),vec,clr,linewidth=lw)

#
def plotnew (dstr='Adend3',dsec=cell.Adend3,xl=((0,tstop/1e3)),useyl=True):
  naxbins=5; tx,ty=-.07,1.04; fsz=20; lw=2
  global vt
  vtime = vt.as_numpy() / 1e3
  figure(); 
  # synaptic current i
  ax=subplot(2,4,1); ax.locator_params(nbins=naxbins); text(tx,ty,'a',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
  llvsy = cell.vsy[dsec]
  idx = int(len(llvsy) / 2)
  clr = ['r','g','b']
  if not dconf['cvodeactive']:
    for i in xrange(3): plotfixeddt(llvsy[idx][i],clr[i],lw)
  else:
    for i in xrange(3): plot(vtime,llvsy[idx][i],clr[i],linewidth=lw)
  #legend([r'$i_{AMPA}$',r'$i_{NMDA}$',r'$i_{GABAA}$'],loc='best')
  #xlabel('Time (s)'); 
  title(r'$syn_i$ (nA)',fontsize=fsz); xlim(xl); 
  if useyl: ylim((-1.22,.06)); #ylim((-.64,.06))
  # cai
  ax=subplot(2,4,2); ax.locator_params(nbins=naxbins);
  if not dconf['cvodeactive']:
    plotfixeddt(cell.soma_cai,'b',lw)
    plotfixeddt(cell.__dict__[dstr+'_cai'],'r',lw);
  else:
    plot(vtime,cell.soma_cai,'b',linewidth=lw); 
    plot(vtime,cell.__dict__[dstr+'_cai'],'r',linewidth=lw);
  title(r'$Ca^{2+}_{cyt}$ (mM)',fontsize=fsz); xlim(xl); 
  if useyl: ylim((0,0.0094))  
  #legend(['soma','dend'],loc='best'); 
  text(tx,ty,'b',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
  # ip3
  ax=subplot(2,4,3); ax.locator_params(nbins=naxbins);
  if not dconf['cvodeactive']:
    plotfixeddt(cell.soma_ip3cyt,'b',lw);
    plotfixeddt(cell.__dict__[dstr+'_ip3cyt'],'r',lw);
  else:
    plot(vtime,cell.soma_ip3cyt,'b',linewidth=lw);
    plot(vtime,cell.__dict__[dstr+'_ip3cyt'],'r',linewidth=lw);
  title(r'$IP_3$ (mM)',fontsize=fsz); xlim(xl); 
  if useyl: ylim((0,0.214))  
  text(tx,ty,'c',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
  # caer
  ax=subplot(2,4,4); ax.locator_params(nbins=naxbins);
  if not dconf['cvodeactive']:
    plotfixeddt(cell.soma_caer,'b',lw); 
    plotfixeddt(cell.__dict__[dstr+'_caer'],'r',lw);
  else:
    plot(vtime,cell.soma_caer,'b',linewidth=lw); 
    plot(vtime,cell.__dict__[dstr+'_caer'],'r',linewidth=lw);
  title(r'$Ca^{2+}_{ER}$ (mM)',fontsize=fsz); xlim(xl); 
  if useyl: ylim((0,1.3))
  text(tx,ty,'d',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
  # cAMP
  ax=subplot(2,4,5); ax.locator_params(nbins=naxbins);
  if not dconf['cvodeactive']:
    plotfixeddt(cell.soma_Ihp1,'b',lw); 
    plotfixeddt(cell.__dict__[dstr+'_Ihp1'],'r',lw);
  else:
    plot(vtime,cell.soma_Ihp1,'b',linewidth=lw); 
    plot(vtime,cell.__dict__[dstr+'_Ihp1'],'r',linewidth=lw);
  title(r'$cAMP$',fontsize=fsz); xlim(xl); 
  if useyl: ylim((0,0.025))
  text(tx,ty,'e',fontsize=fsz,fontweight='bold',transform=ax.transAxes); xlabel('Time (s)',fontsize=fsz);
  # gh
  ax=subplot(2,4,6); ax.locator_params(nbins=naxbins);
  if not dconf['cvodeactive']:
    plotfixeddt(cell.soma_Ihm,'b',lw); 
    plotfixeddt(cell.__dict__[dstr+'_Ihm'],'r',lw);
  else:
    plot(vtime,cell.soma_Ihm,'b',linewidth=lw); 
    plot(vtime,cell.__dict__[dstr+'_Ihm'],'r',linewidth=lw);
  title(r'$g_h$',fontsize=fsz); xlim(xl); 
  if useyl: ylim((0.011,0.168))
  text(tx,ty,'f',fontsize=fsz,fontweight='bold',transform=ax.transAxes); xlabel('Time (s)',fontsize=fsz);
  # dend ions
  ax=subplot(2,4,7); ax.locator_params(nbins=naxbins); 
  text(tx,ty,'g',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
  #li = [cell.__dict__[dstr + s] for s in ['_ik','_ina','_ica','_ih']]
  #for vec,clr in zip(li,['b','g','r','c']): plot(tv/1e3,vec.as_numpy(),clr)
  li = [cell.__dict__[dstr + s] for s in ['_ih']]
  if not dconf['cvodeactive']:
    for vec,clr in zip(li,['r']): plotfixeddt(vec,clr,lw)
  else:
    for vec,clr in zip(li,['r']): plot(vtime,vec,clr,linewidth=lw)
  #legend([r'$i_K$',r'$i_{Na}$',r'$i_{Ca}$',r'$i_h$'],loc='best'); 
  xlabel('Time (s)',fontsize=fsz); title(r'$i_h$ (nA)',fontsize=fsz); xlim(xl); 
  if useyl: ylim((-.102,-.022))
  if False:
    # soma ions
    ax=subplot(2,4,7); ax.locator_params(nbins=naxbins); 
    text(tx,ty,'h',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
    for vec,clr in zip([cell.soma_ik,cell.soma_ina],['b','r']): plot(vtime,vec,clr,linewidth=lw)
    #legend([r'$i_K$',r'$i_{Na}$'],loc='best'); 
    xlabel('Time (s)',fontsize=fsz); title(r'$soma_i$ (nA)',fontsize=fsz); xlim(xl); 
    if useyl: ylim((-15.5,11.5))
  # voltage
  ax=subplot(2,4,8); ax.locator_params(nbins=naxbins); text(tx,ty,'h',fontsize=fsz,fontweight='bold',transform=ax.transAxes);
  if not dconf['cvodeactive']:
    plotfixeddt(cell.soma_volt,'b',lw); 
    plotfixeddt(cell.__dict__[dstr+'_volt'],'r',lw); 
  else:
    plot(vtime,cell.soma_volt,'b',linewidth=lw); 
    plot(vtime,cell.__dict__[dstr+'_volt'],'r',linewidth=lw); 
  title('Vm',fontsize=fsz); xlim(xl)
  xlabel('Time (s)',fontsize=fsz); 
  if useyl: ylim((-85,40))

if dconf['cvodeactive']: h.cvode.active(1) # turn on variable time-step integration

lspks,lids=array([]),array([])

# interpolate voltage recorded in simulation to a fixed grid (dt millisecond spacing)
# output time,voltage returned
def interpvolt (tsrc,vsrc,dt):
  tdest = h.Vector(); tdest.indgen(0,tstop,dt)
  vdest = h.Vector(); vdest.interpolate(tdest,tsrc,vsrc)
  return tdest,vdest

#
def myrun ():
  global fit # pointer to FinitializeHandler needs to be global so doesn't disappear
  dastart = datetime.now()
  print 'started at:',dastart
  if dconf['useWM']: setupWMInputs()
  h.stdinit()
  ca[er].concentration = dconf['caerinit'] # 100e-6
  ca[cyt].concentration = dconf['cacytinit'] # 100e-6
  if dconf['cvodeactive']: h.cvode.re_init()
  h.continuerun(h.t+tstop) # run for tstop
  daend = datetime.now()
  print 'finished ',tstop,' ms sim at:',daend
  dadiff = daend - dastart;
  print 'runtime:',dadiff, '(',tstop/1e3,' s)'
  # concatenate the results so can view/save all at once
  global lspks,lids
  lspks,lids=array([]),array([])
  for i in xrange(len(ltimevec)):
    lspks=concatenate((lspks,ltimevec[i]))
    lids=concatenate((lids,lidvec[i]))    
  ion(); plotnew(xl=(9.9,11))

myrun()

# save data - output path based on simstr 
def savedata (simstr):
  safemkdir(outdir) # make sure the output directory exists
  fn = outdir + '/' + simstr + '.npz'
  print ' saving to ' , fn
  ne,ni,szslow,szfast = 0,0,0,0
  lE,lI=[],[]
  for c in ce:
    if ice(c.ty):
      lI.append(c.ID)
      ni += 1
    else:
      szslow = int(c.soma_cai.size()) # only E cells
      lE.append(c.ID)
      ne += 1
    szfast = int(c.soma_volt.size())
  lE=array(lE) # lE is list of E cell IDs from this host
  lI=array(lI) # Li is list of I cell IDs from this host
  soma_volt = zeros((ne,szfast)); Adend3_volt = zeros((ne,szfast)); Bdend_volt=zeros((ne,szfast));
  soma_cai = zeros((ne,szslow)); Adend3_cai = zeros((ne,szslow))
  soma_Ihm = zeros((ne,szslow)); Adend3_Ihm = zeros((ne,szslow))
  soma_Ihp1 = zeros((ne,szslow)); Adend3_Ihp1 = zeros((ne,szslow))
  soma_voltI = zeros((ni,szfast));
  soma_caer = zeros((ne,szslow));  Adend3_caer = zeros((ne,szslow));
  soma_caCB=zeros((ne,szslow)); Adend3_caCB=zeros((ne,szslow));
  soma_ip3cyt=zeros((ne,szslow)); Adend3_ip3cyt=zeros((ne,szslow));
  saveExtra = dconf['saveExtra']
  if saveExtra:
    Adend3_iAM = zeros((ne,szfast)); Adend3_iNM = zeros((ne,szfast))
    Adend3_iGB = zeros((ne,szfast)); Adend3_iGA = zeros((ne,szfast))
    soma_iGA = zeros((ne,szfast))
    [Adend3_ina, Adend3_ik, Adend3_ica, Adend3_ih] = [zeros((ne,szfast)) for i in xrange(4)]
    [soma_ina, soma_ik, soma_ica, soma_ih] = [zeros((ne,szfast)) for i in xrange(4)]
  cdx = 0; idx = 0;
  for c in ce:
    if ice(c.ty):
      soma_voltI[idx,:] = c.soma_volt.to_python()
      idx += 1
      continue
    soma_volt[cdx,:] = c.soma_volt.to_python()
    Adend3_volt[cdx,:] = c.Adend3_volt.to_python()
    Bdend_volt[cdx,:] = c.Bdend_volt.to_python()
    soma_cai[cdx,:] = c.soma_cai.to_python()
    Adend3_cai[cdx,:] = c.Adend3_cai.to_python()
    soma_Ihm[cdx,:] = c.soma_Ihm.to_python()
    Adend3_Ihm[cdx,:] = c.Adend3_Ihm.to_python()
    soma_Ihp1[cdx,:] = c.soma_Ihp1.to_python()
    Adend3_Ihp1[cdx,:] = c.Adend3_Ihp1.to_python()
    soma_caer[cdx,:] = c.soma_caer.to_python()
    Adend3_caer[cdx,:] = c.Adend3_caer.to_python()
    soma_caCB[cdx,:] = c.soma_caCB.to_python()
    Adend3_caCB[cdx,:] = c.Adend3_caCB.to_python()
    soma_ip3cyt[cdx,:] = c.soma_ip3cyt.to_python()
    Adend3_ip3cyt[cdx,:] = c.Adend3_ip3cyt.to_python()        
    if saveExtra:
      Adend3_iAM[cdx,:] = c.Adend3_iAM.to_python()
      Adend3_iNM[cdx,:] = c.Adend3_iNM.to_python()
      Adend3_iGB[cdx,:] = c.Adend3_iGB.to_python()
      Adend3_iGA[cdx,:] = c.Adend3_iGA.to_python()
      soma_iGA[cdx,:] = c.soma_iGA.to_python()
      Adend3_ina[cdx,:] = c.Adend3_ina.to_python()
      Adend3_ik[cdx,:] = c.Adend3_ik.to_python()
      Adend3_ica[cdx,:] = c.Adend3_ica.to_python()
      Adend3_ih[cdx,:] = c.Adend3_ih.to_python()
      soma_ina[cdx,:] = c.soma_ina.to_python()
      soma_ik[cdx,:] = c.soma_ik.to_python()
      soma_ica[cdx,:] = c.soma_ica.to_python()
      soma_ih[cdx,:] = c.soma_ih.to_python()    
    cdx += 1
  if saveExtra: # with somaVolt, synaptic currents
    numpy.savez(fn,lctyID=array(lctyID),lX=array(lX),lY=array(lY),lZ=array(lZ),vt=numpy.array(vt.to_python()),lspks=lspks,lids=lids,lE=lE,lI=lI,soma_volt=soma_volt,soma_voltI=soma_voltI,Adend3_volt=Adend3_volt,Bdend_volt=Bdend_volt,soma_cai=soma_cai,Adend3_cai=Adend3_cai,soma_Ihm=soma_Ihm,Adend3_Ihm=Adend3_Ihm,soma_Ihp1=soma_Ihp1,Adend3_Ihp1=Adend3_Ihp1,soma_caer=soma_caer,Adend3_caer=Adend3_caer,soma_caCB=soma_caCB,Adend3_caCB=Adend3_caCB,Adend3_ip3cyt=Adend3_ip3cyt,soma_ip3cyt=soma_ip3cyt,Adend3_iAM=Adend3_iAM,Adend3_iNM=Adend3_iNM,Adend3_iGB=Adend3_iGB,Adend3_iGA=Adend3_iGA,soma_iGA=soma_iGA,Adend3_ina=Adend3_ina,Adend3_ik=Adend3_ik,Adend3_ica=Adend3_ica,Adend3_ih=Adend3_ih,soma_ina=soma_ina,soma_ik=soma_ik,soma_ica=soma_ica,soma_ih=soma_ih)
  else: # no somaVolt, synaptic currents
    numpy.savez(fn,lctyID=array(lctyID),lX=array(lX),lY=array(lY),lZ=array(lZ),vt=numpy.array(vt.to_python()),lspks=lspks,lids=lids,lE=lE,lI=lI,Adend3_volt=Adend3_volt,Bdend_volt=Bdend_volt,soma_cai=soma_cai,Adend3_cai=Adend3_cai,soma_Ihm=soma_Ihm,Adend3_Ihm=Adend3_Ihm,soma_Ihp1=soma_Ihp1,Adend3_Ihp1=Adend3_Ihp1,soma_caer=soma_caer,Adend3_caer=Adend3_caer,soma_caCB=soma_caCB,Adend3_caCB=Adend3_caCB,Adend3_ip3cyt=Adend3_ip3cyt,soma_ip3cyt=soma_ip3cyt)

if saveout:
  print 'saving data'
  savedata(simstr)

# get an array of times of interest
def gettimes (nqm,rampt=3e3):
  nqm.tog("DB")
  nMem = int(nqm.size())
  baselinest = dconf['baset']
  times = [] # periods of interest
  vstartt,vendt = Vector(),Vector()
  vstartt.copy(nqm.getcol("startt"))
  vendt.copy(nqm.getcol("endt"))
  #vstartt.add(tstart); vendt.add(tstart)
  if nMem > 0:
    times = [(baselinest,vstartt[0])]
    times.append((vstartt[0],vendt[0]))
    times.append((vendt[0],vendt[0]+rampt))
    for row in xrange(1,nMem):
      times.append((vendt[row-1]+rampt,vstartt[row]))
      times.append((vstartt[row],vendt[row]))
      times.append((vendt[row],vendt[row]+rampt))
    lastt = vendt[-1]+rampt
    inct = 1e3; wint = 1e3
    while lastt < tstop:
      times.append((lastt,lastt+wint))
      lastt += inct
  return times

# get an NQS with spikes (snq)
def getsnq ():
  snq = NQS('id','t','ty','ice')
  snq.v[0].from_python(lids)
  snq.v[1].from_python(lspks)
  for i in xrange(len(lids)):
    snq.v[2].append(lctyID[int(lids[i])])
    snq.v[3].append(h.ice(lctyID[int(lids[i])]))
  return snq

# return x with only two decimal places after x
def twop (x): return round(x,2)

# print spike rates in each population in the periods of interest (baseline, signal, recall)
def prspkcount (snq,nqm,rampt=1e3,times=None,quiet=False):
  snq.verbose=0
  if times is None: times = gettimes(nqm,rampt) # periods of interest
  lfa=[]
  nqm.tog("DB")
  nMem = int(nqm.size())
  for row in xrange(nMem):
    lfa.append([]); 
    vact = Vector(1); vact.x[0] = 0
    for (startt,endt) in times: 
      fa = twop(1e3*snq.select("id","EQW",vact,"t","[]",startt,endt) / ((endt-startt)*vact.size()))
      lfa[-1].append(fa)
      if not quiet:
        print "t=",startt,"-",endt,". ActE:",fa,"Hz."
  snq.verbose=1
  return lfa

# print spike rates in each population in the periods of interest (baseline, signal, recall)
def getspkcount (snq,binsz,quiet=False):
  snq.verbose=0
  times = arange(0,h.tstop,binsz)
  lfa=[]
  for startt in times:
    lfa.append([]); 
    vact = Vector(1); vact.x[0] = 0
    fa = twop(1e3*snq.select("id","EQW",vact,"t","[]",startt,startt+binsz) / (binsz*vact.size()))
    lfa[-1].append(fa)
    if not quiet: print "t=",startt,"-",startt+binsz,". ActE:",fa,"Hz."
  snq.verbose=1
  return times,lfa

#
def getonecellfrd (lfa,binsz,t1,t2):
  a1 = mean(lfa[int(t1[0]/binsz):int(t1[1]/binsz)])
  a2 = mean(lfa[int(t2[0]/binsz):int(t2[1]/binsz)])
  if a1 > 0.: frd = (a2 - a1) / a1
  else: frd = 0.
  return a1,a2,frd

# print spike counts & draw activated and non-activated rates
def drawrat (ld,incsz=500.0,winsz=1e3,quiet=False):
  snq=getsnq(ld);
  times=[]
  for i in xrange(2*int(tstop/winsz)):times.append( (i*incsz,(i*incsz)+winsz) )
  lfa,lfn = prspkcount(snq,nqm,rampt=1e3,times=times,quiet=quiet)
  rat = []; tt=[];
  for i in xrange(len(lfa[0])): rat.append(lfa[0][i]/lfn[0][i])
  for tup in times: tt.append(0.5*(tup[0]+tup[1])) # time midpoints
  plot(tt,rat,'b',linewidth=2);plot(tt,rat,'bo',markersize=10);tight_layout();show() #

if dconf['doquit']: quit()