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




numaxon=1
numsoma=1
numbasal=52
numapical=70
numtrunk=49

xopen("geo5038804.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 =  .025
AXONM = 5
gkdr = 0.01
celsius = 35.0  
KMULT =  0.03
KMULTP = 0.03
ghd=0.00005

objref g, b,c, stim0, stim1, stim2, stim, p,  nc0, nc1, nc2, ncstim, syn


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}

access soma

freq=50
geom_nseg()
tot=0
forall {tot=tot+nseg}
distance()

tstop=40000

b = new VBox()
b.intercept(1)
g = new Graph()
g.size(0,tstop,-65.1,-63.5)
g.addvar("apical_dendrite[65].v(0.5)",3,1,2*tstop,0,2)
g.label(0.3,0.1,"dendritic membrane potential")
g.exec_menu("10% Zoom out")
g.xaxis(1)
c = new Graph()
c.size(0,tstop,0,4)
c.xaxis(1)
c.addexpr("syn.Np", 2, 1)
c.addexpr("syn.Nd", 4, 1)
c.exec_menu("10% Zoom out")
b.intercept(0)
b.map()


rel=0.5
soma {

stim0= new NetStim(.5)
stim0.number=1
stim0.interval=10
stim0.start=250


stim= new NetStim(.5)
stim.number=50
stim.interval=200
stim.start=1100

stim1= new NetStim(rel)
stim1.number=4
stim1.interval=10
stim1.start=1100

ncstim = new NetCon(stim, stim1, 0, 0, 0.001)

stim2 = new NetStim(.5)
stim2.number =1
stim2.interval =10
stim2.start = 35000

}

forsec "axon" {   
                insert nax gbar_nax=gna * AXONM
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
}

forsec "soma" {   
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
}

for i=0, numbasal-1 dendrite[i] {
		insert hd ghdbar_hd=ghd vhalfl_hd=-73
                insert na3 gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { xdist = distance(x)
                	ghdbar_hd(x) = ghd*(1+3*xdist/100)
                		if (xdist > 100){
					vhalfl_hd=-81
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-73
                        		gkabar_kap(x) = KMULTP*(1+xdist/100)
               				}
		}
}
                
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

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

numap=65
apical_dendrite[numap]{ 
syn= new ltpltd(.5)
}



nc0= new NetCon(stim0, syn, 0, 0, 0.8e-3)
nc1 = new NetCon(stim1, syn, 0, 0, 0.8e-3)
nc2= new NetCon(stim2, syn, 0, 0, 0.8e-3)




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

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

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


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

xpanel(" ")
xlabel("  ")
xradiobutton("TBS LTP","setTBSLTP()",1)
xradiobutton("TBS LTD","setTBSLTD()")
xradiobutton("HFS LTP (50 Hz)","setHFSLTP()")
xradiobutton("LFS LTD (3 Hz)","setLFSLTD()")
xlabel("  ")
xbutton("run", "run()")
xpanel()

 

proc setTBSLTP() {
ncstim.weight=0.001
stim.number=50
stim.interval=200
stim1.number=4
stim1.interval=10
}

proc setTBSLTD() {
ncstim.weight=0.001
stim.number=10
stim.interval=1000
stim1.number=2
stim1.interval=10
}


proc setHFSLTP() {
ncstim.weight=0
stim1.number =500
stim1.interval =20	
}

proc setLFSLTD() {
ncstim.weight=0
stim1.number =30
stim1.interval =333	
}