load_file("nrngui.hoc")
secondorder=2 
tstep=0
period=2
dt=0.1
tstop=2000
nbcell = 200
TC=0.4
objref Bcell[nbcell]

	begintemplate BasketCell
ndend1=4
ndend2=4
ndend3=4
ndend4=4

public  pre_list, connect_pre, subsets, is_connected, is_art, stim
public soma, bcdend1, bcdend2, bcdend3, bcdend4
public all, adend, bdend, cdend, ddend
create soma, bcdend1[ndend1], bcdend2[ndend2], bcdend3[ndend3], bcdend4[ndend4]

objref syn, pre_list, stim


objref syn
proc init() {
	pre_list = new List()
	subsets()
	temp()
	synapse()
	}

objref all, adend, bdend, cdend, ddend

proc subsets() { local i

	objref all, adend, bdend, cdend, ddend
	all = new SectionList()

     		soma all.append()
		bcdend1[0] all.append()
		bcdend2[0] all.append()
		bcdend3[0] all.append()
		bcdend4[0] all.append()

		bcdend1[1] all.append()
		bcdend2[1] all.append()
		bcdend3[1] all.append()
		bcdend4[1] all.append()

		bcdend1[2] all.append()
		bcdend2[2] all.append()
		bcdend3[2] all.append()
		bcdend4[2] all.append()

		bcdend1[3] all.append()
		bcdend2[3] all.append()
		bcdend3[3] all.append()
		bcdend4[3] all.append()


	adend  = new SectionList()
		for i=0,3{
		bcdend1 [i] adend.append()}

	bdend  = new SectionList()
		for i=0,3{
		bcdend2 [i] adend.append()}

	cdend  = new SectionList()
		for i=0,3{
		bcdend3 [i] adend.append()}

	ddend  = new SectionList()
		for i=0,3{
		bcdend4 [i] adend.append()}

}

