#after drawfig_fitting_composite_late_rxdrel6_varstim_no1stgluEE_old2.py: removed the I->E part (later restored them)
#drawfig_fitting_composite_late_rxdrel6_varstim.py drawfig_fitting_composite_late_rxdrel6_varstim_no1stgluEE.py: use no1stglu for E->E synapses
#cp drawpresynGABAB_actpasaxon.py drawfig_fitting_composite.py
#cp presynAP_test_gCaNtype_axonspine_train_withrelease_testgamma_withreluctant_varygbars.py presynGABAB.py
from neuron import h,rxd
from pylab import *
import mytools
import scipy.io

Nstims = 2
decay = 1.0
minCa = 50e-6
axonPasLen = 100
axonActLen = 20

Vesmax = 0.001    #Concentration of vesicles (mM)
Fmax = 0.001 #Concentration of fusion factor (mM)
initfile = 'None'
tolstochange = 'Gi,GiaGDP,Gibg'.split(',')
tolschange = [float(x) for x in '2,2,2'.split(',')]
gaba_input_dur = 3.0
blockedStr = 'None'
blockeds = []
block_factors = []
isyntype = 1 #1 for E->E connection (OnlyExp1_RGS used), 2 for I->E connection (OnlyExp2_RGS used)

spineLen = 0.8
spineDiam = 0.4
doInterpolate = 1
doDraw = 0
fileOut = 'None'
stimAmp = 0.3
kcoeff_f = 5e-6
kcoeff_b = 1.0
kcoeff_2f = 100.0
kcoeff_2b = 1.0
kcoeff_3 = 1.0
kcoeff_4 = 1.0

gnabar = 3.9  #L5PCbiophysics4.hoc
gkbar = 0.077 #L5PCbiophysics4.hoc
gl = 0.00003  #L5PCbiophysics4.hoc
el = -90.0    #L5PCbiophysics4.hoc
gCaNtypebar_single = 2.7e-6 #uS (https://www.nature.com/articles/nn.2657) #13e-6 #uS (https://www.ncbi.nlm.nih.gov/books/NBK6181/)
minCa = 50e-6
params_all = ['107.034,514.651,167.633,11.026,66.931,3.009,1.904,2.241,1507.740','117.497,596.380,166.400,10.972,56.312,2.905,1.798,2.486,2569.532','133.092,528.522,185.581,8.574,25.868,3.186,2.125,2.729,2815.426','91.979,819.933,169.073,5.239,22.576,4.547,2.153,2.962,525.594','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','127.868,614.980,360.316,2.912,15.822,3.778,2.492,3.119,1872.786','119.681,591.014,91.960,89.687,88.948,3.453,2.300,2.755,2950.763','170.128,194.903,421.922,2.138,22.729,2.763,1.791,2.151,2061.804']

params_str = 'p0'

data_inppr_ydata_nonscaled = [0.4978165938864628, 0.6462882096069869, 0.6965065502183405, 0.7467248908296943, 0.8275109170305677, 0.8668122270742359, 0.8973799126637554,  ] #https://link.springer.com/content/pdf/10.1007/BF00169135.pdf (Olpe et al. 1994)
data_inppr_ydata_GABABblocked = [0.8078602620087335, 0.7991266375545851, 0.7991266375545851, 0.8296943231441048, 0.8362445414847161, 0.868995633187773, 0.9301310043668121 ] #

data_exppr_scaled = [[147.51475085724883, 0.7605217633066633], #Local and diffuse synaptic actions of GABA in the hippocampus. J.S.Isaacson, J.M.Solis, R.A.Nicoll, 1993
                         [198.375088814062, 0.7229294739118345],
                         [299.67563559976526, 0.664576318309598],
                         [497.7047357202435, 0.7310354947329398],
                         [693.0153532482777, 0.8064038800160637],
                         [892.7743968366749, 0.903557196255908],
                         [1295.6040900806277, 0.9651926106700441],
                         [1591.7704117883293, 1.0000293472552593]]
data_exppr_nonscaled = [[153.11004784688993, 97.40181268882175],
                        [202.10526315789474, 100],
                        [303.157894736842, 100.12084592145015],
                        [500.6698564593301, 100.96676737160121],
                        [900.287081339713, 101.99395770392749],
                        [1601.531100478469, 98.18731117824774],
                        [150.0478468899521, 76.13293051359517],
                        [199.0430622009569, 72.44712990936556],
                        [297.0334928229666, 66.58610271903322],
                        [499.13875598086116, 73.05135951661632],
                        [898.7559808612441, 90.33232628398791],
                        [1600, 100]]
