import sys print sys.path import os import string from neuron import * from datetime import datetime 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.0" 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") CTYPi = 60.0; STYPi = 18.0 from pyinit import * from neuron import h, gui from vector import * from nqs import * from labels import * import random from pylab import * import ConfigParser tl = tight_layout ion() rcParams['lines.markersize'] = 15 rcParams['lines.linewidth'] = 4 rcParams['font.size'] = 25 defbinsz=10 stimdel = 500 # delay from stim to measured firing rate peak response defnCPU=8 print 'Using ' , defnCPU, ' cpus by default.' config = ConfigParser.ConfigParser() # determine config file name def setfcfg (): fcfg = "physiol.cfg" # default config file name for i in xrange(len(sys.argv)): if sys.argv[i].endswith(".cfg") and os.path.exists(sys.argv[i]): fcfg = sys.argv[i] print "config file is " , fcfg return fcfg fcfg=setfcfg() # config file name config.read(fcfg) def conffloat (base,var): return float(config.get(base,var)) def confint (base,var): return int(config.get(base,var)) def confstr (base,var): return config.get(base,var) tstart = 0 tstop = conffloat('run','tstop') + tstart binsz = conffloat("run","binsz") # bin size (in milliseconds) for MUA recdt = conffloat('run','recdt') recvdt = conffloat('run','recvdt') vtslow=Vector(); vtslow.indgen(tstart,tstop-recdt,recdt); vtslow=numpy.array(vtslow.to_python()) vtfast=Vector(); vtfast.indgen(tstart,tstop-recvdt,recvdt); vtfast=numpy.array(vtfast.to_python()) simstr = config.get('run','simstr') def makefname (simstr,pcid): return 'data/' + simstr + '_pc_' + str(pcid) + '.npz' ix,ixe=[1e9 for i in xrange(CTYPi)],[-1e9 for i in xrange(CTYPi)] numc=[0 for i in xrange(CTYPi)] allcells,icells,ecells=0,0,0 # reload the variables from config file def reloadvars (): global tstart,loadstate,tstop,binsz,recdt,recvdt,vtslow,vtfast,simstr,numc,allcells,ecells,icells,startt tstart = conffloat('run','loadtstop') loadstate = confint('run','loadstate') if loadstate == 0: tstart = 0 # only use previous end time if loading state tstop = conffloat('run','tstop') + tstart binsz = conffloat("run","binsz") # bin size (in milliseconds) for MUA recdt = conffloat('run','recdt') recvdt = conffloat('run','recvdt') vtslow=Vector(); vtslow.indgen(tstart,tstop-recdt,recdt); vtslow=numpy.array(vtslow.to_python()) vtfast=Vector(); vtfast.indgen(tstart,tstop-recvdt,recvdt); vtfast=numpy.array(vtfast.to_python()) simstr = config.get('run','simstr') numc=[0 for i in xrange(CTYPi)] allcells,ecells,icells=0,0,0 # setup start/end indices for cell types def makeix (lctyID): global ix,ixe,allcells,ecells,icells allcells,ecells,icells=0,0,0 for i in xrange(CTYPi): ix[i]=1e9 ixe[i]=-1e9 numc[i]=0 for i in xrange(len(lctyID)): ty = lctyID[i] numc[ty]+=1 allcells+=1 ix[ty] = min(ix[ty],i) ixe[ty] = max(ixe[ty],i) if h.ice(ty): icells+=1 else: ecells+=1 lfastvar = ['volt', 'iAM', 'iNM', 'iGB', 'iGA','ina','ik','ica','ih'] # reads data from files saved by mpisim.py - must get nhost correct def readdat (simstr,nhost): ld = {} # dict for data from a single host for pcid in xrange(nhost): fn = makefname(simstr,pcid) ld[pcid]=numpy.load(fn) return ld # load/concat data from mpisim.py (which saves data from each host separately) - must get nhost correct def loaddat (simstr,nhost,quiet=False): ld = readdat(simstr,nhost) # dict for data from a single host lids,lspks=array([]),array([]) # stitch together spike times for pcid in xrange(nhost): lids=concatenate((lids,ld[pcid]['lids'])) lspks=concatenate((lspks,ld[pcid]['lspks'])) ldout={} # concatenated output from diff hosts ldout['lids']=lids; ldout['lspks']=lspks; for k in ['lctyID','vt','lX', 'lY', 'lZ']: ldout[k]=ld[0][k][:] makeix(ldout['lctyID']) ncells = len(ld[0]['lctyID']) vt = ld[0]['vt'][:] locvarlist = [] for loc in ['soma','Adend3']: for var in ['volt']: locvarlist.append(loc + '_' + var) locvarlist.append('Bdend_volt') for k in locvarlist: if not k in ld[0]: continue # skip? if not quiet: print k if k.endswith('volt') or lfastvar.count(k.split('_')[1])>0: sz = len(vtfast) else: sz = len(vtslow) ldout[k] = zeros((ncells,sz)) ldk = ldout[k] for pcid in xrange(nhost): d = ld[pcid][k] lE = ld[pcid]['lE']; for row,cdx in enumerate(lE): ldk[cdx,:] = d[row,0:sz] if k == 'soma_volt': # special case for soma voltage of I cells lI=ld[pcid]['lI']; dI=ld[pcid]['soma_voltI']; row=0 for row,cdx in enumerate(lI): ldk[cdx,:] = dI[row,0:sz] for pcid in xrange(nhost): del ld[pcid].f ld[pcid].close() del ld[pcid] return ldout # colors for spikes for raster -- based on cell type def getcolors (lspks,lids,lctyID): lclrs = [] for i in xrange(len(lids)): ty = lctyID[int(lids[i])] if IsLTS(ty): lclrs.append('b') elif ice(ty): lclrs.append('g') else: lclrs.append('r') return lclrs # def drawraster (ld,sz=2): lspks,lids,lctyID=ld['lspks']/1e3,ld['lids'],ld['lctyID'] try: lclrs = ld['lclrs'] except: lclrs = getcolors(lspks,lids,lctyID) ld['lclrs'] = lclrs scatter(lspks,lids,s=sz**2,c=lclrs,marker='+') allcells = len(lctyID) xlim((tstart/1e3,tstop/1e3)); ylim((0,allcells)); tight_layout(); xlabel('Time (s)',fontsize=30); ax = gca() ax.set_yticks([]) # get cells in vid if they're of specified type (in lty) def getall (ix,ixe,lty=[E2,E4,E5R,E5B,E5P,E6,E6C]): vact = Vector() for ct in lty: vtmp = Vector(); vtmp.indgen(ix[ct],ixe[ct],1); vact.append(vtmp) return vact # get an array of times of interest def gettimes (tstart,tstop,binsz): vt = h.Vector() vt.indgen(tstart,tstop,binsz) return vt # get an NQS with spikes (snq) def getsnq (ld): lspks,lids,lctyID=ld['lspks'],ld['lids'],ld['lctyID'] 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])])) ld['snq']=snq return snq # print spike rates in each population in the periods of interest (baseline, signal, recall) def prspkcount (snq,tstart,tstop,binsz,times=None,quiet=False): snq.verbose=0 if times is None: times = gettimes(tstart,tstop,binsz) # periods of interest lfa=[] vact = getall(ix,ixe,lty=[E2,E4,E5R,E5B,E5P,E6,E6C]) for startt in times: endt = startt + binsz if endt > tstop: break fa = round(1e3*snq.select("id","EQW",vact,"t","[]",startt,endt) / ((endt-startt)*vact.size()),2) lfa.append(fa) if not quiet: stro = "t=" + str(startt) + "-" + str(endt) + ". E:" + str(fa) + " Hz." for ct in [I2, I2L, I5, I5L, I6, I6L]: fb = round(1e3*snq.select("id","[]",ix[ct],ixe[ct],"t","[]",startt,endt) / ((endt-startt)*numc[ct]),2) stro += ' ' + CTYP[ct] + ':' + str(fb) fi = round(1e3*snq.select("ice",1,"t","[]",startt,endt) / ((endt-startt)*icells),2) stro += ' I:' + str(fi) for ct in [E2, E4, E5R, E5B, E5P, E6, E6C]: fb = round(1e3*snq.select("id","[]",ix[ct],ixe[ct],"t","[]",startt,endt) / ((endt-startt)*numc[ct]),2) stro += ' ' + CTYP[ct] + ':' + str(fb) fe = round(1e3*snq.select("ice",0,"t","[]",startt,endt) / ((endt-startt)*ecells),2) stro += ' E:' + str(fe) print stro snq.verbose=1 return lfa # getfnq - make an NQS with ids, firing rates, types def getfnq(snq,lctyID,skipms=200): snq.verbose=0; snq.tog("DB"); fnq = h.NQS("id","freq","ty") tf = tstop - skipms # duration we're considering for frequency calc for i in xrange(allcells): n = float( snq.select("t",">",skipms,"id",i) ) # number of spikes fnq.append(i, n*1e3/tf, lctyID[i]) snq.verbose=1 return fnq # pravgrates - print average firing rates using fnq def pravgrates(snq,skipms=500): snq.verbose=0 for ty in xrange(CTYPi): if numc[ty] < 1: continue print CTYP[ty], ' avg rate = ', getrate(snq,ty,500.0)[2] snq.verbose=1 ###########################################################3 # reload variables, reset config file name def myreload (Simstr,ncpu=defnCPU, quiet=False): global fcfg,simstr if not quiet: print 'simstr was ' , simstr simstr = Simstr print 'simstr: ' , simstr fcfg = 'backupcfg/' + simstr + '.cfg' if len(config.read(fcfg))==0: print 'myreload ERRA: could not read ', fcfg, '!!!' return None reloadvars() if not quiet: print fcfg try: ld=loaddat(simstr,ncpu,quiet); return ld except: print 'could not reload data from' , Simstr return None ld=None # def mydrawrast (lld): for i,ld in enumerate(lld): figure(); drawraster(ld,sz=4); title(lsimstr[i]) xlabel('') xlim((0.5,2)) fsz= 20 ax = gca(); ax.set_yticks([]); ax.set_xticks([]); ax.set_xlabel(''); tx=-0.03 lty = [(ix[ct]+(ixe[ct]-ix[ct])*0.3)/allcells for ct in [E2,E4,E5R,E5B,E5P,E6C,E6]] for ty,lbl in zip(lty,['E2','E4','E5a','E5b','E5P','E6C','E6']): addtext(1,1,[1],[lbl],tx=tx,ty=ty,c='r',fsz=fsz+5) lty = [(ix[ct]+(ixe[ct]-ix[ct])*0.075)/allcells for ct in [I2,I5,I6]] for ty,lbl in zip(lty,['I2','I5','I6']): addtext(1,1,[1],[lbl],tx=tx,ty=ty,c='g',fsz=fsz) lty = [(ix[ct]+(ixe[ct]-ix[ct])*0.25)/allcells for ct in [I2L,I5L,I6L]] for ty,lbl in zip(lty,['I2L','I5L','I6L']): addtext(1,1,[1],[lbl],tx=tx,ty=ty,c='b',fsz=fsz) # def addtext (row,col,lgn,ltxt,tx=-0.025,ty=1.03,c='k',fsz=30): for gn,txt in zip(lgn,ltxt): if row == col and col == 1: ax = gca() else: ax = subplot(row,col,gn) text(tx,ty,txt,fontweight='bold',transform=ax.transAxes,fontsize=fsz,color=c); # def Gaddtext (G,row,col,txt,tx=-0.025,ty=1.03,c='k',fsz=30): ax = subplot(G[row,col]) text(tx,ty,txt,fontweight='bold',transform=ax.transAxes,fontsize=fsz,color=c); def naxbin (ax,nb): ax.locator_params(nbins=nb); lsimstr = ['d2', 'dystonia', 'latchup'] lld,lsnq = [],[] for simstr in lsimstr: ld = loaddat(simstr,defnCPU) lld.append(ld) snq = getsnq(ld) lsnq.append(snq) mydrawrast(lld)