// High frequency (50 Hz) synaptic stimulation of two stimuli (A, C) located
// on two different tip sections individually (A, C) and together (A+C).
// This experiment is used to show how A+B type stimualtion can be sublinear, 
// linear or supralinear depending on the strength and location of stimuli A,B. 

// Inhibition: One GABA_a/GABA_b synapse in each stimulated branch
// The variable "times" is used to select among a set of blockade cases to test
// times: 0 = control , 1 = block_A, 2 =  block_NMDA, 3 = block_A_NMDA, 4 = block_Na, 5 =  block_Ca

times = 0

//load_proc("nrnmainmenu")          
//load_template("ExperimentControl")  // load needed templates
//load_template("EPSPTuning")  
//load_template("RangeRef")         
load_file("nrngui.hoc")
load_file("../../template/load_templates.hoc")

objref econ                         // initialize template parameters
show_errs=1
debug_lev=1
econ=new ExperimentControl(show_errs,debug_lev)
econ.self_define(econ)
econ.morphology_dir = "../../morphology/n123"                                // set location for morphology files
econ.add_lib_dir("Terrence","../../lib")                                     // set location for library files
econ.generic_dir    = "../../experiment/"                                    // set location for cell-setup file

econ.data_dir  = "data"                                                     // set directory to store data
sprint(econ.syscmd, "mkdir -p %s", econ.data_dir)
system(econ.syscmd)
                                   // maximum nseg number
actual_resolution=75               // used in ..lib/choose-secs.hoc (not currently) to count 
desired_resolution=1               // how many copies  of a given synapse to put in the band

econ.xopen_geometry_dependent("cell")                                        // load raw cell morphology	
econ.xopen_geometry_dependent("cell-analysis")                               // load user-defined semantics on morphology
cell_analysis(econ)

printf("Opening cell setup\n")                                               // load cell-setup to	
econ.xopen_generic("cell-setup")                                             // specify all mechanisms,	
printf("Opened. Setting up cell\n")                                          // membrane properties etc
maximum_segment_length=actual_resolution
cell_setup(econ)

// Set simulation parameters for the experiment 

econ.defvar("Simulation Control", "tstop", "250", "Defines when the simulation stops.")
econ.defvar("Simulation Control", "dt", "0.1", "Timestep")
econ.defvar("Simulation Control", "steps_per_ms", "10", "How many points are plotted per ms")
setdt()

// open files with NMDA/AMPA, GABA_A/AMPA and GABA_B/GABA_A ratios
econ.xopen_geometry_dependent("nmda-ampa-ratio")
econ.xopen_geometry_dependent("gabab-gabaa-uniform-ratio")

// Open file with tuned AMPA conductance values for all sections
strdef tunings_file, select, temp, accstr
objref tune_epsp_list
tune_epsp_list=new List()
sprint(tunings_file, "%s", "tunings")
xopen("../tune-synapses/tunings.dat")

// Open library functions that will be needed

econ.xopen_library("Terrence","choose-secs")    // used to randomly select sections from a list
econ.xopen_library("Terrence","salloc")         // used to allocate synapses on sections
econ.xopen_library("Terrence","deduce-ratio")   // used to extract NMDA/AMPA, GABA_A/AMPA and GABA_B/GABA_B ratios
econ.xopen_library("Terrence","basic-graphics") // used to plot graphics 
econ.xopen_library("Terrence","spikecount")     // used to count spikes


temporal_offset=10  // synapses are stimulated simultaneously after 10ms

all_synapses = 100     // Maximum number of AMPA/NMDA synapses
gaba_synapses = all_synapses
objref ampa[all_synapses], nmda[all_synapses], gabaa[gaba_synapses], gabab[gaba_synapses], splot
objref somavrec, vf, vtip[2], somarecf, dendrecf
objref tip_list, vsoma, rpid, sl, Branch_ref[all_synapses]
strdef recordsec

double synapses[2], maxdendrec[2], meandendrec[2]
max_cases = 26
double all_syn[2*max_cases]
min_syn = 2
max_syn = 15

//create all synapse pairs for stimuli A, B

all_syn[0] = 2
all_syn[1] = 2
all_syn[2] = 2
all_syn[3] = 4
all_syn[4] = 4
all_syn[5] = 2
all_syn[6] = 3
all_syn[7] = 3
all_syn[8] = 3
all_syn[9] = 2
all_syn[10] = 2
all_syn[11] = 3
all_syn[12] = 4
all_syn[13] = 4
all_syn[14] = 5
all_syn[15] = 5
all_syn[16] = 5
all_syn[17] = 2
all_syn[18] = 2
all_syn[19] = 5
all_syn[20] = 5
all_syn[21] = 4
all_syn[22] = 4
all_syn[23] = 5
all_syn[24] = 6
all_syn[25] = 6
all_syn[26] = 6
all_syn[27] = 4
all_syn[28] = 4
all_syn[29] = 6
all_syn[30] = 6
all_syn[31] = 5
all_syn[32] = 5
all_syn[33] = 6
all_syn[34] = 4
all_syn[35] = 8
all_syn[36] = 8
all_syn[37] = 4
all_syn[38] = 8
all_syn[39] = 5
all_syn[40] = 5
all_syn[41] = 8
all_syn[42] = 7
all_syn[43] = 7
all_syn[44] = 8
all_syn[45] = 6
all_syn[46] = 6
all_syn[47] = 8
all_syn[48] = 7
all_syn[49] = 8
all_syn[50] = 9
all_syn[51] = 6

