load_file("nrngui.hoc")
cvode_active(1)

Vrest = -64
dt = 0.05
celsius = 35.0  
freq=20

numaxon=1
numsoma=1
numbasal=52
numapical=81

Rm = 25370
Cm    = 1.41
RaAll= 150 
AXONM = 5

gna =  .0
gkdr = 0.0
KMULT =  0.0
gkm=0
gkd=0.0
gc=0
gcal=gc
gcan=gc
gcat=gc
gKc=0
gahp=0
ghd=0.00001

objref rc, rf
use_mcell_ran4()
lowindex = mcell_ran4_init(1)
rc = new Random()
rc.MCellRan4(12345)
rc.uniform(0,80)
rf = new Random()
rf.MCellRan4(10946353)
rf.uniform(0,1)


tstop=1000

xopen("geo-cell1zr.hoc")
xopen("fixnseg.hoc")           

nsyn=30

objref stim, time, y, y2, infile, ifile, currt, curr, syngaba[nsyn], synglu[nsyn]
objref ncgaba[nsyn], ncglu[nsyn], s, apc, spikes

spikes = new Vector()
access soma

apc= new APCount(.5)
apc.thresh=0

distance()

s = new NetStim(.5)
s.interval =50
s.number=10
s.noise=1
s.start=20

low=0
high=100

soma {

for k=0, nsyn-1 {
	syngaba[k] = new Exp2Syn(1)
	syngaba[k].e=-44
	syngaba[k].tau1=1.4
	syngaba[k].tau2=27

/*
	syngaba[k].e=-35
	syngaba[k].tau1=1
	syngaba[k].tau2=10
*/

	synglu[k]= new Exp2Syn(1)
	synglu[k].e=0
	synglu[k].tau1=.4
	synglu[k].tau2=4.8


}
}




for i=0, nsyn-1 { 
    flag=0
		while (flag==0) {
		comp=int(rc.repick()+0.5)
		tmp=rf.repick()
		apical_dendrite[comp] {if (distance(tmp)<low || distance(tmp)>high) {flag=0} else{flag=1}}

}

    apical_dendrite[comp] {
	syngaba[i].loc(tmp)
	print i, comp, tmp, distance(tmp)
//	ppp.printf("exc %d %d %g %g %g\n", 1, i, comp, tmp, distance(tmp))	
	}
}

soma {
	stim = new IClamp(0.5)
	stim.amp=0
	stim.dur=400
	stim.del=50
}


forall {insert pas area(.5)}

forsec "dendrite" { 
	insert ds
	insert hd 
        insert na3 
        insert kdr 
	insert kap 
	insert cacum depth_cacum=diam/2
        insert cal 
        insert can 
        insert cat 
	insert cagk  
	insert KahpM95 
}

forsec "soma" { 
	insert ds
	insert hd 
        insert na3 
        insert kdr 
	insert kap 
	insert km
	insert kd
	insert cacum depth_cacum=diam/2
        insert cal 
        insert can 
        insert cat 
	insert cagk  
	insert KahpM95 
}

forsec "axon" {   
	insert na3 
        insert kdr 
        insert kap 
}

        forall {v=Vrest e_pas=Vrest g_pas = 1/Rm Ra=RaAll cm=Cm ek=-90 ena=55}
	geom_nseg()
	distance()
	tot=0
	forall {tot=tot+nseg}
	maxdist=0
	forall for(x) {if (distance(x)>maxdist) {maxdist=distance(x)}}

forsec "axon" Ra=RaAll/3
forall if(ismembrane("hd")) {ehd_hd=-30}

load_file("ec.ses")

proc init() {
	access soma
        forall {
		v=Vrest e_pas=Vrest
	        if (ismembrane("cal")) {
                gcalbar_cal=gc
                gcanbar_can=gc
                gcatbar_cat=gc
		gbar_cagk= gKc 
		gbar_KahpM95 = gahp 
		}
	}

forsec "axon" {   
	gbar_na3=gna*AXONM 
        gkdrbar_kdr=gkdr 
        gkabar_kap = KMULTP  sh_kap=0
}

forsec "soma" {   
	ghdbar_hd=ghd
        gbar_na3=gna  
        gkdrbar_kdr=gkdr 
        gkabar_kap = KMULTP 
	gbar_km= gkm
	gkdbar_kd = gkd
}

for i=0, numbasal-1 dendrite[i] {
	ghdbar_hd=ghd
        gbar_na3=gna 
        gkdrbar_kdr=gkdr
	gkabar_kap=KMULTP
}
                
forsec "apical_dendrite" {
	ghdbar_hd=ghd
        gbar_na3=gna 
        gkdrbar_kdr=gkdr
	gkabar_kap=KMULTP

}


	finitialize(v)
        fcurrent()
	finitialize(v)
        forall for(x) {
	if (ismembrane("cal")) {e_pas(x)=v(x)+(i_hd(x)+ina(x)+ik(x)+ica(x))/g_pas(x)
		} else {
		e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)
		}
	}
	cvode.re_init()
}

//proc fig9c() {
gna =  .022
gkdr = 0.01
KMULT =  0.02
gc=1.e-5
gKc=0
gkm=0.017
gkd=0.0
gahp=0.0
gcal=gc
gcan=gc
gcat=gc
stim.amp=0
for k=0, nsyn-1 {
	ncglu[k] = new NetCon(s, synglu[k], 0, 0, 5.5e-3)
	ncgaba[k] = new NetCon(s, syngaba[k], 0, 0, 0e-3)
}
for nsim=1,10 {
/*
for i=0, nsyn-1 { 
    flag=0
		while (flag==0) {
		comp=int(rc.repick()+0.5)
		tmp=rf.repick()
		apical_dendrite[comp] {if (distance(tmp)<low || distance(tmp)>high) {flag=0} else{flag=1}}

}

    apical_dendrite[comp] {
	syngaba[i].loc(tmp)
	print i, comp, tmp, distance(tmp)
//	ppp.printf("exc %d %d %g %g %g\n", 1, i, comp, tmp, distance(tmp))	
	}
}
*/
	run()
	print " sim ",nsim, apc.n, " spikes"
	spikes.append(apc.n)
}
spikes.printf()
print " final glu: ", spikes.mean(), spikes.stdev(), spikes.stderr()