//--------------------------------------------------------------------
// Simulation of  a single GDP
// Simulation of AMPA and GABA events with defined correlation
//---------------------------------------------------------------------

// ------------Definition of Parameters -------------------------------
// --------------------------------------------------------------------


// Define Seed Values

  seed_AMPA = 8
  seed_GABA = seed_AMPA +2

// Define Cl--Concentration
   Cl_Steps = 6  // Number of Different [Cl-]i  
   objref Cl_List
   Cl_List = new Vector(Cl_Steps)
   Cl_List.x[0] = 5 //mM
   Cl_List.x[1] = 10
   Cl_List.x[2] = 15
   Cl_List.x[3] = 25
   Cl_List.x[4] = 35
   Cl_List.x[5] = 50


// Determination Parameters AMPA ------------------------------------
G_AMPA =  0.000305 * 1
DECAY_AMPA = 11
nampasyn = 100 //Number of AMPA-Synapses  // calculated was 107
aninputs = 1 // number of inputs per synapse


// Determining Parameters GABA allone ---------------------------
G_GABA = 0.000789   // synaptic weight according to miniature events 
DECAY_GABA = 37
P_GABA = 0.18
ngabasyn = 100  // from experiments = 534
gninputs = 1 // number of inputs per synapse


// Definition of various runtime parameters --------------------------

lenghtoutputvec = 6000  // Number of Lines for output (< 32000 for Excel-Figures)

ndend=56  // Number of dendrites

tstop = 1500  // Duration
v_init = -60   // Initial voltage
dt = 0.01       // Step Interval in ms

// definition of boolean variables
TRUE = 1
FALSE = 0
AmpaCtrl_Steps = 2   // 0 = without AMPA, 1 = With AMPA


// ------------Procedures and Functions -------------------------------
// --------------------------------------------------------------------

// Function MakeShort ---------------------------------------//
// Inputs: $1 Objref to Inputvector                          //
//         $2 Objref to Outoutvector                         //
//         $3 lenoutvec  desired lendth of Outputvector      //
//                                                           //
// Reduce Inputvec to Outputvev by averaging n elements      //
// n (reducing factor) = floor(Inputvec.size() / lenoutvec)  //
// ----------------------------------------------------------//

obfunc MakeShort() {local i, n

  n = int($o1.size()/$3)
  $o2.resize($3)
  for i=0, $3-1 {
    $o2.x[i] = $o1.mean(i*n, (i+1)*n-1)
  }
  return $o2
}   //  End of function





// ---------Definition of objects -------------------------------------
// --------------------------------------------------------------------


// Objects for Synapses ---------------------------------------------------------
objref gabasyn[ngabasyn]                            // Definition of synapse objects
objref ampasyn[nampasyn]


// random function for localization of synapses
objref rand_ampa_loc, rand_gaba_loc

// random function for localization of synapses in which dendrite
objref rand_ampa_dend, rand_gaba_dend

// random function for synapses parameters
objref rand_gaba_t, rand_gaba_g
objref rand_ampa_t, rand_ampa_g

// definition of Vectors for Gaba-Stimulation (t_vec = timestamps t_vecr = sorted timestamps, g_vec = rel conductance)
objref gabastim[ngabasyn], gaba_t_vec[ngabasyn], gaba_t_vecr[ngabasyn],synpulsegaba[ngabasyn], gaba_g_vec[ngabasyn]
// definition of Vectors for AMPA-Stimulation
objref ampastim[nampasyn], ampa_t_vec[nampasyn], ampa_t_vecr[nampasyn], synpulseampa[nampasyn], ampa_g_vec[nampasyn]


// Define vectors to link modelled parameter output ---------------------------------
objref timevec, voltvec, clivec[ndend]                                    // vectors linked to parameter-pointers
objref clivec_aver                                                        // vectors for Averagees over all Dendrites
objref shorttimevec, shortvoltvec, shortclivec, spacevec                  // shorter Vectors for output

// Matrix for output 0 = time, 1 = Voltage, 2 = average Cli , 3 to 3+ndend, Cli n Dend[i]
objref Outmatrix, TimeStampMatrix
objref ampa_timestamp, gaba_timestamp

// Define Name of Output-File
strdef OutFileName, OutTimestampFileName

// Define Output File
objref OutFile, OutTimestampFile


// Generate vectors and matrices -------------------------------------
voltvec = new Vector()
timevec = new Vector()
clivec_aver = new Vector()
shortvoltvec = new Vector()
shorttimevec = new Vector()
shortclivec = new Vector()
spacevec = new Vector()
for i=0, ndend-1 {
  clivec[i] = new Vector()
}
Outmatrix = new Matrix(lenghtoutputvec+1, Cl_Steps*AmpaCtrl_Steps*3+1)
TimeStampMatrix = new Matrix(ngabasyn+1, Cl_Steps*3+1)


