/////////////////////////////////
//
// Figure 5 - varying peak calcium concentration by varying beta parameter
// High calcium configuration
// Jaffe, Wang, and Brenner (2011)
//
// Original code based on Aradi and Holmes (1999)
//
/////////////////////////////////

load_file("nrngui.hoc")

v_init = -70
tstart = 0
tstop = 30
dt = 0.025
celsius = 35

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				// Ca Channesl
      gtcabar_Ca = 0
      gncabar_Ca = 0
      glcabar_Ca = 0.0018

   insert aabBK				// BK Channels
      gakbar_aabBK = 0			   
      gabkbar_aabBK = 0
      cascale_aabBK = .01		// Beta factorg
      tau_aabBK = 1   			
      base_aabBK = 4			

    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 = 5
stim.dur = 1
stim.amp = 2

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

init()

access soma

proc advance() {		

     	// output is in four columns: time, soma.v, calcium, I-BK

	fprint("%g %g %g %g\n",t,soma.v,cal_aabBK, (gak_aabBK+gabk_aabBK)*(soma.v+85))

	fadvance()
}

strdef tname



proc DualSimulate() {	// Perform two simulations, one with alphabeta4 and one with alpha

      // Calcium channels
     
      gtcabar_Ca = 0
      gncabar_Ca = 0
      glcabar_Ca = 0.0045
 
      // BK Channels - alphabeta4 on
 
      gakbar_aabBK = 0
      gabkbar_aabBK = .05

      sprint(tname,"trace_ab_%d.dat",$1)	// output file name
      wopen(tname)
      run()					

      // BK Channels - alpha on

      gakbar_aabBK = .05
      gabkbar_aabBK = 0

      sprint(tname,"trace_a_%d.dat",$1)		// output file name
      wopen(tname)
      run()
      wopen()
}

// Loop through parameter file systematically varying beta

ropen("varybeta.parm")		// beta parameter file

for i=0,7 {
    cascale_aabBK = fscan()
    DualSimulate(i)
}

cascale_aabBK = 0
DualSimulate(8)