PIR gamma oscillations in network of resonators (Tikidji-Hamburyan et al. 2015)

 Download zip file 
Help downloading and running models
Accession:183718
" ... The coupled oscillator model implemented with Wang–Buzsaki model neurons is not sufficiently robust to heterogeneity in excitatory drive, and therefore intrinsic frequency, to account for in vitro models of ING. Similarly, in a tightly synchronized regime, the stochastic population oscillator model is often characterized by sparse firing, whereas interneurons both in vivo and in vitro do not fire sparsely during gamma,but rather on average every other cycle. We substituted so-called resonator neural models, which exhibit class 2 excitability and postinhibitory rebound (PIR), for the integrators that are typically used. This results in much greater robustness to heterogeneity that actually increases as the average participation in spikes per cycle approximates physiological levels. Moreover, dynamic clamp experiments that show autapse-induced firing in entorhinal cortical interneurons support the idea that PIR can serve as a network gamma mechanism. ..."
Reference:
1 . Tikidji-Hamburyan RA, Martínez JJ, White JA, Canavier CC (2015) Resonant Interneurons Can Increase Robustness of Gamma Oscillations. J Neurosci 35:15682-95 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Entorhinal cortex;
Cell Type(s): Wide dynamic range neuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; Python;
Model Concept(s): Gamma oscillations;
Implementer(s): Tikidji-Hamburyan, Ruben [ruben.tikidji.hamburyan at gmail.com] ;
# -*- coding: utf-8 -*-
"""
/***********************************************************************************************************\

 This NEURON + Python script associated with paper:                                                        
 Ruben A. Tikidji-Hamburyan, Joan José Martínez, John A. White, Carmen C. Canavier                         
    Resonant interneurons can increase robustness of gamma oscillations                                    
                                                                                                           
 Network of 300 Izhikevich's resonators connected by double-exponential synapses is modeled                
 All parameters may be set up by command line arguments listed AFTER script name. The command should be:   
  nrngui -nogui -python network.py [PARAMETERS]                                                            
 A list of avalible parameters can be printed out by command:                                              
  nrngui -nogui -python network.py --help
                                                                                                           
 To replicate Figure 4 run:                                                                                
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=3e-7 -nstart='(900.,0.21e-5,1000)' -Istd=0 -gui -view 
                                                                                                           
 To replicate points from Figure 8A run:                                                                   
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=3e-7 -Istd=Xe-5 -gui=off
    X is std dev in nA                                                                                     
                                                                                                           
 To replicate points from Figure 8B run:                                                                   
  nrngui -nogui -python network.py -Iapp=0.2e-5,X -gsyn=3e-7 -Istd=0.0 -gui=off
    X is CV multiple on 0.2e-5, mean value                                                                 
                                                                                                           
 To replicate points from Figure 8C run:                                                                   
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=3e-7 -Istd=0.0 -gsynscale=1,X -gui=off
    X is CV of total synaptic conductance distributed within population.                                   
                                                                                                           
 To replicate points from Figure 8D run:                                                                   
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=3e-7 -Istd=0.0 -delay=3,X -gui=off
    X is delay standard deviation in ms                                                                    
                                                                                                           
 === NOTES FOR ALL Figures 8 simulations ===                                                               
 ===  1. the results of simulations is saved into network.csv file.                                        
 ===  2. to see traces and rasterplot remove -gui=off and add -gui -view for any command above              
 ===  3. when SPC felt below 0.15 the network is very sensitive to initial conditions and result may be    
         slightly different from the paper, specifically for large Iapp CVs.                                   
                                                                                                           
 To replicate points on Figure 7:
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=XXe-5 -tau2=YY -Istd=0 -gui -view 
    XX and YY is gsyn in nS and tau_fall in ms                                                             
 
 To replicate Figure 9C3:                                                                                                          
  nrngui -nogui -python network.py -Istd=0 -Iapp=0.2e-5 -gsyn=0.005e-5 -delay=1 -view -gui -tstop=5000 -F=1,0.1125 -c=299 -Vinit=-65 -sort=F

 To replicate Figure 10C1:
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=0.05e-5 -Istd=26.1e-7 -Vinit=-65 -c=299 -gui -view -tstop=5000
 
 To replicate Figure 10C2:
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=0.03e-5 -Istd=3.915e-5 -Vinit=-65 -c=299 -gui -view -tstop=5000

 To replicate Figure 10C3:
  nrngui -nogui -python network.py -Iapp=0.15e-5 -gsyn=0.03e-5 -Istd=10.455e-5 -Vinit=-65 -c=299 -gui -view -tstop=5000

 Copyright: Ruben Tikidji-Hamburyan <rth@nisms.krinc.ru> Feb.2013 - Aug.2015                               

\************************************************************************************************************/
"""

import numpy as np
import sys,os,csv,threading
import random as rnd
from neuron import h

###### Paramters:
ncell		= 300			# number of neurons in population
ncon		= 40#(20,60)#40	# number of input connections per neuron
#ntype		= "TypeI"		# use: 'TypeI' for Class 1 neuron
#ntype		= "RS-moca"			# 	'RS' - resonator
ntype		= "RS"			# 	'RS' - resonator
#ntype		= "FS"			# 	'FS' - fast spiking
methods		= {
	'R2'		: True,
	'maxFreq'	: 200.0,		# max frequency
	'peakDetec' : True,			# Turn on/off peak detector
	'gkernel'	: (10.0,50.0),	# Kernel and size (5,25),#
	"netFFT"	: False,		# Turn on/off network FFT
	"nrnFFT"	: False,		# Turn on/off neuron FFT
	'netISI'	: 30001,		# max net ISI
	'nrnISI'	: 30001,		# max neuron ISI
	'cliptrn'	: False,#1000,#False alse,#500,#False,	# Clip transience for first n ms or False
	'traceView'	: False,		# Trace and nulcline view
	'traceWidth': 55.0,			#
	'tracetail'	: 'conductance',#'total current',#'firing rate',#'conductance', #'current' 'total current',
	'save'		: False,		# To save all data from model. Don't turn on!
	'patview'	: True,			# Turn on/off Pattern vew
	'gui'		: True,
	'git'		: False,		# Turn on/off git core (Never turn-on at the head node!!!!!!!)
	'gif'		: False,		# Generate gif instead pop up on a screan.
	'fftkernel'	: 5.0,
	'isikernel'	: 5.0,
	'corefunc'	: (4,8,64),
	'coreindex'	: False,		# Turn on/off Core indexing
	'corelog'	: 'network',
	'noexit'	: False,
	'FPcurve'	: False,
	'GPcurve'	: False,
	'IGcurve'	: False,
	'Connectom'	: False,
	'G-Spikerate'	: False,
	'F-Spikerate'	: False,
	'Gtot-dist'	: False,
	'Gtot-rec'	: False,#True,	#record all gtotal in neurons
	'Gtot-stat'	: True, #False, #record gtot statistic
	'sycleprop'	: False,#True,
	'external'	: False,
	## Min
	'extprop'	: 0.5,				# Calculate probability to fire after external input

	'timestep'	: 0.01,#0.005,
	'sortbysk'	: False, #'ST',#'F',#'GT',#'NC',#'GT','G',#'I', #'F',#'I',#'F',#False,			#Do not use
	'taunorm'	: False,#True,#False,#True,
	'nstart'	: False, #
	'cliprst'	: False,#10,#False,#20,
	'TaS'		: True,
	'lastspktrg': True,
	'fullrast'	: True,
	'bias'		: False, #0.0000015
	'gtot-dist'	: 'NORM',#'LOGN', #LOGN - lognormal, 'NORM' - normal
	'gsyn-dist' : 'LOGN',#'NORM',#'LOGN', #same
	'initset'	: None,
	'cycling'	: False, #4,False
	'popfr'		: False,		#calculate population firing rate
}

#Neuron paramters:
if ntype == "RS":
	#RS neuron:
	#a,b,c,d,U,V	= 0.1,0.26,-65.0,-1.0,-16.083619212816,-61.86007190636312
#	a,b,c,d,U,V	= 0.1,0.26,-65.0,-1.0,-13.96214168703548,-65.81070843517298# Seady Satate at Iapp=0.15
#	a,b,c,d,U,V	= 0.1,0.26,-65.0,-1.0,(-13.96214168703548,2),(-65.81070843517298,5)
#	a,b,c,d,U,V	= 0.1,0.26,-65.0,-1.0,(-16.083619212816,5),(-61.86007190636312,20)
	a,b,c,d,U,V	= 0.1,0.26,-65.0,-1.0,(-15.083619212816,5),(-51.86007190636312,20)
	F 			= 1.#1.25 #0.25
	Iapp		= 0.0000015	# Iapp. it's equal to 0.1955 in XPP or Matlab. 0.000001955 in the rheobase for RS
#	Iapp		= (0.000001955,0.0000008)	#Don't use it. It's fake!
	Istdev		= 0.000001	# SD for 0.01ms step. it is equal to 0.09 in XPP
#	Istdev		= 0.0000008
#	Istdev		= 0.00000125
#	Istdev		= 0.000008
##Exp2Syn
	weight		= (3e-7, 0.0)	# Synaptic conductance. It is equal to 0.05 in XPP and Matlab
##e2ssn
#	weight		= (1e-7*7.62224658625,0.0) #e2ssn
#	weight		= (3e-7*7.62224658625,0.0) #e2ssn
#	weight		= (5e-7*7.62224658625,0.0) #e2ssn
	delay		= (0.1,0)			# delay in a system
#	Iapp		= 0.00002	# Iapp. it's equal to 0.1955 in XPP or Matlab. 0.000001955 in the rheobase for RS
#	Istdev		= 0.000005		# SD for 0.01ms step. it is equal to 0.2 in XPP
#	delay		= (5.0,0)
	gsynscale	= 1.0

elif ntype == "RS-moca":
	#RS neuron:
	a,b,c,d,U,V	= 0.1,0.26,-65.0,2.0,-16.083619212816,-61.86007190636312
	#a,b,c,d,U,V	= 0.1,0.26,-65.0,-1.0,(-16.083619212816,5),(-61.86007190636312,20)
	F 			= 1.#1.25 #0.25
	Iapp		= 0.0000015	# Iapp. it's equal to 0.1955 in XPP or Matlab. 0.000001955 in the rheobase for RS
#	Iapp		= (0.000001955,0.0000008)	#Don't use it. It's fake!
	Istdev		= 0.00000#1	# SD for 0.01ms step. it is equal to 0.09 in XPP
