load_file("nrngui.hoc")

v_init = -70
tstart = 0
tstop = 200
dt = 0.025

load_file("channels.ses")
objref p
p = new PWManager()
for i = 2,5 {p.hide(i)}

proc celldef() {
  topol()
  subsets()
  geom()
  biophys()
  geom_nseg()
}

create soma, axon[4], GCL[2], prox[2], middle[2], distal[2]

proc topol() { local i
  connect axon(0), soma(1)
  for i = 1, 3 connect axon[i](0), axon[i-1](1)
  for i = 0, 1 connect GCL[i](0), soma(0)
  for i = 0, 1 connect prox[i](0), GCL[i](1)
  for i = 0, 1 connect middle[i](0), prox[i](1)
  for i = 0, 1 connect distal[i](0), middle[i](1)
  basic_shape()
}
proc basic_shape() {
  soma {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(15, 0, 0, 1)}
  axon {pt3dclear() pt3dadd(15, 0, 0, 1) pt3dadd(30, 0, 0, 1)}
  axon[1] {pt3dclear() pt3dadd(30, 0, 0, 1) pt3dadd(45, 0, 0, 1)}
  axon[2] {pt3dclear() pt3dadd(45, 0, 0, 1) pt3dadd(60, 0, 0, 1)}
  axon[3] {pt3dclear() pt3dadd(60, 0, 0, 1) pt3dadd(120, 0, 0, 1)}
  GCL {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(-29, 30, 0, 1)}
  GCL[1] {pt3dclear() pt3dadd(0, 0, 0, 1) pt3dadd(-29, -29, 0, 1)}
  prox {pt3dclear() pt3dadd(-29, 30, 0, 1) pt3dadd(-44, 30, 0, 1)}
  prox[1] {pt3dclear() pt3dadd(-29, -29, 0, 1) pt3dadd(-44, -29, 0, 1)}
  middle {pt3dclear() pt3dadd(-44, 30, 0, 1) pt3dadd(-59, 30, 0, 1)}
  middle[1] {pt3dclear() pt3dadd(-44, -29, 0, 1) pt3dadd(-59, -29, 0, 1)}
  distal {pt3dclear() pt3dadd(-59, 30, 0, 1) pt3dadd(-119, 30, 0, 1)}
  distal[1] {pt3dclear() pt3dadd(-59, -29, 0, 1) pt3dadd(-119, -29, 0, 1)}
}

objref all, dendrites, GCLs, proxs, middles, distals, allaxon
proc subsets() { local i
  objref all, dendrites, GCLs, proxs, middles, distals, allaxon
  all = new SectionList()
    soma all.append()
    for i=0, 3 axon[i] all.append()
    for i=0, 1 GCL[i] all.append()
    for i=0, 1 prox[i] all.append()
    for i=0, 1 middle[i] all.append()
    for i=0, 1 distal[i] all.append()

  dendrites = new SectionList()
    for i=0, 1 GCL[i] dendrites.append()
    for i=0, 1 prox[i] dendrites.append()
    for i=0, 1 middle[i] dendrites.append()
    for i=0, 1 distal[i] dendrites.append()

  GCLs = new SectionList()
    for i=0, 1 GCL[i] GCLs.append()

  proxs = new SectionList()
    for i=0, 1 prox[i] proxs.append()

  middles = new SectionList()
    for i=0, 1 middle[i] middles.append()

  distals = new SectionList()
    for i=0, 1 distal[i] distals.append()

  allaxon = new SectionList()
    for i=0, 3 axon[i] allaxon.append()

}
proc geom() {
  forsec dendrites {  L = 150  diam = 3  }
  forsec GCLs {  L = 50  }
  forsec allaxon {  L = 50  }
   axon.diam = 0.9
   axon[1].diam = 0.7
   axon[2].diam = 0.5
   axon[3].diam = 0.4
  soma {  L = 16.8  diam = 16.8  }
  axon[3] {  L = 1400  }
}
proc geom_nseg() {
   forsec dendrites { nseg = 4  }
   forsec GCLs { nseg = 2  }
   forsec allaxon { nseg = 1  }
   soma { nseg = 2  }
   axon[3] { nseg = 28  }
}
proc biophys() {
  forsec all {
    Ra = 210
    cm = 1
    insert pas
      g_pas = 2.5e-05
      e_pas = -70
  }
  forsec dendrites {
    cm = 1.6
    insert pas
      g_pas = 4e-05
      e_pas = -70
  }
  forsec GCLs {
    cm = 1
    insert pas
      g_pas = 2.5e-05
      e_pas = -70
    insert Na
      gmax_Na = 0.018
    insert fKDR
      gmax_fKDR = 0.004
    insert sKDR
      gmax_sKDR = 0.003
  }
  forsec proxs {
    insert Na
      gmax_Na = 0.013
    insert fKDR
      gmax_fKDR = 0.004
    insert sKDR
      gmax_sKDR = 0.003
  }
  forsec middles {
    insert Na
      gmax_Na = 0.008
    insert fKDR
      gmax_fKDR = 0.001
    insert sKDR
      gmax_sKDR = 0.003
  }
  forsec distals {
    insert fKDR
      gmax_fKDR = 0.001
    insert sKDR
      gmax_sKDR = 0.004
  }
  forsec allaxon {
    insert Na
      gmax_Na = 0.21
    insert fKDR
      gmax_fKDR = 0.028
  }
  soma {
    insert Ca
      gtcabar_Ca = 0.0018
      gncabar_Ca = 0
      glcabar_Ca = 0
    insert CadepK
      gbkbar_CadepK = 0
      gskbar_CadepK = 0
    insert Na
      gmax_Na = 0.12
    insert fKDR
      gmax_fKDR = 0.016
    insert sKDR
      gmax_sKDR = 0.003
  }
}
access soma

celldef()

forsec all {
	if (ismembrane("Na")) {ena = 45}
	ek = -85
}

// insert current clamp at the soma
objref stim
stim = new IClamp(0.5)
stim.del = 20
stim.dur = 0.5
stim.amp = 2

proc init() {local i
	finitialize(v_init)
	fcurrent()
	t = tstart
}

init()

load_file("graph.ses")

proc advance() {
	fadvance()
}

// declaration of simulation procedures
proc control() {
	soma.gbkbar_CadepK = 0
	soma.gskbar_CadepK = 0
	run()
}

proc bk() {
	soma.gbkbar_CadepK = 0.00226
	soma.gskbar_CadepK = 0
	run()
}

proc sk() {
	soma.gbkbar_CadepK = 0
	soma.gskbar_CadepK = 0.00135
	run()
}