// mitral.hoc
//
// An olfactory bulb mitral cell model displaying mixed-mode oscillations of 
// 	the membrane potential.
// From Rubin DB and Cleland TA (2006) "Dynamical mechanisms of odor processing 
// 	in olfactory bulb mitral cells."  J. Neurophysiology. 
// April 2006
//
// Based on the model of Davison AP, Feng J, Brown D (2000) Brain Res Bulletin 51(5):393-399
 

xopen("$(NEURONHOME)/lib/hoc/noload.hoc") // standard run tools
xopen("tabchannels.hoc")

// Set morphological parameters

Atotal		= 100000	// um2
Len		= 100		// um
RM		= 100000	// ohm.cm2
Erest		= -65		// mV
p 		= 0.051
q 		= 0.084
r 		= 0.328
gpg 		= 5.86e-5	// S.cm-2
gsp 		= 5.47e-5	// S.cm-2
gsd 		= 1.94e-4	// S.cm-2
alphas 		= 1.37
alphag		= 1.85

// Create cell

create soma, glom, prim, dend, s2d, s2p, p2g
access soma

// Connect the cell

soma connect s2p(0),0
s2p connect prim(0),1
prim connect p2g(0),1
p2g connect glom(0),1
soma connect s2d(0),1
s2d connect dend(0),1

// Create subsets 

objref real
proc subsets() {
	real = new SectionList()				// "real" sections that have currents
		soma real.append				
		prim real.append
		dend real.append
		glom real.append
}
subsets()

objref cvode
cvode  = new CVode(1)

// Morphology procedures

proc set_ra() {  			// the argument is the conductance (S.cm-2)
   Ra = (PI*1e4)/(4*Atotal) * ( 1/$1 )	// ohm.cm
}

proc set_size() {			// the argument is the membrane area (um2)
   diam = $1/(PI*Len)			// um
}

proc change_params() {

      Asoma = p*Atotal
      Aglom = q*Atotal
      Aprim = r*Atotal
      Adend = Atotal - Asoma - Aglom - Aprim
      soma {
        set_size(Asoma)
      }
      glom {
        set_size(Aglom)
	}
      prim {
        set_size(Aprim)
      }  
      dend {
        set_size(Adend)
      }
      s2d { set_ra(gsd) }
      s2p { set_ra(gsp) }
      p2g { set_ra(gpg) }
}

// Set cell properties

forsec "*2*" {
   L = 1
   diam = 1
}

soma {
   insert pas
   insert nafast
   insert kfasttab
   insert kslowtab
   insert kA
   insert kca3
   insert lcafixed
   insert cad
   insert kO
   depth_cad = 8
   L			= Len
   Ra			= 1e-7		
   e_pas 		= Erest		
   g_pas 		= 1/RM		
   gnabar_nafast 		= 1.3532	
   gkbar_kfasttab 	= 0.1956
   gkbar_kslowtab 	= 0.0028
   gkbar_kA 		= 0.012
   gkbar_kca3 		= 0.12
   gcabar_lcafixed 	= 0.0040
   gkbar_kO			= 0
   eO 			= -90
}

glom {
   insert pas
   insert kslowtab
   insert lcafixed
   insert cad
   depth_cad = 8
   L			= Len
   Ra			= 1e-7
   e_pas 		= Erest
   g_pas		= 1/RM
   gkbar_kslowtab	= 0.020
   gcabar_lcafixed	= 0.0095
}

prim {
   insert pas
   insert nafast
   insert kfasttab
   insert kslowtab
   insert kA
   insert kca3
   insert lcafixed
   insert cad
   depth_cad = 8
   insert NaP
   insert Ih
   insert kO
   L			= Len
   Ra			= 1e-7
   e_pas 		= Erest
   g_pas		= 1/RM
   gkbar_kfasttab	= 0.00123
   gnabar_nafast	= 0.00534
   gkbar_kslowtab	= 0.00174
   gcabar_lcafixed	= 0.004
   gkbar_kA 		= 0.00735
   gkbar_kca3 		= 0.1
   gcabar_lcafixed 	= 0.0040
   gbar_NaP             = 0.00042
   gkhbar_Ih            = 0.0024
   gkbar_kO 		= 0
   eO				= -90
}

dend {
   insert pas
   insert kfasttab
   insert nafast
   L			= Len
   Ra			= 1e-7
   e_pas 		= Erest
   g_pas		= 1/RM
   gkbar_kfasttab	= 0.0330
   gnabar_nafast	= 0.0226
}

// Set areas and link resistances

change_params()				

// Set reversal potentials

forall if (ismembrane("ca_ion")) {
      eca = 70	
	cai = 0.00001	
	cao = 2		
	ion_style("ca_ion",3,2,0,0,1)
}

forall if (ismembrane("na_ion")) {
	ena = 55	
}

forall if (ismembrane("k_ion")) {
	ek  = -90	
}

load_file("nrngui.hoc")

