// The pyramidal cell with full dendritic tree and axon arbor for similating action potential generation process 
//of paramidal cell in layer V of prefrontal cortex. 

//This code is used to calculate the decay spacial constant of main axon. 

//This neuronal model was used for the following papers:
/*
Shu Y, Duque A, Yu YG, McCormick DA.  Properties of action potential
initiation in neocortical pyramidal cells: evidence from whole cell
axon recordings.  Journal of Neurophysiology. 97: 746-760, (2007).

Shu Y, Hasenstaub A, Duque A, Yu Y, McCormick DA (2006) Modulation of
intracortical synaptic potentials by presynaptic somatic membrane
potential.  Nature 441: 761-765.
*/

//----Yuguo Yu, 2007


// ew = 0.01  this determines the coupling strengh
//  soma set_stim(.5,0.25,0,10000) this determines the DC input to the soma

objref sh, st, st1, axonal, dendritic, dendritic_only

// load_proc("nrnmainmenu")
load_file("nrngui.hoc")

create soma
access soma

tstop = 300
//dt =  0.005
dt =  0.03
steps_per_ms = 1/dt


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

ra        = 100  //before is 150
global_ra = ra
rm        = 40000
c_m       = 0.7
cm_myelin = 0.04
//g_pas_node = 0.002
g_pas_node = 0.02

vsoma_init    = -65
v_init    = -70
celsius   = 37

Ek = -90
Ena = 60



gna_dend = 20
gna_node = 7500
gna_soma = gna_dend*37.5

gkv_axon = 1000
gkv_soma = 400

gca = .3
gkm = .1
gkca = 3

gca_soma = gca
gkm_soma = gkm
gkca_soma = gkca
 
//  forsec "node" gbar_kv = gkv_axon  

// --------------------------------------------------------------
// Axon geometry
// --------------------------------------------------------------
n_axon_seg = 15

create soma,iseg,hill,myelin[n_axon_seg],node[n_axon_seg],nod, cnod[8]
create ccnod[23]
create c11[11],c12[4], c13[30], c14[41], c15[9], c110[2], c111[2]
create c21[5], c22[19], c23[2]
create c31[5], c32