cols = ['#99AA00','#6600FF']
redcols = ['#FF6600','#CC0044']

if len(sys.argv) > 1:
  gnabar = float(sys.argv[1])
if len(sys.argv) > 2:
  gkbar = float(sys.argv[2])
if len(sys.argv) > 3:
  gl = float(sys.argv[3])
if len(sys.argv) > 4:
  el = float(sys.argv[4])
if len(sys.argv) > 5:
  decay = float(sys.argv[5])
if len(sys.argv) > 6:
  minCa = float(sys.argv[6])
if len(sys.argv) > 7:
  params_str = sys.argv[7]
if len(sys.argv) > 8:
  axonPasLen = int(float(sys.argv[8]))
if len(sys.argv) > 9:
  axonActLen = int(float(sys.argv[9]))
if len(sys.argv) > 10:
  stimAmp = float(sys.argv[10])
if len(sys.argv) > 11:
  kcoeff_f = float(sys.argv[11])
if len(sys.argv) > 12:
  kcoeff_b = float(sys.argv[12])
if len(sys.argv) > 13:
  kcoeff_2f = float(sys.argv[13])
if len(sys.argv) > 14:
  kcoeff_2b = float(sys.argv[14])
if len(sys.argv) > 15:
  kcoeff_3 = float(sys.argv[15])
if len(sys.argv) > 16:
  kcoeff_4 = float(sys.argv[16])
if len(sys.argv) > 17:
  Fmax = float(sys.argv[17])
if len(sys.argv) > 18:
  Vesmax = float(sys.argv[18])
if len(sys.argv) > 19:
  doInterpolate = int(float(sys.argv[19]))

  
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])}

gaba_input_flux = params['gaba_flux']

styles = ['-','--']
widths = [0.3,0.6]

def boxoff(ax):
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.get_xaxis().tick_bottom()
    ax.get_yaxis().tick_left()

f,axarr = subplots(5,4)
axarr_right = []
for iax in range(0,4):
  for iay in range(0,5):
    axarr[iay,iax].tick_params(axis='both', which='major', labelsize=5)
    #for tick in axarr[iay,iax].xaxis.get_major_ticks() + axarr[iay,iax].yaxis.get_major_ticks():
    #  tick.label.set_fontsize(5)
    boxoff(axarr[iay,iax])

axarr_right.append(axarr[4,0].twinx())
axarr_right.append(axarr[4,2].twinx())
for iax in [0,1]:
  axarr_right[iax].spines['right'].set_color('red')
  axarr_right[iax].tick_params(colors='red', which='both')  # 'both' refers to minor and major axes
  axarr_right[iax].spines['top'].set_visible(False)
  axarr_right[iax].spines['left'].set_visible(False)

for iax in range(0,2):
  for iax2 in range(0,2):
    axarr[0,iax*2+iax2].set_position([0.1+0.48*iax+0.185*iax2,0.715,0.145,0.205])
    axarr[1,iax*2+iax2].set_position([0.1+0.48*iax+0.185*iax2,0.60,0.145,0.06])
    axarr[2,iax*2+iax2].set_position([0.1+0.48*iax+0.185*iax2,0.485,0.145,0.06])
    axarr[3,iax*2+iax2].set_position([0.1+0.48*iax+0.185*iax2,0.37,0.145,0.06])
    #axarr[3,iax*2+iax2].set_position([0.07+0.48*iax+0.22*iax2,0.785-0.20*3-0.04,0.16,0.16])
    axarr[4,iax*2+iax2].set_position([0.09+0.48*iax+0.24*iax2,0.065,0.12,0.23])
  #axarr[3,iax].set_position([0.07+0.52*iax,0.73-0.22*3,0.34,0.18])
