/* --------------------------------------------------------------
   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  
   
   modified by K. Diba for stochastic simulations and 
   enhanced BPAP for trailing spikes in spike trains
 
   -------------------------------------------------------------- */

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 = 14 
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_snap   *=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 = 0

// --------------------------------------------------------------
// 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 deterministic_sna = 0
  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 {
    uninsert sna

    insert snap   eta_snap = gna_dend/gamma_snap
    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("sna")) {
    ena = Ena
    vshift_sna = -5
  }
  forall if(ismembrane("snap")) {
    ena = Ena
    vshift_snap = -5
  }  
  forall if(ismembrane("ska")) {
    ek = Ek
    lmin_ska = 16.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 =  start_at + 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 = 2.0
}

// 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 = 235809314 // seed for Random number generator

init()
random.init_seed = -1       //required to continue using same seed