// Start of Input generation -------------------------------------------

// Initialize Random Functions -----------
  rand_ampa_loc = new Random(seed_AMPA)
  rand_gaba_loc = new Random(seed_GABA)
  rand_ampa_dend = new Random(seed_AMPA)
  rand_gaba_dend = new Random(seed_GABA)
  rand_ampa_t = new Random(seed_AMPA)
  rand_ampa_g = new Random(seed_AMPA)
  rand_gaba_t = new Random(seed_GABA)
  rand_gaba_g = new Random(seed_GABA)

//Define properties of random Function
  rand_gaba_t.normal(600, 20000) //9000 was original from experimental results
  rand_gaba_g.normal(1, 0.28) //rel variance of GABA according to experimental results
  rand_ampa_t.normal(600, 20000) //8500 was original from experimental results
  rand_ampa_g.normal(1, 0.26) //rel variance of AMPA according to experimental results

// generate Vectors --- (gniputs, aninputs defines number of inputs per synapse) ------
  for i = 0, ngabasyn- 1 {
    gaba_t_vec[i] = new Vector(gninputs)
    gaba_t_vecr[i] = new Vector(gninputs)
    gaba_g_vec[i] = new Vector(gninputs)
  }
  gaba_timestamp = new Vector(ngabasyn)

  for i = 0, nampasyn-1 {
    ampa_t_vec[i] = new Vector(aninputs)
    ampa_t_vecr[i] = new Vector(aninputs)
    ampa_g_vec[i] = new Vector(aninputs)
  }
  ampa_timestamp = new Vector(nampasyn)


  // Distribute GABA synapses -----------------------------------------------------------
  for k=0, ngabasyn - 1 {
      pos = rand_gaba_loc.uniform(0,ndend-1)
      pos2 = rand_gaba_dend.uniform (0.001, 0.999)
      dend_0[pos]{
        gabasyn[k] = new gaba(pos2)
        gabasyn[k].tau1 = 0.1
        gabasyn[k].tau2 = DECAY_GABA
        gabasyn[k].P = P_GABA
      }
   }

  // distribute AMPA synapses ------------------------------------------------------------------------
  for k=0, nampasyn-1 {
    pos = rand_ampa_loc.uniform(0,ndend-1)
    pos3 = rand_ampa_dend.uniform (0.001, 0.999)
    dend_0[pos]{
      ampasyn[k] = new Exp2Syn(pos3)
      ampasyn[k].tau1 = 0.1
      ampasyn[k].tau2 = DECAY_AMPA
    }                                              
  }

//-- Simulation starts here -----------------------------------------------------------------
//-------------------------------------------------------------------------------------------