#	Istdev		= 0.0000008
#	Istdev		= 0.00000125
#	Istdev		= 0.000008
##Exp2Syn
	weight		= (3e-7, 0.0)	# Synaptic conductance. It is equal to 0.05 in XPP and Matlab
##e2ssn
#	weight		= (1e-7*7.62224658625,0.0) #e2ssn
#	weight		= (3e-7*7.62224658625,0.0) #e2ssn
#	weight		= (5e-7*7.62224658625,0.0) #e2ssn
	delay		= (0.1,0)			# delay in a system
#	Iapp		= 0.00002	# Iapp. it's equal to 0.1955 in XPP or Matlab. 0.000001955 in the rheobase for RS
#	Istdev		= 0.000005		# SD for 0.01ms step. it is equal to 0.2 in XPP
#	delay		= (5.0,0)
	gsynscale	= 1.0

elif ntype == "FS":
	#FS neuron:
	a,b,c,d,U,V	= 0.02,0.2,-65.0,2.0,(-12.836,2),(-64.183,10)
#	a,b,c,d,U,V	= 0.02,0.2,-65.0,2.0,-12.836,-64.183
#	a,b,c,d,U,V	= 0.02,0.2,-65.0,2.0,-13.836,-12.183
	F 			= 1.#1.25 #0.25
	Iapp		= 0.0000015	# Iapp. it's equal to 0.1955 in XPP or Matlab. 0.000001955 in the rheobase for RS
#	Iapp		= (0.000001955,0.0000008)	#Don't use it. It's fake!
	Istdev		= 0.00000#1	# SD for 0.01ms step. it is equal to 0.09 in XPP
#	Istdev		= 0.0000008
#	Istdev		= 0.00000125
#	Istdev		= 0.000008
##Exp2Syn
	weight		= (3e-7, 0.0)	# Synaptic conductance. It is equal to 0.05 in XPP and Matlab
##e2ssn
#	weight		= (1e-7*7.62224658625,0.0) #e2ssn
#	weight		= (3e-7*7.62224658625,0.0) #e2ssn
#	weight		= (5e-7*7.62224658625,0.0) #e2ssn
	delay		= (0.1,0)			# delay in a system
#	Iapp		= 0.00002	# Iapp. it's equal to 0.1955 in XPP or Matlab. 0.000001955 in the rheobase for RS
#	Istdev		= 0.000005		# SD for 0.01ms step. it is equal to 0.2 in XPP
#	delay		= (5.0,0)
	gsynscale	= 1.0

elif ntype == "TypeI":
	#Type I neuron:
	a,b,c,d,U,V	= 0.02,-0.1,-55,6,6.3908,-63.9
	F = 1.0
	Iapp		= 0.00010		# Iapp. it's equal to 10 in XPP or Matlab# 0.00022562 in the rheobase for TypeI
	Istdev		= 0.00019		# SD for 0.01ms step. it is equal to 19 in XPP
##Exp2Syn
#	weight		= (2.5e-5,0.0)	# Synaptic conductance. It is equal to 2.5 in XPP and Matlab
##e3ssn
	weight		= (2.5e-5*7.62224658625,0.0)	# Synaptic conductance. It is equal to 2.5 in XPP and Matlab
#	delay		= (7.5,0)		# delay in a system
	delay		= (0.1,0)		# delay in a system
#	Iapp		= 0.00032		# Iapp. it's equal to 10 in XPP or Matlab# 0.00022562 in the rheobase for TypeI
#	Istdev		= 0.00005		# SD for 0.01ms step. it is equal to 19 in XPP

#Synaptic paramters:
ST1,ST2,SE	= 2.0,5.0,-70.0
#ST1,ST2,SE	= 2.0,10.0,-70.0
synscaler	= None

#Simulation paramters:
tstop		= 10001			#10000
tvl,tvr		= 0, 1500 #1000 #1000 # 2000
trustbrak	= 5
nsig		= 3
unlock		= 12 #24, 64, 78, 9959 12(?) ###12 - doesn't work!!! try somwthing different.......
#####

class neuron:
	def __init__(self):
		self.soma = h.Section()
		self.soma.L = 1.
		self.soma.diam=1/np.pi
		self.soma.nseg=1
		self.soma.cm=1
		self.izh	= h.izhcur(0.5, sec=self.soma)
		self.innp	= h.InNp(0.5, sec=self.soma)
		self.rnd	= h.Random(np.random.randint(0,32562))
		self.innp.noiseFromRandom(self.rnd)
		self.innp.dur	= h.tstop
		self.innp.delay	= 0
		self.innp.per	= 0.1
		self.innp.mean	= 0.0
		self.innp.stdev	= 0.0
		self.syn	= h.Exp2Syn(0.5, sec=self.soma)
		#self.syn	= h.e2ssn(0.5, sec=self.soma)
		self.syn.e		= -75.0
		self.syn.tau1	= 2.0
		self.syn.tau2	= 10.0
		######## Recorders ##########
		self.spks	= h.Vector()
		self.sptr	= h.APCount(.5, sec=self.soma)
		self.sptr.thresh = 25
		self.sptr.record(self.spks)
		if methods['gui']:
			self.inoise	= h.Vector()
			self.inoise.record(self.innp._ref_i)
			self.volt	= h.Vector()
			self.volt.record(self.soma(0.5)._ref_v)
			self.uolt	= h.Vector()
			self.uolt.record(self.izh._ref_u)
			self.isyn	= h.Vector()
			self.isyn.record(self.syn._ref_i)
			self.gsyn	= h.Vector()
			self.gsyn.record(self.syn._ref_g)
		######## Registrations ###### 
		self.gsynscale	= 0.0
		self.concnt		= 0.0
		self.gtotal		= 0.0
		self.tsynscale	= 1.0
	def setparams(self,a = None, b = None, c = None, d = None, e = None, f = None, g = None, V=None, U = None, F = None, Iapp = None, Insd = None, SynE = None, SynT1 = None, SynT2 = None):
		if a != None : self.izh.a = a
		if b != None : self.izh.b = b
		if c != None : self.izh.c = c
		if d != None : self.izh.d = d
		if e != None : self.izh.e = e
		if f != None : self.izh.f = f
		if g != None : self.izh.g = g
		if V != None : self.soma(0.5).v = V
		if U != None : self.izh.uinit = U
		if F != None : self.izh.F = F
		if methods['bias']: self.izh.I = methods['bias']
		########
		if Iapp != None : self.innp.mean  = Iapp
		if Insd != None : self.innp.stdev = Insd
		self.innp.dur = h.tstop
		self.innp.delay = 0
		self.innp.per = 0.1
		########
		if SynE != None:  self.syn.e	= SynE
		if SynT1 != None: self.syn.tau1	= SynT1
		if SynT2 != None: self.syn.tau2	= SynT2
	def addnoise(self,Iapp=0.,Insd=0.,delay=0.,dur=0.,per=0.1):
		self.andnoise = h.InNp(0.5, sec=self.soma)
		self.andrnd	= h.Random(np.random.randint(0,32562))
		self.andnoise.noiseFromRandom(self.andrnd)
		self.andnoise.dur	= h.tstop
		self.andnoise.mean  = Iapp
		self.andnoise.stdev = Insd
		self.andnoise.delay	= delay
		self.andnoise.per	= per
		self.andnoise.dur	= dur

#class symulation:
	#def __init___(self,params):
		#if params.get("a",False):

def onclick1(event):
	if not hasattr(onclick1,"aix"):
		aix=zooly.add_subplot(111)
	onclick1.et = event.xdata
	
	### BUG
	onclick1.tl, onclick1.tr = onclick1.et-methods['traceWidth'], onclick1.et+methods['traceWidth']
	onclick1.idx, = np.where( (t > onclick1.tl) * (t < onclick1.tr))
	
	if not hasattr(onclick1,"marks"):
		onclick1.marks = []
		onclick1.marks.append( p.plot([onclick1.tl,onclick1.tl],[-80,30],"r--",lw=2)[0] )
		onclick1.marks.append( p.plot([onclick1.tr,onclick1.tr],[-80,30],"r--",lw=2)[0] )
	else:
		onclick1.marks[0].set_xdata([onclick1.tl,onclick1.tl])
		onclick1.marks[1].set_xdata([onclick1.tr,onclick1.tr])
	
	if not hasattr(onclick1,"lines"):
		onclick1.lines = []
		for n in neurons:
			volt = np.array(n.volt)
			onclick1.lines.append(aix.plot(t[onclick1.idx],volt[onclick1.idx])[0])
	else:
		vmin,vmax = 1000,-1000
		for ind,n in map(None,xrange(ncell),neurons):
			volt = np.array(n.volt)
			if vmin > volt[onclick1.idx].min():vmin = volt[onclick1.idx].min()
			if vmax < volt[onclick1.idx].max():vmax = volt[onclick1.idx].max()
			onclick1.lines[ind].set_xdata(t[onclick1.idx])
			onclick1.lines[ind].set_ydata(volt[onclick1.idx])
			onclick1.lines[ind].set_linewidth(1)
			onclick1.lines[ind].set_ls("-")

		aix.set_xlim(onclick1.tl,onclick1.tr)
		#print vmin,"---",vmax
		aix.set_ylim(vmin,vmax)
	if hasattr(zoolyclickevent,"lines"):
		del zoolyclickevent.lines
	moddy.canvas.draw()
	zooly.canvas.draw()
	mainfig.canvas.draw()

def getnulls(vmin,vmax,gsyn,inoise,ibias):
	vx=np.linspace(vmin,vmax,200)
	u0=b*vx
#	v0 =0.04*vx*vx+5.*vx+140.-(ibias+gsyn*(vx-SE))*100000.
	v0 =0.04*vx*vx+5.*vx+140.-ibias*100000.
	v0n=0.04*vx*vx+5.*vx+140.-(inoise+gsyn*(vx-SE))*100000.
	#print "DB>\n vmin=%g, vmax=%g\n b=%g\n u0min=%g, u0max=%g\n<DB"%(vmin,vmax,b,u0[0],u0[-1])
	#print "gsyn=%g, biase Iapp=%g, total Iapp = %g"%(gsyn*100000.,ibias*100000.,inoise*100000.)
	return vx,u0,v0,v0n
	
def numsptk(postidx,idxrange):
	prespikes = np.array([])
	trange=t[idxrange]
	sptk = np.zeros(trange.size)
	for nidx in OUTList[postidx]:
		sptime = np.array(neurons[nidx].spks)
		sptime = sptime[ np.where( (sptime > trange[0]) * (sptime < trange[-1]) ) ]
		prespikes = np.append(prespikes,sptime)
	
	prespikes = np.sort(prespikes)
	#print prespikes
	accumulator = 0
	for tm in trange:
		mp = np.where(prespikes < tm)[0]
		sptk[np.where( trange == tm )] = mp.size
	return sptk
	
