#cp calcFIplanes_combe.py sim_baclofenwashout_combe.py
#cp calcFIplanes_withbaclofen_blocked_givenst2_givenT_scalesomaticI_varconc.py  calcFIplanes_combe.py 
#cp calcifcurves_hay_withGABAinputs.py calcifcurves_hay_withGABAinputs_inline.py # don't use the functions of the template, try to make them inline instead
#cp ../modulhcn_hay/calcifcurves.py calcifcurves_hay.py
#cp ../haymod3e/runcontrol_check.py runcontrol.py
# runcontrols
# A script for determining the control neuron F-I curve and limit cycle.
#
# The input code for the hoc-interface is based on BAC_firing.hoc by Etay Hay (2011)
#
# Tuomo Maki-Marttunen, Oct 2014
# (CC BY)

from neuron import h
from neuron import rxd
import mytools
import pickle
import numpy as np
import sys
import time
from pylab import *


GIRKcond = 0.008
tstop = 60000
T = 1000
baclofen_max = 100.0 #uM; the maximal agonist concentration
baclofen_dur = 25000.0
somaticIscale = 1.0
dodraw = 1
blockedStr = ''
params_all = ['107.034,514.651,167.633,11.026,66.931,3.009,1.904,2.241,1507.740','91.805,846.811,67.856,9.858,71.217,5.000,1.817,2.170,3000.000','133.092,528.522,185.581,8.574,25.868,3.186,2.125,2.729,2815.426','108.557,729.183,292.710,2.473,16.923,4.406,2.127,2.973,696.729','111.531,708.687,165.616,10.797,93.218,2.903,1.728,2.235,2797.567','109.869,469.365,114.137,59.053,73.263,3.573,2.385,3.028,2852.993','84.387,57.393,126.013,4.253,57.544,5.000,1.722,2.122,929.204','126.993,135.047,348.945,3.733,18.257,3.800,2.647,3.237,2804.842','105.438,642.913,142.107,19.602,40.482,3.599,2.391,2.974,1386.488','91.403,459.570,97.658,3.443,65.833,5.000,1.540,1.942,2447.114','71.410,380.140,143.521,2.991,45.796,10.000,2.044,2.836,2126.204','146.055,702.846,310.725,3.041,13.076,3.334,2.351,3.166,2112.405','90.417,993.718,98.403,6.879,95.287,3.944,1.572,2.204,1413.654','119.681,591.014,91.960,89.687,88.948,3.453,2.300,2.755,2950.763','123.170,616.560,117.913,74.842,46.625,3.952,2.609,3.283,3000.000']

params_str = 'p1'
GIRKcond_single = 33e-3 #nS #https://www.jneurosci.org/content/jneuro/25/15/3787.full.pdf

if len(sys.argv) > 1:
  GIRKcond_single = float(sys.argv[1])
if len(sys.argv) > 2:
  tstop = int(float(sys.argv[2]))
if len(sys.argv) > 3:
  T = int(float(sys.argv[3]))
if len(sys.argv) > 4:
  baclofen_max = float(sys.argv[4])
if len(sys.argv) > 5:
  baclofen_dur = float(sys.argv[5])
if len(sys.argv) > 6:
  blockedStr = sys.argv[6]
if len(sys.argv) > 7:
 params_str = sys.argv[7]
if len(sys.argv) > 8:
  somaticIscale = float(sys.argv[8])
if len(sys.argv) > 9:
  dodraw = int(float(sys.argv[9]))

if blockedStr == 'None':
  blockedStr = ''
if params_str[0] == 'p':
  iparam = int(params_str[1:])
  params = {'k[0]': float(params_all[iparam].split(',')[0]), 'k[1,9]': float(params_all[iparam].split(',')[1]), 'k[3,4,5,6,7,8]': float(params_all[iparam].split(',')[2]), 'k[15,17,19,21,23]': float(params_all[iparam].split(',')[3]), 'k[16,18,20,22,24]': float(params_all[iparam].split(',')[4]), 'OnlyExp0_RGS': float(params_all[iparam].split(',')[5]), 'OnlyExp1_RGS': float(params_all[iparam].split(',')[6]), 'OnlyExp2_RGS': float(params_all[iparam].split(',')[7]), 'gaba_flux': float(params_all[iparam].split(',')[8])}
