load_file("nrngui.hoc")
cvode.active(0)
mid=13
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

gsyn=0
tstim=1.2
strength = 1 /*namps*/
tstop=90
npulses=10

objref g, b, syn[npulses]

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

access soma[5]
distance()

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] {
              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)
                	}
              				}
							}
				}

apical[mid] {
        for j=0,npulses-1 {
        syn[j] = new AlphaSynapse(0.5)
        syn[j].onset=10+10*j
        syn[j].gmax=gsyn
        syn[j].tau=3
        syn[j].e = 0
        }
}


b = new VBox()
b.intercept(1)
g = new Graph()
g.size(0,tstop,-70,30)
g.color(2)
g.label(0.6,0.7,"40*KA inactiv. at 250um")
g.color(3)
g.label(0.6,0.05,"v apical at 250um")
g.xaxis(1)
g.begin()
xpanel("")
xbutton("run Fig.1C", "runu()")
xpanel()
b.intercept(0)
b.map()


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() {

	init()
	t=0

access soma[5]

        nstim=3
        fstim(nstim)
	for i=0,nstim-1 {
        fstim(i, 0.5, 2+i*40, tstim, strength)
	}

apical[mid] {
        for j=0,npulses-1 {
        syn[j].gmax=gsyn
        }
}

        access soma[5]

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


proc runu() {

gsyn=0
g.addexpr("apical[mid].l_kad(0.5)*40",2,1, 1.2*tstop,0,2)
g.addvar("apical[mid].v(0.5)",3,1,1.2*tstop,0,2)
gsyn=4000e-6
//cvode = new CVode(1)
run()
g.exec_menu("Keep Lines")
g.addexpr("apical[mid].l_kad(0.5)*40",2,2, 1.2*tstop,0,2)
g.addvar("apical[mid].v(0.5)",3,2,1.2*tstop,0,2)
g.begin()
gsyn=0
run()
}