def getprespikes(postidx,tl,tr):
	postspk = []
	for nidx in OUTList[postidx]:
		#DB>>
		#print nidx,":",
		#<<DB
		for nspk in neurons[nidx].spks[ np.where( (neurons[nidx].spks >= tl)*(neurons[nidx].spks < tr) ) ]:
			#DB>>
			#print nspk,
			#<<DB
			postspk.append([nspk,nidx] )
		#DB>>
		#print 
		#<<DB
		
	return np.array( postspk )
	
def zoolyclickevent(event):
	zoolyclickevent.spikesymbol = "."
	if not hasattr(onclick1,"lines"): return
	et = event.xdata
	ev = event.ydata
	idx = np.where( np.abs(t-et)<h.dt)[0][0]
	#DB>>
	#print idx, et,ev
	#<<DB
	vmax = abs(neurons[0].volt.x[idx] - ev)
	zoolyclickevent.imax = 0
	for ind,n in map(None,xrange(ncell),neurons):
		onclick1.lines[ind].set_linewidth(1)
		onclick1.lines[ind].set_ls("-")
		if vmax > abs(n.volt.x[idx] - ev) :
			vmax = abs(n.volt.x[idx] - ev)
			zoolyclickevent.imax = ind
		#print vmax,n.volt.x[idx],ev
	onclick1.lines[zoolyclickevent.imax].set_linewidth(4)
	onclick1.lines[zoolyclickevent.imax].set_ls("--")
	zooly.canvas.draw()
	#print "#",zoolyclickevent.imax,"\t v:",neurons[zoolyclickevent.imax].volt.x[idx],"\t g:",neurons[zoolyclickevent.imax].gsyn.x[idx]
	
	zoolyclickevent.v,zoolyclickevent.u = np.array(neurons[zoolyclickevent.imax].volt),np.array(neurons[zoolyclickevent.imax].uolt)
	zoolyclickevent.g,zoolyclickevent.i = np.array(neurons[zoolyclickevent.imax].gsyn),np.array(neurons[zoolyclickevent.imax].inoise)
	vmin,vmax =zoolyclickevent.v[onclick1.idx].min(),zoolyclickevent.v[onclick1.idx].max()
	vx,u0,v0,v0n = getnulls(vmin,vmax,zoolyclickevent.g[idx], zoolyclickevent.i[idx],neurons[zoolyclickevent.imax].innp.mean)
	zoolyclickevent.sptk = numsptk(zoolyclickevent.imax,onclick1.idx)
	zoolyclickevent.rst = getprespikes(zoolyclickevent.imax,onclick1.tl, onclick1.tr)
	#print "size sptk %d, idx %d time %d"%(sptk.size,onclick1.idx.size,tprin[onclick1.idx].size)
	if not hasattr(zoolyclickevent,"lines"):
		zoolyclickevent.lines = []
		#zoolyclickevent.lines.append(faxi.plot(zoolyclickevent.rst[:,0],zoolyclickevent.rst[:,1],"k|",ms=9,lw=5)[0])
		zoolyclickevent.lines.append(faxi.plot(zoolyclickevent.rst[:,0],zoolyclickevent.rst[:,1],"k"+zoolyclickevent.spikesymbol,ms=9,lw=5)[0])
		zoolyclickevent.lines.append(vaxi.plot(tprin[onclick1.idx],zoolyclickevent.v[onclick1.idx],"k-")[0])
		zoolyclickevent.lines.append(uaxi.plot(tprin[onclick1.idx],zoolyclickevent.u[onclick1.idx],"k-")[0])
		zoolyclickevent.lines.append(gaxi.plot(tprin[onclick1.idx],zoolyclickevent.g[onclick1.idx],"k-")[0])
		zoolyclickevent.lines.append(naxi.plot(zoolyclickevent.v[onclick1.idx],zoolyclickevent.u[onclick1.idx],"k-")[0])
		zoolyclickevent.lines.append(naxi.plot(vx,u0,"k--")[0])
		zoolyclickevent.lines.append(naxi.plot(vx,v0,"k-.")[0])
		zoolyclickevent.lines.append(naxi.plot(vx,v0n,"k.")[0])
		zoolyclickevent.lines.append(iaxi.plot(tprin[onclick1.idx],zoolyclickevent.i[onclick1.idx],"k-")[0])
		#zoolyclickevent.lines.append(saxi.plot(tprin[onclick1.idx],zoolyclickevent.sptk,"k-")[0])
	else:
		zoolyclickevent.lines[0].set_xdata(zoolyclickevent.rst[:,0])
		zoolyclickevent.lines[0].set_ydata(zoolyclickevent.rst[:,1])
		zoolyclickevent.lines[1].set_xdata(tprin[onclick1.idx])
		zoolyclickevent.lines[1].set_ydata(zoolyclickevent.v[onclick1.idx])
		zoolyclickevent.lines[2].set_xdata(tprin[onclick1.idx])
		zoolyclickevent.lines[2].set_ydata(zoolyclickevent.u[onclick1.idx])
		zoolyclickevent.lines[3].set_xdata(tprin[onclick1.idx])
		zoolyclickevent.lines[3].set_ydata(zoolyclickevent.g[onclick1.idx])
		zoolyclickevent.lines[4].set_xdata(zoolyclickevent.v[onclick1.idx])
		zoolyclickevent.lines[4].set_ydata(zoolyclickevent.u[onclick1.idx])
		zoolyclickevent.lines[5].set_xdata(vx)
		zoolyclickevent.lines[5].set_ydata(u0)
		zoolyclickevent.lines[6].set_xdata(vx)
		zoolyclickevent.lines[6].set_ydata(v0)
		zoolyclickevent.lines[7].set_xdata(vx)
		zoolyclickevent.lines[7].set_ydata(v0n)
		zoolyclickevent.lines[8].set_xdata(tprin[onclick1.idx])
		zoolyclickevent.lines[8].set_ydata(zoolyclickevent.i[onclick1.idx])
		#zoolyclickevent.lines[9].set_xdata(tprin[onclick1.idx])
		#zoolyclickevent.lines[9].set_ydata(zoolyclickevent.sptk)
	faxi.set_ylim(0,ncell)
	vaxi.set_ylim(zoolyclickevent.v[onclick1.idx].min(),zoolyclickevent.v[onclick1.idx].max())
	uaxi.set_ylim(zoolyclickevent.u[onclick1.idx].min(),zoolyclickevent.u[onclick1.idx].max())
	gaxi.set_ylim(zoolyclickevent.g[onclick1.idx].min(),zoolyclickevent.g[onclick1.idx].max())
	iaxi.set_ylim(zoolyclickevent.i[onclick1.idx].min(),zoolyclickevent.i[onclick1.idx].max())
	#saxi.set_ylim(0,zoolyclickevent.sptk[-1]+2)
	faxi.set_xlim(onclick1.tl, onclick1.tr)
	naxi.set_xlim(vmin,vmax)
	naxi.set_ylim(zoolyclickevent.u[onclick1.idx].min(),zoolyclickevent.u[onclick1.idx].max())
	moddy.canvas.draw()
	
def zoolykeyevent(event):
	if not hasattr(zoolyclickevent,"lines"): return
	if event.key == "K":
		v,u,g,i = np.array(neurons[zoolyclickevent.imax].volt),np.array(neurons[zoolyclickevent.imax].uolt),np.array(neurons[zoolyclickevent.imax].gsyn),np.array(neurons[zoolyclickevent.imax].inoise)
		sptk = numsptk(zoolyclickevent.imax,onclick1.idx)
		rst = getprespikes(zoolyclickevent.imax,onclick1.tl, onclick1.tr)
		zoolyclickevent.lines.append(faxi.plot(zoolyclickevent.rst[:,0],zoolyclickevent.rst[:,1],zoolyclickevent.spikesymbol,ms=9,lw=5)[0])
		zoolyclickevent.lines.append(vaxi.plot(tprin[onclick1.idx],v[onclick1.idx])[0])
		zoolyclickevent.lines.append(uaxi.plot(tprin[onclick1.idx],u[onclick1.idx])[0])
		zoolyclickevent.lines.append(gaxi.plot(tprin[onclick1.idx],g[onclick1.idx])[0])
		zoolyclickevent.lines.append(naxi.plot(v[onclick1.idx],u[onclick1.idx])[0])
		zoolyclickevent.lines.append(iaxi.plot(tprin[onclick1.idx],i[onclick1.idx])[0])
		#zoolyclickevent.lines.append(saxi.plot(tprin[onclick1.idx],sptk)[0])
	elif event.key == "X":
		for lin in zoolyclickevent.lines:
			lin.remove()
		del zoolyclickevent.lines
	moddy.canvas.draw()	

def moddykeyevent(event):
	zoolykeyevent(event)


def neuronsoverview(event):
	global vindex
	if event.key == "up":
		vindex += 1
		if vindex >= ncell : vindex = ncell -1
	elif event.key == "down":
		vindex -= 1
		if vindex < 0 : vindex = 0
	elif event.key == "home": vindex = 0
	elif event.key == "end" : vindex = ncell -1
	if event.key == "pageup":
		vindex += 10
		if vindex >= ncell : vindex = ncell -1
	elif event.key == "pagedown":
		vindex -= 10
		if vindex < 0 : vindex = 0
	vtrace.set_ydata( np.array(neurons[vindex].volt)[:tprin.size])
	mainfig.canvas.draw()
		