else:
  params = {'k[0]': float(params_str.split(',')[0]), 'k[1,9]': float(params_str.split(',')[1]), 'k[3,4,5,6,7,8]': float(params_str.split(',')[2]), 'k[15,17,19,21,23]': float(params_str.split(',')[3]), 'k[16,18,20,22,24]': float(params_str.split(',')[4]), 'OnlyExp0_RGS': float(params_str.split(',')[5]), 'OnlyExp1_RGS': float(params_str.split(',')[6]), 'OnlyExp2_RGS': float(params_str.split(',')[7]), 'gaba_flux': float(params_str.split(',')[8])}

  
for icell in range(0,1):
  v0 = -80
  ca0 = 0.0001

  distalpoint = 300

  h("""
load_file("stdlib.hoc")
load_file("stdrun.hoc")
objref cvode
cvode = new CVode()
cvode.active(1)
cvode.atol(1e-7)
cvode.maxstep(5)
load_file("loadcell.hoc")
objref st1, vclamp, vclamp_irec
soma st1 = new IClamp(0.5)
soma vclamp = new SEClamp(0.5)

objref vsoma, vdend, cadend, casoma
vsoma = new Vector()
casoma = new Vector()
vdend = new Vector()
cadend = new Vector()
objref sl,ns,tvec
tvec = new Vector()
sl = new List()
double siteVec[2]
sl = locateSites("apic","""+str(distalpoint)+""")
maxdiam = 0
for(i=0;i<sl.count();i+=1){
  dd1 = sl.o[i].x[1]
  dd = apic[sl.o[i].x[0]].diam(dd1)
  if (dd > maxdiam) {
    j = i
    maxdiam = dd
  }
}
siteVec[0] = sl.o[j].x[0]
siteVec[1] = sl.o[j].x[1]
apic[siteVec[0]] cvode.record(&v(siteVec[1]),vdend,tvec)
apic[siteVec[0]] cvode.record(&cai(siteVec[1]),cadend,tvec)
soma cvode.record(&v(0.5),vsoma,tvec)
soma cvode.record(&cai(0.5),casoma,tvec)
""")
  #additionalblockeds = ['gCa_HVAbar_Ca_HVA','gCa_LVAstbar_Ca_LVAst','gIhbar_Ih','gImbar_Im','gK_Pstbar_K_Pst','gK_Tstbar_K_Tst','gNaTa_tbar_NaTa_t','gNaTs2_tbar_NaTs2_t','gNap_Et2bar_Nap_Et2','gSK_E2bar_SK_E2','gSKv3_1bar_SKv3_1']
  myMechs =  ['Ca_HVA','Ca_LVAst','Ih','Im','K_Pst','K_Tst','NaTa_t','Nap_Et2','SK_E2','SKv3_1','decay_CaDynamics_E2','gamma_CaDynamics_E2','minCai_CaDynamics_E2','']

  species = ['gaba', 'gabaOut', 'GABABR', 'gabaGABABR', 'Gi', 'GABABRGi', 'gabaGABABRGi', 'gabaGABABRGibg', 'GiaGTP', 'GiaGDP', 'RGS', 'GiaGTPRGS', 'Gibg', 'GIRK', 'VGCC', 'GIRKGibg', 'GIRKGibg2', 'GIRKGibg3', 'GIRKGibg4', 'VGCCGibg']
  tolstochange = 'Gi,GiaGDP,Gibg'.split(',')
  tolschange = [float(x) for x in '2,2,2'.split(',')]
  tolscales = [1.0 for i in range(0,len(species))]
  for itol in range(0,len(tolstochange)):
    for ispec in range(0,len(species)):
      if tolstochange[itol] == species[ispec]:
        tolscales[ispec] = tolscales[ispec]*10**(-tolschange[itol])

  Napic = len(h.apic)
  Nbasal = len(h.dend)
  #Insert biochemistry:
  cyts = []
  seclist = []
  dist2s = []
  for iiapic in range(0,Napic):
    mysec = h.apic[iiapic]
    seclist.append(mysec)
    #print(mysec.name()+""" insert GIRK_Yim""")
    h(mysec.name()+""" insert GIRK_Yim""")
    cyts.append(rxd.Region(mysec, name='cyt_apic'+str(iiapic), nrn_region='i'))
    coords = [[mysec.x3d(i),mysec.y3d(i)] for i in range(0,mysec.n3d())]
    dist2s.append((mean(coords,axis=0)[0]-0)**2 + (mean(coords,axis=0)[1]-650)**2)
  for iibasal in range(0,Nbasal):
    mysec = h.dend[iibasal]
    seclist.append(mysec)
    #print(mysec.name()+""" insert GIRK_Yim""")
    h(mysec.name()+""" insert GIRK_Yim""")
    cyts.append(rxd.Region(mysec, name='cyt_basal'+str(iibasal), nrn_region='i'))
    coords = [[mysec.x3d(i),mysec.y3d(i)] for i in range(0,mysec.n3d())]
    dist2s.append((mean(coords,axis=0)[0]-0)**2 + (mean(coords,axis=0)[1]-650)**2)
  #cyt = rxd.Region(seclist, name='cyt', nrn_region='i')
  print("cyt OK")
  iGIRK = [i for i in range(0,len(species)) if species[i] == 'GIRK'][0] #13    
  
  h("""
  tstop = """+str(tstop)+"""
  """)
  #//distributeSyn()
  
  ks = [1.0]*25
  ks[0]   = 0.0005            # gaba --> gabaOut (forward)
  ks[1]   = 5.555000000000001 # gaba + GABABR <-> gabaGABABR (forward)
  ks[2]   = 0.005             # gaba + GABABR <-> gabaGABABR (backward)
  ks[3]   = 150.0             # gabaGABABR + Gi <-> gabaGABABRGi (forward)
  ks[4]   = 0.00025           # gabaGABABR + Gi <-> gabaGABABRGi (backward)
  ks[5]   = 0.000125          # gabaGABABRGi --> gabaGABABRGibg + GiaGTP (forward)
  ks[6]   = 0.001             # gabaGABABRGibg --> gabaGABABR + Gibg (forward)
  ks[7]   = 75.0              # GABABR + Gi <-> GABABRGi (forward)
  ks[8]   = 0.000125          # GABABR + Gi <-> GABABRGi (backward)
  ks[9]   = 5.555000000000001 # gaba + GABABRGi <-> gabaGABABRGi (forward)
  ks[10]  = 0.005             # gaba + GABABRGi <-> gabaGABABRGi (backward)
  ks[11]  = 2.0               # GiaGTP + RGS <-> GiaGTPRGS (forward)
  ks[12]  = 0.002             # GiaGTP + RGS <-> GiaGTPRGS (backward)
  ks[13]  = 0.03              # GiaGTPRGS --> GiaGDP + RGS (forward)
  ks[14]  = 1250.0            # GiaGDP + Gibg --> Gi (forward)
  ks[15]  = 14.0              # GIRK + Gibg <-> GIRKGibg (forward)
  ks[16]  = 0.001             # GIRK + Gibg <-> GIRKGibg (backward)
  ks[17]  = 14.0              # GIRKGibg + Gibg <-> GIRKGibg2 (forward)
  ks[18]  = 0.001             # GIRKGibg + Gibg <-> GIRKGibg2 (backward)
  ks[19]  = 14.0              # GIRKGibg2 + Gibg <-> GIRKGibg3 (forward)
  ks[20]  = 0.001             # GIRKGibg2 + Gibg <-> GIRKGibg3 (backward)
  ks[21]  = 14.0              # GIRKGibg3 + Gibg <-> GIRKGibg4 (forward)
  ks[22]  = 0.001             # GIRKGibg3 + Gibg <-> GIRKGibg4 (backward)
  ks[23]  = 14.0              # VGCC + Gibg <-> VGCCGibg (forward)
  ks[24]  = 0.001             # VGCC + Gibg <-> VGCCGibg (backward)
    
  initvalues = [0.0, 0.0, 0.00039999999999999996, 0.0, 0.0026, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.001, 9.999999999999999e-05, 0.0, 0.0, 0.0, 0.0, 0.0]

  #run_fitc_seed5_N2000_3_obj.pdf : initvalues and ks from this fit
  initvalues[10] = initvalues[10]*params['OnlyExp0_RGS'] #Postsynaptic
  ks[0] = ks[0]*params['k[0]']
  for i in [1,9]:
    ks[i] = ks[i]*params['k[1,9]']
  for i in [2,10]:
    ks[i] = ks[i]*5.555*params['k[1,9]']*0.00011/0.005
  for i in [3,4,5,6,7,8]:
    ks[i] = ks[i]*params['k[3,4,5,6,7,8]']
  for i in [15,17,19,21,23]:
    ks[i] = ks[i]*params['k[15,17,19,21,23]']
  for i in [16,18,20,22,24]:
    ks[i] = ks[i]*params['k[16,18,20,22,24]']
  initvalues[14] = 0 # No VGCC postsynaptically

  blockedStrs = blockedStr.split(',')
  for blockedStrThis in blockedStrs:
    if blockedStrThis.find('x') > -1:
      if blockedStrThis[0:3] == 'IC_':
        xind = blockedStrThis.rfind('x')
        i_species = [i for i in range(0,len(species)) if species[i] == blockedStrThis[3:xind]]
        if len(i_species) == 1:
          initvalues[i_species[0]] = initvalues[i_species[0]]*float(blockedStrThis[xind+1:])
          print("initvalues["+str(i_species[0])+"] = initvalues["+str(i_species[0])+"]*"+str(float(blockedStrThis[xind+1:])))
        else:
          qwe
      else:
        iblocked = [i for i in range(0,len(myMechs)) if myMechs[i] in blockedStrThis]
        if len(iblocked) == 0:
          print('blockedStrThis not recognized')
        else:
          xind = blockedStrThis.rfind('x')
          mysuff = myMechs[iblocked[0]]
          print("""forall if(ismembrane(\""""+mysuff+"""\")) """+blockedStrThis[0:xind]+""" = """+blockedStrThis[0:xind]+""" * """+blockedStrThis[xind+1:])
          h("""forall if(ismembrane(\""""+mysuff+"""\")) """+blockedStrThis[0:xind]+""" = """+blockedStrThis[0:xind]+""" * """+blockedStrThis[xind+1:])
  
  conc_GIRK = initvalues[iGIRK]
  h("""forall if(ismembrane("GIRK_Yim")) {
  myarea = 0
  myvolume = 0
  for (iseg=0; iseg<nseg; iseg+=1) {
    myarea = myarea + area((0.5+iseg)/nseg)
    myvolume = myvolume + area((0.5+iseg)/nseg)*L/nseg
  }
  myarea_cm2 = myarea*1e-8
  myvolume_l = myvolume*1e-15
  N_GIRK = """+str(conc_GIRK)+"""*1e-3*myvolume_l*6.022e23 // mM -> M -> mol -> number of molecules
  gGIRKbar_tot = N_GIRK*"""+str(GIRKcond_single/1000)+""" // 0.033 nS = 33e-6 uS
  gbar_GIRK_Yim = gGIRKbar_tot/myarea_cm2*1e-6 // S/cm2
 } """)
  print("GIRKs OK")

  
  def reg2dist2(cytname):
    #print(str(cytname))
    secname = cytname[4:]
    isecs = [i for i in range(0,len(seclist)) if seclist[i].name() == secname]
    if len(isecs) > 1:
      print("Error: many with same name")
      error("Error: many with same name")
    if len(isecs) == 0:
      secname = cytname[4:]
      isecs = [i for i in range(0,len(seclist)) if seclist[i].name() == secname]
    isec = isecs[0]
    return dist2s[isec]

  def sec2dist2(mysec):
    #print(str(mysec))
    isec = [i for i in range(0,len(seclist)) if seclist[i] == mysec][0]
    return dist2s[isec]

  #gaba = rxd.Species(cyts, name='gaba', charge=0, initial=lambda nd: baclofen_max, atolscale=tolscales[0])
  #gabaOut = rxd.Species(cyts, name='gabaOut', charge=0, initial=lambda nd: baclofen_max, atolscale=tolscales[1])
  gaba = rxd.Species(cyts, name='gaba', charge=0, initial=lambda nd: 0, atolscale=tolscales[0])
  gabaOut = rxd.Species(cyts, name='gabaOut', charge=0, initial=lambda nd: 0, atolscale=tolscales[1])
  GABABR = rxd.Species(cyts, name='GABABR', charge=0, initial=initvalues[2], atolscale=tolscales[2])
  gabaGABABR = rxd.Species(cyts, name='gabaGABABR', charge=0, initial=initvalues[3], atolscale=tolscales[3])
  Gi = rxd.Species(cyts, name='Gi', charge=0, initial=initvalues[4], atolscale=tolscales[4])
  GABABRGi = rxd.Species(cyts, name='GABABRGi', charge=0, initial=initvalues[5], atolscale=tolscales[5])
  gabaGABABRGi = rxd.Species(cyts, name='gabaGABABRGi', charge=0, initial=initvalues[6], atolscale=tolscales[6])
  gabaGABABRGibg = rxd.Species(cyts, name='gabaGABABRGibg', charge=0, initial=initvalues[7], atolscale=tolscales[7])
  GiaGTP = rxd.Species(cyts, name='GiaGTP', charge=0, initial=initvalues[8], atolscale=tolscales[8])
  GiaGDP = rxd.Species(cyts, name='GiaGDP', charge=0, initial=initvalues[9], atolscale=tolscales[9])
  RGS = rxd.Species(cyts, name='RGS', charge=0, initial=initvalues[10], atolscale=tolscales[10])
  GiaGTPRGS = rxd.Species(cyts, name='GiaGTPRGS', charge=0, initial=initvalues[11], atolscale=tolscales[11])
  Gibg = rxd.Species(cyts, name='Gibg', charge=0, initial=initvalues[12], atolscale=tolscales[12])
  GIRK = rxd.Species(cyts, name='GIRK', charge=0, initial=initvalues[13], atolscale=tolscales[13])
  VGCC = rxd.Species(cyts, name='VGCC', charge=0, initial=initvalues[14], atolscale=tolscales[14])
  GIRKGibg = rxd.Species(cyts, name='GIRKGibg', charge=0, initial=initvalues[15], atolscale=tolscales[15])
  GIRKGibg2 = rxd.Species(cyts, name='GIRKGibg2', charge=0, initial=initvalues[16], atolscale=tolscales[16])
  GIRKGibg3 = rxd.Species(cyts, name='GIRKGibg3', charge=0, initial=initvalues[17], atolscale=tolscales[17])
  GIRKGibg4 = rxd.Species(cyts, name='GIRKGibg4', charge=0, initial=initvalues[18], atolscale=tolscales[18])
  VGCCGibg = rxd.Species(cyts, name='VGCCGibg', charge=0, initial=initvalues[19], atolscale=tolscales[19])
  #gaba_flux_rates = []
  #gabaOut_decay_rates = []
  #reaction_gaba_fluxes = []
  #reaction_gabaOut_decays = []
  #for icyt in range(0,len(seclist)):
  #  gaba_flux_rates.append(rxd.Parameter(cyts[icyt], initial=0))
  #  reaction_gaba_fluxes.append(rxd.Rate(gaba, gaba_flux_rates[icyt])) # gaba
  #for icyt in range(0,len(seclist)):
  #  gabaOut_decay_rates.append(rxd.Parameter(cyts[icyt], initial=0))
  #  reaction_gabaOut_decays.append(rxd.Rate(gabaOut, -gabaOut*gabaOut_decay_rates[icyt])) # gaba

  #print('working on cyt '+str(cyt))
  #reaction000 = rxd.Reaction(gaba, gabaOut, ks[0], ks[0]) #Backward rate added to model bath-application of baclofen
  reaction000 = rxd.Reaction(gaba, gabaOut, ks[0]) #Backward rate set zero. Add another reaction where this rate is ks[0] by default but 0 once gaba/baclofen is washed out
  reaction001 = rxd.Reaction(gaba + GABABR, gabaGABABR, ks[1], ks[2])
  reaction002 = rxd.Reaction(gabaGABABR + Gi, gabaGABABRGi, ks[3], ks[4])
  reaction003 = rxd.Reaction(gabaGABABRGi, gabaGABABRGibg + GiaGTP, ks[5])
  reaction004 = rxd.Reaction(gabaGABABRGibg, gabaGABABR + Gibg, ks[6])
  reaction005 = rxd.Reaction(GABABR + Gi, GABABRGi, ks[7], ks[8])
  reaction006 = rxd.Reaction(gaba + GABABRGi, gabaGABABRGi, ks[9], ks[10])
  reaction007 = rxd.Reaction(GiaGTP + RGS, GiaGTPRGS, ks[11], ks[12])
  reaction008 = rxd.Reaction(GiaGTPRGS, GiaGDP + RGS, ks[13])
  reaction009 = rxd.Reaction(GiaGDP + Gibg, Gi, ks[14])
  reaction010 = rxd.Reaction(GIRK + Gibg, GIRKGibg, ks[15], ks[16])
  reaction011 = rxd.Reaction(GIRKGibg + Gibg, GIRKGibg2, ks[17], ks[18])
  reaction012 = rxd.Reaction(GIRKGibg2 + Gibg, GIRKGibg3, ks[19], ks[20])
  reaction013 = rxd.Reaction(GIRKGibg3 + Gibg, GIRKGibg4, ks[21], ks[22])
  reaction014 = rxd.Reaction(VGCC + Gibg, VGCCGibg, ks[23], ks[24])
  #print('cyt = '+str(cyt)+' done')

  gabaOut_backward_rates = []
  rates_gabaOut_to_gaba_effect_on_gabaOut = []
  rates_gabaOut_to_gaba_effect_on_gaba = []
  for icyt in range(0,len(seclist)):
    gabaOut_backward_rates.append(rxd.Parameter(cyts[icyt], initial=ks[0]))
    rates_gabaOut_to_gaba_effect_on_gabaOut.append(rxd.Rate(gabaOut, -gabaOut*gabaOut_backward_rates[icyt])) # gaba to gabaOut, rate of subtraction from gabaOut
    rates_gabaOut_to_gaba_effect_on_gaba.append(rxd.Rate(gaba, gabaOut*gabaOut_backward_rates[icyt])) # gaba to gabaOut, rate of adding to gaba

  #Recordings:
  vec_t = h.Vector()
  vec_t.record(h._ref_t)
  vecs_all = []
  for icyt in range(0,min(2,len(seclist))):
    vecs = []
    mysec = seclist[icyt]
    for ispec in range(0,len(species)):
      vecs.append(h.Vector())
    if icyt == 0:
      print('First recording, first species...')
      timenow = time.time()
    vecs[0].record(gaba.nodes(mysec)(0.5)[0]._ref_concentration)
    if icyt == 0:
      print('First recording set in '+str(time.time()-timenow)+' seconds')
    vecs[1].record(gabaOut.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[2].record(GABABR.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[3].record(gabaGABABR.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[4].record(Gi.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[5].record(GABABRGi.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[6].record(gabaGABABRGi.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[7].record(gabaGABABRGibg.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[8].record(GiaGTP.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[9].record(GiaGDP.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[10].record(RGS.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[11].record(GiaGTPRGS.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[12].record(Gibg.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[13].record(GIRK.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[14].record(VGCC.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[15].record(GIRKGibg.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[16].record(GIRKGibg2.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[17].record(GIRKGibg3.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs[18].record(GIRKGibg4.nodes(mysec)(0.5)[0]._ref_concentration)
    #vecs[19].record(VGCCGibg.nodes(mysec)(0.5)[0]._ref_concentration)
    vecs_all.append(vecs[:])
  print('Recordings OK')
  
  spikes_thisI1 = []
  timenow = time.time()
  h("""
v_init = """+str(v0)+"""
cai0_ca_ion = """+str(ca0)+"""
st1.amp = """+str(0)+"""
st1.del = """+str(tstop-T)+"""
st1.dur = """+str(tstop)+"""
vclamp.dur1 = 10
vclamp.dur2 = 10
vclamp.dur3 = """+str(tstop)+"""
vclamp.amp1 = -60
vclamp.amp2 = -60
vclamp.amp3 = -60
vclamp_irec = new Vector()
vclamp_irec.record(&vclamp.i,1.0)
""")
  h.init()

  def introduce_gaba_concentration():
    for cyt in cyts:
      gaba[cyt].concentration = baclofen_max/1000
      gabaOut[cyt].concentration = baclofen_max/1000
    h.cvode.re_init()

  def reset_gaba_concentration():
    for icyt in range(0,len(cyts)):
      gaba[cyts[icyt]].concentration = [0 for x in gaba[cyts[icyt]].concentration]
      gabaOut[cyts[icyt]].concentration = [0 for x in gabaOut[cyts[icyt]].concentration]
      GABABR[cyts[icyt]].concentration = [GABABR[cyts[icyt]].concentration[j] + gabaGABABR[cyts[icyt]].concentration[j] for j in range(0,len(GABABR[cyts[icyt]].concentration))]
      GABABRGi[cyts[icyt]].concentration = [GABABRGi[cyts[icyt]].concentration[j] + gabaGABABRGi[cyts[icyt]].concentration[j] for j in range(0,len(GABABRGi[cyts[icyt]].concentration))]
      gabaGABABR[cyts[icyt]].concentration = [0 for x in gabaGABABR[cyts[icyt]].concentration]
      gabaGABABRGi[cyts[icyt]].concentration = [0 for x in gabaGABABRGi[cyts[icyt]].concentration]
      gabaGABABRGibg[cyts[icyt]].concentration = [0 for x in gabaGABABRGibg[cyts[icyt]].concentration]
      gabaOut_backward_rates[icyt].value = 0 #Setting this makes sure that gabaOut becomes a sink, drains the rest of gaba which might be bound to GABABRs at the time of the event
    h.cvode.re_init()

  def print_stuff():
    print('t = '+str(h.t))
    print('[gaba] = '+str(gaba[cyts[0]].concentration))
    print('[gabaOut] = '+str(gabaOut[cyts[0]].concentration))
    print('[gabaGABABR]+[gabaGABABRGi]+[gabaGABABRGibg] = '+str(gabaGABABR[cyts[0]].concentration[0]+gabaGABABRGi[cyts[0]].concentration[0]+gabaGABABRGibg[cyts[0]].concentration[0]))
    print('gabaOut_backward_rates[0] = '+str(gabaOut_backward_rates[0].value))

  h.cvode.event(1000, lambda: introduce_gaba_concentration())
  h.cvode.event(1000-0.01, lambda: print_stuff())
  h.cvode.event(1000+0.01, lambda: print_stuff())
  h.cvode.event(1000+baclofen_dur, lambda: reset_gaba_concentration())
  h.cvode.event(1000+baclofen_dur-0.01, lambda: print_stuff())
  h.cvode.event(1000+baclofen_dur+0.01, lambda: print_stuff())
  
  h.continuerun(tstop)
  print("Simulation done in "+str(time.time()-timenow)+" seconds")
  times=np.array(h.tvec)
  Vsoma=np.array(h.vsoma)
  I_Vclamp=np.array(h.vclamp_irec)
  Vdend=np.array(h.vdend)
  Casoma=np.array(h.casoma)
  Cadend=np.array(h.cadend)
  spikes = mytools.spike_times(times,Vsoma,-35,100)

  if dodraw:
    f,axs = subplots(3,2)
    axarr = axs.reshape(prod(axs.shape),).tolist()
    axarr[0].plot(times,Vsoma,lw=0.3); axarr[0].set_ylabel('V soma',fontsize=6)
    axarr[1].plot(list(range(0,len(I_Vclamp))),I_Vclamp,lw=0.3); axarr[1].set_ylabel('I VClamp',fontsize=6)
    axarr[2].plot(times,vecs_all[0][0],lw=0.3)
    axarr[2].plot(times,vecs_all[-1][0],'r--',lw=0.3); axarr[2].set_ylabel('Gaba',fontsize=6)
    axarr[3].plot(times,vecs_all[0][1],lw=0.3)
    axarr[3].plot(times,vecs_all[-1][1],'r--',lw=0.3); axarr[3].set_ylabel('GabaOut',fontsize=6)
    axarr[4].plot(times,vecs_all[0][7],lw=0.3)
    axarr[4].plot(times,vecs_all[-1][7],'r--',lw=0.3); axarr[4].set_ylabel('gabaGABABRGibg',fontsize=6)
    axarr[5].plot(times,vecs_all[0][12],lw=0.3)
    axarr[5].plot(times,vecs_all[-1][12],'r--',lw=0.3); axarr[5].set_ylabel('Gibg',fontsize=6)
    f.savefig('combe_'+params_str+('' if len(blockedStr)==0 else '_'+blockedStr)+'_'+str(GIRKcond_single)+'_tstop'+str(tstop)+'_T'+str(T)+'_dur'+str(baclofen_dur)+'_somIsc'+str(somaticIscale)+'_max'+str(baclofen_max)+'.pdf')

  picklelist = [spikes,times,I_Vclamp]
  file = open('combe_'+params_str+('' if len(blockedStr)==0 else '_'+blockedStr)+'_'+str(GIRKcond_single)+'_tstop'+str(tstop)+'_T'+str(T)+'_dur'+str(baclofen_dur)+'_somIsc'+str(somaticIscale)+'_max'+str(baclofen_max)+'.sav', 'wb')
  pickle.dump(picklelist,file)
  file.close()

  h("""forall if(ismembrane("GIRK_Yim")) print secname()," ",gbar_GIRK_Yim""")