from neuron import *
# MPI
pc = h.ParallelContext() # MPI: Initialize the ParallelContext class
nhosts = int(pc.nhost()) # Find number of hosts
if pc.id()==0: print('\nSetting up network...')
import sys
import os
import string
h("strdef simname, allfiles, simfiles, output_file, datestr, uname, osname, comment")
h.simname=simname = "dystdemo"
h.allfiles=allfiles = "pyinit.py geom.py mpisim.py"
h.simfiles=simfiles = "pyinit.py geom.py mpisim.py"
h("runnum=1")
runnum = 1.0
h.datestr=datestr = "16jun21"
h.output_file=output_file = "data/16jun21.14"
h.uname=uname = "x86_64"
h.osname=osname="linux"
h("templates_loaded=0")
templates_loaded=0
h("xwindows=1.0")
xwindows = 1.0
h.xopen("nrnoc.hoc")
h.xopen("init.hoc")
h("proc setMemb () { }") # so e_pas will not get modified
CTYPi = 60.0
STYPi = 20.0
from pyinit import *
from labels import *
delm = numpy.zeros( (CTYPi, CTYPi) )
deld = numpy.zeros( (CTYPi, CTYPi) )
pmat = numpy.zeros( (CTYPi, CTYPi) )
synloc = numpy.zeros( (CTYPi, CTYPi) )
from geom import *
from vector import *
from nqs import *
import random
from pylab import *
from datetime import datetime
#########################################################################
# global params
verbose = dconf['verbose']
ISEED = dconf['iseed']
WSEED = dconf['wseed']
PSEED = dconf['pseed']
scale = dconf['scale']
gGID = 0 # global ID for cells
pmatscale = 1.0 # 1.0 / scale
spiketh = -15 # spike threshold, 10 mV is NetCon default, lower it for all cells
simstr = dconf['simstr']
saveout = dconf['saveout']
recdt = dconf['recdt']
recvdt = dconf['recvdt']
saveconns = dconf['saveconns']
indir = dconf['indir']
outdir = dconf['outdir']
# 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
# backup the config file
def backupcfg (simstr):
safemkdir('backupcfg')
fout = 'backupcfg/' + simstr + '.cfg'
if os.path.exists(fout):
print 'removing prior cfg file' , fout
os.system('rm ' + fout)
os.system('cp ' + fcfg + ' ' + fout) # fcfg created in geom.py via conf.py
if pc.id()==0:
backupcfg(simstr) # backup the config file
print "config file is " , fcfg
h.tstop = tstop = dconf['tstop']
tstart = 0.0
h.dt = dconf['dt']
h.steps_per_ms = 1/h.dt
h.v_init = -65
h.celsius = 37
h.fracca_MyExp2SynNMDABB = dconf['nmfracca'] # fraction of NMDA current that is from calcium
rdmsec = dconf['rdmsec']
EEGain = dconf['EEGain']
EIGainFS = dconf['EIGainFS']
EIGainLTS = dconf['EIGainLTS']
IEGain = dconf['IEGain']
IIGain = dconf['IIGain']
IIGainLTSFS = IIGain
IIGainFSLTS = IIGain
IIGainLTSLTS = IIGain
IIGainFSFS = IIGain
GB2R = dconf['GB2R'];
NMAMREE = dconf['NMAMREE']
NMAMREI = dconf['NMAMREI']
mGLURR = dconf['mGLURR'] # ratio of mGLUR weights to AM2 weights
cpernet = [] # cells of a given type for network
wmat = numpy.zeros( (CTYPi, CTYPi, STYPi) ) # internal weights
wmatex = numpy.zeros( (CTYPi, STYPi) ) # external weights
ratex = numpy.zeros( (CTYPi, STYPi) ) # external rates
EXGain = dconf['EXGain']
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
sgrhzNMI = dconf['sgrhzNMI'] # external NM inputs to I cells
sgrhzMGLURE = dconf['sgrhzMGLURE'] # external mGLUR inputs to E cells
sgrhzGB2 = dconf['sgrhzGB2'] # external inputs onto E cell GB2 synapses
# params for swire
colside = 120 # 120 microns squared
slambda = 100 # space constant for wiring falloff (used for I -> X)
axonalvelocity = 10000 # axonal velocity in um/ms -- this is 10 mm/s
#########################################################################
# setwmatex - set weights of external inputs to cells
def setwmatex ():
for ct in xrange(CTYPi):
for sy in xrange(STYPi):
ratex[ct][sy]=0
wmatex[ct][sy]=0
for ct in xrange(CTYPi):
if cpernet[ct] <= 0: continue
if IsLTS(ct): # dendrite-targeting interneurons (LTS cells)
ratex[ct][AM2]=sgrhzEI
ratex[ct][NM2] = sgrhzNMI
ratex[ct][GA]=sgrhzII
ratex[ct][GA2]=sgrhzII
wmatex[ct][AM2] = 0.02e-3
wmatex[ct][NM2] = 0.02e-3
wmatex[ct][GA]= 0
wmatex[ct][GA2]= 0.2e-3 # * 0
elif ice(ct): # soma-targeting interneurons (basket/FS cells)
ratex[ct][AM2]=sgrhzEI
ratex[ct][NM2] = sgrhzNMI
ratex[ct][GA]=sgrhzII
ratex[ct][GA2]=sgrhzII
wmatex[ct][AM2] = 0.02e-3 * 5.0
wmatex[ct][NM2] = 0.02e-3 * 5.0
wmatex[ct][GA]= 0
wmatex[ct][GA2]= 0.2e-3
else: # E cells
ratex[ct][MG]=sgrhzMGLURE
ratex[ct][AM2]=sgrhzEE
ratex[ct][NM2]=sgrhzNME
ratex[ct][GA]=sgrhzIE
ratex[ct][GA2]=sgrhzIE
ratex[ct][GB2]=sgrhzGB2
wmatex[ct][MG] = 1.5
wmatex[ct][AM2] = 0.02e-3
wmatex[ct][NM2] = 0.02e-3
wmatex[ct][GA] = 0.2e-3
wmatex[ct][GA2] = 0.2e-3
wmatex[ct][GB2] = 5e-3
for sy in xrange(STYPi): wmatex[ct][sy] *= EXGain # apply gain control
# set number of cells of a type in the network at scale==1
def setcpernet ():
global cpernet
cpernet = []
for i in xrange(CTYPi): cpernet.append(0)
cpernet[E2] = 300
cpernet[I2] = 37
cpernet[I2L] = 37
cpernet[E4] = 173
cpernet[E5R] = 150
cpernet[E5B] = 196
cpernet[E5P] = 196
cpernet[I5] = 89
cpernet[I5L] = 89
cpernet[E6] = 179
cpernet[E6C] = 179
cpernet[I6] = 45
cpernet[I6L] = 45
# synapse locations DEND SOMA AXON
def setsynloc ():
for ty1 in xrange(CTYPi):
for ty2 in xrange(CTYPi):
if ice(ty1):
if IsLTS(ty1):
synloc[ty1][ty2]=DEND # distal [GA2] - from LTS
else:
synloc[ty1][ty2]=SOMA # proximal [GA] - from FS
else:
synloc[ty1][ty2]=DEND # E always distal. use AM2,NM2
# setdelmats -- setup delm,deld
def setdelmats ():
for ty1 in xrange(CTYPi):
for ty2 in xrange(CTYPi):
if synloc[ty1][ty2]==DEND and ice(ty2):
# longer delays at dendrites of interneurons since they are currently single compartment
delm[ty1][ty2]=2.0
deld[ty1][ty2]=0.2
else:
delm[ty1][ty2]=2.0
deld[ty1][ty2]=0.2
# weight params
def setwmat ():
for ty1 in xrange(CTYPi):
for ty2 in xrange(CTYPi):
for sy in xrange(STYPi): wmat[ty1][ty2][sy]=0
# E2 -> I weight
wmat[E2][I2][AM2] = wmat[E2][I2L][AM2] = 0.78
wmat[E2][I5][AM2] = 0.11
wmat[E2][I5L][AM2] = 1.01
# E4 -> I weight
wmat[E4][I2][AM2] = wmat[E4][I2L][AM2] = 0.3625
wmat[E4][I5][AM2] = 1.0775
wmat[E4][I5L][AM2] = 0.1225
wmat[E4][I6][AM2] = wmat[E4][I6L][AM2] = 0.4375
# E5R (IT E5a) -> I weight
wmat[E5R][I2][AM2] = wmat[E5R][I2L][AM2] = 0.3625
wmat[E5R][I5][AM2] = 1.0775
wmat[E5R][I5L][AM2] = 0.1225
wmat[E5R][I6][AM2] = wmat[E5R][I6L][AM2] = 0.4375
# E5B (IT E5b) -> I weight
wmat[E5B][I2][AM2] = wmat[E5B][I2L][AM2] = 0.3625
wmat[E5B][I5][AM2] = 1.0775
wmat[E5B][I5L][AM2] = 0.1225
wmat[E5B][I6][AM2] = wmat[E5B][I6L][AM2] = 0.4375
# E5P (PT E5b) -> I weight
wmat[E5P][I2][AM2] = wmat[E5P][I2L][AM2] = 0.3625
wmat[E5P][I5][AM2] = 1.0775
wmat[E5P][I5L][AM2] = 0.1225
wmat[E5P][I6][AM2] = wmat[E5P][I6L][AM2] = 0.4375
# E6 (IT) -> I weight
wmat[E6][I5][AM2] = wmat[E6][I5L][AM2] = 0.24786
wmat[E6][I6][AM2] = wmat[E6][I6L][AM2] = 0.53
# E6C (CT) -> I weight
wmat[E6C][I5][AM2] = wmat[E6C][I5L][AM2] = 0.24786
wmat[E6C][I6][AM2] = wmat[E6C][I6L][AM2] = 0.53
epops = [E2, E4, E5R, E5B, E5P, E6, E6C]
prty = E2; # E2 -> E weight
for poty,wght in zip(epops,[0.6416, 0.3700, 0.5049, 0.4446, 0.4446, 0.0000, 0.0000]): wmat[prty][poty][AM2] = wght
prty = E4 # E4 -> E weight
for poty,wght in zip(epops,[0.7368, 0.6416, 0.6383, 0.8981, 0.8981, 1.9082, 1.9082]): wmat[prty][poty][AM2] = wght
prty = E5R # E5R -> E weight
for poty,wght in zip(epops,[0.5243, 0.4180, 0.5716, 0.8320, 0.8320, 0.3220, 0.3220]): wmat[prty][poty][AM2] = wght
prty = E5B # E5B -> E weight
for poty,wght in zip(epops,[0.2342, 0.1700, 0.3194, 0.6855, 0.6855, 0.4900, 0.4900]): wmat[prty][poty][AM2] = wght
# E5P -> E weight
wmat[E5P][E5P][AM2] = 0.6855
prty = E6 # E6 -> E weight
for poty,wght in zip(epops,[0.0000, 0.0000, 0.1365, 0.3092, 0.3092, 0.5300, 0.5300]): wmat[prty][poty][AM2] = wght
prty = E6C # E6C -> E weight
for poty,wght in zip(epops,[0.0000, 0.0000, 0.1365, 0.3092, 0.3092, 0.5300, 0.5300]): wmat[prty][poty][AM2] = wght
# I -> X
for prty,sy in zip([I2,I2L],[GA,GA2]):
for poty in [E2,I2,I2L]:
wmat[prty][poty][sy] = 1.5
if IsLTS(prty) and not ice(poty): wmat[prty][poty][GB2] = 1.5 * GB2R
for prty,sy in zip([I5,I5L],[GA,GA2]):
for poty in [E4,E5R,E5B,E5P,I5,I5L]:
wmat[prty][poty][sy] = 1.5
if IsLTS(prty) and not ice(poty): wmat[prty][poty][GB2] = 1.5 * GB2R
for prty,sy in zip([I6,I6L],[GA,GA2]):
for poty in [E6,E6C,I6,I6L]:
wmat[prty][poty][sy] = 1.5
if IsLTS(prty) and not ice(poty): wmat[prty][poty][GB2] = 1.5 * GB2R
# gain control
for ty1 in xrange(CTYPi):
for ty2 in xrange(CTYPi):
for sy in xrange(STYPi):
if wmat[ty1][ty2][sy] > 0:
if ice(ty1): # I -> X
if ice(ty2):
if IsLTS(ty1): # LTS -> I
if IsLTS(ty2): # LTS -> LTS
gn = IIGainLTSLTS
else: # LTS -> FS
gn = IIGainLTSFS
else: # FS -> I
if IsLTS(ty2): # FS -> LTS
gn = IIGainFSLTS
else: # FS -> FS
gn = IIGainFSFS
else: # I -> E
gn = IEGain
else: # E -> X
if ice(ty2): # E -> I
if IsLTS(ty2): # E -> LTS
gn = EIGainLTS
else: # E -> FS
gn = EIGainFS
else: # E -> E
gn = EEGain
if sy==AM2:
wmat[ty1][ty2][MG] = wmat[ty1][ty2][AM2] * mGLURR
if verbose: print 'AM2:',wmat[ty1][ty2][AM2],'mGLURR:',wmat[ty1][ty2][MG]
if sy==AM2:
if ice(ty2): # E -> I
wmat[ty1][ty2][NM2] = wmat[ty1][ty2][AM2] * NMAMREI
else: # E -> E
wmat[ty1][ty2][NM2] = wmat[ty1][ty2][AM2] * NMAMREE
wmat[ty1][ty2][sy] *= gn
def setpmat ():
for ii in xrange(CTYPi):
for jj in xrange(CTYPi): pmat[ii][jj]=0
# E2 -> I wiring
pmat[E2][I2] = pmat[E2][I2L] = 0.1871
pmat[E2][I5] = 0.02
pmat[E2][I5L] = 0.217
# E4 -> I wiring
pmat[E4][I2] = pmat[E4][I2L] = 0.0222
pmat[E4][I5] = 0.1906
pmat[E4][I5L] = 0.0349
pmat[E4][I6] = pmat[E4][I6L] = 0.0155
# E5R (IT E5a) -> I wiring
pmat[E5R][I2] = pmat[E5R][I2L] = 0.0222
pmat[E5R][I5] = 0.1906
pmat[E5R][I5L] = 0.0349
pmat[E5R][I6] = pmat[E5R][I6L] = 0.0155
# E5B (IT E5b) -> I wiring
pmat[E5B][I2] = pmat[E5B][I2L] = 0.0222
pmat[E5B][I5] = 0.1906
pmat[E5B][I5L] = 0.0349
pmat[E5B][I6] = pmat[E5B][I6L] = 0.0155
# E5P (PT E5b) -> I wiring
pmat[E5P][I2] = pmat[E5P][I2L] = 0.0222
pmat[E5P][I5] = 0.1906
pmat[E5P][I5L] = 0.0349
pmat[E5P][I6] = pmat[E5P][I6L] = 0.0155
# E6 (IT) -> I wiring
pmat[E6][I5] = pmat[E6][I5L] = 0.0249
pmat[E6][I6] = pmat[E6][I6L] = 0.0234
# E6C (CT) -> I wiring
pmat[E6C][I5] = pmat[E6C][I5L] = 0.0249
pmat[E6C][I6] = pmat[E6C][I6L] = 0.0234
epops = [E2, E4, E5R, E5B, E5P, E6, E6C]
prty = E2; # E2 -> E wiring
for poty,prob in zip(epops,[0.1503, 0.1100, 0.0509, 0.0124, 0.0690, 0,0]): pmat[prty][poty] = prob
prty = E4 # E4 -> E wiring
for poty,prob in zip(epops,[0.0523, 0.1503, 0.0413, 0.0143, 0.0120, 0.0030, 0.0030]): pmat[prty][poty] = prob
prty = E5R # E5R -> E wiring
for poty,prob in zip(epops,[0.0355, 0.0275, 0.1810, 0.0142, 0.0205, 0.0136, 0.0136]): pmat[prty][poty] = prob
prty = E5B # E5B -> E wiring
for poty,prob in zip(epops,[0.0167, 0.0293, 0.0536, 0.1810, 0.0448, 0.0160, 0.0160]): pmat[prty][poty] = prob
pmat[E5P][E5P] = 0.1810 # E5P -> E5P
prty = E6 # E6 -> E wiring
for poty,prob in zip(epops,[0, 0, 0.0334, 0.0271, 0.0277, 0.0282, 0.0234]): pmat[prty][poty] = prob
prty = E6C # E6C -> E wiring
for poty,prob in zip(epops,[0, 0, 0.0334, 0.0271, 0.0277, 0.0234, 0.0282]): pmat[prty][poty] = prob
# I -> X
for prty in [I2,I2L]:
for poty in [E2,I2,I2L]: pmat[prty][poty] = 1.0
for prty in [I5,I5L]:
for poty in [E4,E5R,E5B,E5P,I5,I5L]: pmat[prty][poty] = 1.0
for prty in [I6,I6L]:
for poty in [E6,E6C,I6,I6L]: pmat[prty][poty] = 1.0
for ii in xrange(CTYPi):
for jj in xrange(CTYPi): pmat[ii][jj]*=pmatscale
numc = [0 for i in xrange(CTYPi)]; # number of cells of a type
ix = [0 for i in xrange(CTYPi)]; #starting index of a cell type (into self.ce list)
ixe = [0 for i in xrange(CTYPi)]; #ending index of a cell type
allcells,ecells,icells = 0,0,0
div = zeros( (CTYPi, CTYPi) )
conv = zeros( (CTYPi, CTYPi) )
syty1 = zeros( (CTYPi, CTYPi) ) # stores synapse codes (from labels.py)
syty2 = zeros( (CTYPi, CTYPi) ) # stores synapse code (from labels.py)
syty3 = zeros( (CTYPi, CTYPi) ) # stores synapse code (from labels.py)
sytys1 = {} # dictionary of synapse names
sytys2 = {} # dictionary of synapse names
sytys3 = {} # dictionary of synapse names
SOMA = 0; BDEND = 1; ADEND1 = 2; ADEND2 = 3; ADEND3 = 4;
dsecnames = ['soma','Bdend','Adend1','Adend2','Adend3']
def setdivmat ():
import math
for ty1 in xrange(CTYPi):
for ty2 in xrange(CTYPi):
if pmat[ty1][ty2] > 0.0:
div[ty1][ty2] = math.ceil(pmat[ty1][ty2]*numc[ty2])
conv[ty1][ty2] = int(0.5 + pmat[ty1][ty2]*numc[ty1])
# setup cell-type-to-cell-type synapse-type information
def setsyty ():
for ty1 in xrange(CTYPi): # go thru presynaptic types
for ty2 in xrange(CTYPi): # go thru postsynaptic types
syty1[ty1][ty2] = syty2[ty1][ty2] = syty3[ty1][ty2] = -1 # initialize to invalid
if numc[ty1] <= 0 or numc[ty2] <= 0: continue
if ice(ty1): # is presynaptic type inhibitory?
if IsLTS(ty1): # LTS -> X
syty1[ty1][ty2] = GA2 # code for dendritic gabaa synapse
if ice(ty2): # LTS -> Io
sytys1[(ty1,ty2)] = "GABAss"
else: # LTS -> E
syty2[ty1][ty2] = GB2 # code for denritic gabab synapse
sytys1[(ty1,ty2)] = "GABAs"
sytys2[(ty1,ty2)] = "GABAB"
else: # BAS -> X
syty1[ty1][ty2] = GA # code for somatic gabaa synapse
sytys1[(ty1,ty2)] = "GABAf"
else: # E -> X
syty1[ty1][ty2] = AM2 # code for dendritic ampa synapse
syty2[ty1][ty2] = NM2 # code for dendritic nmda synapse
if ice(ty2): # E -> I
sytys1[(ty1,ty2)] = "AMPA"
sytys2[(ty1,ty2)] = "NMDA"
else: # E -> E
sytys1[(ty1,ty2)] = "AMPA"
sytys2[(ty1,ty2)] = "NMDA"
sytys3[(ty1,ty2)] = "mGLUR"
syty3[ty1][ty2] = MG # use MG -- for mGluR
lctyID,lctyClass = [],[]
# setup some convenient data structures
def setix (scale):
import math
global allcells,ecells,icells
for i in xrange(CTYPi):
numc[i] = int(math.ceil(cpernet[i]*scale))
if numc[i] > 0:
ty = PyrAdr
if ice(i):
if IsLTS(i): ty = LTS
else: ty = FS
for j in xrange(numc[i]):
lctyClass.append(ty)
lctyID.append(i)
allcells += numc[i]
if ice(i): icells += numc[i]
else: ecells += numc[i]
sidx = 0
for i in xrange(CTYPi):
if numc[i] > 0:
ix[i] = sidx
ixe[i] = ix[i] + numc[i] - 1
sidx = ixe[i] + 1
setdivmat()
setsyty()
# setcellpos([pseed,network diameter in microns])
def setcellpos (pseed=4321,cside=colside):
rdm=Random(); rdm.ACG(pseed)
cellsnq = NQS("id","ty","ice","xloc","yloc","zloc")
cellsnq.clear(allcells) # alloc space
lX,lY,lZ=[],[],[]
ldy = {}
ldy[E2] = (160, 420)
ldy[E4] = (420, 570)
ldy[E5R] = (570, 700)
ldy[E5B] = (700, 1040)
ldy[E5P] = (700, 1040)
ldy[E6] = (1040, 1350)
ldy[E6C] = (1040, 1350)
ldy[I2] = (160, 420)
ldy[I2L] = (160, 420)
ldy[I5] = (420, 1040)
ldy[I5L] = (420, 1040)
ldy[I6] = (1040, 1350)
ldy[I6L] = (1040, 1350)
for i in xrange(allcells):
ctyp = lctyID[i]
[x,y,z] = [rdm.uniform(0,cside), rdm.uniform(ldy[ctyp][0],ldy[ctyp][1]), rdm.uniform(0,cside)]
cellsnq.append(i,ctyp,ice(ctyp),x,y,z); lX.append(x); lY.append(y); lZ.append(z);
return cellsnq,lX,lY,lZ
setcpernet() # setup number of cells per network
setwmatex() # setup matrices of external inputs
setsynloc() # setup synapse location matrices
setdelmats() # setup delay matrices
setwmat() # setup weight matrix
setpmat() # setup connectivity matrix
setix(scale)
cellsnq,lX,lY,lZ=setcellpos()
ce = [] # cells on the host
gidvec = [] # gids of cells on the host
lncrec,ltimevec,lidvec=[],[],[] # spike recorders
dlids = {} # map from gid back to ce index
# create the cells
pcID = int(pc.id());
maxcells=0
cperhost = int(allcells/nhosts)
maxcells = cperhost
extra = allcells - cperhost*nhosts
if extra > 0: # check if any remainder cells
if pcID < extra: # first hosts get extra cell
maxcells += 1 # assign an extra cell if any remainders
gid = pcID * (cperhost + 1)
else: # rest of hosts do not
gid = extra*(cperhost+1) + (pcID-extra) * cperhost
else: # even division? all hosts get equal cells
gid = pcID * cperhost
for i in xrange(maxcells):
ct = lctyID[gid]
cell = lctyClass[gid](0+i*50,0,0,gid,ct)
cell.x,cell.y,cell.z = lX[gid],lY[gid],lZ[gid]
dlids[gid] = len(ce) # map from gid back to ce index
ce.append(cell)
gidvec.append(gid)
pc.set_gid2node(gid,pcID)
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,gid)
ncrec.threshold = spiketh # 10 mV is default, lower it for FS cells
ltimevec.append(timevec); lidvec.append(idvec); lncrec.append(ncrec)
pc.cell(gid,lncrec[-1],1) # 1 as 3rd arg means this cell can be source for events
gid += 1
print(' Number of cells on node %i: %i' % (pcID,len(ce)))
pc.barrier()
# wire the network using a NQS table
nccl = []
def wirenq (cnq):
global nccl
nccl = [] # NetCon list for connections between cells
cnq.tog("DB")
vid1,vid2,vwt1,vwt2,vdel,vsec=cnq.getcol("id1"),cnq.getcol("id2"),cnq.getcol("wt1"),cnq.getcol("wt2"),cnq.getcol("del"),cnq.getcol("sec")
vwt3 = cnq.getcol("wt3")
for i in xrange(int(cnq.v[0].size())):
prid = int(vid1[i])
poid = int(vid2[i])
if not pc.gid_exists(poid): continue # only make the connection on a node that has the target
ty1 = lctyID[prid]
ty2 = lctyID[poid]
sname = dsecnames[int(vsec[i])] # which section is the synapse on?
syn = sname + sytys1[(ty1,ty2)]
wt1 = vwt1[i]
delay = vdel[i]
targ = ce[dlids[poid]]
nc1 = pc.gid_connect(prid, targ.__dict__[syn].syn)
nc1.delay = delay; nc1.weight[0] = wt1; nc1.threshold = spiketh; nccl.append(nc1)
wt2 = vwt2[i]
if wt2 > 0: # two synapses? (i.e., AMPA and NMDA)
syn = sname + sytys2[(ty1,ty2)]
if syn in targ.__dict__: # since GABAB not in Bdend
nc2 = pc.gid_connect(prid, targ.__dict__[syn].syn)
nc2.delay = delay; nc2.weight[0] = wt2; nc2.threshold = spiketh; nccl.append(nc2)
wt3 = vwt3[i]
if wt3 > 0: # three synapses? (i.e., AMPA and NMDA and mGLUR)
if verbose: print 'mGLUR synapse wt3 > 0:',wt3
syn = sname + sytys3[(ty1,ty2)]
if syn in targ.__dict__: # make sure target has this synapse (needed since Bdend does not have mGLUR)
nc3 = pc.gid_connect(prid, targ.__dict__[syn].syn)
nc3.delay = delay; nc3.weight[0] = wt3; nc3.threshold = spiketh; nccl.append(nc3)
#
def picksec (prty, poty, rdm):
if ice(poty): return SOMA
if ice(prty): # I -> E
if IsLTS(prty): # LTS -> E
if rdmsec: return rdm.discunif(ADEND1,ADEND3)
else: return ADEND3
else:
return SOMA
else: # E -> E
if rdmsec: return rdm.discunif(BDEND,ADEND3)
else: return ADEND3
# swire - spatial wiring: wires the network using pmat and cell positions
# (wiring probability affected by distance btwn cells)
# slambda (global) specifies length-constant for spatially-dependent fall-off in wiring probability
# slambda is only used for I->X wiring; for E->X wiring uses fixed probability
def swire (wseed):
global slambda
from math import sqrt,exp
[vidx,vdel,vtmp,vwt1,vwt2,vwt3,vprob] = [Vector() for x in xrange(7)]
z = 0
if slambda <= 0:
print "swire WARN: invalid slambda=", slambda, "setting slambda to ", colside/3
slambda=colside/3
slambdasq = slambda**2 # using squared distance
h.vrsz(1e4,vidx,vdel,vtmp)
rdm=Random(); rdm.ACG(wseed) #initialize random # generator
rdm.uniform(0,1)
vprob.resize(allcells**2); vprob.setrand(rdm)
pdx=0 # index into vprob
connsnq=NQS("id1","id2","del","wt1","wt2","wt3","sec")
connsnq.clear(1e3*allcells)
for prid in xrange(allcells):
vrsz(0,vidx,vdel,vwt1,vwt2,vwt3)
prty=lctyID[prid]
ic1=ice(prty)
for poty in xrange(0,CTYPi):
if numc[poty] > 0 and pmat[prty][poty]>0:
pbase = pmat[prty][poty]
for poid in xrange(ix[poty],ixe[poty]+1): # go thru postsynaptic cells
if prid==poid: continue # no self-connects
ic2=ice(lctyID[poid])
dx = lX[prid] - lX[poid]
dy = lY[prid] - lY[poid]
dz = lZ[prid] - lZ[poid]
ds = sqrt(dx**2 + dy**2 + dz**2) # Connectivity fall-off depends on 3D distance
if ic1: prob = exp(-ds/slambda) # probability of connect falls off only for I->X
else: prob = pbase
if prob >= vprob[pdx]: # rdm.uniform(0,1)
mindelay = delm[prty][poty]-deld[prty][poty]
maxdelay = delm[prty][poty]+deld[prty][poty]
delay=rdm.uniform(mindelay,maxdelay) # synaptic delay
delay += ds/axonalvelocity # add axonal delay
vidx.append(poid); vdel.append(delay)
if syty1[prty][poty]>=0: vwt1.append(wmat[prty][poty][int(syty1[prty][poty])])
else: vwt1.append(0)
if syty2[prty][poty]>=0: vwt2.append(wmat[prty][poty][int(syty2[prty][poty])])
else: vwt2.append(0)
if syty3[prty][poty]>=0: vwt3.append(wmat[prty][poty][int(syty3[prty][poty])])
else: vwt3.append(0)
pdx += 1
for ii in xrange(int(vidx.size())): connsnq.append(prid,vidx[ii],vdel[ii],vwt1[ii],vwt2[ii],vwt3[ii],picksec(prty , lctyID[int(vidx[ii])], rdm))
wirenq(connsnq) # do the actual wiring based on self.connsnq
return connsnq
connsnq=swire(WSEED)
pc.barrier() # wait for wiring to get completed
# setup rxd for E cells
# get list of all Sections associated with an excitatory cell
def getesec ():
esec = []
for cell in ce:
if ice(cell.ty): continue
for s in cell.all_sec: esec.append(s)
return esec
def pcidpr (s):
global pcID
print 'host',pcID,':',s
### RXD ###
[cyt,er,cyt_er_membrane,ca,caextrude,serca,leak,CB,caCB,buffering]=[None for i in xrange(10)]
rxdsec=getesec() # Section list for use with rxd
pc.barrier()
if len(rxdsec) > 0: # only use rxd if there are viable Sections
from neuron import rxd
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)
### RYR - based on Sneyd et al, 2003
# constants
k_a_pos = 1500000000000.0 # mM^-4/ms
k_a_neg = 0.0288 # /ms
k_b_pos = 1500000000.0 # mM^-3/ms
k_b_neg = 0.3859 # /ms
k_c_pos = 0.00175 # /ms
k_c_neg = 0.0001 # /ms
v1ryr = dconf['v1ryr'] # /ms
Ka_4 = k_a_neg / k_a_pos # Ka**4
Kb_3 = k_b_neg / k_b_pos # Kb**3
Kc = k_c_neg / k_c_pos
# w_state is fraction of RYR not in C2 state (closed state), ie fraction of RYR that is open
# w_infinity - equ: 29
c3ryr = (ca[cyt]**3)/Kb_3;
c4ryr = Ka_4/(ca[cyt]**4);
w_inf = (1.0 + c4ryr + c3ryr) / (1.0+(1.0/Kc)+ c4ryr + c3ryr)
w_state = rxd.State(cyt_er_membrane, initial=0.9999) # 0)#0.9999) # check if initial == 0???? or can put it as w_inf
# equ: 8 (which is the same as equ 22)
w_rate = rxd.Rate(w_state, 1.0 - w_state[cyt_er_membrane] / w_inf)
# P_ryr - gating variable - equ: 7 (which is the same as 27)
# (open probability)
ryr_gate = w_state[cyt_er_membrane] * (1.0 + c3ryr) / (1.0 + c4ryr + c3ryr)
# the following is extracted from equ 9 and 15
k_ryr = v1ryr*ryr_gate
ryr = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], k_ryr, k_ryr, membrane=cyt_er_membrane)
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, recdt )
c.Adend3_ip3cyt = Vector(tstop/h.dt)
c.Adend3_ip3cyt.record( ip3[cyt].nodes(c.Adend3)(0.5)[0]._ref_concentration, recdt )
# 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, recdt )
c.Adend3_caer = Vector(tstop/h.dt)
c.Adend3_caer.record( ca[er].nodes(c.Adend3)(0.5)[0]._ref_concentration, recdt )
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, recdt )
c.Adend3_caCB = Vector(tstop/h.dt)
c.Adend3_caCB.record( caCB.nodes(c.Adend3)(0.5)[0]._ref_concentration, recdt )
setreccaer() # NB: only record from RXD variables after ALL species setup!
setreccacb() # otherwise, the pointers get messed up.
setrecip3()
setmGLURflux()
# setup inputs - first noise inputs
def getsyns ():
syns = {} # mapping of synapse names, first index is ice, second is synapse code
syns[ (0,MG) ] = ["Adend3mGLUR","Adend2mGLUR","Adend1mGLUR"]
syns[ (0,AM2) ] = ["Adend3AMPA","Adend2AMPA","Adend1AMPA","BdendAMPA"]
syns[ (1,AM2) ] = "somaAMPA"
syns[ (0,NM2) ] = ["Adend3NMDA","Adend2NMDA","Adend1NMDA","BdendNMDA"]
syns[ (1,NM2) ] = "somaNMDA"
syns[ (0,GB2) ] = ["Adend3GABAB","Adend2GABAB","Adend1GABAB"]
syns[ (0,GA2) ] = ["Adend3GABAs","Adend2GABAs","Adend1GABAs"]
syns[ (1,GA2) ] = "somaGABAss"
syns[ (0,GA) ] = "somaGABAf"
syns[ (1,GA) ] = "somaGABAf"
return syns
dsstr = ['AMPA', 'NMDA', 'GABAS', 'mGLUR', 'GABAB']
# adds synapses across dendritic fields for the E cells
def addsyns ():
for cell in ce:
cell.dsy = {}; cell.vsy = {}
if ice(cell.ty): continue
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
#print 'loc:',loc
lsy.append(Synapse(sect=sec,loc=loc,tau1=0.05,tau2=5.3,e=0)); lvsy.append(h.Vector())#AMPA
lsy.append(SynapseNMDA(sect=sec,loc=loc,tau1NMDA=15,tau2NMDA=150,r=1,e=0)) # NMDA
lvsy.append(h.Vector())
lsy.append(Synapse(sect=sec,loc=loc,tau1=0.2,tau2=20,e=-80)) # GABAA_slow
lvsy.append(h.Vector())
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(h.Vector())
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(h.Vector())
lsy.append(Synapse(sect=sec,loc=loc,tau1=0.05,tau2=5.3,e=0) ) # AMPA
lvsy.append(h.Vector())
lsy.append(SynapseNMDA(sect=sec,loc=loc,tau1NMDA=15,tau2NMDA=150,r=1,e=0)) # NMDA
lvsy.append(h.Vector())
llsy.append(lsy); llvsy.append(lvsy);
cell.dsy[sec] = llsy; cell.vsy[sec] = llvsy;
addsyns()
#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,100)
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
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 = tstart # 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)
def makeNoiseNetStims (simdur,rdmseed):
nsl = [] #NetStim List
ncl = [] #NetCon List
nrl = [] #Random List for NetStims
nrlsead = [] #List of seeds for NetStim randoms
syns = getsyns() ;
for cell in ce: # go through cell types, check weights,rates of inputs
ct = cell.ty # get cell type code
if ice(ct): # only has 1 compartment
for sy in xrange(STYPi):
if wmatex[ct][sy] <= 0.0 or ratex[ct][sy] <= 0: continue
syn = syns[(ice(ct),sy)]
if type(syn) == list:
for idx,SYN in enumerate(syn):
makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,SYN,wmatex[ct][sy],1e3/ratex[ct][sy],simdur,rdmseed*(cell.ID+1)*(idx+1))
else:
makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,syn,wmatex[ct][sy],1e3/ratex[ct][sy],simdur,rdmseed*(cell.ID+1))
else: # E cells - need to distribute noise over all sections
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,MG,GB2]):
if ratex[ct][sy] > 0. and wmatex[ct][sy] > 0.:
makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,lsy[i].syn,wmatex[ct][sy],(1e3/ratex[ct][sy]),simdur,rdmseed*(cell.ID+1)*(i+1));
sec = cell.Bdend; llsy = cell.dsy[sec];
for lsy in llsy:
for i,sy in enumerate([AM2,NM2,GA2]):
if ratex[ct][sy] > 0. and wmatex[ct][sy] > 0.:
makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,lsy[i].syn,wmatex[ct][sy],(1e3/ratex[ct][sy]),simdur,rdmseed*(cell.ID+1)*(i+4));
sec = cell.soma; llsy = cell.dsy[sec];
for i,sy in enumerate([GA,AM,NM]):
if ratex[ct][sy] > 0. and wmatex[ct][sy] > 0.:
for lsy in llsy:
makeNoiseNetStim(cell,nsl,ncl,nrl,nrlsead,lsy[i].syn,wmatex[ct][sy],(1e3/ratex[ct][sy]),simdur,rdmseed*(cell.ID+1)*(i+7)); rdmseed+=1
return nsl,ncl,nrl,nrlsead
nsl,ncl,nrl,nrlsead = makeNoiseNetStims(tstart+tstop,ISEED)
pc.barrier() # wait for completion of NetStim creation
#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)
# handler for printing out time during simulation run
def fi():
for i in xrange(int(tstart),int(tstart+tstop),100): h.cvode.event(i, "print " + str(i))
if pc.id() == 0: fih = h.FInitializeHandler(1, fi)
vt=Vector(); vt.record(h._ref_t); # record time
pc.barrier() # wait for NetStims to get setup
####################################################################################
### simulation run here
def myrun ():
pc.set_maxstep(10)
dastart,daend=None,None
if pc.id()==0:
dastart = datetime.now()
print 'started at:',dastart
h.stdinit()
if len(rxdsec)>0: # any sections with rxd?
ca[er].concentration = dconf['caerinit'] # 100e-6
ca[cyt].concentration = dconf['cacytinit'] # 100e-6
pc.psolve(h.t+tstop) # run for tstop
pc.barrier() # Wait for all hosts to get to this point
if pc.id()==0:
daend = datetime.now()
print 'finished ',tstop,' ms sim at:',daend
dadiff = daend - dastart;
print 'runtime:',dadiff, '(',tstop/1e3,' s)'
if dconf['dorun']: myrun()
# concatenate the results so can view/save all at once
lspks,lids=array([]),array([])
for host in xrange(nhosts): # is this loop required? can't just post messages from given host?
if host == pc.id():
for i in xrange(len(ltimevec)):
lspks=concatenate((lspks,ltimevec[i]))
lids=concatenate((lids,lidvec[i]))
# save data - output path based on simstr and pcid
def savedata (simstr,pcid):
safemkdir(outdir)
fn = outdir + '/' + simstr + '_pc_' + str(pcid) + '.npz'
print 'host ' , pcid, ' saving to ' , fn
ne,ni,szfast = 0,0,0
lE,lI=[],[]
for c in ce:
if ice(c.ty):
lI.append(c.ID)
ni += 1
else:
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_voltI = zeros((ni,szfast));
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()
cdx += 1
numpy.savez(fn,lctyID=array(lctyID),lX=array(lX),lY=array(lY),lZ=array(lZ),vt=vt.as_numpy(),lspks=lspks,lids=lids,lE=lE,lI=lI,Adend3_volt=Adend3_volt,Bdend_volt=Bdend_volt)
pc.barrier()
####################################################################################
if saveout: # save the sim data
if pcID == 0: print 'saving data'
savedata(simstr,pcID)
if saveconns and pcID == 0: # save connectivity data
print 'saving connections'
h.batch_flag = 1 # in case file already exists
connsnq.tog("DB")
connsnq.sv(outdir + '/' + simstr + 'connsnq.nqs')
pc.runworker()
pc.done()
if nhosts > 1: h.quit() # this means was likely running in batch mode