#===============================================#
#               MAIN PROGRAMM                   #
#===============================================#
if __name__ == "__main__":
	if len(sys.argv) > 1:
		def patternmatch(pattern,arg):
			if arg[:len(pattern)] != pattern: return None
			return arg[len(pattern):]
		for arg in sys.argv:
			redarg = patternmatch("-corelog=",arg)
			if redarg != None: methods["corelog"] = redarg
			redarg = patternmatch("-gsyn=",arg)
			if redarg != None:
				redarg = redarg.split(",")
				if type(redarg) is list and  len(redarg) >= 2:
					weight = (float(redarg[0]),float(redarg[1]))
				else:
					weight = float(redarg[0])
			redarg = patternmatch("-gif=",arg)
			if redarg != None: methods['gif'] = redarg
			redarg = patternmatch("-Iapp=",arg)
			if redarg != None:
				redarg = redarg.split(",")
				if type(redarg) is list and  len(redarg) >= 2:
					Iapp = (float(redarg[0]),float(redarg[1]))
				else:
					Iapp = float(redarg[0])
			
			redarg = patternmatch("-Istd=",arg)
			if redarg != None: Istdev = float(redarg)
			redarg = patternmatch("-gui=",arg)
			if redarg != None:
				if redarg == "False" or redarg =="off" or redarg =="OFF" or redarg == '0':
					methods["gui"] = False
			redarg = patternmatch("-git",arg)
			if redarg != None: methods["git"] = True
			redarg = patternmatch("-core=",arg)
			if redarg != None: methods['corefunc'] = int(redarg)
			redarg = patternmatch("-F=",arg)
			if redarg != None:
				redarg = redarg.split(",")
				if type(redarg) is list and len(redarg) >= 2:
					F = (float(redarg[0]),float(redarg[1]))
				else:
					F = float(redarg[0])
			redarg = patternmatch("-noexit",arg)
			if redarg != None: methods["noexit"] = True
			redarg = patternmatch("-gsynscale=",arg)
			if redarg != None:
				redarg = redarg.split(",")
				if type(redarg) is list and  len(redarg) >= 2:
					gsynscale = (float(redarg[0]),float(redarg[1]))
				else:
					gsynscale = float(redarg[0])
			redarg = patternmatch("-FPcurve",arg)
			if redarg != None: methods["FPcurve"] = True
			redarg = patternmatch("-GPcurve",arg)
			if redarg != None: methods["GPcurve"] = True
			redarg = patternmatch("-dt=",arg)
			if redarg != None: methods["timestep"] = float(redarg)
#			redarg = patternmatch("-tstop=",arg)
#			if redarg != None: tstop = float(redarg)
			redarg = patternmatch("-n=",arg)
			if redarg != None: ncell = int(redarg)
			redarg = patternmatch("-c=",arg)
			if redarg != None:
				redarg = redarg.split(",")
				if len(redarg) == 1: ncon = int(redarg[0])
				if len(redarg) == 2: ncon = ( int(redarg[0]), int(redarg[1]) )
				if len(redarg) == 3: ncon = ( int(redarg[0]), int(redarg[1]) , int(redarg[2]))
			### Synapses
			redarg = patternmatch("-tau1=",arg);
			if redarg != None: ST1					= float(redarg)
			redarg = patternmatch("-tau2=",arg);	
			if redarg != None: ST2					= float(redarg)
			redarg = patternmatch("-Esyn=",arg);	
			if redarg != None: SE					= float(redarg)
			redarg = patternmatch("-taunorm=",arg);
			if redarg != None: methods['taunorm']	= bool(int(redarg))
			redarg = patternmatch("-tsynscaler=",arg);	
			if redarg != None: exec "synscaler = "+redarg

			redarg = patternmatch("-external=",arg);
			if redarg != None:
				exec( "methods['external'] = "+redarg )

			redarg = patternmatch("-delay=",arg)
			if redarg != None:
				redarg = redarg.split(",")
				if type(redarg) is list  and len(redarg) >= 2:
					delay = (float(redarg[0]),float(redarg[1]))
				else:
					delay = (float(redarg[0]), 0.0)
			
			redarg = patternmatch("-view",arg)
			if redarg != None: tstop			= tvr + 1
			redarg = patternmatch("-tstop=",arg)
			if redarg != None:
				tvr			= float(redarg)
				tstop		= tvr + 1
			redarg = patternmatch("-bias=",arg)
			if redarg != None: methods['bias'] = float(redarg)
			redarg = patternmatch("-sort=",arg)
			if redarg != None: methods['sortbysk'] = redarg
			redarg = patternmatch("-gtot-dist=",arg)
			if redarg != None: methods['gtot-dist'] = redarg
			redarg = patternmatch("-gsyn-dist=",arg)
			if redarg != None: methods['gsyn-dist'] = redarg
			redarg = patternmatch("-initset=",arg)
			if redarg != None: exec "methods['initset'] ="+redarg
			
			redarg = patternmatch("-Vinit=",arg)
			if redarg != None:
				methods['initset'] = None
				exec( "V = "+redarg)
			redarg = patternmatch("-Uinit=",arg)
			if redarg != None:
				methods['initset'] = None
				exec( "U = "+redarg)
			redarg = patternmatch("-nstart=",arg)
			if redarg != None:
				exec( "methods['nstart'] = "+redarg)
			if arg == '-h' or arg == '-help' or arg == '--h' or arg == '--help':
				print __doc__
				print 
				print "USAGE: nrngui -nogui -python network.py [parameters]"
				print "\nPARAMETERS:"
				print "-n=          number of neurons in population"
				print "-c=          number of connections per neuron"
				print "-Iapp=       apply current. Use scaling factor 1e-5 to get nA"
				print "             Iapp may be a constant or mean,standard deviation across population."
				print "-Istd=       amplitude of noise. Should be scaled by 1e-5 to get nA"
				print "-gui=ON|OFF  Turn on/off gui and graphs"
				print "-F=          Set up neuron dynamics scale factor"
				print "             F may be a constant or mean,standard deviation across population."
				print "-gsyn=       conductance of single synapse. Use scaling factor 1e-5 to get nS"
				print "             gsyn may be a constant or mean,standard deviation for all synapses in model"
				print "-gsynscale=  total synaptic conductance.  Use scaling factor 1e-5 to get nS"
				print "             gsynscale may be a constant or mean,standard deviation for all neurons within population"
				print "-tau1=       rising time constant in ms"
				print "-tau2=       falling time constant in ms"
				print "-Esyn=       synaptic reversal potential in mV"
				print "-taunorm=0|1 On or Off normalization by space under the curve"
				print "-tsynscaler= scaling coefficient for synaptic time constants"
				print "-delay=      axonal delay in ms"
				print "             delay may be a constant or mean,standard deviation for all synapses in model"
				print "-view        limits simulation and save memory"
				
				exit(0)
			
			

	
	with open("network.start","w") as fd:
		for arg in sys.argv: fd.write("%s "%arg)
	if (ST1 != 2.0 or ST2 != 5.0) and methods['taunorm']:
		from norm_translation import getscale
		nFactor = getscale(2.0,5.0,ST1,ST2)
		if type(weight) == tuple or type(weight) == list:
			weight		= (weight[0]*nFactor, weight[1]*nFactor)
		else:
			weight		= (weight*nFactor, 0.)

###DB>
	print "=================================="
	print "===         :PARAMETERS:       ==="
	print "=================================="
	print "  > a=",a,"b=",b,"c=",c,"d=",d
	print "  > U0=",U,"V0=",V
	print "  > F           = ",F
	print "  > Bias        = ",methods['bias']
	print "  > Iapp        = ",Iapp
	print "  > Istdev      =  %g"%Istdev
	print "  > gsyn        = ",weight
	print "  > gsyn dist.  = ", methods['gsyn-dist'] 
	print "  > delay       = ",delay
	print "  > gsynscale   = ", gsynscale
	print "  > gtotal dis. = ", methods['gtot-dist'] 
	print "  > dt          = ", methods["timestep"]
	print "  > ncell       = ", ncell
	print "  > ncon        = ", ncon
	print "  > Syn tau1/2  = ", ST1,"/",ST2
	print "  > Syn E       = ", SE
	print "  > Syn taunorm = ", methods['taunorm']
	print "  > Syn t scaler= ", synscaler
	print "  > External    = ", methods['external']
	print "  > tstop       = ", tstop
	print "  > Sorted by   : ", methods['sortbysk']
	print "  > Init setup  : ", methods['initset']
	print "  > Noise Start : ", methods['nstart']
	print "==================================\n"
