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

dist=1
rel=0.1

numaxon=1
numsoma=1
numbasal=51
numapical=42
numtrunk=49

xopen("pc2b.hoc")             // geometry file
xopen("fixnseg.hoc")           

Rm = 14600
RmDend = Rm/1
RmSoma = Rm
RmAx = Rm

Cm    = 3.2
CmSoma= Cm
CmAx  = Cm
CmDend = Cm*1

RaAll= 50
RaSoma=50  
RaAx = 50

Vrest = -62
dt = 0.1
gna =  .024
AXONM = 5
gkdr = 0.008
celsius = 35.0  
KMULT =  0.004
ghd=20e-6

objref g, b,c, f,time, stim, distrx, distry, outfile, cdistry
objref distri, p, s, rsyn, sref, blist[numtrunk], str1, str2, aplist
strdef filename, dends, trunk

for i=0, numtrunk-1 {
blist[i] = new Vector()
}

aplist = new Vector(numapical)

forsec "axon" {insert pas e_pas=Vrest g_pas = 1/RmAx Ra=RaAx cm=CmAx}
forsec "soma" {insert pas e_pas=Vrest g_pas = 1/RmSoma Ra=RaSoma cm=CmSoma}
forsec "dend"{insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}
forsec "dend_5" {insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}
forsec "apic" {insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}

access soma

freq=100
geom_nseg()
access soma
distance()

//*********mapping bifurcations******************


for i=0, numapical-1 apic[i] {
	while (!issection("dend_5.*")) {
	sref = new SectionRef()
	access sref.parent
	sprint(dends, secname())
	}
	for k=0, numtrunk-1 dend_5[k] {
	sprint(trunk,secname())
	x=strcmp(dends, trunk)
	if (x==0) {blist[k].append(i)  aplist.x[i]=k}
	}
}
//************************************************


tstim=1
strength = 1 /*namps*/
tstop=550

dend_5[10] {
stim= new IClamp(.5)
stim.amp=-1
stim.dur=200
stim.del=50
}

forsec "axon" {   
                insert nax gbar_nax=gna*AXONM
                insert kdr gkdrbar_kdr=gkdr*AXONM
                insert kap gkabar_kap = KMULT*AXONM
}

forsec "soma" {   
		insert hd ghdbar_hd=ghd	vhalfl_hd=-75
                insert na3 gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULT
}

for i=0, numbasal-1 dend[i] {
		insert hd ghdbar_hd=ghd vhalfl_hd=-75 
                insert na3 gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x,0) { xdist = distance(x)
                		if (xdist > 100){
					vhalfl_hd=-80
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-75
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}
}
                
forsec "dend_5" {
		insert hd ghdbar_hd=ghd
                insert na3 gbar_na3=gna 
                insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x,0) { xdist = distance(x)
                	ghdbar_hd(x) = ghd*(1+100/(1+exp((256-xdist)/26)))
                	if (xdist>500) {xdist=500}
                		if (xdist > 100){
					vhalfl_hd=-80
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-75
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}
}

for i=0, numapical-1 apic[i] {
	insert hd 
        insert na3 gbar_na3=gna 
        insert kdr gkdrbar_kdr=gkdr
	insert kap
	insert kad
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
                		if (xdist > 100){
					vhalfl_hd=-80
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-75
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}

	vhalfl_hd = dend_5[aplist.x[i]].vhalfl_hd
       	ghdbar_hd = dend_5[aplist.x[i]].ghdbar_hd(1)

}

proc init() {
	t=0

forsec "axon" {gkabar_kap = KMULT*AXONM}

forsec "soma" {   
	ghdbar_hd=ghd
        gkabar_kap = KMULT
}

for i=0, numbasal-1 dend[i] {
	ghdbar_hd=ghd 
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
                		if (xdist > 100){
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}

}
                
forsec "dend_5" {
	ghdbar_hd=ghd
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
                	ghdbar_hd(x) = ghd*(1+100/(1+exp((256-xdist)/26)))
	               	if (xdist>500) {xdist=500}
                		if (xdist > 100){
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}
}

for i=0, numapical-1 apic[i] {
        gbar_na3=gna
        gkdrbar_kdr=gkdr
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
	               	if (xdist>500) {xdist=500}
                		if (xdist > 100){
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}

       	ghdbar_hd = dend_5[aplist.x[i]].ghdbar_hd(1)
}

        forall {
        v=Vrest
        if (ismembrane("nax") || ismembrane("na3")) {ena=55}
        if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {ek=-90}
        if (ismembrane("hd") ) {ehd_hd=-30}
	}
	finitialize(Vrest)
        fcurrent()

        forall {
	for (x) {
	if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)}
	if (ismembrane("hd")) {e_pas(x)=e_pas(x)+i_hd(x)/g_pas(x)}
		}
	}
	cvode.re_init()
	cvode.event(tstop)
	access soma
}

load_file("rebound.ses")

proc runc() {
KMULT =  0.004
ghd=20e-6
def()
run()
}

proc run4ap() {
KMULT =  0.002
ghd=20e-6
def()
run()
}

proc run4apzd() {
KMULT =  0.002
ghd=0
def()
run()
}


proc def() {
forsec "axon" {gkabar_kap = KMULT*AXONM}

forsec "soma" {   
	ghdbar_hd=ghd
        gkabar_kap = KMULT
}

for i=0, numbasal-1 dend[i] {
	ghdbar_hd=ghd 
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
                		if (xdist > 100){
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}

}
                
forsec "dend_5" {
	ghdbar_hd=ghd
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
                	ghdbar_hd(x) = ghd*(1+100/(1+exp((256-xdist)/26)))
	               	if (xdist>500) {xdist=500}
                		if (xdist > 100){
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}
}

for i=0, numapical-1 apic[i] {
        gbar_na3=gna
        gkdrbar_kdr=gkdr
	gkabar_kap=0
	gkabar_kad=0

		for (x,0) { xdist = distance(x)
	               	if (xdist>500) {xdist=500}
                		if (xdist > 100){
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
                        		gkabar_kap(x) = KMULT*(1+xdist/100)
               				}
		}

       	ghdbar_hd = dend_5[aplist.x[i]].ghdbar_hd(1)
}

}