//-- Outer Loop Variation of Cl- --------------------------------------------------
  Cl_Step = 0
  while (Cl_Step < Cl_Steps){

    // 1. define Cl- concentration in all elements of the cell ------------------------------------
    forsec all {
      cli0_cldif_CA3_NKCC1_HCO3 = Cl_List.x[Cl_Step]
      cli_Start_cldif_CA3_NKCC1_HCO3 = Cl_List.x[Cl_Step]
      cli_cldif_CA3_NKCC1_HCO3 = Cl_List.x[Cl_Step]
    }

   // 2a. Generate timestamps/conductances for AMPA and GABA synapses (at same timepoint) --------------------------------------
    for f=0, nampasyn-1 {
      for i=0, aninputs-1 {
          t = abs(rand_ampa_t.repick())
          ampa_timestamp.x[f] = t
          gaba_timestamp.x[f] = t 
          ampa_t_vec[f].x[i]= t
          gaba_t_vec[f].x[i]=t         
      }
    }
    ampa_timestamp.sort()
    for f=0, nampasyn-1 {                     
      ampa_t_vecr[f] = ampa_t_vec[f].sort()
    }     
    gaba_timestamp.sort()
    for f=0, ngabasyn - 1 {
      gaba_t_vecr[f] = gaba_t_vec[f].sort()
    }

  // 2d. Generate conductances for GABA synapses --------------------------------------
    for f=0, ngabasyn-1 {
      for i=0, gninputs-1 {
          g = rand_gaba_g.repick()
          gaba_t_vec[f].x[i]= gaba_timestamp.x[f]
          gaba_g_vec[f].x[i]= abs(G_GABA*g)
      }
    }

   // 2d. Generate conductances for AMPA synapses --------------------------------------
    for f=0, nampasyn-1 {
      for i=0, aninputs-1 {
          g = rand_ampa_g.repick()
          ampa_g_vec[f].x[i]=abs(G_AMPA*g)
      }
    }

  // Save the vectors to analyze the temporal correlation
    ActCol = (Cl_Step)*3
    TimeStampMatrix.x[0][ActCol] = Cl_List.x[Cl_Step]
    TimeStampMatrix.x[1][ActCol] = seed_AMPA
    TimeStampMatrix.setcol(ActCol+1,  gaba_timestamp)
    TimeStampMatrix.setcol(ActCol+2,  ampa_timestamp)

  // 3. generate Vecstim-vectors from the sorted timestamp-vectors -------------------------------
    for i=0, ngabasyn -1 {
      gabastim[i] = new VecStim()
      gabastim[i].play(gaba_t_vecr[i])    // GABA stimulator
    }                                               

    for i=0, nampasyn-1{
      ampastim[i] = new VecStim()
      ampastim[i].play(ampa_t_vecr[i])    // AMPA stimulator
    }    
  
    // Run the simulation for each stimulation pattern once with and once without AMPA co-stimulation
  
    for AmpaCtrl_Step = 0, AmpaCtrl_Steps-1 {
      // 4. Play the Vecstim objects to the synapses ---------------------------------------------
      for i=0, ngabasyn-1 {
        synpulsegaba[i] = new NetCon(gabastim[i], gabasyn[i], 0, 0, gaba_g_vec[i].x[0])
      }                                                  // GABA NetCon
     
      for i=0, nampasyn-1 {
        synpulseampa[i] = new NetCon(ampastim[i], ampasyn[i], 0, 0, ampa_g_vec[i].x[0]*AmpaCtrl_Step)
      }    
                                             // Ampa NetCon

      printf("%g of %g;",(Cl_Step*AmpaCtrl_Steps+AmpaCtrl_Step+1), (Cl_Steps*AmpaCtrl_Steps))
      printf("seed- %g;[Cl-]i=%gmM;AMPA=%g:",seed_AMPA, Cl_List.x[Cl_Step], AmpaCtrl_Step)
      // 5. Link Objects to Output-Vectors -----------------------------------
      timevec.record(&t)   // Time vector
      voltvec.record(&v(.5)) // Volt vector in soma
      for i=0, ndend-1 {
        clivec[i].record(&dend_0[i].cli(0.5))
      }

      // 6. Run Simulation --------------------------------------------------------
      run()
 
      // 7. Put Data in Output Vector ------------------------------------------------------  
      MakeShort(timevec, shorttimevec, lenghtoutputvec)
      Outmatrix.setcol(0, shorttimevec)
      spacevec.resize(shorttimevec.size()+1) //spavevec caries information about the properties during a loop)
      spacevec.fill(0) 
      spacevec.x[0] = Cl_List.x[Cl_Step]
      spacevec.x[1] = seed_AMPA
      spacevec.x[2] = G_AMPA*1000*AmpaCtrl_Step
      spacevec.x[3] = 777
      ActCol = (Cl_Step*AmpaCtrl_Steps+AmpaCtrl_Step)*3
      Outmatrix.setcol(ActCol+1, spacevec)
      MakeShort(voltvec, shortvoltvec, lenghtoutputvec)
      Outmatrix.setcol(ActCol+2, shortvoltvec)

      // Calculate average [Cli] in dendrites -------
        clivec_aver.mul(0) // empty vector
        for i=0, ndend-1 { 
          clivec_aver.resize(clivec[i].size())
          clivec_aver.add(clivec[i])
        } 
        clivec_aver.div(ndend)
      // - end of averaging

      MakeShort(clivec_aver, shortclivec, lenghtoutputvec)
      Outmatrix.setcol(ActCol+3, shortclivec)
      printf(", max [Cl-]i %g, min [Cl-]i %g \n", clivec_aver.max, clivec_aver.min)
    } // End of forr loop with/without AMPA
     
    Cl_Step+=1  // Goto next Cl- Concentration
  }  // End of outer loop --------------------------------- 


  // Save the Data --------------------------------------------------------------------
  OutTimestampFile = new File()
  sprint(OutTimestampFileName, "TimeStamps_gGABA-0.789_nGABA-100_gAMPA-0.305_Div-Cl_100percCorrel_seed-%g.asc",seed_AMPA)
  OutTimestampFile.wopen(OutTimestampFileName)
  TimeStampMatrix.fprint(OutTimestampFile, "\t%.12g") 
  OutTimestampFile.close

  OutFile = new File()
  sprint(OutFileName, "Result_Single_GDP_gGABA-0.789_nGABA-100_gAMPA-0.305_Div-Cl_100percCorrel_seed-%g.asc",seed_AMPA)
  OutFile.wopen(OutFileName) 
  Outmatrix.fprint(OutFile, "\t%.12g")     
  OutFile.close
// - end of save data ---------------------------------------------------------------