###<DB

	
	if methods["gui"]:
		import matplotlib.pyplot as plt
		print "=================================="
		print "===        GUI turned ON       ==="
		print "==================================\n"
	
	h.init()
	h.tstop 	= tstop
	#h.v_init 	= V
	h.dt		= methods["timestep"]
	

	#### Create Neurons and setup noise and Iapp
	print "=================================="
	print "===        Create Neurons      ==="
	print "==================================\n"
	neurons = [ neuron() for x in xrange(ncell) ]
	
	if type(V) is float or type(V) is int:
		xV = V*np.ones(ncell)
	elif type(V) is tuple:
		xV = V[0]+np.random.randn(ncell)*V[1]
	elif type(V) is str:
		xV = np.genfromtxt(V)
	if type(U) is float or type(U) is int:
		xU = U*np.ones(ncell)
	elif type(U) is tuple:
		xU = U[0]+np.random.randn(ncell)*U[1]
	elif type(U) is str:
		xU = np.genfromtxt(U)
			
	for n,i in zip(neurons,xrange(ncell)):
		if type(F) is float or type(F) is int:
			Fx = F
		elif type(F) is tuple:
			Fx = F[0] + np.random.randn()*F[1]
			while Fx <= 0.0: Fx = F[0] + np.random.randn()*F[1]
		if type(Iapp) is float or type(Iapp) is int:
			xIapp = Fx*Iapp
		elif type(Iapp) is tuple:
			xIapp = (Fx*Iapp[0]+Fx*np.random.randn()*Iapp[1])
		if synscaler != None:
			if type(synscaler) is float or type(synscaler) is int:
				n.tsynscale = float(synscaler)
				xST1,xST2 = ST1 * n.tsynscale, ST2 * n.tsynscale
			elif type(synscaler) is list or type(synscaler) is tuple:
				if len(synscaler) >= 2:
					n.tsynscale = -1.0
					while( n.tsynscale < 0.0 ):
						n.tsynscale = synscaler[0] + np.random.randn()*synscaler[1]
					xST1,xST2 = ST1 * n.tsynscale, ST2 * n.tsynscale
				else :
					n.tsynscale = float(synscaler[0])
					xST1,xST2 = ST1 * n.tsynscale, ST2 * n.tsynscale
			else:
				xST1,xST2 = ST1,ST2
		else:
			xST1,xST2 = ST1,ST2


		n.setparams(a=a, b=b, c=c, d=d, V=xV[i], U=xU[i], F=Fx, SynT1=xST1, SynT2=xST2, SynE=SE, Iapp=-xIapp, Insd=Istdev)
	
	if methods['initset'] != None:
		for setting in methods['initset']:
			if type(setting[0]) is int:
				itr=[setting[0]]
			elif type(setting[0]) is tuple or type(setting[0]) is list:
				if len(setting[0]) == 1: itr=setting[0][0]
				if len(setting[0]) == 2: itr=range(setting[0][0],setting[0][1])
				if len(setting[0]) == 3: itr=range(setting[0][0],setting[0][1],setting[0][2])
			for n in itr:
				if 'V' in setting[1]:
					if type(setting[1]['V']) is float or type(setting[1]['V']) is int:
						neurons[n].soma(0.5).v = setting[1]['V']
					elif type(setting[1]['V']) is tuple or type(setting[1]['V']) is list:
						neurons[n].soma(0.5).v = setting[1]['V'][0]+setting[1]['V'][1]*np.random.randn()
				if 'U' in setting[1]:
					if type(setting[1]['U']) is float or type(setting[1]['U']) is int:
						neurons[n].izh.uinit  = setting[1]['U']
					elif type(setting[1]['U']) is tuple or type(setting[1]['U']) is list:
						neurons[n].izh.uinit  = setting[1]['U'][0]+setting[1]['U'][1]*np.random.randn()
						
	#DB>>
	#for n in neurons:
	#	print n.soma.v, n.izh.uinit
	#exit(0)
	#<<DB
		

	if methods['nstart']:
		for n in neurons:
			n.addnoise(Iapp=0.,Insd=methods['nstart'][1],delay=methods['nstart'][0],dur=methods['nstart'][2],per=0.1)

	h.finitialize()
	h.fcurrent()
	h.frecord_init()

	t = h.Vector()
	t.record(h._ref_t)


	#### Create Connection List:
	OUTList = [ [] for x in xrange(ncell)]
	if ncon > 0:
		print "=================================="
		print "===        Map Connections     ==="
		print "==================================\n"
		for toid in xrange(ncell):
			from_excaption = [ 0 for x in xrange(ncell) ]
			from_excaption[toid] = 1
			upcnt = 0
			total = 0
			if type(ncon) is tuple or type(ncon) is list:
				if len(ncon) > 1:
					neurons[toid].concnt = int(np.random.random()*(ncon[1]-ncon[0]) + ncon[0])
				else:
					neurons[toid].concnt = ncon[0]
			else:
				neurons[toid].concnt = ncon
			while  upcnt < 10000*ncell:
				upcnt += 1
				fromid = rnd.randint(0, ncell-1)
				if from_excaption[fromid] == 1 : continue
				upcnt  = 0
				total += 1
				from_excaption[fromid] = 1
				OUTList[toid].append(fromid)
				if total >= neurons[toid].concnt :break
			else:
				sys.stderr.write("Couldn't obey connections conditions\nNeuron:%d TOTLA:%d CURRENT:%d\n"%(toid,ncon,total))
				for x in OUTList:
					sys.stderr.write("ID:%d #%d\n"%(x[0],x[1]))
				sys.exit(1)

	#DB>
		#for i in OUTList:
			#print len(i)
			#for j in i:	print "%03d"%(j),
			#print
		#sys.exit(0)
	#<DB
	if methods['cycling']:
		print "=================================="
		print "===      Cycles counting       ==="
		print "=================================="
		mat=np.matrix( np.zeros((ncell,ncell)) )
		for i,vec in map(None,xrange(ncell),OUTList):
			mat[i,vec]=1
		kx = []
		for cnt in xrange(methods['cycling']):
			kx.append(np.trace(mat)/ncell)
			mat = mat.dot(mat)
		print "  > Cyclopedic numbers : ",kx
		print "==================================\n"
		#methods['cycling'] = kx
		del mat
		
		
	

	#### Create Conneactions:
	print "=================================="
	print "===    Make the Connections    ==="
	print "==================================\n"
	connections = []
	for x in map(None,xrange(ncell),OUTList):
		if type(gsynscale) is tuple :
			if gsynscale[1] <= 0: gx = gsynscale[0]
			elif methods["gtot-dist"] == "NORM":
				### Trancated normal
				gx = gsynscale[1]*np.random.randn()+gsynscale[0]
				while gx < 0.0 : gx = gsynscale[1]*np.random.randn()+gsynscale[0]
			elif methods["gtot-dist"] == "LOGN":
				### Lognormal
				gx = np.random.lognormal(mean=mp.log(gsynscale[0])-gsynscale[1]**2/2., sigma=gsynscale[1])
		else:
			gx = float(gsynscale)
		neurons[x[0]].gsynscale = gx
		for pre in x[1]:
			if type(delay) is tuple :
				if delay[1] > 0:
					dx = -1
					while dx < 0.1: dx = delay[1]*np.random.randn()+delay[0]
				else:
					dx = float(delay[0])
			else:
				dx = float(delay)
			if type(weight) is tuple :
				if weight[1] <= 0: wx = weight[0]
				elif methods["gsyn-dist"] == "NORM":
					#### Trancated normal
					wx = weight[1]*np.random.randn()+weight[0]
					while wx < 0.0 : wx = weight[1]*np.random.randn()+weight[0]
				elif methods["gsyn-dist"] == "LOGN":
					### Lognormal
					wx = np.random.lognormal(mean=np.log(weight[0])-weight[1]**2/2., sigma=weight[1])
			else:
				wx = float(weight)
			if type(ncon) is tuple or type(ncon) is list:
				if len(ncon) > 2:
					wx *= float(ncon[2])/float(neurons[x[0]].concnt)
			if methods['taunorm'] and synscaler != None:
				#DB print "norm by factor",1./neurons[x[0]].tsynscale
				wx /= neurons[x[0]].tsynscale
			#####DB>>
			#print "DB:gx=",gx,"dx=",dx,"wx=",wx
			#####<<DB
			#connections.append( h.NetCon(neurons[pre].soma(0.5)._ref_v,neurons[x[0]].syn,
					#25, dx, gx*wx,
					#sec=neurons[pre].soma) )
			connections.append( (h.NetCon(neurons[pre].soma(0.5)._ref_v,neurons[x[0]].syn,
					25, dx, gx*wx,
					sec=neurons[pre].soma),pre,x[0]) )
			neurons[x[0]].gtotal += gx*wx
	#DB>>
	#plt.figure(0)
	#w=np.array([c[0].weight[0] for c in connections])
	#print np.mean(w), np.std(w)
	#plt.hist(w,bins=50,range=(0,1e-6))
	#plt.show()
	#exit(0)
	#<<DB

	
	if methods['external']:
		ex_netstim	= h.NetStim(.5, sec = neurons[0].soma)
		ex_netstim.start	= methods['external'][0]
		ex_netstim.noise	= 0
		ex_netstim.interval	= methods['external'][1]
		ex_netstim.number	= methods['external'][2]
		ex_syn,ex_netcon = [],[]
		for n in neurons:
			if len(methods['external']) >=9:
				if rnd.random() > methods['external'][8]: continue
			ex_syn_new = h.Exp2Syn(0.5, sec=n.soma)
			ex_syn_new.e	= methods['external'][4]
			ex_syn_new.tau1	= methods['external'][5]
			ex_syn_new.tau2	= methods['external'][6]
			ex_syn.append(ex_syn_new)
			exdelay = -1.0
			if len(methods['external']) < 8:
				exdelay = 0.1
			elif type(methods['external'][7]) is tuple :
				while exdelay <= 0.0 : exdelay = methods['external'][7][0]+np.random.randn()*methods['external'][7][1]
			elif type(methods['external'][7]) is float or type(methods['external'][7]) is int :
				exdelay = methods['external'][7]
			
			ex_netcon_new	= h.NetCon(ex_netstim, ex_syn_new, 1,exdelay ,methods['external'][3], sec = n.soma)
			ex_netcon.append(ex_netcon_new)

	print "=================================="
	print "===           RUN              ==="
	print "==================================\n"
	npc = h.ParallelContext()
	if type(methods['corefunc']) is int:
		npc.nthread(methods['corefunc'])
		sys.stderr.write( "Setup %g core\n"%(methods['corefunc']) )
	elif delay[0] > h.dt*2 or delay[1] > h.dt*2 :
		#### Setup parallel context if there are delays
		if not os.path.exists("/etc/beowulf") and os.path.exists("/sysini/bin/busybox"):
			#I'm not on head node. I can use all cores
			methods['corefunc'] = methods['corefunc'][2]
			npc.nthread(methods['corefunc'])
			sys.stderr.write( "Setup %g core\n"%(methods['corefunc']) )
		elif os.path.exists("/etc/beowulf"):
			#I'm on head node. I grub only half :)
			methods['corefunc'] = methods['corefunc'][1]
			npc.nthread(methods['corefunc'])
			sys.stderr.write( "Setup %g core\n"%(methods['corefunc']) )
		else:
			#I'm on Desktop :(
			methods['corefunc'] = methods['corefunc'][0]
			npc.nthread(methods['corefunc'])
			sys.stderr.write( "Setup %g core\n"%(methods['corefunc']) )
	h.run()

	if methods["save"]:
		print "=================================="
		print "===     Saving the Results     ==="
		print "==================================\n"
		fd = open("network-param.txt","w")
		fd.write("ncell=%d\nncon=%d\nntype=%s\n"%(ncell,ncon,ntype))
		fd.write("a=%g\nb=%g\nc=%g\nd=%g\nU=%g\nV=%g\n"%(a,b,c,d,U,V))
		fd.write("Iapp=%g\nIstdev=%g\nweight=(%g,%g)\ndelay=(%g,%g)\n"%(Iapp,Istdev,weight[0],weight[1],delay[0],delay[1]))
		fd.write("ST1=%g\nST2=%g\nSE=%g\n"%(ST1,ST2,SE))
		fd.write("tstop=%g\ndt=%g\nnsampls=%d"%(tstop,h.dt,t.size()))
		fd.close()
		fdv = open("network-volt.csv","w")
		fdn = open("network-noise.csv","w")
		fds = open("network-syn.csv","w")
		for idx in xrange(int(np.round(t.size()))):
			vstr = "%g"%t.x[idx]
			nstr = "%g"%t.x[idx]
			sstr = "%g"%t.x[idx]
			for n in neurons:
				vstr +=",%g"%n.volt.x[idx]
				nstr +=",%g"%n.inoise.x[idx]
				sstr +=",%g"%n.isyn.x[idx]
			fdv.write(vstr+"\n")
			fdn.write(nstr+"\n")
			fds.write(sstr+"\n")
		fdv.close()
		fdn.close()
		fds.close()

	print "=================================="
	print "===          Analysis          ==="
	print "==================================\n"
	
	t = np.array(t)
	if methods['gui']:
		plt.rc('text', usetex = True )
		plt.rc('font', family = 'serif')
		plt.rc('svg', fonttype = 'none')
		mainfig = plt.figure(1)
		if methods['traceView']:
			cid = mainfig.canvas.mpl_connect('button_press_event', onclick1)
		p=plt.subplot(411)
		tprin=np.array(t)
		tprin = tprin[ np.where( tprin < tvr ) ]
		vindex = (ncell-1)/2
		vtrace, = plt.plot(tprin,np.array(neurons[vindex].volt)[:tprin.size],"k")
		plt.ylim(ymax=40.)
		mainfig.canvas.mpl_connect('key_press_event',neuronsoverview)
		plt.ylabel("Voltage (mV)", fontsize=16)
		if methods["external"]:
			ex0 = methods["external"][0]
			ex1 = methods["external"][1]
			for ex2 in xrange(methods["external"][2]):
				plt.plot([ex0+ex1*ex2,ex0+ex1*ex2],[0,30],"r--")
		plt.subplot(412,sharex=p)
		nurch = np.arange(1,ncell+1,1)
		if methods['sortbysk']:
			if methods['sortbysk'] == 'F':
				nindex = [ (neurons[i].izh.F,i) for i in xrange(ncell)]
				nindex.sort()
				#print nindex
				for i in xrange(ncell):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'I':
				nindex = [ (-neurons[i].innp.mean,i) for i in xrange(ncell)]
				nindex.sort()
				#print nindex
				for i in xrange(ncell):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'G':
				nindex = [ (-neurons[i].gsynscale,i) for i in xrange(ncell)]
				nindex.sort()
				#print nindex
				for i in xrange(ncell):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'NC':
				nindex = [ (-neurons[i].concnt,i) for i in xrange(ncell)]
				nindex.sort()
				#print nindex
				for i in xrange(ncell):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'GT':
				nindex = [ (-neurons[i].gtotal,i) for i in xrange(ncell)]
				nindex.sort()
				#print nindex
				for i in xrange(ncell):
					nurch[nindex[i][1]]=i
			if methods['sortbysk'] == 'ST':
				nindex = [ (neurons[i].tsynscale,i) for i in xrange(ncell)]
				nindex.sort()
				#print nindex
				for i in xrange(ncell):
					nurch[nindex[i][1]]=i
				
			#print nurch
	pmean = 0.
	pcnt  = 0


	meancur = np.zeros(t.size)
	spbins  = np.zeros( int(np.floor(tstop))+1 )
	specX	= np.arange(spbins.size, dtype=float)
	specX	*= 1000.0/tstop
	#pnum	= specX.size/2
	pnum 	= 200.*tstop/1000.0
	specX	= specX[:pnum]
	if methods["nrnFFT"]:
		specN	= np.zeros(pnum)