proc temp() {
	soma {nseg=1 L=20 diam=15} // changed L & diam
	bcdend1 [0] {nseg=1 L=75 diam=4}	
	bcdend2 [0] {nseg=1 L=75 diam=3}
	bcdend3 [0] {nseg=1 L=75 diam=2}
 	bcdend4 [0] {nseg=1 L=75 diam=1}

	bcdend1 [1] {nseg=1 L=75 diam=4}
	bcdend2 [1] {nseg=1 L=75 diam=3}
	bcdend3 [1] {nseg=1 L=75 diam=2}
	bcdend4 [1] {nseg=1 L=75 diam=1}
 		 
	bcdend1 [2] {nseg=1 L=50 diam=4} 	
	bcdend2 [2] {nseg=1 L=50 diam=3}
	bcdend3 [2] {nseg=1 L=50 diam=2}
	bcdend4 [2] {nseg=1 L=50 diam=1} 
	
	bcdend1 [3] {nseg=1 L=50 diam=4}
	bcdend2 [3] {nseg=1 L=50 diam=3}
	bcdend3 [3] {nseg=1 L=50 diam=2}
	bcdend4 [3] {nseg=1 L=50 diam=1} 	


	forsec all {
		insert ccanl
	catau_ccanl = 10
	caiinf_ccanl = 5.e-6
		insert borgka
	gkabar_borgka=0.00015
		insert nca  // HAV-N- Ca channel
	gncabar_nca=0.0008   //check to modify- original 0.004
		insert lca 
	glcabar_lca=0.005
		insert gskch
	gskbar_gskch=0.000002
		insert cagk
	gkbar_cagk=0.0002
	insert tonic
    	g_tonic =  0.000008  //0.001
   	e_tonic = -74
	}

	soma {insert ichan2  //ildikos ichan
	gnatbar_ichan2=0.12  
	gkfbar_ichan2=0.013  
	gl_ichan2 = 0.00014
	cm=1.4
	} 

	forsec adend {insert ichan2
	gnatbar_ichan2=0.12  
	gkfbar_ichan2=0.013
	gl_ichan2 = 0.00012
	cm=1.4
	}		
	forsec bdend {insert ichan2
	gnatbar_ichan2=0.0
	gkfbar_ichan2=0.00
	gl_ichan2 = 0.00012
	cm=1.4}
		
	forsec cdend {insert ichan2
	gnatbar_ichan2=0.0
	gkfbar_ichan2=0.00
	gl_ichan2 = 0.00012
	cm=1.4}

	forsec ddend {insert ichan2
	gnatbar_ichan2=0.0
	gkfbar_ichan2=0.00
	gl_ichan2 = 0.00012
	cm=1.4}

	connect bcdend1[0](0), soma(1)
	connect bcdend1[1](0), soma(1)
	connect bcdend1[2](0), soma(0)
	connect bcdend1[3](0), soma(0)
	
	connect bcdend2[0](0), bcdend1[0](1)
	connect bcdend2[1](0), bcdend1[1](1)
	connect bcdend2[2](0), bcdend1[2](1)
	connect bcdend2[3](0), bcdend1[3](1)
	
	connect bcdend3[0](0), bcdend2[0](1)
	connect bcdend3[1](0), bcdend2[1](1)
	connect bcdend3[2](0), bcdend2[2](1)
	connect bcdend3[3](0), bcdend2[3](1)

	connect bcdend4[0](0), bcdend3[0](1)
	connect bcdend4[1](0), bcdend3[1](1)
	connect bcdend4[2](0), bcdend3[2](1)
	connect bcdend4[3](0), bcdend3[3](1)


		forsec all {Ra=100}
	forsec all {enat = 55 ekf = -90  ek=-90  elca=130	esk=-90
		 el_ichan2 =-75
		cao_ccanl=2 }  // make catau slower70e-3 	cao=2 cai=50.e-6 


	}

	objref syn  
	proc synapse() {
	bcdend4 [0] syn = new Exp2Syn(0.5)	
	syn.tau1 = 2	syn.tau2 = 6.3	syn.e = 0    //2 and 6.3         
	pre_list.append(syn)
	
	bcdend2 [1] syn = new Exp2Syn(0.5)	//BC(GABA) syn to apical IML dend Bartos
	syn.tau1 = 0.16		syn.tau2 = 1.8	syn.e = -74    //0.16 and 1.8
	pre_list.append(syn)
//******************************Tonic GABA synapses**************

	bcdend3[0] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum  2
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend3[1] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	3
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend3[2] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	4
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend3[3] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	5
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend2[0] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	6
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend1[2] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	7
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend2[2] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	8
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)
	
	bcdend2[3] syn = new Exp2Syn(0.5) // TONIC GABA synapse based on Rossi evoked data from cerebellum	9
	syn.tau1 = 7	syn.tau2 = 200	syn.e = -74
	pre_list.append(syn)

	}

	proc connect_pre() {  // $o1 target point process, $o2 returned NetCon
	soma $o2 = new NetCon (&v(1), $o1) 
	}

		func is_art()  { return 0 }
	endtemplate BasketCell
//**********************************************************************************
//********************************************************************

// NETWORK SPECIFICATION INTERFACE

	for i=0, nbcell-1 {Bcell[i] = new BasketCell(i)} 
	

//*********************************************************************


objref nclist, netcon, cells, net_c, net_bc
{  cells = new List()
nclist = new List()
}
 func cell_append() {cells.append($o1) 
	return cells.count -1}

func nc_append() {

	if ($3 >= 0 )	{
		cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon)
		netcon.weight = $4	netcon.delay = $5		netcon.threshold = $6
	} 
	nclist.append(netcon)
	return nclist.count-1
		}
func is_connected() {local i, c
	c=0
	for i=0, nclist.count-1 {
	net_c= nclist.object(i)
	if (($o1 == net_c.postcell())  && ($o2 == net_c.precell())) {c=1}
}
return c
}
objref vbc2bc
	{
	vbc2bc = new Vector(nbcell, 0)
	}

//initiating randm number generator

objref rdsyna, rdbc2bc, rdgap, rddend
ropen("/proc/uptime")		// get a seed  that is changing based on the processing time
	 {			
 	rseed = fscan()		// so simulation will not start with the same seed
	ropen()		
	}