proc create_axon() {

 create iseg,hill,myelin[n_axon_seg],node[n_axon_seg],nod, cnod[8]
create ccnod[23]
create c11[11],c12[4], c13[30], c14[41], c15[9], c110[2], c111[2]
create c21[5], c22[19], c23[2]
create c31[5], c32

  soma {
    equiv_diam = sqrt(area(.5)/(4*PI))

    // area = equiv_diam^2*4*PI
  }
  if (numarg()) equiv_diam = $1

  iseg {                // initial segment 
     L = 40  
     nseg = 10
     diam = equiv_diam/15        // see Sloper and Powell 1982, Fig.71
  }

  hill {                
    L = 10
    nseg = 3
    diam(0:1) = 3*iseg.diam:iseg.diam
  }

  // construct myelinated axon with nodes of ranvier

    nod {           // unmyelinated main axon of PFC is long about 300 micronmeters
      nseg = 30
       //L = 1000         
       L = 300         
     diam = iseg.diam*0.9    
    }

    for i=0,7 {
    cnod[i] {           // 1st order collaterals
      nseg = 20
     diam = iseg.diam*0.5       // .5?  nodes are thinner than axon
            }
              }

    for i=0,22 {
    ccnod[i] {            //second order collaterals
      nseg = 20
     diam = iseg.diam*0.4      
            }
              }
    
     cnod[0] {L = 1780}
     cnod[1] {L = 3040}
     cnod[2] {L = 1600}
     cnod[3] {L = 1500}
     cnod[4] {L = 300}
     cnod[5] {L = 60}
     cnod[6] {L = 830}
     cnod[7] {L = 880}
     
    for i=0,10 {            //3rd order collaterals
      c11[i] {
      nseg = 10
     diam = iseg.diam*0.35       
               }
              }
    for i=0,3 {       
      c12[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }
    for i=0,29 {           
      c13[i] {
      nseg = 10
     diam = iseg.diam*0.35  
               }
              }
    for i=0,40 {           
      c14[i] {
      nseg = 10
     diam = iseg.diam*0.35 
               }
              }
    for i=0,8 {          
      c15[i] {
      nseg = 10
     diam = iseg.diam*0.35    
               }
              }
    for i=0,1 {           // 
      c110[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }
    for i=0,1 {           // 
      c111[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }
    
    for i=0,4 {           // 
      c21[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }
    for i=0,18 {           //
      c22[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }
    for i=0,1 {           // 
      c23[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }

    for i=0,4 {           // 
      c31[i] {
      nseg = 10
     diam = iseg.diam*0.35   
               }
              }
              
      c32 {
      nseg = 5
     diam = iseg.diam*0.35     
      L = 130
          }
          
     //12 subcol for 1st col
     ccnod[0] {L = 1000} ccnod[1] {L = 625} ccnod[2] {L = 1310} ccnod[3] {L = 1160}
     ccnod[4] {L = 850} ccnod[5] {L = 460} ccnod[6] {L = 160} ccnod[7] {L = 480}
     ccnod[8] {L = 110} ccnod[9] {L = 480} ccnod[10] {L = 420} ccnod[11] {L = 126}
        //11 subsubcol for ccnod[0]
         c11[0] {L = 760} c11[1] {L = 360} c11[2] {L = 260} c11[3] {L = 105}
         c11[4] {L = 360} c11[5] {L = 280} c11[6] {L = 130} c11[7] {L = 180}
         c11[8] {L = 710} c11[9] {L = 400} c11[10] {L = 435}
        // 4 subsubcol for ccnod[1]
         c12[0] {L = 400} c12[1] {L = 510} c12[2] {L = 90} c11[3] {L = 150}
     
       //30 subsubcol for ccnod[2]
         c13[0] {L = 1060} c13[1] {L = 438} c13[2] {L = 272} c13[3] {L = 264}
         c13[4] {L = 510} c13[5] {L = 85} c13[6] {L = 220} c13[7] {L = 610}
         
         c13[8] {L = 700} c13[9] {L = 40} c13[10] {L = 178}
         c13[11] {L = 890} c13[12] {L = 600} c13[13] {L = 440} c13[14] {L = 615}
         c13[15] {L = 444} c13[16] {L = 330} c13[17] {L = 325} c13[18] {L = 80}
         c13[19] {L = 150} c13[20] {L = 130} c13[21] {L = 560} c13[22] {L = 480}
         
         c13[23] {L = 606} c13[24] {L = 310} c13[25] {L = 290} c13[26] {L = 108} c13[27] {L = 200}
         c13[28] {L = 72} c13[29] {L = 60} 

       //41 subsubcol for ccnod[3]
         c14[0] {L = 1075} c14[1] {L = 50} c14[2] {L = 818} c14[3] {L = 290}  c14[4] {L = 103} 
         c14[5] {L = 657} c14[6] {L = 513} c14[7] {L = 226}
         c14[8] {L = 500} c14[9] {L = 254} 
         c14[10] {L = 700} c14[11] {L = 491} c14[12] {L = 584} c14[13] {L = 205} 
         c14[14] {L = 278} c14[15] {L = 411} c14[16] {L = 286} 
         c14[17] {L = 681} c14[18] {L = 57}  c14[19] {L = 480} 
         
         c14[20] {L = 504} c14[21] {L = 277}
         c14[22] {L = 487} c14[23] {L = 422} c14[24] {L = 326} c14[25] {L = 60} c14[26] {L = 109} c14[27] {L = 366} 
         c14[28] {L = 391} c14[29] {L = 69} 
         c14[30] {L = 544} c14[31] {L = 366} c14[32] {L = 32} c14[33] {L = 120} c14[34] {L = 504} c14[35] {L = 68} 
         c14[36] {L = 512} c14[37] {L = 382} c14[38] {L = 370} c14[39] {L = 50} c14[40] {L = 60} 
                  
                  
       //9 subsubcol for ccnod[4]
         c15[0] {L = 430} c15[1] {L = 658} c15[2] {L = 302} c15[3] {L = 160}  c15[4] {L = 560} 
         c15[5] {L = 160} c15[6] {L = 330} c15[7] {L = 170}
         c15[8] {L = 330} 
         
         //2 subsubcol for ccnod[9]
         c110[0] {L = 300} c110[1] {L = 200}
         
         //2 subsubcol for ccnod[10]
         c111[0] {L = 90} c111[1] {L = 260}
     
     //4 subcol for 2nd col
     ccnod[12] {L = 650} ccnod[13] {L = 1425} ccnod[14] {L = 320} ccnod[15] {L = 930}
        //5 subsubcol for ccnod[12]
        c21[0] {L = 450} c21[1] {L = 225} c21[2] {L = 130}  c21[3] {L = 110}  c21[4] {L = 130} 
        //19 subsubcol for ccnod[13]
        c22[0] {L = 690} c22[1] {L = 615} c22[2] {L = 290}  c22[3] {L = 50}  c22[4] {L = 155} 
        c22[5] {L = 450} c22[6] {L = 190} c22[7] {L = 80}  c22[8] {L = 115}  c22[9] {L = 60} 
        c22[10] {L = 770} c22[11] {L = 260} c22[12] {L = 45}  c22[13] {L = 90}  c22[14] {L = 510} 
        c22[15] {L = 220} c22[16] {L = 250} c22[17] {L = 220}  c22[18] {L = 320} 
       //2 subsubcol for ccnod[14]
        c23[0] {L = 140} c23[1] {L = 175}
        
     //4 subcol for 3rd col
     ccnod[16] {L = 470} ccnod[17] {L = 125} ccnod[18] {L = 335} ccnod[19] {L = 230}
         c31[0] {L = 445} c31[1] {L = 305} c31[2] {L = 155}  c31[3] {L = 360} c31[4] {L = 140} 
     
     //2 subcol for 4th col
      ccnod[20] {L = 450} ccnod[21] {L = 50}
        
     //1 subcol for 7th col
      ccnod[22] {L = 60}
      
  for i=0,n_axon_seg-1 {
    myelin[i] {         // myelin element
      nseg = 20
      L = 200
      diam = iseg.diam         
    }
    
    node[i] {           // nodes of Ranvier
      nseg = 1
      L = 1.0           
      diam = iseg.diam*.5      
    }
  }
}
   
proc create_axon_connection() {

    soma connect hill(0), 0.5
  hill connect iseg(0), 1
  iseg connect nod(0), 1
  nod  connect node[0](0), 1
   for i=0,n_axon_seg-2  { 
      node[i] connect myelin[i](0), 1 
      myelin[i] connect node[i+1](0), 1
  }


  // nod=500, iseg=50, hill=10
  nod connect cnod[0](0), 0.06
  nod connect cnod[1](0), 0.16
  nod connect cnod[2](0), 0.42
  nod connect cnod[3](0), 0.7
  nod connect cnod[4](0), 0.91
  nod connect cnod[5](0), 0.96
  node[6] connect cnod[6](0), 0.5
  node[9] connect cnod[7](0), 0.5
    
  cnod[0] connect ccnod[0](0), 0.018
    ccnod[0] connect c11[0](0), 0.022
      c11[0] connect c11[1](0), 0.08
      c11[0] connect c11[2](0), 0.16
      c11[2] connect c11[3](0), 0.1

    ccnod[0] connect c11[4](0), 0.035
      c11[4] connect c11[5](0), 0.018
      c11[4] connect c11[6](0), 0.036
      c11[5] connect c11[7](0), 0.4
    
    ccnod[0] connect c11[8](0), 0.056
      c11[8] connect c11[9](0), 0.6
   
     ccnod[0] connect c11[10](0), 0.2

     cnod[0] connect ccnod[1](0), 0.066
     ccnod[1] connect c12[0](0), 0.001
       c12[0] connect c12[1](0), 0.12
       c12[0] connect c12[2](0), 0.2
       c12[1] connect c12[3](0), 0.08
    cnod[0] connect ccnod[2](0), 0.13
     ccnod[2] connect c13[0](0), 0.14
       c13[0] connect c13[1](0), 0.4
       c13[0] connect c13[2](0), 0.45
       c13[2] connect c13[3](0), 0.2
       c13[0] connect c13[4](0), 0.5
  	c13[4] connect c13[5](0), 0.16
	c13[4] connect c13[6](0), 0.25
       
     ccnod[2] connect c13[7](0), 0.25
     ccnod[2] connect c13[8](0), 0.28
        c13[8] connect c13[9](0), 0.15
	c13[8] connect c13[10](0), 0.35
     
     ccnod[2] connect c13[11](0), 0.33
        c13[11] connect c13[12](0), 0.6
	c13[12] connect c13[13](0), 0.05
	
    ccnod[2] connect c13[14](0), 0.35
        c13[14] connect c13[15](0), 0.15
	c13[15] connect c13[16](0), 0.2
	c13[15] connect c13[17](0), 0.25
	c13[17] connect c13[18](0), 0.1
	c13[17] connect c13[19](0), 0.3
	c13[15] connect c13[20](0), 0.3
	
    ccnod[2] connect c13[21](0), 0.38
        c13[21] connect c13[22](0), 0.3
	
    ccnod[2] connect c13[23](0), 0.42
        c13[23] connect c13[24](0), 0.08
        c13[24] connect c13[25](0), 0.08
        c13[25] connect c13[26](0), 0.2
        c13[24] connect c13[27](0), 0.25
        
        c13[23] connect c13[28](0), 0.1
        c13[23] connect c13[29](0), 0.15

    cnod[0] connect ccnod[3](0), 0.144
    ccnod[3] connect c14[0](0), 0.035
        c14[0] connect c14[1](0), 0.45
        
    ccnod[3] connect c14[2](0), 0.155
        c14[2] connect c14[3](0), 0.36
        c14[2] connect c14[4](0), 0.71
        
    ccnod[3] connect c14[5](0), 0.264
        c14[5] connect c14[6](0), 0.21
        c14[5] connect c14[7](0), 0.29
        
    ccnod[3] connect c14[8](0), 0.32
        c14[8] connect c14[9](0), 0.354
        
    ccnod[3] connect c14[10](0), 0.35
        c14[10] connect c14[11](0), 0.1
        c14[10] connect c14[12](0), 0.143
        c14[10] connect c14[13](0), 0.236
        c14[10] connect c14[14](0), 0.243
        c14[10] connect c14[15](0), 0.4
        c14[10] connect c14[16](0), 0.57
        
    ccnod[3] connect c14[17](0), 0.4
        c14[17] connect c14[18](0), 0.39
    ccnod[3] connect c14[19](0), 0.45
    ccnod[3] connect c14[20](0), 0.47
        c14[20] connect c14[21](0), 0.08
    
    ccnod[3] connect c14[22](0), 0.49
        c14[22] connect c14[23](0), 0.1
        c14[23] connect c14[24](0), 0.1
        c14[23] connect c14[25](0), 0.2
        c14[22] connect c14[26](0), 0.2
        c14[22] connect c14[27](0), 0.3
        
    ccnod[3] connect c14[28](0), 0.51
        c14[28] connect c14[29](0), 0.45
        
    ccnod[3] connect c14[30](0), 0.55
        c14[30] connect c14[31](0), 0.01
        c14[31] connect c14[32](0), 0.1
        c14[31] connect c14[33](0), 0.2
        c14[30] connect c14[34](0), 0.03
        c14[34] connect c14[35](0), 0.01
        
    ccnod[3] connect c14[36](0), 0.6
        c14[36] connect c14[37](0), 0.05
        c14[36] connect c14[38](0), 0.27
        c14[36] connect c14[39](0), 0.37
        
    ccnod[3] connect c14[40](0), 0.76

  cnod[0] connect ccnod[4](0), 0.486
    ccnod[4] connect c15[0](0), 0.05
    ccnod[4] connect c15[1](0), 0.2
    c15[1] connect c15[2](0), 0.5
    
    ccnod[4] connect c15[3](0), 0.25
    ccnod[4] connect c15[4](0), 0.3
    c15[4] connect c15[5](0), 0.2
    c15[4] connect c15[6](0), 0.4
    c15[4] connect c15[7](0), 0.6
    ccnod[4] connect c15[8](0), 0.45
    
      
  cnod[0] connect ccnod[5](0), 0.494
  cnod[0] connect ccnod[6](0), 0.568
  cnod[0] connect ccnod[7](0), 0.576
  cnod[0] connect ccnod[8](0), 0.606
  cnod[0] connect ccnod[9](0), 0.628
        ccnod[9] connect c110[0](0), 0.1
        ccnod[9] connect c110[1](0), 0.4
  
  cnod[0] connect ccnod[10](0), 0.652
        ccnod[10] connect c111[0](0), 0.1
        ccnod[10] connect c111[1](0), 0.25
  
  cnod[0] connect ccnod[11](0), 0.679
  
  cnod[1] connect ccnod[12](0), 0.0157
    ccnod[12] connect c21[0](0), 0.031
    c21[0] connect c21[1](0), 0.156
    c21[0] connect c21[2](0), 0.333
    c21[0] connect c21[3](0), 0.49
    ccnod[12] connect c21[4](0), 0.29
  
  cnod[1] connect ccnod[13](0), 0.0173
    ccnod[13] connect c22[0](0), 0.042
    c22[0] connect c22[1](0), 0.06
    c22[1] connect c22[2](0), 0.4
    c22[1] connect c22[3](0), 0.7
    c22[0] connect c22[4](0), 0.12
    
    ccnod[13] connect c22[5](0), 0.07
      c22[5] connect c22[6](0), 0.4
      c22[5] connect c22[7](0), 0.5
    
    ccnod[13] connect c22[8](0), 0.2
      c22[8] connect c22[9](0), 0.2
    
    ccnod[13] connect c22[10](0), 0.4
      c22[10] connect c22[11](0), 0.4
      c22[10] connect c22[12](0), 0.5
      c22[10] connect c22[13](0), 0.6
    
    ccnod[13] connect c22[14](0), 0.55
      c22[14] connect c22[15](0), 0.5
    
    ccnod[13] connect c22[16](0), 0.7
      c22[16] connect c22[17](0), 0.4
    ccnod[13] connect c22[18](0), 0.8
    

  cnod[1] connect ccnod[14](0), 0.0535
    ccnod[14] connect c23[0](0), 0.2
    ccnod[14] connect c23[1](0), 0.4
        
  
  cnod[1] connect ccnod[15](0), 0.179
    
  cnod[2] connect ccnod[16](0), 0.025
    ccnod[16] connect c31[0](0), 0.02
    ccnod[16] connect c31[1](0), 0.03
    c31[1] connect c31[2](0), 0.4
    ccnod[16] connect c31[3](0), 0.1
    c31[3] connect c31[4](0), 0.1
       
  cnod[2] connect ccnod[17](0), 0.691
    ccnod[17] connect c32(0), 0.7
  cnod[2] connect ccnod[18](0), 0.697
  cnod[2] connect ccnod[19](0), 0.809

  cnod[3] connect ccnod[20](0), 0.267
  cnod[3] connect ccnod[21](0), 0.867
  cnod[6] connect ccnod[22](0), 0.78
       }

// --------------------------------------------------------------
// Spines
// --------------------------------------------------------------

      // Based on the "Folding factor" described in
      // Jack et al (1989), Major et al (1994)
      // note, this assumes active channels are present in spines 
      // at same density as dendrites

spine_dens = 1
      // just using a simple spine density model due to lack of data on some 
      // neuron types.

spine_area = 0.83 // um^2  -- K Harris

proc add_spines() { local a
  forsec $o1 {
    a =0
    for(x) a=a+area(x)

    F = (L*spine_area*spine_dens + a)/a

    L = L * F^(2/3)
    for(x) diam(x) = diam(x) * F^(1/3)
  }
}



proc init_cell() {

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


  // exceptions along the axon
    forsec "myelin" cm = cm_myelin*1
   forsec "myelin" g_pas = 0.00002
    forsec "node" g_pas = g_pas_node
    
  // na+ channels
  forall insert na

  forsec dendritic gbar_na = gna_dend*1
  forsec "myelin" gbar_na = gna_dend*1
  hill.gbar_na = gna_node
  iseg.gbar_na = gna_node
  //iseg.g_pas = 1/rm
  //nod.g_pas = 1/rm
  forsec "node" gbar_na = gna_node/5
  
  nod.gbar_na = gna_node/3
  forsec "cnod" gbar_na = gna_node/8.333
  forsec "ccnod" gbar_na = gna_node/10
    
  forsec "c11" gbar_na = gna_node/20
  forsec "c12" gbar_na = gna_node/20
  forsec "c13" gbar_na = gna_node/20
  forsec "c14" gbar_na = gna_node/20
  forsec "c15" gbar_na = gna_node/20
  forsec "c110" gbar_na = gna_node/20
  forsec "c111" gbar_na = gna_node/20
  forsec "c21" gbar_na = gna_node/20
  forsec "c22" gbar_na = gna_node/20
  forsec "c23" gbar_na = gna_node/20
  forsec "c31" gbar_na = gna_node/20
  forsec "c32" gbar_na = gna_node/20
  
 
//  forsec "node" gbar_kv = gkv_axon  
  // this is added by Yu at Oct.13

  // kv delayed rectifier channels
  iseg { insert kv  gbar_kv = gkv_axon }
  hill { insert kv  gbar_kv = gkv_axon }
  soma { insert kv  gbar_kv = gkv_soma }
  nod { insert kv  gbar_kv = gkv_axon/2 }
  forsec "cnod" { insert kv  gbar_kv = gkv_axon/5 }
  forsec "ccnod" { insert kv  gbar_kv = gkv_axon/10 }
  
  forsec "c11" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c12" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c13" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c14" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c15" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c110" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c111" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c21" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c22" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c23" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c31" { insert kv  gbar_kv = gkv_axon/20 }
  forsec "c32" { insert kv  gbar_kv = gkv_axon/20 }
      
  // dendritic channels
  forsec dendritic {
    insert km    gbar_km  = gkm
    insert kca   gbar_kca = gkca
    insert ca    gbar_ca = gca
    insert cad
  }

  soma {
    gbar_na = gna_soma
    gbar_km = gkm_soma
    gbar_kca = gkca_soma
    gbar_ca = gca_soma
  }

 
  forall if(ismembrane("k_ion")) ek = Ek
  forall if(ismembrane("na_ion")) {
    ena = Ena
    // seems to be necessary for 3d cells to shift Na kinetics -5 mV
    vshift_na = -5
  }
  forall if(ismembrane("ca_ion")) {
    eca = 140
    ion_style("ca_ion",0,1,0,0,0)
    vshift_ca = 0
  }
}

proc load_3dcell() {

// $s1 filename

  aspiny = 0
  forall delete_section()    //forall delete_section will remove all sections. 
  xopen($s1)
  access soma

  dendritic = new SectionList()

  // make sure no compartments exceed 50 uM length
  forall {
    diam_save = diam
    n = L/50
    nseg = n + 1
    if (n3d() == 0) diam = diam_save
    dendritic.append()
  }    


  // show cell
  sh = new PlotShape()
  sh.size(-300,300,-300,300)

  dendritic_only = new SectionList()
  forsec dendritic dendritic_only.append()
  soma  dendritic_only.remove()

  create_axon()
  create_axon_connection() 

  init_cell()

  if (!aspiny) add_spines(dendritic_only,spine_dens)

}


load_3dcell("j4a.hoc") 
// load_3dcell("cells/j4a.hoc") 
//load_3dcell("cells/j7.hoc") 
//load_3dcell("cells/j8.hoc") 
//load_3dcell("cells/lcAS3.hoc") 

nrnmainmenu()
nrncontrolmenu()
//newPlotV()

// --------------------------------------------------------------
// synapses ----June 5, 2005
// --------------------------------------------------------------
ne=3    // number of excitatory synapses
ni=3    // number of inhibitory synapses
create dummy
dummy {L = 20 diam = 20}
objref espike[ne], ispike[ni], esyn[ne], isyn[ni], esconn[ne], isconn[ni]

// Excitatory synaptic connections
// set synaptic position and properties
ew = 0.001  //for -80mv for 1/rm
  //ew = 0.0001 //for *100
soma esyn[0] = new Exp2Syn(0.5)
soma esyn[1] = new Exp2Syn(0.5)
soma esyn[2] = new Exp2Syn(0.5)

//dend11[1] esyn[0] = new Exp2Syn(0.5)
//dend11[1] esyn[1] = new Exp2Syn(0.5)
//dend11[1] esyn[2] = new Exp2Syn(0.5)
for i=0, ne-1 {
  esyn[i].tau1 = 1
  esyn[i].tau2 = 10
  esyn[i].e = 0
  // spike generator to excitatory connections
  dummy espike[i] = new NetStim(0.5)
//  espike[i].interval = 2
  espike[i].start = 50
//  espike[i].number = 0
//  espike[i].noise = 0
  // connect spike generator to synapse
  esconn[i] = new NetCon(espike[i], esyn[i], -20, 1, ew)
  esconn[0] = new NetCon(espike[0], esyn[0], -20, 1, ew*0.3)
  }
  espike[0].interval = 8
  espike[0].number = 10000
  espike[0].noise = 1
  
  espike[1].interval = 15
//espike[1].interval = 8
// espike[1].interval = 10
  espike[1].number = 10000
  espike[1].noise = 1

  espike[1].interval = 25
//  espike[1].interval = 20
//  espike[2].interval = 20
  espike[2].number = 10000
  espike[2].noise = 1
  

// Inhibitory synaptic connections
// set synaptic position and properties
//iw = 0.0004
iw = 0.00
soma isyn[0] = new Exp2Syn(0.5)
soma isyn[1] = new Exp2Syn(0.5)
soma isyn[2] = new Exp2Syn(0.5)

for i=0, ni-1 {
  isyn[i].tau1 = 10
  isyn[i].tau2 = 50
  isyn[i].e = -80
  // spike generator to inhibitory connections
  dummy ispike[i] = new NetStim(0.5)

  // connect spike generator to synapse
  isconn[i] = new NetCon(ispike[i], isyn[i], -20, 1, iw)
// netcon = new NetCon(source, target, threshould, delay, weight)
//http://www.neuron.yale.edu/neuron/docs/help/neuron/neuron/classes/netcon.html
}
  ispike[0].interval = 5
  ispike[0].number = 10000
  ispike[0].start = 5
  ispike[0].noise = 1
  
  ispike[1].interval = 20
  ispike[1].number = 10000
  ispike[1].start = 5
  ispike[1].noise = 1

  ispike[2].interval = 100
  ispike[2].number = 10000
  ispike[2].start = 5
  ispike[2].noise = 1

//  ispike[i].interval = 100
 // ispike[i].start = 5
 // ispike[i].number = 1000
 // ispike[i].noise = 0.5


// Create shape plot to display synapse positions
// excitatory in red, inhibitory in blue
objref ns
ns = new Shape()
for i=0, ne-1 {
ns.point_mark(esyn[i], 2)
}
//for i=0, ni-1 {
//ns.point_mark(isyn[i], 3)
//}


// Set excitatory weights
proc eweight() {
  for i=0, ne-1 esconn[i].weight = $1
}

// Set inhibitory weights
proc iweight() {
  for i=0, ni-1 isconn[i].weight = $1
}

xpanel("Connection Weights")
xvalue("Exc. Weight", "ew", 1, "eweight(ew)", 0, 0)
xvalue("Inh. Weight", "iw", 1, "iweight(iw)", 0, 0)
xpanel()

// --end of synapse------------------------------------------------------------

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

  st=new IClamp(.5)
  st.dur = 5000
  st.del = 5
  st.amp = .0

  st1=new IClamp(.5)
  st1.dur = 5000
  st1.del = 5
  st1.amp = .0

// --------------------------------------------------------------
// create useful graphs & panels
// --------------------------------------------------------------

nrnmainmenu()
//nrncontrolmenu()

proc newplot() { local wd,ht
    graphItem = new Graph(0)
//    if (numarg()==1) $o1=graphItem // make a pointer
//    if (numarg()==3) $o3=graphItem
//    if (numarg()>1) { wd=$1 ht=$2 } else { wd=500 ht=300 }
   wd=300 ht=200
    graphItem.save_name("graphList[0].")
    graphList[0].append(graphItem)
    graphItem.view(0,-80,tstop,140,900,200,wd,ht)
  }

/*
newplot()
graphItem.addvar("esyn[0].i",3,1)
graphItem.addvar("esyn[1].i",7,1)
graphItem.addvar("esyn[2].i",5,1)

newPlotV()
graphItem.addvar("iseg.v(0)",17,1)
graphItem.addvar("iseg.v(.1)",2,1)
graphItem.addvar("iseg.v(.2)",3,1)
graphItem.addvar("iseg.v(.3)",4,1)
graphItem.addvar("iseg.v(.4)",5,1)
graphItem.addvar("iseg.v(.5)",6,1)
graphItem.addvar("iseg.v(.6)",7,1)
graphItem.addvar("iseg.v(.7)",8,1)
graphItem.addvar("iseg.v(.8)",9,1)
graphItem.addvar("iseg.v(.9)",10,1)
graphItem.addvar("iseg.v(1)",15,1)

graphItem.addvar("node[0].v(.1)",12,1) //plot membrane potentials of different parts of axon
//myelin[1] for (x) print x*L, v(x)

newplot()
graphItem.addvar("soma.ina(0.5)",7,1)
graphItem.addvar("hill.ina(0.5)",1,1)
graphItem.addvar("iseg.ina(0.5)",2,1)
graphItem.addvar("node[0].ina(0.5)",5,1)
*/

newplot()
graphItem.addvar("soma.v(0.5)",1,1)
//graphItem.addvar("hill.v(0.5)",2,1)
graphItem.addvar("iseg.v(0.7)",2,1)
//graphItem.addvar("nod.v(0.2)",4,1)
//graphItem.addvar("nod.v(0.6)",5,1)
//graphItem.addvar("nod.v(0.8)",6,1)
graphItem.addvar("node[0].v(0.5)",3,1)

//graphItem.addvar("nod.v(0.1)",3,1)
//graphItem.addvar("nod.v(0.4)",4,1)
//graphItem.addvar("nod.v(0.6)",5,1)
//graphItem.addvar("nod.v(0.8)",6,1)
//graphItem.addvar("nod.v(0.9)",7,1)
//graphItem.addvar("nod.v(1)",7,1)

objectvar scene_vector_[14]
objectvar save_window_, rvp_
save_window_ = new Graph(0)
save_window_.size(0,350,-80,40)
scene_vector_[5] = save_window_
{save_window_.view(0, -80, 350, 150, 33, 410, 300, 200)}
flush_list.append(save_window_)
save_window_.save_name("flush_list.")
objectvar rvp_
rvp_ = new RangeVarPlot("v")
soma rvp_.begin(0)
nod rvp_.end(1)
//rvp_.origin(11.5382)
save_window_.addobject(rvp_, 2, 1, 0.8, 0.9)  //addobject: Adds the RangeVarPlot to the list of items to be plotted on flush


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

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

//see C:\yale_journal\neuron\nrntuthtml\nrntuthtml\tutorial\tutE.html
objref rect, recv, rec_iseg, rec_hill, rec_axon1, rec_axon3, rec_axon5, rec_axon7, rec_axon9 //used to record values of variables
objref rec_syn1, rec_syn2,rec_syn3, rec_axon0, rec_axon2, rec_axon4, rec_axon6, rec_axon8, rec_axon10, rec_axon11
objref csna, chna, cina, cana1, cana2, cana3, cana4, cana5, cana6, cana7, cana8, cana9, cana10
objref seg0, seg1, seg2, seg3, seg4, seg5, seg6, seg7, seg8, seg9, seg10
objref seg11, seg12, seg13, seg14, seg15, seg16, seg17, seg18, seg19, seg20
objref rec_axon21, rec_axon22, rec_axon23, rec_axon24, rec_axon25, rec_axon26, rec_axon27, rec_axon28, rec_axon29, rec_axon30
objref rec_axon31, rec_axon32, rec_axon33, rec_axon34, rec_axon35, rec_axon36, rec_axon37, rec_axon38, rec_axon39, rec_axon40
objref rec_axon41, rec_axon42, rec_axon43, rec_axon44, rec_axon45, rec_axon46, rec_axon47, rec_axon48, rec_axon49, rec_axon50
objref rec_axon51, rec_axon52, rec_axon53, rec_axon54, rec_axon55, rec_axon56, rec_axon57, rec_axon58, rec_axon59, rec_axon60, rec_axon61
objref rec_caxon0, rec_caxon2, rec_caxon4, rec_caxon6, rec_caxon8, rec_caxon10
objref rec_caxon1, rec_caxon3, rec_caxon5, rec_caxon7, rec_caxon9, rec_caxon11,rec_caxon12,rec_caxon13,rec_caxon14
objref rec_axon20, rec_axon12, rec_axon14, rec_axon16, rec_axon18, rec_axon13, rec_axon11, rec_axon15, rec_axon17, rec_axon19
objref rec_isy1, rec_isy2, rec_isy3,rec_axon31
rec_isy1=new Vector()
rec_isy2=new Vector()
rec_isy3=new Vector()


rect = new Vector()
recv = new Vector()
rec_hill=new Vector()
rec_iseg=new Vector()
rec_axon0=new Vector()
rec_axon1=new Vector()
rec_axon2=new Vector()
rec_axon3=new Vector()
rec_axon4=new Vector()
rec_axon5=new Vector()
rec_axon6=new Vector()
rec_axon7=new Vector()
rec_axon8=new Vector()
rec_axon9=new Vector()
rec_axon10=new Vector()
rec_axon11=new Vector()
rec_axon12=new Vector()
rec_axon13=new Vector()
rec_axon14=new Vector()
rec_axon15=new Vector()
rec_axon16=new Vector()
rec_axon17=new Vector()
rec_axon18=new Vector()
rec_axon19=new Vector()
rec_axon20=new Vector()

seg0=new Vector()
seg1=new Vector()
seg2=new Vector()
seg3=new Vector()
seg4=new Vector()
seg5=new Vector()
seg6=new Vector()
seg7=new Vector()
seg8=new Vector()
seg9=new Vector()
seg10=new Vector()

seg11=new Vector()
seg12=new Vector()
seg13=new Vector()
seg14=new Vector()
seg15=new Vector()
seg16=new Vector()
seg17=new Vector()
seg18=new Vector()
seg19=new Vector()
seg20=new Vector()

rec_axon21=new Vector()
rec_axon22=new Vector()
rec_axon23=new Vector()
rec_axon24=new Vector()
rec_axon25=new Vector()
rec_axon26=new Vector()
rec_axon27=new Vector()
rec_axon28=new Vector()
rec_axon29=new Vector()
rec_axon30=new Vector()

rec_axon31=new Vector()
rec_axon32=new Vector()
rec_axon33=new Vector()
rec_axon34=new Vector()
rec_axon35=new Vector()
rec_axon36=new Vector()
rec_axon37=new Vector()
rec_axon38=new Vector()
rec_axon39=new Vector()
rec_axon40=new Vector()

rec_axon51=new Vector()
rec_axon52=new Vector()
rec_axon53=new Vector()
rec_axon54=new Vector()
rec_axon55=new Vector()
rec_axon56=new Vector()
rec_axon57=new Vector()
rec_axon58=new Vector()
rec_axon59=new Vector()
rec_axon60=new Vector()
rec_axon61=new Vector()

rec_caxon0=new Vector()
rec_caxon1=new Vector()
rec_caxon2=new Vector()
rec_caxon3=new Vector()
rec_caxon4=new Vector()
rec_caxon5=new Vector()
rec_caxon6=new Vector()
rec_caxon7=new Vector()
rec_caxon8=new Vector()
rec_caxon9=new Vector()
rec_caxon10=new Vector()
rec_caxon11=new Vector()
rec_caxon12=new Vector()
rec_caxon13=new Vector()
rec_caxon14=new Vector()
rec_syn1=new Vector()
rec_syn2=new Vector()
rec_syn3=new Vector()
csna=new Vector()
chna=new Vector()
cina=new Vector()
cana1=new Vector()
cana2=new Vector()
cana3=new Vector()
cana4=new Vector()

/*
objref savdata //used for saving file
savdata = new File()  //savdata is a file pointer
savdata.wopen("neuron_soma.dat")  //open a file with a defined name to save data
objref savdata1 //used for saving file
savdata1 = new File()  //savdata is a file pointer
savdata1.wopen("ina.dat")  //open a file with a defined name to save data
*/
objref savdata2 //used for saving file
savdata2 = new File()  //savdata is a file pointer
savdata2.wopen("seg.dat")  //open a file with a defined name to save data

/*
objref savdata3 //used for saving file
savdata3 = new File()  //savdata is a file pointer
savdata3.wopen("cnod1.dat")  //open a file with a defined name to save data

objref savdata4 //used for saving file
savdata4 = new File()  //savdata is a file pointer
savdata4.wopen("input_out.dat")  //open a file with a defined name to save data
*/

/*
//------controlled by Gluct.mod---------------------------------------------------------------
access soma

objref fl
fl = new Gfluct2(0.5)


//fl.std_e = 0.012		// 4 times larger
//fl.std_i = 0.0264

fl.std_e = 0.012/1000		// 4 times larger
fl.std_i =  0.0264/1000
//fl.std_e = 0		// 4 times larger
//fl.std_i = 0
*/

//********************stimulus injection***************************************
soma set_stim1(.5,0,0,200000)
proc soma_inj() {
soma set_stim(.5,0,0,8)
//soma set_stim(.5,2,0,8)
//soma set_stim(.5,0.1,100,300)  //for 1/rm, vinit=-80, 0.1 is the stimulus threshold
//st.loc($1)  st.amp = $2   st.del = $3 (delay)  st.dur = $4
//********************stimulus injection***************************************


recv.record(&soma.v(0.5))   //put value of soma.v to recv
rect.record(&t)             //put value of t to rect 
rec_hill.record(&hill.v(0.5))

seg0.record(&iseg.v(0))
seg1.record(&iseg.v(.125))
seg2.record(&iseg.v(.25))
seg3.record(&iseg.v(.375))
seg4.record(&iseg.v(.5))
seg5.record(&iseg.v(.625))
seg9.record(&iseg.v(0.7)) // this is the initiation point
seg6.record(&iseg.v(.75))
seg7.record(&iseg.v(.875))
seg8.record(&iseg.v(1))


rec_axon1.record(&nod.v(.025))
rec_axon2.record(&nod.v(.05))
rec_axon3.record(&nod.v(0.075))
rec_axon4.record(&nod.v(0.1))
rec_axon5.record(&nod.v(0.125))
rec_axon6.record(&nod.v(0.15))
rec_axon7.record(&nod.v(0.175))
rec_axon8.record(&nod.v(0.2))
rec_axon9.record(&nod.v(0.225))
rec_axon10.record(&nod.v(0.25))

rec_axon11.record(&nod.v(.275))
rec_axon12.record(&nod.v(.3))
rec_axon13.record(&nod.v(.325))
rec_axon14.record(&nod.v(0.35))
rec_axon15.record(&nod.v(0.375))
rec_axon16.record(&nod.v(0.4))
rec_axon17.record(&nod.v(0.425))
rec_axon18.record(&nod.v(0.45))
rec_axon19.record(&nod.v(0.475))
rec_axon20.record(&nod.v(0.5))

rec_axon21.record(&nod.v(.525))
rec_axon22.record(&nod.v(.55))
rec_axon23.record(&nod.v(0.575))
rec_axon24.record(&nod.v(0.6))
rec_axon25.record(&nod.v(0.625))
rec_axon26.record(&nod.v(0.65))
rec_axon27.record(&nod.v(0.675))
rec_axon28.record(&nod.v(0.7))
rec_axon29.record(&nod.v(0.725))
rec_axon30.record(&nod.v(0.75))

rec_axon31.record(&nod.v(.775))
rec_axon32.record(&nod.v(.8))
rec_axon33.record(&nod.v(.825))
rec_axon34.record(&nod.v(0.85))
rec_axon35.record(&nod.v(0.875))
rec_axon36.record(&nod.v(0.9))
rec_axon37.record(&nod.v(0.925))
rec_axon38.record(&nod.v(0.95))
rec_axon39.record(&nod.v(0.975))
rec_axon40.record(&nod.v(1))

rec_axon51.record(&node[0].v(.5))
rec_axon52.record(&node[1].v(.5))
rec_axon53.record(&node[2].v(0.5))
rec_axon54.record(&node[3].v(0.5))
rec_axon55.record(&node[4].v(0.5))
rec_axon56.record(&node[5].v(0.5))
rec_axon57.record(&node[6].v(0.5))
rec_axon58.record(&node[7].v(0.5))
rec_axon59.record(&node[8].v(0.5))
rec_axon60.record(&node[9].v(0.5))
rec_axon61.record(&node[14].v(0.5))

rec_isy1.record(&esyn[0].i)
rec_isy2.record(&esyn[1].i)
rec_isy3.record(&esyn[2].i)


rec_caxon0.record(&cnod[0].v(0))
rec_caxon11.record(&cnod[0].v(.025))
rec_caxon12.record(&cnod[0].v(.05))
rec_caxon13.record(&cnod[0].v(.075))
rec_caxon1.record(&cnod[0].v(.1))
rec_caxon14.record(&cnod[0].v(.15))
rec_caxon2.record(&cnod[0].v(.2))
rec_caxon3.record(&cnod[0].v(0.3))
rec_caxon4.record(&cnod[0].v(0.4))
rec_caxon5.record(&cnod[0].v(0.5))
rec_caxon6.record(&cnod[0].v(0.6))
rec_caxon7.record(&cnod[0].v(0.7))
rec_caxon8.record(&cnod[0].v(0.8))
rec_caxon9.record(&cnod[0].v(0.9))
rec_caxon10.record(&cnod[0].v(1))


csna.record(&soma.ina(0.5))
chna.record(&hill.ina(0.5))
cina.record(&iseg.ina(0.5))
cana1.record(&nod.ina(0.3))
cana2.record(&nod.ina(.6))
cana3.record(&nod.ina(1))


run() // this is the right position

// the following write the data to the file pointer savdata
for i=0,rect.size()-1 {
//    savdata.printf("%g %g\n", rect.x(i), recv.x(i))
savdata2.printf("%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g \n", rect.x(i),recv.x(i),rec_hill.x(i),seg0.x(i),seg1.x(i),seg2.x(i),seg3.x(i),seg4.x(i),seg5.x(i),seg9.x(i),seg6.x(i),seg7.x(i),seg8.x(i),rec_axon1.x(i),rec_axon2.x(i),rec_axon3.x(i),rec_axon4.x(i),rec_axon5.x(i),rec_axon6.x(i),rec_axon7.x(i),rec_axon8.x(i),rec_axon9.x(i),rec_axon10.x(i),rec_axon11.x(i),rec_axon12.x(i),rec_axon13.x(i),rec_axon14.x(i),rec_axon15.x(i),rec_axon16.x(i),rec_axon17.x(i),rec_axon18.x(i),rec_axon19.x(i),rec_axon20.x(i),rec_axon21.x(i),rec_axon22.x(i),rec_axon23.x(i),rec_axon24.x(i),rec_axon25.x(i),rec_axon26.x(i),rec_axon27.x(i),rec_axon28.x(i),rec_axon29.x(i),rec_axon30.x(i),rec_axon31.x(i),rec_axon32.x(i),rec_axon33.x(i),rec_axon34.x(i),rec_axon35.x(i),rec_axon36.x(i),rec_axon37.x(i),rec_axon38.x(i),rec_axon39.x(i),rec_axon40.x(i))

}
savdata2.close()
}
xpanel("Stimulus Injection")
xbutton("inject soma","soma_inj()")
xpanel()