#	specV	= np.zeros(t.size())
	if 10 < methods["nrnISI"] <= 3000:
		isi		= np.zeros(methods["nrnISI"])
	if methods['coreindex']:
		coreindex = [0.0, 0.0]

	if methods['gui']:
		xrast = np.array([])
		yrast = np.array([])
	for (idx,n) in map(None,xrange(ncell),neurons):
		n.spks = np.array(n.spks)
		if methods['gui']:
			spk = n.spks[ np.where (n.spks < tvr) ]
			if not methods['cliprst']:
				xrast = np.append(xrast,spk)
				yrast = np.append(yrast,np.repeat(nurch[idx],spk.size))
			elif idx%methods['cliprst'] == 0:
				xrast = np.append(xrast,spk)
				yrast = np.append(yrast,np.repeat(nurch[idx],spk.size))
		
		if methods['cliptrn']:
			fstidx = np.where(n.spks > methods['cliptrn'] )[0][0]
			aisi = n.spks[fstidx+1:] - n.spks[fstidx:-1]
		else:
			aisi = n.spks[1:] - n.spks[:-1]
		if methods['coreindex']:
			coreindex[0] += np.sum((aisi[1:] - aisi[:-1])/aisi[:-1])
			coreindex[1] += aisi.size - 1
		if 10 < methods["nrnISI"] <= 3000:
			for i in aisi[ np.where(aisi < methods["nrnISI"]) ]:
				isi[ int(np.floor(i)) ] += 1.0
		pmean += np.sum(aisi)
		pcnt  += n.spks.size - 1
		if methods['gui']:
			if methods['tracetail'] == 'total current':
				meancur += np.array(n.isyn.x) + np.array(n.inoise.x)
			elif methods['tracetail'] == 'current':
				meancur += np.array(n.isyn.x)
			elif methods['tracetail'] == 'conductance':
				meancur += np.array(n.gsyn.x)
		spn	= np.zeros(spbins.size)
		for sp in n.spks:
			spbins[ int( np.floor(sp) ) ] +=1
			spn[ int( np.floor(sp) ) ] +=1
		if methods['cliptrn']:
			spn = spn[methods['cliptrn']:]
		if methods["nrnFFT"]:
			fft = np.abs( np.fft.fft(spn)*1.0/tstop )**2
			specN += fft[:pnum]
	if methods['gui']:
		plt.plot(xrast,yrast,"k|",mew=0.75,ms=1)#,ms=10)
		
		if methods['fullrast']	: plt.ylim(ymin=0,ymax=ncell)
		else			: plt.ylim(ymin=0)
	if methods["nrnFFT"]:
		specN /= float(ncell)#	specV /= float(ncell)
	
	if methods['cliptrn']:
		spbins = spbins[methods['cliptrn']:]
	
	if methods['popfr']:
		popfr = np.mean(spbins)
		print "=================================="
		print "===       MEAN FIRING RATE     ==="
		print "  > MFR =           ",popfr
		print "==================================\n"

	if (methods["netFFT"] or methods["nrnFFT"]):
		print "=================================="
		print "===            FFT             ==="
		print "==================================\n"
		fft = np.abs( np.fft.fft(spbins)*1.0/tstop )**2

	##EN
	#probscale = np.zeros(ncell + 1)
	#probscale[0] = 1./float(ncell + 1)
	#for x in range(1,ncell + 1):
		#probscale[x] = probscale[0]*probscale[x-1]
	#pspbin = np.array([ probscale[int(x)] for x in  spbins] )
	#en = np.sum( (-1)*pspbin*np.log(pspbin) )
	
	if methods['coreindex']:
		coreindex = coreindex[0]/coreindex[1]
		print coreindex, 1./(1.+ abs(coreindex))
		sys.exit(0)
		#with open("coreindex.csv","w") as fd:
			#for i in coreindex: fd.write("%g\n"%i)
		#coreindex = np.corrcoef(coreindex[:-1],y=coreindex[1:])[0][1]
	
	#external stimulation index
	if methods['external'] and methods['extprop']:
		print "=================================="
		print "===      Spike Probability     ==="
		print "==================================\n"
		spprop = 0
		for etx in xrange(methods['external'][2]):
			lidx = int( np.floor(methods['external'][0]+methods['external'][1]*etx) )
			ridx = int( np.floor(lidx + methods['external'][1]*methods['extprop']) )
			spprop += float( np.sum(spbins[lidx:ridx]) )
		spprop /= methods['external'][2]*ncell
			

	if methods["peakDetec"] or methods["R2"]:
		print "=================================="
		print "===         Peak Detector      ==="
		print "==================================\n"
		kernel = np.arange(-methods["gkernel"][1],methods["gkernel"][1],1.)
		kernel = np.exp(kernel**2/methods["gkernel"][0]/(-methods["gkernel"][0]))
		module = np.convolve(spbins,kernel)
		module = module[kernel.size/2:1-kernel.size/2]
		#spbinmax = (np.diff(np.sign(np.diff(module))) < 0).nonzero()[0] + 1
		#spbinmin = (np.diff(np.sign(np.diff(module))) > 0).nonzero()[0] + 1
		spbinmark = []
		for idx in (np.diff(np.sign(np.diff(module))) < 0).nonzero()[0] + 1:
			spbinmark.append([idx,1])
		for idx in (np.diff(np.sign(np.diff(module))) > 0).nonzero()[0] + 1:
			spbinmark.append([idx,-1])
		spbinmark.sort()
		spbinmark = np.array(spbinmark)
		peakmark  = []
		spc,ccnt = 0.,0.
		for mx in np.where( spbinmark > 0 )[0]:
			if mx <= 2 or mx >= (spbinmark.size/2 -2):continue
			if spbinmark[mx-1][1] > 0 or spbinmark[mx+1][1] > 0 or spbinmark[mx][1] < 0:continue
			peakmark.append(spbinmark[mx])
			ccnt += 1
			spc += np.sum(spbins[spbinmark[mx-1][0]:spbinmark[mx+1][0]])
		if ccnt > 0:
			spc /= ccnt
		
			

	##R2
	##Per
	if methods["R2"]:
		print "=================================="
		print "===             R2             ==="
		print "==================================\n"
		X,Y,Rcnt,netpermean,netpercnt=0.,0.,0.,0.0,0.0
		phydist  = []
		for mx in np.where( spbinmark > 0 )[0]:
			if mx >= (spbinmark.size/2 - 3):continue
			if spbinmark[mx+1][1] > 0 or spbinmark[mx+2][1] < 0 or spbinmark[mx][1] < 0:continue
			Pnet = float(spbinmark[mx+2][0] - spbinmark[mx][0])
			netpermean += Pnet
			netpercnt  += 1.
			for n,i in map(None,spbins[spbinmark[mx][0]:spbinmark[mx+2][0]],xrange(spbinmark[mx+2][0] - spbinmark[mx][0])):
				phyX = np.cos(np.pi*2.*float(i)/Pnet)
				phyY = np.sin(np.pi*2.*float(i)/Pnet)
				X += n*phyX
				Y += n*phyY
				Rcnt += n
				if methods['sycleprop']:
					#phydist.append( (360.*np.arctan2(phyY,phyX)/2/np.pi,n) )
					phydist.append( (np.arctan2(phyY,phyX),n) )
		if Rcnt > 0.:
			R2 = np.sqrt((X/Rcnt)**2+(Y/Rcnt)**2)
		if netpercnt > 0.:
			netpermean /= ( netpercnt - 1)
		if methods['sycleprop']:
			phydist = np.array(phydist)
			phydist[:,1] /= np.sum(phydist[:,1])
			phyhist,phyhistbins = np.histogram(phydist[:,0], bins=37, weights=phydist[:,1],range=(-np.pi-np.pi/36,np.pi+np.pi/36))
			
		
	if 10 < methods["netISI"] < 3000:
		print "=================================="
		print "===          NET ISI           ==="
		print "==================================\n"
		netisi	= np.zeros(methods["netISI"])
		lock = threading.RLock()
		
		def calcnetisi(ns):
			global netisi, lock
			scans	= np.zeros(ncell,dtype=int)
			localnetisi = np.zeros(methods["netISI"])
			for n in ns:
				for sp in n.spks:
					for (idx,m) in map(None,xrange(ncell),neurons):
						if m.spks.size < 2 : continue
						while m.spks[scans[idx]] <= sp and scans[idx] < m.spks.size - 1 : scans[idx] += 1
						if m.spks[scans[idx]] <= sp : continue
						if m == n and m.spks[scans[idx]] - sp < 1e-6 : continue
						aisi = m.spks[scans[idx]] - sp
						if int(aisi) >= methods["netISI"] : continue
						localnetisi[ int(aisi) ] += 1
			with lock:
				netisi += localnetisi
		pids = [ threading.Thread(target=calcnetisi, args=(neurons[x::methods['corefunc']],)) for x in xrange(methods['corefunc']) ]
		for pidx in pids:
			pidx.start()
			#print pidx, "starts"
		for pidx in pids:
			pidx.join()
			#print pidx,"finishs"
			
			
	if methods["TaS"] or methods['lastspktrg']:
		print "=================================="
		print "===           T & S            ==="
		print "==================================\n"
		allspikes,activeneurons = [],0.
		for n in neurons:
			allspikes += list(n.spks)
			if n.spks.size !=0 :activeneurons += 1.
		allspikes.sort()
		allspikes = np.array(allspikes)
		TaSisi = allspikes[1:]-allspikes[:-1]
		if methods['lastspktrg']:
			lastspktrg = int( np.mean(allspikes) > tstop/4. )
			
		
		del allspikes
		mean1TaSisi = np.mean(TaSisi)
		TaSindex	= (np.sqrt(np.mean(TaSisi**2) - mean1TaSisi**2)/mean1TaSisi - 1.)/np.sqrt(activeneurons) 
		
		
	if methods['Gtot-dist'] :
		print "=================================="
		print "===   G-scaler distribution    ==="
		print "==================================\n"
		gsk = [ n.gsynscale for n in neurons ]
		gskhist,gskbins = np.histogram(gsk, bins=ncell/25, normed=True)#/10)#, normed=True)
		gskhist /= np.sum(gskhist)
		del gsk
		#DB>>
		#print gskhist
		#print gskbins
		#print np.sum(gskhist)
		#<<DB
		
	if methods['Gtot-stat']:
		print "=================================="
		print "===     G-total Statistics     ==="
		agtot = np.array([ n.gtotal/n.concnt for n in neurons ])
		mgtot = np.mean(agtot)
		sgtot = np.std(agtot)
		print "  > mean   gtotal =           ",mgtot
		print "  > stderr gtotal =           ",sgtot
		print "  > CV     gtotal =           ",sgtot/mgtot
		print "==================================\n"

			
	#EN
	#p.set_title("Mean individual Period = %g, Sychrony(Entropy) = %g(%g)"%(pmean/pcnt,1./(1.+en),en))
	
	##R2
	if methods['gui']:
		title = "Mean individual Period = %g"%(pmean/pcnt)
		if methods['popfr']:
			title += 'Mean FR =%g'%popfr
		if methods["R2"]:
			if Rcnt > 0 :
				title += r", $R^2$ = %g, Mean network Period = %g, Spike per cycle = %g"%(R2**2,netpermean,spc)
			else:
				title += ", *Fail to estimate network period*"
		elif methods["peakDetec"]:
			title += ", Spike per cycle = %g. "%(spc)
		if methods['TaS']:
			title += ", TaS = %g"%TaSindex
		if methods['lastspktrg']:
			title += ", LST = %g"%lastspktrg
		p.set_title(title)

		plt.subplot(413,sharex=p)
		if methods['cliptrn']:
			nppoints = np.arange(tvl+methods['cliptrn'],tvr,1.0)
			plt.bar(nppoints,spbins[:tvr-methods['cliptrn']],0.5,color="k")
			hight = spbins[:tvr-methods['cliptrn']].max()
			if methods["peakDetec"] or methods["R2"] :
				for mark in spbinmark:
					if mark[0]+methods['cliptrn'] < tvl or mark[0]+methods['cliptrn'] > tvr: continue
					if mark[1] > 0:
						plt.plot([mark[0]+methods['cliptrn'],mark[0]+methods['cliptrn']],[0,hight],"r--")
					else:
						plt.plot([mark[0]+methods['cliptrn'],mark[0]+methods['cliptrn']],[0,hight],"b--")
		else:
			nppoints = np.arange(tvl,tvr,1.0)
			plt.bar(nppoints,spbins[tvl:tvr],0.5,color="k")
			hight = spbins[tvl:tvr].max()
			if methods["peakDetec"] or methods["R2"] :
				for mark in spbinmark:
					if mark[0] < tvl or mark[0] > tvr: continue
					if mark[1] > 0:
						plt.plot([mark[0],mark[0]],[0,hight],"r--")
					else:
						plt.plot([mark[0],mark[0]],[0,hight],"b--")