// Make a list of all trunk sections to choose sections from
tip_list=new SectionList()
forsec apical_tip_list {
       tip_list.append()
}
forsec apical_tip_list_addendum {
       tip_list.append()
}

// Use these two sections to generate A,C, A+C traces for figure 7 in 
// Poirazi. P, Brannon T. and Mel B.W, 'Arithmetic of Subthreshold Synaptic Summation in a Model'

//tip_list=new SectionList()
//apical_dendrite[82]  tip_list.append()
//apical_dendrite[51]  tip_list.append()

// make a reference list for all sections used
m=0
forsec tip_list {
   Branch_ref[m] = new SectionRef() 
   m = m + 1
}

Deadtime_GLU=dt
Deadtime_NMDA=dt
hertz=50				      // frequency of stimulation for all synapses 
synch=0					      // synapses are stimulated randomly (NOT synchronously)
perio=0					      // spike trains for each synapse are NOT periodic
dendritic_spike_threshold=-25
somatic_spike_threshold=0 

//Proceedures for the different cases to be tested

proc Ca_block(){   // Block all Ca++ channels
  forall if(ismembrane("cat")) { 
     for (x) { gcatbar_cat(x) = 0 }
  }
  forall if(ismembrane("calH")) { 
     for (x) { gcalbar_calH(x) = 0 }
  }
  forall if(ismembrane("cal")) { 
     for (x) { gcalbar_cal(x) = 0 }
  }
  forall if(ismembrane("car")) { 
     for (x) { gcabar_car(x) = 0 }
  }
  forall if(ismembrane("somacar")) { 
     for (x) { gcabar_somacar(x) = 0 }
  }
}

proc Na_block(){  // Block all Na+ channels
  forall if(ismembrane("hha2")) { 
     for (x) { gnabar_hha2(x) = 0 }
  }
  forall if(ismembrane("hha_old")) { 
     for (x) { gnabar_hha_old(x) = 0 }
  }
}

proc NMDA_block(){  // Block NMDA current 
 for i=0, synapses[0] + synapses[1] -1 {
   nmda[i].gmax = 0  
 }
}

proc A_block() {  // Block all A-type K+ channels 
f = 0.2
forall if(ismembrane("kad")) { // distal conductances
     for(x) { gkabar_kad(x)= gkabar_kad(x)*f }
  } else if(ismembrane("kap")) { // proximal conductances
     for(x) { gkabar_kap(x)= gkabar_kap(x)*f }
  }
}

proc A_NMDA_block(){ // block both A-current and NMDA current
 NMDA_block()
 A_block()
}

if (times == 0 || times == 2) {
  temporal_offset = 10
}else {
  temporal_offset = 120
}
tstop = tstop + temporal_offset 

addgraph_2("soma.v(0.5)",0,tstop, -72,-10)