//rseed=12345

rdgap = new Random(rseed)                            
proc new_rdgap() {rdgap.discunif(-20,20)}    // use for gap junc
new_rdgap()
rddend = new Random(rseed)                            
proc new_rddend() {rddend.discunif(0,3)}    
new_rddend()
rdbc2bc = new Random(rseed)			// use for syn.connections 
proc new_rdbc2bc() {rdbc2bc.discunif(-15, 15)}
new_rdbc2bc()
rdsyna = new Random(rseed)		// initialize random distr.
proc new_rdsyna() {rdsyna.discunif(1, 1)}
new_rdsyna()

//	NETWORK INITIATION
	for i = 0, nbcell-1 {cell_append(Bcell[i])} 
	

//*********************************************************************************
//**********************************************************************************


objref randomVector
randomVector = new Vector(nbcell)


proc initBcell() { local i,j

for  i=0, nbcell-1 {

	randomVector.resize(nbcell)
	for k = 0, randomVector.size-1 {randomVector.x[k]=0}
	while (randomVector.sum < 30) {
	Gauz3  = rdbc2bc.repick()
	dist = abs(Gauz3)
	axdel = dist/5*0.2
	if (Gauz3==0){Gauz3  = rdbc2bc.repick()}
	if (i+Gauz3==i){Gauz3  = rdbc2bc.repick()}
	if (Gauz3==0){Gauz3  = rdbc2bc.repick()}
	if (i+Gauz3==i){Gauz3  = rdbc2bc.repick()}
	if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell }
	if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell} 
	if ((i+Gauz3 >-1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3}
	if ((randomVector.x[npost] == 0) && (vbc2bc.x[npost] < 31)) {
	   randomVector.x[npost] += 1
	   vbc2bc.x[npost] += 1
	}
	dbr = rdsyna.repick()
	}
	for k = 0, randomVector.size-1 {
	    if (randomVector.x[k] == 1) {
	    	nc_append(i, k, dbr, 7.6e-3, 0.8+axdel, 0)  
		u=k
		v=u+1
		w=u-1
		r=0.4*TC
		if (v > 199) {v = v-1}
		if (w < 0) {w = w+1}
		nc_append(i, v, 8, r*0.0125e-3, 8, 10)
		nc_append(i, v, 9, r*0.0125e-3, 8, 10)
		nc_append(i, u, 8, TC*0.0125e-3, 8, 10)
		nc_append(i, u, 9, TC*0.0125e-3, 8, 10)
		nc_append(i, w, 8, r*0.0125e-3, 8, 10)
		nc_append(i, w, 9, r*0.0125e-3, 8, 10)

	    }
	}
   }
}


objref gaps[20000]
for i=0,19999{
gaps[i] = new Gap(0.5)
gaps[i].r = 100000}


n=0
for i=0,199{
d = rddend.repick() 
pre=i    //for bcells
for j=1,50{
post = j+i
n +=2 
if (post == pre){post =pre-post}   
if (post>199) {post = post-200}                                                                                   	  
Bcell[pre].bcdend3[d] gaps[n-2].loc(0.5)
Bcell[post].bcdend3[d] gaps[n-1].loc(0.5)                                                      
setpointer gaps[n-2].vgap, Bcell[post].bcdend3[d].v(0.5)
setpointer gaps[n-1].vgap, Bcell[pre].bcdend3[d].v(0.5)
//print "Bcell", pre, "gap pre", n-2, "is connected post Bcell", post, "post gap", n-1
}
}

//******************************************
objref rddelay,rdamp
rddelay= new Random()
rddelay.uniform(-20,50)
rdamp= new Random()
rdamp.normal(0.6, 0.0006)   //0.5% het Iinj and mean Iinj is 0.6nA 