#Ts_plot = [[350.0,1000.0],[300.0,900.0]]
Ts_plot = [[350.0],[300.0]]
for isyntype in [0,1]:
 if isyntype == 1:
  Ts = [150.0, 200.0, 300.0, 500.0, 700.0, 900.0, 1300.0, 1600.0]
  axarr[4,2*isyntype+1].plot(Ts, array(data_exppr_scaled)[:,1], 'k.',zorder=10,label='Exp.') #. data, with GABAB\nrelative to GABAB-blocked')
  #axarr[3,2*isyntype].plot(array(data_exppr_nonscaled)[0:6,0], array(data_exppr_nonscaled)[0:6,1], 'r.',color='#CC6600',zorder=10,label='Experiments, no GABABR') #. data, with GABAB\nrelative to GABAB-blocked')
  #axarr[3,2*isyntype].plot(array(data_exppr_nonscaled)[6:12,0], array(data_exppr_nonscaled)[6:12,1], 'm.',color='#CC00CC',zorder=10,label='Experiments, with GABABR') #. data, with GABAB\nrelative to GABAB-blocked')
  axarr_right[1].plot(array(data_exppr_nonscaled)[0:6,0], array(data_exppr_nonscaled)[0:6,1]/100, 'r.',lw=0.5,color=redcols[0],zorder=10,label='No GABABR (Exp.)') #. data, with GABAB\nrelative to GABAB-blocked')
  axarr_right[1].plot(array(data_exppr_nonscaled)[6:12,0], array(data_exppr_nonscaled)[6:12,1]/100, 'm.',lw=0.5,color=redcols[1],zorder=10,label='With GABABR (Exp.)') #. data, with GABAB\nrelative to GABAB-blocked')
 else:
  Ts = [350.0, 550.0, 750.0, 1000.0, 2000.0, 3000.0, 4000.0]
  axarr[4,2*isyntype+1].plot(Ts, [data_inppr_ydata_nonscaled[i]/data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'k.',zorder=10,label='Exp.') #'Exp. data, with GABAB\nrelative to GABAB-blocked')
  #axarr[3,2*isyntype].plot(Ts, [data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'r.',color='#CC6600',zorder=10,label='Experiments, no GABABR') #'Exp. data, with GABAB\nrelative to GABAB-blocked')
  #axarr[3,2*isyntype].plot(Ts, [data_inppr_ydata_nonscaled[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'m.',color='#CC00CC',zorder=10,label='Experiments, with GABABR') #'Exp. data, with GABAB\nrelative to GABAB-blocked')
  axarr_right[0].plot(Ts, [data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'r.',lw=0.5,color=redcols[0],zorder=10,label='No GABABR (Exp.)') #'Exp. data, with GABAB\nrelative to GABAB-blocked')
  axarr_right[0].plot(Ts, [data_inppr_ydata_nonscaled[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'m.',lw=0.5,color=redcols[1],zorder=10,label='With GABABR (Exp.)') #'Exp. data, with GABAB\nrelative to GABAB-blocked')
 #axarr.plot(Ts, [data_inppr_ydata_nonscaled[i]/data_inppr_ydata_GABABblocked[i] for i in range(0,len(data_inppr_ydata_nonscaled))],'r.')
 NT_amps_all = []
 NT_firstamps_all = []
 times_all_all = []
 maxPoint = 0

 from matplotlib.collections import PatchCollection
 for i in [0,1,2,3]:
  for j in [0,1]:
   for k in [0,2]:
    for startT in [0,350] if k==0 else [0]:
     polygon = Polygon(array([[startT+x for x in [0,3,3,0]],[-100,-100,60,60]]).T)
     p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
     p.set_facecolor('#BBBBBB')
     p.set_edgecolor(None)
     axarr[i,k+j].add_collection(p)


 
 for iGABAB in [0,1]:
  NTs_all = []
  times_all = []
  NT_amps = []
  NT_firstamps = []
  for iT in range(0,len(Ts)):
    addition = ('_no1stglu_recVGCC_EE' if isyntype == 1 else '_recVGCC_IE') + ('' if iGABAB else '_GABABRx0.0')
    A = scipy.io.loadmat("presynGABA_rxdrel6"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_"+str(Ts[iT])+"ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+'.mat')
    print('Loeaded '+"presynGABA_rxdrel6"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_"+str(Ts[iT])+"ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+'.mat')
    times = A['times'][0]
    NT = A['NT'][0]
    Vrec = A['Vrec'][0]
    Carec = A['Carec'][0]
    VGCCrec = A['VGCCrec'][0]
    VGCCGibgrec = A['VGCCGibgrec'][0]
    VGCCfrac = [VGCCGibgrec[i]/(VGCCrec[i]+VGCCGibgrec[i]) for i in range(0,len(VGCCrec))]
    NT_amps.append(max([NT[i] for i in range(0,len(times)) if times[i] >= 1000+Ts[iT]]))
    if isyntype == 1 :
      addition = '_EE' + ('' if iGABAB else '_GABABRx0.0')
      A2 = scipy.io.loadmat("presynGABA_rxdrel6_recVGCC"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_"+str(Ts[iT])+"ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+'.mat')
      print('Loeaded '+"presynGABA_rxdrel6_recVGCC"+addition+"_gnabar"+str(gnabar)+"_gkbar"+str(gkbar)+"_gl"+str(gl)+"_el"+str(el)+"_decay"+str(decay)+'_minCa'+str(minCa)+"_"+params_str+"_Ves"+str(Vesmax)+"_Fmax"+str(Fmax)+"_"+str(Ts[iT])+"ms_N"+str(Nstims)+"_len"+str(axonPasLen)+","+str(axonActLen)+'_amp'+str(stimAmp)+'_kcoeff'+str(kcoeff_f)+"_"+str(kcoeff_b)+"_"+str(kcoeff_2f)+"_"+str(kcoeff_2b)+"_"+str(kcoeff_3)+"_"+str(kcoeff_4)+'.mat for determining the control NT amplitude')
      times2 = A2['times'][0]
      NT2 = A2['NT'][0]
      NT_firstamps.append(max([NT2[i] for i in range(0,len(times)) if times2[i] < 1000+Ts[iT]]))
    else:
      NT_firstamps.append(max([NT[i] for i in range(0,len(times)) if times[i] < 1000+Ts[iT]])) #This doesn't work now that there's no AP for the first input. Instead, we'll have to load the corresponding file without "_no1stglu" addition and check the 1st amplitude

    NTs_all.append(NT[:])
    times_all.append(times[:])
    if Ts[iT] in Ts_plot[isyntype]:
      axarr[0,2*isyntype+iGABAB].plot(times-1000,Vrec,lw=0.5,color=cols[iGABAB],label='Control' if iGABAB else 'GABABR\nblocked')
      axarr[1,2*isyntype+iGABAB].plot(times-1000,VGCCfrac,lw=0.5,color=cols[iGABAB],label='Control' if iGABAB else 'GABABR\nblocked')
      axarr[2,2*isyntype+iGABAB].plot(times-1000,1000*Carec,lw=0.5,color=cols[iGABAB],label='Control' if iGABAB else 'GABABR\nblocked')
      axarr[3,2*isyntype+iGABAB].plot(times-1000,NT,lw=0.5,color=cols[iGABAB],label='Control' if iGABAB else 'GABABR\nblocked')
      #qwe
      for iax in range(0,4):
        axarr[iax,2*isyntype].set_xlim([920-1000,2100-1000])
        axarr[iax,2*isyntype+1].set_xlim([920-1000,2100-1000])
    it = [i for i in range(0,len(times)) if times[i] < min(Ts)][-1]
    print("baseline NT = "+str(NT[it])+", V = "+str(Vrec[it])+", Ca = "+str(Carec[it])+" at time "+str(times[it])+", firstamp="+str(NT_firstamps[-1]))

  #axarr.plot([1+i for i in range(0,len(Ts))],NT_amps)
  NT_amps_all.append(NT_amps[:])
  NT_firstamps_all.append(NT_firstamps[:])
  #axarr.set_xticks([1+i for i in range(0,len(Ts))])
  #axarr.set_xticklabels([str(x)+' ms' for x in Ts])
  maxPoint = max(maxPoint,max(NT_amps))
  #axarr.plot([1+i for i in range(0,len(Ts))], [data_inppr_ydata_nonscaled[i]/data_inppr_ydata_GABABblocked[i]*maxPoint for i in range(0,len(data_inppr_ydata_nonscaled))],'r.')
  axarr[4,2*isyntype].plot(Ts,[NT_amps_all[iGABAB][i] for i in range(0,len(NT_amps_all[0]))],'-',lw=0.5,color=cols[iGABAB],label='With GABABR' if iGABAB else 'No GABABR')
  print("iGABAB="+str(iGABAB)+", min max = "+str(min([NT_amps_all[iGABAB][i] for i in range(0,len(NT_amps_all[0]))]))+" "+str(max([NT_amps_all[iGABAB][i] for i in range(0,len(NT_amps_all[0]))])))
 axarr[4,2*isyntype+1].plot(Ts,[NT_amps_all[1][i]/NT_amps_all[0][i] for i in range(0,len(NT_amps_all[0]))],'k-',lw=0.5,label='Pred.')

 axarr[0,2*isyntype].set_ylabel('Spine $V_m$ (mV)',fontsize=6)
 axarr[1,2*isyntype].set_ylabel('Fraction\nGibg-\nbound',fontsize=6) 
 axarr[0,2*isyntype+1].set_yticklabels([])
 axarr[1,2*isyntype+1].set_yticklabels([])
 #axarr[2,2*isyntype].set_ylabel('Spine\n[Ca$^{2+}$]\n($\mu$M)',fontsize=6)
 #axarr[3,2*isyntype].set_ylabel('[neuro-\ntransmitter]\n(mM)',fontsize=6)
 axarr[2,2*isyntype].set_ylabel('Spine\n[Ca$^{2+}$]',fontsize=6)
 axarr[3,2*isyntype].set_ylabel('[neuro-\ntransmitter]',fontsize=6)
 axarr[2,2*isyntype].plot([950,950],[1,2],'k-',lw=0.8)
 axarr[2,2*isyntype].text(900,1.5,'1 $\mu$M',fontsize=5,ha='right',va='center')
 axarr[3,2*isyntype].plot([950,950],[1,2],'k-',lw=0.8)
 axarr[3,2*isyntype].text(900,1.5,'1 mM',fontsize=5,ha='right',va='center')
 axarr[4,2*isyntype].set_ylabel('peak\n[neurotransmitter]\n(mM)',fontsize=6)
 #axarr_right[isyntype].set_ylim([0,1.05-0.1*(isyntype==0)])
 #axarr[3,2*isyntype].set_ylim([0,1.05*maxPoint])
 axarr[4,2*isyntype].legend(fontsize=4.3,frameon=False,loc=5)
 axarr[4,2*isyntype+1].legend(fontsize=4.3,frameon=False,loc=4)

 axarr_right[isyntype].set_ylabel('PPR',fontsize=6,rotation=270,color='#FF0000')
 axarr_right[isyntype].legend(fontsize=4.3,frameon=False,loc=4)
 axarr[4,2*isyntype].set_ylim([0,1.2*max(max(NT_firstamps_all))])
 axarr_right[isyntype].set_ylim([0,1.2])
 
 for tick in axarr_right[isyntype].xaxis.get_major_ticks() + axarr_right[isyntype].yaxis.get_major_ticks():
   tick.label2.set_fontsize(4)
   tick.label2.set_color('#FF0000')
 f.savefig("fig2.pdf")

axarr[4,2*0+1].set_ylabel('GABABR-mediated PPR',fontsize=6)
axarr[4,2*1+1].set_ylabel('GABABR-mediated PPR',fontsize=6)
#axarr[4,2*1+1].set_ylabel('GABABR-mediated\ndecrease in neuro-\ntransmitter release',fontsize=6)

for iax in range(0,4):
  #axarr[0,iax].set_ylim([-90,25])
  #axarr[1,iax].set_ylim([0,0.00028])
  #axarr[2,iax].set_ylim([0,4.2])
  for iay in range(0,1):
    axarr[iay,iax].legend(fontsize=4.3,frameon=False)
  axarr[3,iax].set_xlabel('time (ms)',fontsize=6)
axarr[0,0].set_title('                                  I->E synapse',fontsize=8)
axarr[0,2].set_title('                                  E->E synapse',fontsize=8)

for iax in range(0,4):
  axarr[3,iax].set_yticks([0,2])
  axarr[0,iax].set_ylim([-90,50])
  axarr[1,iax].set_ylim([0,1])
  axarr[2,iax].set_ylim([0,3])
  axarr[3,iax].set_ylim([0,2.5])

  axarr[2,iax].spines['left'].set_visible(False)
  axarr[2,iax].set_yticks([])
  axarr[3,iax].spines['left'].set_visible(False)
  axarr[3,iax].set_yticks([])

  
for iax in range(0,2):
  axarr[4,2*iax].set_xlabel('ISI (ms)',fontsize=6)
  axarr[4,2*iax+1].set_xlabel('ISI (ms)',fontsize=6)
  for iay in range(0,5):
    pos = axarr[iay,2*iax].get_position()
    f.text(pos.x0 - 0.05 - 0.016*(iay>1) - 0.02*(iay==1), pos.y1 - 0.01 + 0.02*(iay < 4), chr(ord('B')+iay+6*iax), fontsize=12)
  pos = axarr[iay,2*iax+1].get_position()
  f.text(pos.x0 - 0.075, pos.y1 - 0.01, chr(ord('B')+iay+6*iax+1), fontsize=12)
  #pos = axarr[3,iax].get_position()
  #f.text(pos.x0 - 0.06, pos.y1 - 0.01, chr(ord('E')+4*iax), fontsize=12)
f.savefig("fig2.pdf")