// --------------------------------------------------------------

// Investigating  ANTIDROMIC ACTION POTENTIALS

// Using Diana Smetters p19 rat L5 pyramid

// Huguenard/McCormick HH kinetics 
//  na3h5.mod, kd3h5.mod

// Set up to reproduce Stuart and Sakmann (1994), Fig. 2a

// --------------------------------------------------------------


load_proc("nrnmainmenu")


// --------------------------------------------------------------
// redefine some things in stdrun.hoc
// --------------------------------------------------------------

tstop = 17.7
steps_per_ms = 40
dt = 0.025
secondorder = 1

// --------------------------------------------------------------
// cell geometry
// --------------------------------------------------------------


xopen("dks577a.hoc")
xopen("dks577.append.hoc")


// --------------------------------------------------------------
// passive & active membrane 
// --------------------------------------------------------------

ra        = 200
global_ra = ra
rm        = 40000
c_m       = 0.75
v_init    = -70
celsius = 23

Ek = -90
Ena = 60
            

gna_dend = 30
gna_node = 30000
gk = 100
frac_gk_apical = 0
gk_node = 0


proc install_channels() {
  
  // passive
  forall {
    insert pas
    Ra = ra
    cm = c_m 
    g_pas = 1/rm
    e_pas = v_init
  }
  // exceptions along the axon
  forsec "myelin" cm = 0.04
  forsec "node" g_pas = 0.02

//  add_spines(apical,sd_apical)
//  add_spines(oblique,sd_oblique)
//  add_spines(distal,sd_distal)
//  add_spines(basal,sd_basal)


  // Na+ channels - Stochastic
  forall insert sna
  forsec dend gbar_sna = gna_dend
  forsec "myelin" gbar_sna = gna_dend
  soma.gbar_sna = gna_dend
  hill.gbar_sna = gna_node
  iseg.gbar_sna = gna_node
  forsec "node" gbar_sna = gna_node

  forall if (ismembrane("sna")) { 
          ena = Ena     // per S&S ena = 90
              for (x) { if (x!=0 &&  x!=1) {
                      D_sna(x) = int(gbar_sna(x)/gama_sna(x))
              }}
  }

  // K+ channels - Stochastic
  forsec basal {
    insert sk
    gbar_sk = gk
    for (x) { if (x!=0 &&  x!=1) {
      D_sk(x) = int(gbar_sk(x)/gama_sk(x))
    }
   }
  }

  soma insert sk
  soma.gbar_sk = gk
  soma { D_sk = int(gbar_sk/gama_sk) } 
    // no K+ in apical dend or axon
  forall if (ismembrane("sk")) ek = Ek         // per S&S ek = -100
}

install_channels()

//RNG
objref rnd
soma rnd = new RNG(0.5)
rnd.seed = 23658859 

// --------------------------------------------------------------
// stimulus
// --------------------------------------------------------------

objectvar st
dsite.sec st = new IClamp(site_loc)
st.amp = 0.25
st.del = 3
st.dur = 40


// --------------------------------------------------------------
// create useful graphs & panels
// --------------------------------------------------------------
load_file("nrngui.hoc")
nrncontrolmenu()

newPlotV()
graphItem.addvar("dsite.sec.v(site_loc)",4,1)
graphItem.addvar("iseg.v(.5)",7,1)

proc set_stim() {
  st.loc($1)  st.amp = $2   st.del = $3   st.dur = $4
}

proc soma_inj() {
  install_channels()
  soma set_stim(.5,.16,3,40)
  run()
}

proc dend_inj() {
  install_channels()
  dsite.sec set_stim(site_loc,.21,3,40)
  run()
}

xpanel("Stuart & Sakmann")
xvalue("gna dend (pS/um2)","gna_dend")
xvalue("gna axon (pS/um2)","gna_node")
xvalue("stim amp","st.amp")
xbutton("inject soma","soma_inj()")
xbutton("inject dend","dend_inj()")
xpanel()
//forall nseg = 1