// Insert stimulus mechanisms 

objref stim1,stim2,syn

prim stim1 = new IClamp(0.5)
stim1.amp = 0
stim1.del = 0
stim1.dur = 0

prim stim2 = new IClamp(0.5)
stim2.amp = 0
stim2.del = 0
stim2.dur = 0

glom syn = new AlphaSynapse(1)
syn.onset = 0
syn.tau = 0
syn.gmax = 0
syn.e = 0

// Generate figures from the paper (GUI)

xpanel("Figures")
xbutton("Figure 2A Top Trace","fig2Atop()")
xbutton("Figure 2A Middle Trace","fig2Amiddle()")
xbutton("Figure 2A Bottom Trace","fig2Abottom()")
xlabel("")
xbutton("Figure 2B Top Trace","fig2Btop()")
xbutton("Figure 2B Second Trace","fig2Bsecond()")
xbutton("Figure 2B Third Trace","fig2Bthird()")
xbutton("Figure 2B Bottom Trace","fig2Bbottom()")
xlabel("")
xbutton("Figure 3A","fig3A()")
xbutton("Figure 3B","fig3B()")
xbutton("Figure 3C Top","fig3Ctop()")
xbutton("Figure 3C Bottom","fig3Cbottom()")
xlabel("")
xbutton("Figure 4A","fig4A()")
xlabel("")
xbutton("Figure 5B","fig5B()")
xlabel("")
xbutton("Figure 6A","fig6A()")
xbutton("Figure 6B","fig6B()")
xbutton("Figure 6C","fig6C()")
xbutton("Figure 6D","fig6D()")
xbutton("Figure 6E","fig6E()")
xlabel("")
xbutton("Figure 7A","fig7A()")
xbutton("Figure 7B","fig7B()")
xlabel("")
xbutton("Figure 8A","fig8A()")
xbutton("Figure 8B","fig8B()")
xbutton("Figure 8C","fig8C()")
xpanel(1,1)

objref g, scene_vector_[4]

proc grapher(){
	g = new Graph(0)
	scene_vector_[3] = g
	{g.view($1, $2, $3, $4, 400, 1, 700, 350)}
	graphList[0].append(g)
	g.save_name("graphList[0].")
	g.addexpr("Membrane Potential","soma.v(0.5)", 1, 1, 0.8, 0.9, 2)
}
tstop = 5000
steps_per_ms = 100
dt = 0.01
grapher(0,-80,tstop,140)

