/* --------------------------------------------------------------
   NEURON code to simulate BAC firing in compartmental models
   of Layer 5 Pyramidal cells
   for detailed description see Schaefer, Larkum, Sakmann, Roth
   "Coincidence detection in pyramidal neurons is tuned by their
   dendritic branching pattern" J. Neurophys. in press
   -------------------------------------------------------------- */

strdef modelName, loadProgram, cellName, outputFile, cellPath, loadProgram
objref trunc, secR, fi,vC, mbSec, random
objref sh, axonal, dendritic, dendritic_only, stimSec
objref st, st2, syn,trunc,mbSec, middleSec

create somaA,iseg,hill,myelin[2],node[2],dendA1_0
access somaA 
sprint(cellPath,"")

//load_proc("nrnmainmenu")

load_file("nrngui.hoc")

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

spA  = 2        // increase in membrane area due to spines

ra        = 80          
global_ra = ra
rm        = 30000           
c_m       = 0.6
cm_myelin = 0.04
g_pas_node = 0.02

v_init    = -73.65            
celsius   = 34          

Ek = -90
Ena = 60

gna_dend = 27 
gna_node = 30000   
gna_soma = 54
gna_myelin = 80  
gkv_axon = 3000 
gkv_soma = 600 
gkv_dend = 30 

gca_dend =1.5
gkm_dend =.1
gkca_dend=3.25

gca_soma = 3 
gkm_soma = 0.2 
gkca_soma = 6.5 
 
gka_soma  = 600        
gka_dend  =  300
gka_slope = 0   // no gradient

tauR   = 80
caiExp = 4
rA     = 0.05
rB     = 0.1

// --------------------------------------------------------------
// Low Threshold Ca Channel to reproduce frequency effect (Larkum, Kaiser, Sakmann, PNAS,1999)
// --------------------------------------------------------------
vh1_sit2=56
vh2_sit2=415
ahc_sit2=30               
v12m_sit2=45
v12h_sit2=65  
amc_sit2=3
vshift_sit2=10
vm1_sit2=50
vm2_sit2=125

it2_init=5   
gca_init=4.5

// --------------------------------------------------------------
// initiation zone in the dendrite 
// with slightly elevated Ca conductance densities
// --------------------------------------------------------------

proc InitZone() {

  mbSec.sec distance(0,1)
  forall for(x) if(distance(x) <100) {     
  eta_ssca = gca_init/gamma_ssca
    eta_sit2 = it2_init/gamma_sit2   
    }
    printf("InitZone for calcium initiated\n",spA)  
}
// --------------------------------------------------------------
// Axon geometry
//
// Similar to Mainen et al (Neuron, 1995)
// --------------------------------------------------------------

n_axon_seg = 5

proc create_axon() {local somaArea
  create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]

  somaA {
    somaArea=0
    for(x) somaArea+=area(x)
    equiv_diam = sqrt(somaArea/(4*PI))
  }

  iseg {                
     pt3dclear() pt3dadd(0,0,0,diam) pt3dadd(0,-1000,0,diam) 
     L = 15
     nseg = 5
     diam = equiv_diam/10    
  }

  hill {                
    pt3dclear() pt3dadd(0,0,0,diam) pt3dadd(0,-1000,0,diam) 
    L = 10
    nseg = 5
    diam(0:1) = 4*iseg.diam:iseg.diam
  }
  for i=0,n_axon_seg-1 {
    myelin[i] {         // myelin element
      pt3dclear() pt3dadd(0,0,0,diam) pt3dadd(0,-1000,0,diam)   
      nseg = 5
      L = 100
      diam = iseg.diam         
    }
    node[i] {           // nodes of Ranvier
      pt3dclear() pt3dadd(0,0,0,diam) pt3dadd(0,-1000,0,diam) 
      nseg = 1
      L = 1.0           
      diam = iseg.diam*.75       // nodes are thinner than axon
    }
  }
  somaA connect hill(0), 0.5
  hill connect iseg(0), 1
  iseg connect myelin[0](0), 1
  myelin[0] connect node[0](0), 1
  for i=0,n_axon_seg-2  { 
      node[i] connect myelin[i+1](0), 1 
      myelin[i+1] connect node[i+1](0), 1
  }
}

// --------------------------------------------------------------
// Spines
// --------------------------------------------------------------
proc add_spines() { 

  // increase all dendritic conductances by factor spA
  // increase dendritic cm and g_pas by same
  // to account for increase in membrane area without changing distances etc

  forsec dendritic_only {
    cm        *=spA 
    g_pas     *=spA 
    eta_sna   *=spA 
    eta_sk   *=spA 
    eta_skm   *=spA 
    eta_skca  *=spA 
    eta_ssca  *=spA 
    eta_ska  *=spA
    eta_sit2 *=spA   } 
    printf("spine factor of %d is now incorporated\n",spA)  
}

