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

dist=7
rel=0.5

numbasal=70

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

Rm = 28000
Cm    = 1
RaAll= 150
RaAx = 50

Vrest = -65
dt = 0.1
gna =  .045
AXONM = 3
gkdr = 0.06
celsius = 35.0  
ka =  0.04
ghd=0.00005
gkm=0.02
gcat=0.0001
gahp=0.00001

objref g, b,c, stim, s, rsyn, nc, m, gt, gk

forsec "axon" {insert pas e_pas=Vrest g_pas = 1/Rm Ra=RaAx cm=Cm}
forsec "soma" {insert pas e_pas=Vrest g_pas = 1/Rm Ra=RaAll cm=Cm}
forsec "dendrite"{insert pas e_pas=Vrest g_pas = 1/Rm Ra=RaAll cm=Cm}
forsec "user5" {insert pas e_pas=Vrest g_pas = 1/Rm Ra=RaAll cm=Cm}

access soma

freq=50
geom_nseg()
distance()

tstop=500

b = new VBox()
b.intercept(1)
g = new Graph()
g.size(0,tstop,-70,20)
g.xaxis(1)
g.exec_menu("10% Zoom out")
g.addexpr("v_soma","soma.v(rel)",1,1, 0.2,0.9,2)
xpanel("",1)
xbutton("fig.4(i)", "runi()")
xbutton("fig.4(ii)", "runii()")
xbutton("fig.4(iii)", "runiii()")
xbutton("fig.4(iv)", "runiv()")
xpanel()
b.intercept(0)
b.map()

forsec "axon" {   
		insert km gbar_km=gkm
                insert nax gbar_nax=gna*AXONM	
                insert kdr gkdrbar_kdr=gkdr*AXONM
                insert kap gkabar_kap = ka
}

forsec "soma" {  
		insert km gbar_km=gkm
		insert hd ghdbar_hd=0	
                insert na3 gbar_na3=gna*AXONM
                insert kdr gkdrbar_kdr=gkdr*AXONM
                insert kap gkabar_kap = ka
                insert cat  gcatbar_cat=gcat
		insert cacum tau_cacum=100 depth_cacum=diam/2
		insert KahpM95 gbar_KahpM95 = gahp 
}

for i=0, numbasal-1 dendrite[i] {
                insert na3 gbar_na3=0
                insert cat  gcatbar_cat=gcat
		insert cacum tau_cacum=100 depth_cacum=diam/2
		insert KahpM95 gbar_KahpM95 = gahp 
}

                
forsec "apical_dendrite" {
		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
		insert cacum tau_cacum=100 depth_cacum=diam/2
		insert KahpM95 gbar_KahpM95 = gahp 
                insert cat  gcatbar_cat=gcat

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

forsec "user5" {
		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
		insert cacum tau_cacum=100 depth_cacum=diam/2
                insert cat  gcatbar_cat=gcat
		insert KahpM95 gbar_KahpM95 = gahp 

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

proc init() {
	t=0
        forall {
        v=Vrest
        if (ismembrane("nax") || ismembrane("na3")) {ena=55}
        if (ismembrane("KahpM95") || 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)}
	if (ismembrane("cat")){e_pas(x)=e_pas(x)+ica(x)/g_pas(x)}
		}
	}
	cvode.re_init()
	cvode.event(tstop)
	access soma
	g.begin()
}

objref playvec1, playvec2
playvec1 = new Vector()
playvec1.append(0, 50, 120)
playvec2 = playvec1.c.indgen()
playvec1.play("f($1)", playvec1, playvec2, 1)
proc f() {
	if (t < 50) {
		axon.gbar_km=factor*gkm
		soma.gbar_km=gkm
	}else if (t>=50 && t<120) {
		axon.gbar_km=factor*gkm-gkm*factor*(t-50)/70 
		soma.gbar_km=gkm-gkm*(t-50)/70 
	}else{
		axon.gbar_km=0
		soma.gbar_km=0
	}
}

proc advance() {
	fadvance()
	g.plot(t)
	g.flush()
	doNotify()
}

proc runi() {
factor=2
soma.gbar_km=gkm
axon.gbar_km=gkm*factor
run()
}

proc runii() {
factor=3
soma.gbar_km=gkm
axon.gbar_km=gkm*factor
run()
}

proc runiii() {
factor=4
soma.gbar_km=gkm
axon.gbar_km=gkm*factor
run()
}

proc runiv() {
factor=6
soma.gbar_km=gkm
axon.gbar_km=gkm*factor
run()
}