proc fig2Atop(){
	reset()
	tstop = 2610
	grapher(2510,-80,100,60)
	stim1.amp = 1.6
	stim1.dur = tstop
	run()
}
proc fig2Amiddle(){
	reset()
	tstop = 1365
	grapher(1265,-80,100,60)
	stim1.amp = 1.0
	stim1.dur = tstop
	run()
}
proc fig2Abottom(){
	reset()	
	tstop = 1200
	grapher(1100,-80,100,60)
	stim1.amp = 0.2
	stim1.dur = tstop
	run()
}
proc fig2Btop(){
	reset()
	tstop = 3100
	grapher(2600,-80,500,140)
	stim1.amp = 1.6
	stim1.dur = tstop
	run()
}
proc fig2Bsecond(){
	reset()
	tstop = 3100
	grapher(2600,-80,500,140)
	stim1.amp = 1.35
	stim1.dur = tstop
	run()
}
proc fig2Bthird(){
	reset()
	tstop = 3100
	grapher(2600,-80,500,140)
	stim1.amp = 1.2
	stim1.dur = tstop
	run()
}
proc fig2Bbottom(){
	reset()
	tstop = 3200
	grapher(2600,-80,500,140)
	stim1.amp = 0.2
	stim1.dur = tstop
	run()
}
proc fig3A(){
	reset()
	tstop = 1800
	grapher(1600,-80,200,60)
	stim1.amp = 0.25
	stim1.dur = tstop
	stim2.del = 1634
	stim2.amp = 0
	stim2.dur = 5
	run()
	g.exec_menu("Keep Lines")
	g.addexpr("Membrane Potential w/Input","soma.v(0.5)", 2, 1, 0.8, 0.9, 2)
	stim2.amp = -3.5
	run()
}
proc fig3B(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	stim1.amp = 0.25
	stim1.dur = tstop
	stim2.del = 2100
	stim2.amp = -5
	stim2.dur = 10
	run()
}
proc fig3Ctop(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	stim1.amp = 0.25
	stim1.dur = tstop
	syn.onset = 2100
	syn.e = -65
	syn.tau = 2.5
	syn.gmax = 5
	g.addexpr("Synaptic Current","(syn.i)/(glom.v(0.5) + 65.00001)", 2, 1, 0.8, 0.9, 2)
	run()
}
proc fig3Cbottom(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	stim1.amp = 0.25
	stim1.dur = tstop
	syn.onset = 2100
	syn.e = -65
	syn.tau = 2.5
	syn.gmax = 25
	g.addexpr("Synaptic Current","(syn.i)/(glom.v(0.5) + 65.00001)", 2, 1, 0.8, 0.9, 2)
	run()
}
proc fig4A(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	syn.onset = 2100
	syn.e = 0
	syn.tau = 50
	syn.gmax = 0.01
	g.addexpr("Synaptic Current","1000*(syn.i/glom.v(0.5))", 2, 1, 0.8, 0.9, 2)
	run()
}
proc fig5B(){
	reset()
	tstop = 1750
	grapher(1600,-100,150,160)
	stim1.amp = 1.15
	stim1.dur = tstop
	run()
	g.exec_menu("Keep Lines")
	g.addexpr("Membrane Potential w/Reduced K","soma.v(0.5)", 2, 1, 0.75, 0.9, 2)
	forall if (ismembrane("k_ion")) {
		ek  = -100	
		}
	stim1.amp = 2.5
	run()
}
proc fig6A(){
	reset()
	tstop = 4000
	grapher(3000,-80,1000,140)
	prim.gkbar_kA   = 0
	soma.gkbar_kA   = 0
	prim.gkbar_kO   = 0.000114
	prim.eO	    = -90
	soma.gkbar_kO   = 0.000219
	soma.eO	    = -90
	run()
}
proc fig6B(){
	reset()
	tstop = 4000
	grapher(3000,-80,1000,140)
	prim.gkbar_kca3 = 0
	soma.gkbar_kca3 = 0
	prim.gkbar_kO   = 0.000166
	prim.eO	    = -90
	soma.gkbar_kO   = 0.0001992
	soma.eO	    = -90
	run()
}
proc fig6C(){
	reset()
	tstop = 4000
	grapher(3000,-80,1000,140)
	prim.gkbar_kA   = 0
	soma.gkbar_kA   = 0
	prim.gkbar_kca3 = 0
	soma.gkbar_kca3 = 0
	prim.gkbar_kO   = 0.00021
	prim.eO	    = -90
	soma.gkbar_kO   = 0.000261
	soma.eO	    = -90
	run()
}
proc fig6D(){
	reset()
	tstop = 4000
	grapher(3000,-80,1000,140)
	prim.gbar_NaP   = 0
	prim.gkbar_kO   = 0.000045
	prim.eO	    = 55
	run()
}
proc fig6E(){
	reset()
	tstop = 5000
	grapher(4000,-80,1000,140)
	prim.gkhbar_Ih  = 0
	prim.gkbar_kO   = 0.00013715
	prim.eO	    = 0
	run()
}
proc fig7A(){
	reset()
	tstop = 4000
	grapher(2000,-80,2000,140)
	stim1.amp = 1.58
	stim1.dur = tstop
	run()
}
proc fig7B(){
	reset()
	tstop = 4000
	grapher(2000,-80,2000,140)
	stim1.amp = 1.58
	stim1.dur = tstop
	soma.gkbar_kA = soma.gkbar_kA/2
	prim.gkbar_kA = prim.gkbar_kA/2
	run()
}
proc fig8A(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	syn.onset = 2100
	syn.e = 0
	syn.tau = 50
	syn.gmax = 0.01
	g.addexpr("Synaptic Current","1000*(syn.i/glom.v(0.5))", 2, 1, 0.8, 0.9, 2)
	run()
}
proc fig8B(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	syn.onset = 2100
	syn.e = 0
	syn.tau = 50
	syn.gmax = 0.01
	prim.gkbar_kO = 0.0001
	prim.eO = -90
	g.addexpr("Synaptic Current","1000*(syn.i/glom.v(0.5))", 2, 1, 0.8, 0.9, 2)
	run()
}
proc fig8C(){
	reset()
	tstop = 2200
	grapher(2000,-80,200,140)
	syn.onset = 2100
	syn.e = 0
	syn.tau = 50
	syn.gmax = 0.04
	prim.gkbar_kO = 0.0001
	prim.eO = -90
	g.addexpr("Synaptic Current","1000*(syn.i/glom.v(0.5))", 2, 1, 0.8, 0.9, 2)
	run()
}
proc reset(){
	stim1.amp = 0
	stim1.dur = 0
	stim1.del = 0
	stim2.amp = 0
	stim2.del = 0
	stim2.dur = 0
	syn.gmax = 0
	syn.e = 0
	syn.tau = 0
	syn.onset = 0
	soma.gkbar_kA   = 0.012
	soma.gkbar_kca3 = 0.12
	soma.gkbar_kO   = 0
	soma.eO	    = 0
	prim.gkbar_kA   = 0.00735
	prim.gkbar_kca3 = 0.1
   	prim.gbar_NaP   = 0.00042
   	prim.gkhbar_Ih  = 0.0024
	prim.gkbar_kO   = 0
	prim.eO	    = 0
	forall if (ismembrane("k_ion")) {
		ek  = -90	
		}
	g.erase_all
}