initzoneflag = 1

// --------------------------------------------------------------
// Initialization routines
// --------------------------------------------------------------
proc init_cell() {

  // passive
  forall {
    insert pas
    Ra = ra 
    cm = c_m 
    g_pas = 1/rm
    e_pas = -70 //v_init
  }

  // exceptions along the axon
  forsec "myelin" cm = cm_myelin
  forsec "node" g_pas = g_pas_node

  // active 
  // axon
  forall insert sna   forsec "myelin" eta_sna = gna_myelin/gamma_sna
  forsec "hill" eta_sna = gna_node/gamma_sna
  forsec "iseg" eta_sna = gna_node/gamma_sna   
  forsec "node" eta_sna = gna_node/gamma_sna
  forsec "iseg" { insert sk  eta_sk = gkv_axon/gamma_sk}
  forsec "hill" { insert sk  eta_sk = gkv_axon/gamma_sk}

  // dendrites
  forsec dendritic_only {
                eta_sna = gna_dend/20
    insert sk    eta_sk = gkv_dend/gamma_sk
    insert skm    eta_skm  = gkm_dend/gamma_skm 
    insert skca   eta_skca = gkca_dend/gamma_skca
    insert ska   eta_ska = gka_dend/gamma_ska 
    insert ssca   eta_ssca = gca_dend/gamma_ssca
    insert sit2   eta_sit2=0      insert cad2
  }

  // soma
  somaA {
                 eta_sna = gna_soma/gamma_sna
    insert sk    eta_sk = gkv_soma/gamma_sk 
    insert skm    eta_skm = gkm_soma/gamma_skm
    insert skca   eta_skca = gkca_soma/gamma_skca 
    insert ska   eta_ska = gka_soma/gamma_ska
    insert ssca   eta_ssca = gca_soma/gamma_ssca
    insert cad2
  }

  forall if(ismembrane("k_ion")) ek = Ek
  forall if(ismembrane("na_ion")) {
    ena = Ena
    vshift_sna = -5
  }
  forall if(ismembrane("ca_ion")) {
    eca = 140
    ion_style("ca_ion",0,1,0,0,0)
    //vshift_ca = 10
  }

// ca diffusion  and kca parameters
  taur_cad2 = tauR 
  caix_skca  = caiExp
  Ra_skca    = rA 
  Rb_skca    = rB 

if (initzoneflag) { InitZone() }
  add_spines()
}

// --------------------------------------------------------------
// loading cell
// --------------------------------------------------------------

proc load_3dcell() {
  // $s1 filename   aspiny = 0
  forall delete_section()
  xopen($s1)
  access somaA
  forsec "axon" delete_section()
  dendritic = new SectionList()   // make sure no compartments exceed 20 uM length
  forall {
//    if(nseg < L/20) { print secname(), " not accurate" nseg=L/20+1 }
    dendritic.append()
  }    
  dendritic_only = new SectionList()
  forsec dendritic dendritic_only.append()
  somaA  dendritic_only.remove()
}

// --------------------------------------------------------------
// Main Loading procedure
// --------------------------------------------------------------

proc LoadNInit() {
  sprint(loadProgram,"%s.nrn",$s1)
  load_3dcell(loadProgram) 
  create_axon()
  init_cell()
}

// --------------------------------------------------------------
// setting stimuli
// --------------------------------------------------------------
proc Inumber2() {
  st2=new IClamp(.5)
  st2.del =  21.1 // 1005.1
  st2.dur = 5 
  st2.amp = 2.2
}

proc IatSoma() {
  st=new IClamp(.5)
  st.del =  1.1 // 1005.1
  st.dur = 5 
  st.amp = 1.2
}

// EPSP : f(t) = (1-exp(-t/chi1)) * exp(-t/chi2) 
//       mit chi1 = 0.5 -2 ms und chi2 = 2-8 ms

proc EPSPAtDend() {
 mbSec.sec  {
   syn = new epsp(1)
   syn.tau0 = 0.8       
   syn.tau1 = 4         
   syn.onset = 3.1 //1007.1  
   syn.imax = 0.5
  }    
}

// ----------------------------------------------------------------
// Run Routine
// ----------------------------------------------------------------

proc DoIt() {local i, j
  init()
  for i=0,999 {
    for j=0,0.1/dt-1 fadvance()
  }
} 
forall delete_section()
cellName="d980329a-1"
printf("%s is the current cell\n",cellName)

LoadNInit(cellName)
IatSoma()
EPSPAtDend()

// --------------------------------------------------------------
// finally opening all graphs and panels
// --------------------------------------------------------------

xopen("Display2.ses")    

random = new RNG(0.5)
random.init_seed = 928102051 // seed for random number generator

init()
random.init_seed = -1   // this step is necessary to continue seed generation after
                        // the initial trial, insteading of reusing the same seed