#			plt.plot(nppoints,module[tvl:tvr]/np.sum(kernel),"k--")
#			plt.plot(nppoints,module[tvl:tvr],"k--")
		plt.ylabel("Rate (ms$^{-1}$)", fontsize=16)


	meancur = meancur / float(-ncell)
#	meancur = meancur / float(ncell)
	if methods['gui']:
		plt.subplot(414,sharex=p)
		if methods['tracetail'] == 'total current' or methods['tracetail'] == 'current':
			if ntype == "RS":
				#Hysteresis of type II
				plt.ylabel("Current (nA)", fontsize=16)
				plt.plot(tprin,meancur[:tprin.size]*10e6)
				plt.plot([tvl,tvr],[0.000001821*10e6,0.000001821*10e6],"r--")
				plt.plot([tvl,tvr],[0.000002625*10e6,0.000002625*10e6],"r--")
				#plt.plot([0,2000],[0.00000,0.0],"r--")
			elif ntype == "TypeI":
				# Sidle-node on type I
				plt.ylabel(r"Current ($\mu$A)", fontsize=16)
				plt.plot(tprin,meancur[:tprin.size]*10e3)
				plt.plot([0,2000],[0.00022562*10e3,0.00022562*10e3],"r--")
			else:
				plt.ylabel("Current (nA)", fontsize=16)
				plt.plot(tprin,meancur[:tprin.size]*10e6)
		elif methods['tracetail'] == 'conductance':
			plt.ylabel("Conductance (nS)", fontsize=16)
			plt.plot(tprin,meancur[:tprin.size]*(-10e6))
		elif methods['tracetail'] == 'firing rate' and ( methods["peakDetec"] or methods["R2"] ):
			plt.ylabel("Firing Rate (ms$^{-1}$)", fontsize=16)
			plt.plot(nppoints,module[tvl:tvr]/np.sum(kernel),"k--")
			hight = np.max(module[tvl:tvr]/np.sum(kernel))
			for mark in spbinmark:
				if mark[0] < tvl or mark[0] > tvr: continue
				if mark[1] > 0:
					plt.plot([mark[0],mark[0]],[0,hight],"k--")
		plt.xlabel("time (ms)", fontsize=16)


	
	if (methods["netFFT"] or methods["nrnFFT"]) and methods['gui']:
		plt.figure(2)
		if methods["netFFT"] and methods["nrnFFT"]:
			pl=plt.subplot(211)
		elif methods["netFFT"]:
			pl=plt.subplot(111)
		if methods["netFFT"]:
			plt.title("Network spectrum")
			plt.bar(specX[1:],fft[1:pnum],0.75)
		if methods["netFFT"] and methods["nrnFFT"]:
			plt.subplot(212,sharex=pl)
		elif methods["nrnFFT"]:
			plt.subplot(111)
		if methods["nrnFFT"]:
			plt.title("Neuronal spectrum")
			plt.bar(specX[1:],specN[1:],0.75)

	#plt.subplot(313,sharex=p)
	#specX =np.arange(0.0,tstop+h.dt,h.dt)
	#specX *= 1000.0/tstop/h.dt
	#pnum = specX.size/2
	#plt.title("Voltage spectrum")
	#plt.plot(specX[1:pnum],specV[1:pnum])
	#plt.xlim(0,200)
	
	if 10 < methods["netISI"] <= 3000 and sum(netisi) > 0: netisi /= sum(netisi)
	if 10 < methods["nrnISI"] <= 3000 and sum(isi) > 0: isi /= sum(isi)
	if (10 < methods["netISI"] <= 3000 or 10 < methods["nrnISI"] <= 3000) and methods['gui']:
		plt.figure(3)
		if 10 < methods["netISI"] <= 3000 and 10 < methods["nrnISI"] <= 3000:
			pl=plt.subplot(211)
		elif 10 < methods["netISI"] <= 3000 :
			plt.subplot(111)
		if 10 < methods["netISI"] <= 3000: 
			plt.title("Network ISI")
			plt.ylabel("Normalized number of interspike intervals", fontsize=16)
			plt.bar(np.arange(methods["netISI"]),netisi,0.75)
		if 10 < methods["netISI"] <= 3000 and 10 < methods["nrnISI"] <= 3000:
			plt.subplot(212)#,sharex=pl)
		elif 10 < methods["nrnISI"] <= 3000:
			plt.subplot(111)
		if 10 < methods["nrnISI"] <= 3000:
			plt.ylabel("Normalized number of interspike intervals", fontsize=16)
			plt.title("Neuronal ISI")
			plt.bar(np.arange(methods["nrnISI"]),isi,0.75,color='k')
			plt.xlim(0,methods["nrnISI"])
			plt.xlabel("Interspike intervals (ms)", fontsize=16)
	
	if methods['traceView'] and methods['gui']:
		zooly = plt.figure(4)
		zooly.canvas.mpl_connect('button_press_event', zoolyclickevent)
		zooly.canvas.mpl_connect('key_press_event', zoolykeyevent)
		moddy = plt.figure(5)
		faxi = plt.subplot2grid((6,10),(0,0),colspan=4,rowspan=2)
		vaxi = plt.subplot2grid((6,10),(2,0),colspan=4,sharex=faxi)
		uaxi = plt.subplot2grid((6,10),(3,0),colspan=4,sharex=faxi)
		gaxi = plt.subplot2grid((6,10),(4,0),colspan=4,sharex=faxi)
		iaxi = plt.subplot2grid((6,10),(5,0),colspan=4,sharex=faxi)
		#saxi = plt.subplot2grid((6,10),(5,0),colspan=4,sharex=faxi)
		naxi = plt.subplot2grid((6,10),(0,5),colspan=6,rowspan=6)
		moddy.canvas.mpl_connect('key_press_event', moddykeyevent)

	if methods['FPcurve'] and methods['gui']:
		plt.figure(6)
		f  = np.array([ [n.izh.F,n.spks.size]       for n in neurons])/tstop*1000. #because ms
		#f = np.sort(f, axis=0)
		plt.plot(f[:,0] , f[:,1],"k^",ms=9)
		xrng=plt.xlim()
		plt.plot(np.array(xrng),[1000./netpermean,1000./netpermean],"k--")

	if methods['GPcurve'] and methods['gui']:
		plt.figure(7)
		f  = np.array([ [n.gsynscale,n.spks.size]       for n in neurons])
		#f = np.sort(f, axis=0)
		plt.plot(f[:,0] ,f[:,1],"k+")
			
	if methods['sycleprop'] and methods['gui']:
		plt.figure(8)
		polarax = plt.subplot(111, polar=True)
		#bars = polarax.bar(phydist[:,1], phydist[:,0], width=0.25, bottom=0.0)
		#np.histogram(phydist[:,0], bins=180, weights=phydist[:,1])
		#polarax.hist(phydist[:,0], bins=36, weights=phydist[:,1])
		polarax.bar(phyhistbins[:-1],phyhist,width=phyhistbins[1]-phyhistbins[0],bottom=0)
		#DB>>
		plt.figure(9)
		plt.bar(phyhistbins[:-1],phyhist,width=phyhistbins[1]-phyhistbins[0],bottom=0)
		#<<DB
	if methods['Gtot-dist'] and methods['gui']:
		plt.figure(10)
		plt.bar(gskbins[:-1],gskhist,width=gskbins[1]-gskbins[0],color="b")
		#plt.hist(gsk,bins=ncell/50)

	if methods['corelog']:
		with open(methods['corelog']+".csv","a") as fd:
			
			#0:Type,1:a,2:b,3:c,4:d,5:(U),6:(V),7:(F),8:(Gscale),9:(Iapp),
			#10:(weight),11:(delay),12:Istdev,13:mean np,14:spc,15:R2,16:mean netp
			fd.write("%s,%g,%g,%g,%g,"%(ntype,a,b,c,d))
			if type(U) is float or type(U) is int:
				fd.write("%g,"%(U))
			elif type(U) is tuple:
				fd.write("%g:%g,"%U)
			elif type(U) is str:
				fd.write("%s,"%U)
			if type(V) is float or type(V) is int:
				fd.write("%g,"%(V))
			elif type(V) is tuple:
				fd.write("%g:%g,"%V)
			elif type(V) is str:
				fd.write("%s,"%V)

			if type(F) is float or type(F) is int:
				fd.write("%g,"%(F))
			else:
				fd.write("%g:%g,"%F)
			if type(gsynscale) is float or type(gsynscale) is int:
				fd.write("%g,"%(gsynscale))
			else:
				fd.write("%g:%g,"%gsynscale)
			if type(Iapp) is float:
				fd.write("%g,"%(Iapp))
			else:
				fd.write("%g:%g,"%Iapp)
			if type(weight) is float:
				fd.write("%g,"%(weight))
			else:
				fd.write("%g:%g,"%weight)
			if type(delay) is float:
				fd.write("%g,"%(delay))
			else:
				fd.write("%g:%g,"%delay)
			if pcnt > 1:
				fd.write("%g,%g,"%(Istdev,pmean/pcnt))
			else:
				fd.write("%g,x,"%(Istdev))
			if methods["R2"]:
				if Rcnt > 0 :
					fd.write("%g,%g,%g,"%(spc,R2,netpermean))
				else:
					fd.write("%g,x,%g,"%(spc,netpermean))
			elif methods["peakDetec"]:
				fd.write("%g,x,x,"%(spc))
			else:
				fd.write("x,x,x,")
			if methods['IGcurve']:
				fd.write("IGF,%d,"%(ncell))
				for n in neurons:
					fd.write("%g:%g:%g,"%(n.innp.mean,n.gsynscale,n.izh.F))
			if methods['Connectom']:
				fd.write("CONNECTOM,%d,"%(len(connections)))
				for n in connections:
					fd.write("%d:%d:%g:%g,"%(n[1],n[2],n[0].weight[0],n[0].delay))
			if methods["netFFT"] or methods["nrnFFT"]:
				fftkernel = np.arange(-methods['fftkernel']*3,methods['fftkernel']*3,1.)
				fftkernel = np.exp(fftkernel**2/(-methods['fftkernel']**2))
			if methods["netFFT"]:
				fftmodule = np.convolve(fft[1:pnum],fftkernel)
				fftmodule = fftmodule[fftkernel.size/2:1-fftkernel.size/2]
				fftmax = (np.diff(np.sign(np.diff(fftmodule))) < 0).nonzero()[0] + 1
				fftmin = (np.diff(np.sign(np.diff(fftmodule))) > 0).nonzero()[0] + 1
				fd.write("netFFTmax,%d,"%fftmax.size)
				for fftm in fftmax:
					fd.write("%g:%g,"%(specX[fftm+1],fft[fftm+1]))
				fd.write("netFFTmin,%d,"%fftmin.size)
				for fftm in fftmin:
					fd.write("%g:%g,"%(specX[fftm+1],fft[fftm+1]))
			if methods["nrnFFT"]:
				fftmodule = np.convolve(specN[1:],fftkernel)
				fftmodule = fftmodule[fftkernel.size/2:1-fftkernel.size/2]
				fftmax = (np.diff(np.sign(np.diff(fftmodule))) < 0).nonzero()[0] + 1
				fftmin = (np.diff(np.sign(np.diff(fftmodule))) > 0).nonzero()[0] + 1
				fd.write("nrnFFTmax,%d,"%fftmax.size)
				for fftm in fftmax:
					fd.write("%g:%g,"%(specX[fftm+1],specN[fftm+1]))
				fd.write("nrnFFTmin,%d,"%fftmin.size)
				for fftm in fftmin:
					fd.write("%g:%g,"%(specX[fftm+1],specN[fftm+1]))
			if 10 < methods["netISI"] <= 3000 or 10 < methods["nrnISI"] <= 3000:
				isikernel = np.arange(-methods['isikernel']*3,methods['isikernel']*3,1.)
				isikernel = np.exp(isikernel**2/(-methods['isikernel']**2))
			if 10 < methods["netISI"] <= 3000 :
				if sum(netisi) > 0:
					isimodule = np.convolve(netisi,isikernel)
					isimodule = isimodule[isikernel.size/2:1-isikernel.size/2]
					isimax = (np.diff(np.sign(np.diff(isimodule))) < 0).nonzero()[0] + 1
					isimin = (np.diff(np.sign(np.diff(isimodule))) > 0).nonzero()[0] + 1
					fd.write("netISImax,%d,"%isimax.size)
					for isim in isimax:
						fd.write("%g:%g,"%(isim-2,netisi[isim-2]))
					fd.write("netISImin,%d,"%isimin.size)
					for isim in isimin:
						fd.write("%g:%g,"%(isim-2,netisi[isim-2]))
				else:
					fd.write("netISImax,0,")
					fd.write("netISImin,0,")
			if 10 < methods["nrnISI"] <= 3000 :
				if sum(isi) > 0:
					isimodule = np.convolve(isi,isikernel)
					isimodule = isimodule[isikernel.size/2:1-isikernel.size/2]
					isimax = (np.diff(np.sign(np.diff(isimodule))) < 0).nonzero()[0] + 1
					isimin = (np.diff(np.sign(np.diff(isimodule))) > 0).nonzero()[0] + 1
					fd.write("nrnISImax,%d,"%isimax.size)
					for isim in isimax:
						fd.write("%g:%g,"%(isim-2,isi[isim-2]))
					fd.write("nrnISImin,%d,"%isimin.size)
					for isim in isimin:
						fd.write("%g:%g,"%(isim-2,isi[isim-2]))
				else:
					fd.write("nrnISImax,0,")
					fd.write("nrnISImin,0,")
			if methods['coreindex']:
				fd.write("%g,"%coreindex)
			if methods['F-Spikerate']:
				fd.write("FNC,%d,%g"%(ncell,tstop))
				for n in neurons:
					fd.write(",%f:%d"%(n.izh.F,n.spks.size))
			if methods['G-Spikerate']:
				fd.write(",GSNC,%d"%(ncell))
				for n in neurons:
					fd.write(",%f:%d"%(n.gsynscale,n.spks.size))
			if methods['TaS']:
				fd.write(",TaS,%g"%TaSindex)
			if methods['lastspktrg']:
				fd.write(",LST,%d"%lastspktrg)
			fd.write(",SYN:%g:%g:%g:%d"%(ST1,ST2,SE,int(methods['taunorm'])) )
			if methods['external'] and methods['extprop']:
				fd.write(",EXTPROP,%g,%g,%g,%d"%(spprop,methods['external'][0],methods['external'][1],methods['external'][2]))
			#Syn stat
			allconv = np.array( [ c[0].weight[0] for c in connections ] )
			fd.write(",STAT:G,%g,%g,%g,%g"%(np.mean(allconv),np.std(allconv),np.min(allconv),np.max(allconv)) )
			allconv = np.array( [ c[0].delay for c in connections ] )
			fd.write(",STAT:D,%g,%g,%g,%g"%(np.mean(allconv),np.std(allconv),np.min(allconv),np.max(allconv)) )
			# TSscaler
			if synscaler != None:
				if type(synscaler) is float or type(synscaler) is int:
					fd.write(",GScaler,{}".format(synscaler))
				elif type(synscaler) is list or type(synscaler) is tuple:
					if len(synscaler) >= 2:
						fd.write(",GScaler,{}:{}".format(synscaler[0],synscaler[1]))
					else:
						fd.write(",GScaler,{}".format(synscaler[0]))
				else:
					fd.write(",GScaler,None")
			else: fd.write(",GScaler,None")
			if methods['popfr']:fd.write(",MFR,%g"%popfr)
			if methods['sycleprop']:
				fd.write(",Cyledist,%d,%g"%(phyhist.shape[0],phyhistbins[1]-phyhistbins[0]))
				for p,n in zip(phyhistbins[:-1],phyhist):
					fd.write(",%g:%g"%(p,n))
			if methods['Gtot-rec']:
				fd.write(",Gtot-rec,%d"%(ncell))
				for n in neurons:
					fd.write(",%g:%g"%(n.gsynscale,n.gtotal))
			if methods['Gtot-stat']:
				agtot = np.array([ n.gtotal/n.concnt for n in neurons ])
				mgtot = np.mean(agtot)
				sgtot = np.std(agtot)
				fd.write(",Gtot-stat,%g,%g,%g"%(mgtot,sgtot,sgtot/mgtot))
			fd.write("\n")

	if methods['gui']:
		if methods['gif']:
			plt.savefig(methods['gif'])
		else:
			plt.show()
	if not methods['noexit']:
		sys.exit(0)