objref stim [200]
for i=0,199{
r=rddelay.repick()
a=abs(rdamp.repick())
//print i, "delay this time is   ", r, "amplitude is  ",a
Bcell[i].soma stim[i] = new IClamp(0.5)
stim[i].del = r
stim[i].dur = 1000
stim[i].amp = a
}
//*******************************************************************************************************
//*********************************Print out Net cons***************************************************
strdef strvar
objref nfile
nfile = new File()

proc saveNet(){ local i
nfile.wopen("N_connect_50GJ_30syn_74mV tonicspill.txt")
	nfile.printf("Precell \tps.tcell \t Synapse \n")
	for i=0, nclist.count-1 {
	nfile.printf("%s\t%s\t%s\n", nclist.object[i].precell, nclist.object[i].postcell, nclist.object[i].syn)}
nfile.printf("TO BC\tBC")
for i= 0, nbcell-1 {nfile.printf("%d\t \n",  vbc2bc.x[i])}
//print "BCout"
nfile.close("N_connect_50GJ_30syn_74mV tonicspill.txt")
}

//**************************************************************************

objref  VmT
objref VmMat[nbcell]
VmT = new Vector()
for i=0, nbcell-1 {
	VmMat[i] = new Vector()
	}

proc VecMx() { local i
	VmT.append(t)
	for i=0, nbcell-1 {
		VmMat[i].append( cells.object[i].soma.v(0.5))
		}
	}
objref Spike[nbcell]
for i=0, nbcell-1 {
	Spike[i] = new Vector()
	}
strdef Spkstr
objref dfile
dfile = new File()

proc SpkMx() { local i, j
	k = 0
     	for i=0, nbcell-1 {
		Spike[i].spikebin(VmMat[i], 0)
		}
dfile.wopen("Spike50gaps_30syn_74mV tonicspill.txt")

	while(k <  VmT.size) {
	for j = 0, nbcell-1 {
	if(Spike[j].x[k] != 0) {
	dfile.printf("%f\t%d\n", VmT.x[k], j)}
	}
	k +=1 }
dfile.close("Spike50gaps_30syn_74mV tonicspill.txt")
	}



strdef strmat
objref efile
efile = new File()

	efile.wopen("Mempot50gaps_30syn_74mV tonicspill.txt")
	efile.printf("t\t")
for i = 0, 199 {
	efile.printf("%s\t", Bcell[i])}
	efile.printf("\n")
efile.close("Mempot50gaps_30syn_74mV tonicspill.txt")

proc sMatrix(){ local  i
	efile.aopen("Mempot50gaps_30syn_74mV tonicspill.txt")
	efile.printf("%f\t", t)
	for i = 0, 199 {
	efile.printf("%f\t", Bcell[i].soma.v(0.5))}
	efile.printf("\n")
	efile.close("Mempot50gaps_30syn_74mV tonicspill.txt")

}


strdef strmat
objref gfile
gfile = new File()

gfile.wopen("gapctesthz_50GJ_30syn_74mV tonicspill.txt")
	gfile.printf("t\t")
for i = 50, 100 {
	gfile.printf("%s\t", Gap[i])}
	gfile.printf("\n")
gfile.close("gapctesthz_50GJ_30syn_74mV tonicspill.txt")

proc gMatrix(){ local  i
	gfile.aopen("gapctesthz_50GJ_30syn_74mV tonicspill.txt")
	gfile.printf("%f\t", t)
	for i=50, 100 {
	gfile.printf("%f\t", Gap[i].i)}
	gfile.printf("\n")
	gfile.close("gapctesthz_50GJ_30syn_74mV tonicspill.txt")
	}

objref r_plt
proc initrPlt() {
	r_plt = new Graph(0)
	r_plt.size(0, tstop,0, cells.count)
	r_plt.label(0.95, 0.02, "ms")
	r_plt.label(0.01, 0.82, "neu")
	r_plt.view(0,0, tstop, cells.count,320,20,300,230)
}
 initrPlt()

proc plotAP() { local i, a
	a=1
 	r_plt.erase()
	while(j < cells.count-1) {
	for i = 0, VmT.size-1 {
	if (Spike[j].x[i] == 1) {
	r_plt.mark(VmT.x[i], j, "T", 5, a, 1)}}
	j+=1}
	r_plt.flush()
	}


