load_file("nrngui.hoc")
cvode.active(1)

mid=12
dist=27
secondorder=2
FARADAY=96520
PI=3.14159

Rm = 28000
RmDend = Rm/2
RmSoma = Rm
RmAx = Rm

Cm    = 1
CmSoma= Cm
CmAx  = Cm
CmDend = Cm*2

RaAll= 150
RaSoma=150  
RaAx = 50

Vrest = -65
dt = 0.025
gna =  .032
AXONM = 2
gkdr = 0.01
celsius = 34.0  
KMULT =  0.048
KMULTP = 0.048

xopen("n160_mod.nrn")             // geometry file

tstim=90
strength = 0.15 /*namps*/
tstop=90

objref g, b
b = new VBox()
b.intercept(1)
g = new Graph()
g.size(0,tstop,-70,30)
g.addvar("soma[5].v(0.5)",2,1, 2*tstop,0,2)
g.addvar("apical[mid].v(0.5)",3,2,2*tstop,0,2)
g.addvar("apical[dist].v(0.5)",4,2,2*tstop,0,2)
g.color(2)
g.label(0.1,0.05,"soma")
g.color(3)
g.label(0.24,0.05,"mid apical (200um)")
g.color(4)
g.label(0.6,0.05,"dist apical (400um)")
g.xaxis(1)
g.begin()
xpanel("")
xbutton("run Fig.1A", "runu()")
xpanel()
b.intercept(0)
b.map()



soma[5] distance()

access soma[5]

axon[1] {
                insert nax    gbar_nax=gna * AXONM ena=55
                insert kdr gkdrbar_kdr=gkdr ek=-90
                insert pas e_pas=Vrest g_pas = 1/RmAx Ra=RaAx cm=CmAx
                insert kap gkabar_kap = KMULTP*0.2 ek=-90
}

axon[0] {   
                insert nax  gbar_nax=gna * AXONM ena=55
                insert kdr  gkdrbar_kdr=gkdr ek=-90
                insert pas e_pas=Vrest g_pas = 1/RmAx Ra=RaAx cm=CmAx
                insert kap gkabar_kap = KMULTP*0.2 ek=-90
}

for i=0,NumSoma-1 soma[i] {   
                insert na3  gbar_na3=gna  ena=55
                insert kdr gkdrbar_kdr=gkdr ek=-90
		ar2_na3=0.8
                insert kap gkabar_kap = KMULTP ek=-90
                insert pas e_pas=Vrest g_pas = 1/RmSoma Ra=RaSoma cm=CmSoma
}

for i=0,NumBasal-1 basal[i] {
                insert na3    gbar_na3=gna ena=55
                insert kdr gkdrbar_kdr=gkdr ek=-90
		ar2_na3=1
                insert kap gkabar_kap=KMULTP ek=-90
                insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend
}
                
for i=0,NumApical-1 apical[i] {
//if (i<40) {printf(" %d, %g \n",i,distance(0.5))}
              insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll  cm=CmDend
              for (x){ xdist = distance(x)
			}
			if (diam>0.5 && distance(0.5)<500) {
                insert na3 ena=55
		ar2_na3=0.5
		gbar_na3=gna
                insert kdr ek=-90
		gkdrbar_kdr=gkdr
		insert kap ek=-90
		insert kad ek=-90
		gkabar_kap=0
		gkabar_kad=0

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

proc init() {
	t=0
        forall {v=Vrest}
	finitialize(v)
        fcurrent()

        forall {
	  for (x) {if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)
		   } else {
			e_pas(x)=v(x)
		   }
		  }
	}
	cvode.re_init()
}

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

proc run() {

	t=0

	access soma[5]
        nstim=1
        fstim(nstim)
        fstim(0, 0.5, 2, tstim, strength)
	init()
        access soma[5]

	while(t<tstop) { step()}
}


proc runu() {
run()
}