load_file("nrngui.hoc")
cvode_active(1)
numbasal=49
numapical=52
numtrunk=53

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

Rm = 28000
RmDend = Rm
RmSoma = Rm
RmAx = Rm

Cm    = 1
CmSoma= Cm
CmAx  = Cm
CmDend = Cm

RaAll= 150
RaSoma=150  
RaAx = 50

Vrest = -65
dt = 0.1
gna =  .03
gnad = gna/2
AXONM = 5
gkdr = 0.01
celsius = 35.0  
KMULT =  0.03
KMULTD = KMULT
ghd=0.00005
nash=0

objref g, b,c, stim, distrx, distry, outfile, cdistry, vsoma,vdend,xgamma,x2gamma
objref distri, p, s, rsyn, nc, sref, blist[numtrunk], str1, str2

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 "dendrite"{insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}
forsec "user5" {insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}

freq=50
geom_nseg()

access soma
distance()
tstop=20

b = new VBox()
b.intercept(1)
g = new Graph()
g.size(0,tstop,-70,30)
g.addvar("soma.v(0.5)",1,1,2*tstop,0,2)
g.xaxis(1)
c = new Graph()
c.size(-300,740,0,100)
c.xaxis(1)
c.exec_menu("10% Zoom out")
c.color(1)
c.label(0.7,0.9," peak AP (trunk)")
c.color(2)
c.label(0.7,0.85," basal dendrites")
c.color(3)
c.label(0.7,0.8," oblique dendrites")
xpanel("")
xbutton("run ", "runp()")
xvalue(" KA dend", "KMULTD")
xvalue(" Na dend", "gnad")
xpanel()
b.intercept(0)
b.map()

p = new PlotShape()
p.exec_menu("Shape Plot")

soma {
stim= new IClamp(0.5)
stim.amp=0.8
stim.dur=1
stim.del=1
}
rel=0.5
forsec "axon" {   
                insert nax
                insert kdr
                insert kap
}

forsec "soma" {
		insert ds
		insert hd
        insert na3
        insert kdr
        insert kap
}

for i=0, numbasal-1 dendrite[i] {
		insert ds
		insert hd
        insert na3
        insert kdr
		insert kap
		insert kad
}
                
forsec "apical_dendrite" {
		insert ds
		insert hd
        insert na3
        insert kdr
		insert kap
		insert kad
}

forsec "user5" {
		insert ds
		insert hd
        insert na3
        insert kdr
		insert kap
		insert kad
}

proc init() {
	t=0
	conductances()
	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
	g.begin()
}


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

proc runp() {
run()
	c.erase()

	distrx=new Vector()
	distry=new Vector()
	cdistry=new Vector()
	distri=new Vector()
	forsec "apical_dendrite" {
		for (x,0) {
			distrx.append(distance(x)) 
			distry.append(vmax_ds(x)-Vrest)
		}
	}
	distry.mark(c,distrx,"O",3,3,2)

	distrx=new Vector()
	distry=new Vector()
	cdistry=new Vector()
	distri=new Vector()
	for i=0, numbasal-1 dendrite[i] {
		for (x,0) {
			distrx.append(-distance(x)) 
			distry.append(vmax_ds(x)-Vrest)
		}
	}
	distry.mark(c,distrx,"O",3,2,2)

	distrx=new Vector()
	distry=new Vector()
	cdistry=new Vector()
	distri=new Vector()
	forsec "user5" {
		for (x,0) {
			distrx.append(distance(x)) 
			distry.append(vmax_ds(x)-Vrest)
		}
	}
	forsec "soma" {
		for (x,0) {
			distrx.append(distance(x)) 
			distry.append(vmax_ds(x)-Vrest)
		}
	}
	distry.mark(c,distrx,"t",5,1,1)

	c.flush()
	doNotify()
}

proc conductances() {
forsec "axon" {   
    gbar_nax=gna * AXONM	
	sh_nax=nash
    gkdrbar_kdr=gkdr
    gkabar_kap = KMULT*1
}

forsec "soma" {   
	ghdbar_hd=ghd
	vhalfl_hd=-73
	ar_na3=1 
	sh_na3=nash 
	gbar_na3=gna
    gkdrbar_kdr=gkdr
    gkabar_kap = KMULT
}

for i=0, numbasal-1 dendrite[i] {
	ghdbar_hd=ghd
	vhalfl_hd=-73
	ar_na3=1 
	gbar_na3=gnad
	sh_na3=nash
    gkdrbar_kdr=gkdr
	gkabar_kap=0
	gkabar_kad=0

	for (x,0) { xdist = distance(x)
    	ghdbar_hd(x) = ghd*(1+3*xdist/100)
            if (xdist > 100){
				vhalfl_hd=-81
            	gkabar_kad(x) = KMULTD*(1+xdist/100)
                } else {
				vhalfl_hd=-73
                gkabar_kap(x) = KMULTD*(1+xdist/100)
               	}
		}
}
                
forsec "apical_dendrite" {
	ghdbar_hd=ghd
	vhalfl_hd=-73
	ar_na3=1 
	gbar_na3=gnad
	sh_na3=nash
    gkdrbar_kdr=gkdr
	gkabar_kap=0
	gkabar_kad=0

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

forsec "user5" {
	ghdbar_hd=ghd
	vhalfl_hd=-73
	ar_na3=1 
	gbar_na3=gnad
	sh_na3=nash
    gkdrbar_kdr=gkdr
	gkabar_kap=0
	gkabar_kad=0

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