proc init() { local dtsav, temp, secsav
finitialize(v_init)
t = -1000
dtsav = dt
secondorder =0
dt= 10
	// if cvode is on, turn it off to do large fixed step
temp= cvode.active()
if (temp!=0) {cvode.active(0)}
while(t<-100) { fadvance() print t}
	//restore cvode if reqd
if (temp!=0) {cvode.active(1)}
dt = dtsav
secondorder =2
t = 0
if (cvode.active()){
cvode.re_init()
}else{
fcurrent()
}
//frecord_init()
}
proc continuerun() {local rt
	eventcount =0
	eventslow =1
	stoprun =0
	if (using_cvode_) {
	cvode.event($1)
	}
	while(t < $1 && stoprun == 0) {
	step()
	sMatrix()
	gMatrix()
	VecMx()
	rt = stopsw()
	if (rt > realtime) {
		realtime = rt
		if (!stdrun_quiet) fastflushPlot()
		doNotify()
		if (realtime == 2 && eventcount > 50) {
			eventslow = int(eventcount/50)+1
		}
		eventcount = 0
	}else{
		eventcount = eventcount +1
		if ((eventcount%eventslow) == 0) {
			doEvents()
		}
	}
	}
	flushPlot()
}

objectvar save_window_, rvp_
objectvar scene_vector_[4]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}

{
xpanel("RunControl", 0)
v_init = -60
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 0
xvalue("t","t", 2 )
tstop = 1000	//1500
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.1
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 10	//40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
xpanel(544,121)
}
{
save_window_ = new Graph(0)
save_window_.size(0,tstop,-80,40)
scene_vector_[2] = save_window_
{save_window_.view(0, -80, tstop, 120, 290, 470, 579.84, 208)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addexpr("Bcell[1].soma.v(0.5)",1,1)
save_window_.addexpr("Bcell[10].soma.v(0.5)",2,1)
save_window_.addexpr("Bcell[20].soma.v(0.5)",3,1)
save_window_.addexpr("Bcell[30].soma.v(0.5)",4,1)
save_window_.addexpr("Bcell[40].soma.v(0.5)",5,1)
save_window_.addexpr("Bcell[50].soma.v(0.5)",6,1)
save_window_.addexpr("Bcell[60].soma.v(0.5)",7,1)
save_window_.addexpr("Bcell[70].soma.v(0.5)",8,1)
save_window_.addexpr("Bcell[80].soma.v(0.5)",9,1)
save_window_.addexpr("Bcell[90].soma.v(0.5)",2,1)
save_window_.addexpr("Bcell[100].soma.v(0.5)",3,1)
save_window_.addexpr("Bcell[120].soma.v(0.5)",4,1)
save_window_.addexpr("Bcell[140].soma.v(0.5)",5,1)
save_window_.addexpr("Bcell[160].soma.v(0.5)",6,1)
save_window_.addexpr("Bcell[180].soma.v(0.5)",7,1)
save_window_.addexpr("Bcell[199].soma.v(0.5)",8,1)
}
save_window_ = new Graph(0)
save_window_.size(0,1000,-2,2)
scene_vector_[3] = save_window_
{save_window_.view(0, -2, tstop, 120, 290, 470, 579.84, 208)}
graphList[1].append(save_window_)
save_window_.save_name("graphList[1].")
save_window_.addexpr("Gap[10].i",2,1)

save_window_ = new Graph(0)
save_window_.size(0,1000,-2,2)
scene_vector_[3] = save_window_
{save_window_.view(0, -2, tstop, 120, 290, 470, 579.84, 208)}
graphList[1].append(save_window_)
save_window_.save_name("graphList[1].")
save_window_.addexpr("Bcell[10].bcdend4[0].igaba_tonic",2,1)
save_window_.addexpr("Bcell[10].soma.igaba_tonic",3,1)


proc rrun(){
initBcell()
saveNet()
run()
SpkMx()
}
rrun()
plotAP()
objectvar scene_vector_[1]
{doNotify()}