r = 0
for try = 1, 30  {   // try selecting a pair of trunk sections for 30 times
   r = r + 1 
   sl = new SectionList()
   rpid = new Random(r-1)
   PID = rpid.discunif(0,m) 
   access Branch_ref[PID].sec    
   //apical_dendrite[51] sl.append() //for figure 7
   secname() sl.append
   sprint(recordsec, "%s.v(0.5)",secname()) 
   addgraph_2(recordsec,0,tstop,-72,-10) // plot trace of current branch
  
   rpid = new Random(r+1) 
   PID = rpid.discunif(0,m)
   access Branch_ref[PID].sec  
   //apical_dendrite[82] sl.append() //for figure 7
   secname() sl.append
   sprint(recordsec, "%s.v(0.5)",secname()) 
   addgraph_2(recordsec,0,tstop,-72,-10) // plot trace of current branch

   
   for c = 0, max_cases-1 {   // for the sections selected, create all possible    
     k=-1                     // A, B stimulation cases
     splot=new Shape()        // make a shape graph 
     all_synapses = 0         // initialize  
     gaba_synapses = 0  
     strdef dir_str
     sprint(dir_str, "data/%s/Apical_Tips/Freq/Both_Tips/",tunings_file) // define data directory
     forsec sl {              // for both selected sections
    
       k = k+1
       synapses[k]= all_syn[2*c+k]  // get number of synapses to put in section 
       //synapses[k]= 6    // for figure 7
       
       print secname(), "is where we assess potency."
       COLOR=k+2
       t = 0
       if ( k > 0) { t = k*synapses[k-1] }
       nseg = synapses[k]
       for si=1,synapses[k] {       // uniformly distribute synapses on sections
            posn = (2*si -1)/(2*synapses[k]) 
            ampa[t+si-1] = new GLU(posn)  
            nmda[t+si-1] = new NMDA(posn)  
            salloc(ampa[t+si-1],nmda[t+si-1],posn)
            splot.point_mark(ampa[t+si-1],COLOR)    // mark synapses on shape graph 
                       
       }
       splot.flush()
       splot.show(1)
       sprint(dir_str, "%s%s",  dir_str, secname()) // set directory name to sections name 
       all_synapses = all_synapses + synapses[k]    // update number of synapses allocated 


      //Add one GABAa/GABAb synapse in the center of the current branch
        posa = 0.5
        posb = 0.5
	gabaa[k] = new GABAa(posa)
    	SALLOC_GABAa(gabaa[k],posa,1,0)        
        gabab[k] = new GABAb(posb)
        SALLOC_GABAb(gabab[k],posb,1,0)
        gaba_synapses = gaba_synapses + 1

      }
    
    // Plot gaba synapses on the same shape graph
       for i = 0, gaba_synapses -1 {
          splot.point_mark(gabaa[i],COLOR+2)
          splot.point_mark(gabab[i],COLOR-2)
       }
       splot.flush()
       splot.show(1)
    
     GABA_flag = 0    // Don't make both AMPA/NMDA and GABA trains in shiftsyn_init        
     // create the stimulation trains for AMPA & NMDA synapses   
     econ.xopen_library("Terrence","shiftsyn-initA")
     shiftsyn_init(all_synapses,tstop,dt,hertz,synch,perio,PID,temporal_offset,GABA_flag,"ampa","nmda")
     
     // Create the stimulation trains for GABA synapses
     econ.xopen_library("Terrence","GABA_shiftsyn")  
     gaba_shift(gaba_synapses,tstop,dt,hertz,synch,perio,PID,temporal_offset, "gabaa","gabab")
    
     // Execute current-blockade proceedure specified 
     
     if (times == 0) {
             sprint(select, "%s", "control")}  
        if (times == 1) {
             A_block()
             sprint(select, "%s", "block_A")}
        if (times == 2) {
             NMDA_block() 
             sprint(select, "%s", "block_NMDA")}
        if (times == 3) {
             A_NMDA_block() 
             sprint(select, "%s", "block_A_NMDA")}
        if (times == 4) {
             Na_block()
             sprint(select, "%s", "block_Na")}
        if (times == 5) {
             Ca_block()
             sprint(select, "%s", "block_Ca")}

     print "secname = ", secname(), " synapses[1] = ", synapses[0], " synapses[2] = ", synapses[1], "case = ", select
  
     somavrec=new Vector(tstop/dt)
     somavrec.record(&soma.v(0.5)) // prepare to record at soma

     ind=0
     forsec sl {                                        // prepare to record at both dendrites
       vtip[ind]=new Vector(tstop/dt)
       sprint(temp, "vtip[%d].record(&%s.v(0.5))", ind, secname())
       execute1(temp)
       ind=ind+1
     }
   
     finitialize(v_init)           // Initialize and run experiment
     fcurrent()
     run()
  
      meanvrec = somavrec.mean(temporal_offset/dt,tstop/dt -1)      // mean somatic depolarization
      maxvrec = somavrec.max(temporal_offset/dt,tstop/dt -1)        // max somatic depolarization
      ind = 0
      forsec sl {
         meandendrec[ind] = vtip[ind].mean(temporal_offset/dt,tstop/dt -1) // mean dendritic depolarization
         maxdendrec[ind] = vtip[ind].max(temporal_offset/dt,tstop/dt -1)   // max dendritic depolarization
         ind = ind + 1      
      }
      soma_spikes =  spikecount(somavrec,somatic_spike_threshold)     // count current number of somatic spikes
        
     sprint(econ.syscmd,  "mkdir -p %s", dir_str)                     // make the output directory
     system(econ.syscmd) 

     vf=new File()
     sprint(econ.tmp_str2, "%s/%s_freq_%d_%d", dir_str, select, min_syn, max_syn)             // make the output file
     
     // Print number of synapses in each section, mean and max somatic voltage and number of spikes
     vf.aopen(econ.tmp_str2)
     vf.printf("%d %g %g %g %g %g %g %g %g\n", synapses[0], synapses[1], meanvrec, maxvrec, meandendrec[0], maxdendrec[0], meandendrec[1], maxdendrec[1], soma_spikes)
     vf.close()
        
    //Print somatic traces for stimulus conditions A+C

     sprint(econ.tmp_str, "data/%s/Apical_Tips/Freq/Both_Tips/%s/",tunings_file, secname())
     somarecf=new File()
     sprint(temp, "%s/Vsoma_%d_%d", dir_str, synapses[0], synapses[1])
     somarecf.wopen(temp)
     somavrec.printf(somarecf, "%g\n")
     somarecf.close()
        
   //Print graphics to eps files
     for i=0,windex {
       sprint(econ.tmp_str, "%s", dir_str)
       sprint(econ.tmp_str2, "%s/graph-%d.eps",econ.tmp_str, i)      
       win[i].printfile(econ.tmp_str2)
     }
  }
}