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.025
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=80
rel=48
soma {
stim= new IClamp(0.5)
stim.amp=0.8
stim.dur=1
stim.del=20
}

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
}

objref ncond, netd, synd
apical_dendrite[rel] {
	netd = new NetStim(0.5)
	netd.number=1
	netd.interval=10
	netd.start=30
	print distance(1)
	synd = new Exp2Syn(1)
	synd.e=0
	synd.tau1=3
	synd.tau2=3
}

ncond= new NetCon(netd,synd,0,0,0.25e-3) 
delay=5
b = new VBox()
b.intercept(1)
g = new Graph()
g.size(0,tstop,-70,-10)
g.color(1)
g.exec_menu("Keep Lines")
g.color(1)
g.label(0.6,0.5,"AP->syn")
g.color(2)
g.label(0.6,0.45,"syn->AP")
xpanel("Fig 4",1)
xbutton("run", "runc()")
xvalue("delay", "delay")
xpanel()
b.intercept(0)
b.map("stdp", 750,270,300,250)


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

proc runc() {
	g.erase()
	netd.start=20+delay	
	color=1
	g.addvar("apical_dendrite[rel].v(.5)",color,2,2*tstop,0,2)
	run()
	color=2
	netd.start=20-delay	
	g.addvar("apical_dendrite[rel].v(.5)",color,2,2*tstop,0,2)
	run()
}