drawfigcortfit.py

from pylab import *
import scipy.io
import scipy.stats
from os.path import exists
from matplotlib.collections import PatchCollection
import pickle

target_standard = 7 #previously 8. Check Fig. 2g from Ulanovsky et al. 2003
target_deviant = 15 #previously 20. Check Fig. 2g from Ulanovsky et al. 2003

targetAMPA = 5000  #Check that these are
targetNoise = 1750 #really the right ones

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

def mybar(ax,x,y,facecolor=[],linewidth=0.3,w=0.4):
  qs = quantile(y, [0,0.25,0.5,0.75,1])
  polygon = Polygon(array([[x-w,x+w,x+w,x-w],[qs[1],qs[1],qs[3],qs[3]]]).T)
  p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
  if type(facecolor) is not list or len(facecolor) > 0:
    p.set_facecolor(facecolor)
  p.set_edgecolor('#000000')
  p.set_linewidth(0.3)
  ax.add_collection(p)
  a2 = ax.plot([x-w,x+w,x,x,x-w,x+w,x,x,x-w,x+w],[qs[0],qs[0],qs[0],qs[2],qs[2],qs[2],qs[2],qs[4],qs[4],qs[4]],'k-',lw=linewidth)
  return [p,a2]

def myhorbarstd(ax,x,y,facecolor=[],linewidth=0.3,w=0.4):
  polygon = Polygon(array([[0,0,mean(x),mean(x)],[y-w,y+w,y+w,y-w]]).T)
  p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
  if type(facecolor) is not list or len(facecolor) > 0:
    p.set_facecolor(facecolor)
  p.set_edgecolor('#000000')
  p.set_linewidth(0.3)
  ax.add_collection(p)
  a2 = ax.plot([mean(x)-std(x),mean(x)+std(x)],[y,y],'k-',lw=linewidth)
  return [p,a2]

MMNtypes = ['Frequency deviant','Omission','Duration deviant','Inv. duration deviant']
  
f,axarr = subplots(5,1)
axarr[0].set_position([0.08, 0.84, 0.5, 0.14])
axarr[1].set_position([0.68, 0.53, 0.3, 0.44])
axarr[2].set_position([0.08, 0.6, 0.5, 0.16])
for iax in [3,4]:
  axarr[iax].set_position([0.08, 0.33-0.16*(iax==4), 0.9, 0.16])
#axarr[4].set_position([0.39, 0.21, 0.26, 0.12])
#axarr[5].set_position([0.70, 0.21, 0.26, 0.12])
#axarr[6].set_position([0.70, 0.21, 0.26, 0.12])

#Dim areas indicating the targets and nontargets
for iax in [2]:
  for itarget in [0,1,2,3,4]:
    x = 1310+500*itarget
    polygon = Polygon(array([[x,x+480,x+480,x],[-1e5,-1e5,1e5,1e5]]).T)
    p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
    p.set_facecolor('#EEEEEE' if itarget != 2 else '#DDDDFF')
    p.set_edgecolor(None)
    axarr[iax].add_collection(p)
for iax in [3]:
  for iMMN in [0,1,2,3]:
    for itarget in [0,1]:
      x = 1400*iMMN+600*itarget
      polygon = Polygon(array([[x,x+500,x+500,x],[-1e5,-1e5,1e5,1e5]]).T)
      print('iMMN='+str(iMMN)+', x = '+str(x))
      p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
      p.set_facecolor('#EEEEEE' if itarget == 1 else '#DDDDFF')
      p.set_edgecolor(None)
      axarr[iax].add_collection(p)


for iax in range(0,len(axarr)):
  axarr[iax].tick_params(axis='both', which='major', labelsize=4, direction='out', width=0.4, length=2)
  boxoff(axarr[iax])

tauNeur = 580.0
gLeak = 4.0
thresh = 0.0

gLeak_addition = '_gLeak'+str(gLeak)

amps = [100.0*i for i in range(0,11)]
nSps = []

its = [800,1300,1800,2300,2800,3300]
i_target = 3
i_nontarget = [1,2,4,5]

threshadd = '' if thresh == -40 else '_thresh'+str(thresh)
dtadd = '_dt0.25'
for iamp in range(0,len(amps)):
  stimAmp = amps[iamp]
  if exists('fIs/fI_highres_tau'+str(tauNeur)+gLeak_addition+'_amp'+str(stimAmp)+threshadd+dtadd+'.mat'):
    A = scipy.io.loadmat('fIs/fI_highres_tau'+str(tauNeur)+gLeak_addition+'_amp'+str(stimAmp)+threshadd+dtadd+'.mat')
  else:
    from brian2 import *
    simdt = 0.25
    defaultclock.dt = simdt*ms
    print('fIs/fI_highres_tau'+str(tauNeur)+gLeak_addition+'_amp'+str(stimAmp)+threshadd+dtadd+'.mat does not exist')
    start_scope()

    eqs = '''
dv/dt = (stimulus(t) + g_leak*(e_leak-v))/tau : 1 (unless refractory)
I : 1
tau : second
e_leak : 1
g_leak : 1
'''

    stimList = [0]*10+[stimAmp]*100
    stimulus = TimedArray(array(stimList), dt=100*ms)

    Population = NeuronGroup(1, eqs, threshold='v>'+str(thresh), reset='v = -80', refractory=2*ms, method='rk4')

    # Configuration for Population
    Population.v = -80
    Population.tau = [tauNeur]*ms #[10]*(2*Nperpop)*ms
    Population.g_leak = [gLeak]
    Population.e_leak = [-80]

    tstop = (10+100)*100
    print("tstop = "+str(tstop))

    V_outputPopulation = StateMonitor(Population, 'v', record=True)
    PopulationSpikeMonitor = SpikeMonitor(Population)

    timenow = time.time()
    run(tstop*ms)
    print("Simulation run in "+str(int(time.time()-timenow))+" seconds, amp="+str(stimAmp))

    scipy.io.savemat('fIs/fI_highres_tau'+str(tauNeur)+gLeak_addition+'_amp'+str(stimAmp)+threshadd+dtadd+'.mat',  {'spikes': [PopulationSpikeMonitor.t/msecond, array(PopulationSpikeMonitor.i)]})
    A = scipy.io.loadmat('fIs/fI_highres_tau'+str(tauNeur)+gLeak_addition+'_amp'+str(stimAmp)+threshadd+dtadd+'.mat')
    #axarr[1].plot(array(V_outputPopulation.t),array(V_outputPopulation.v)[0]+100*iamp,lw=0.1)
  nSps.append(len(A['spikes'][0]))

axarr[0].plot(amps,[x/10.0 for x in nSps],lw=0.6,label='IAF')
axarr[0].set_xlabel('$I$ (pA)',fontsize=6.5)
axarr[0].set_ylabel('$f$ (spikes/sec)',fontsize=6.5)



unpicklefile = open('hay/CTRL.sav','rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
spikfreqsAll = unpickledlist[0]
Is = unpickledlist[-1]

axarr[0].plot([1000*x for x in Is],spikfreqsAll[0],'k--',lw=0.6,label='Hay')
axarr[0].legend(fontsize=6.5,loc='lower right')

axinset = f.add_axes([0.14,0.91,0.07,0.05])
axinset.tick_params(axis='both', which='major', labelsize=4, direction='out', width=0.4, length=2)
boxoff(axinset)
taus = [560.0,570.0,575.0,580.0,585.0,590.0,595.0,600.0,605.0,610.0,620.0,630.0,640.0] #200 250 300 320 340 360 380 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 565 570 575 580 585 590 595 600 605 610 615 620 625 630 635 640 645 650 655 660 665 670 680 690 700 710 720 730 740 750 760 770 780 790 800 820 840 860 880 900 920 940 960 980 1000 1050 1100 1150 1200 1250 1300 1400 1500 1600 1700 1800 1900 2000
errs = []
AUCs = []
dtadd = ''
for itau in range(0,len(taus)):
  nSps = []
  for iamp in range(0,len(amps)):
    stimAmp = amps[iamp]
    A = scipy.io.loadmat('fIs/fI_highamps_tau'+str(taus[itau])+gLeak_addition+threshadd+'_amp'+str(stimAmp)+dtadd+'.mat')
    nSps.append(len(A['spikes'][0]))
  errs.append(sum([abs(nSps[i]/10.0-spikfreqsAll[0][i]) for i in range(0,len(Is))]))
  AUCs.append(sum([nSps[i]/10.0 for i in range(0,len(Is))]))
#axinset.plot(taus,errs,'k-',lw=0.4)
axinset.plot(taus,AUCs,'k.-',lw=0.4,ms=0.8,mew=0.8)
axinset.plot([min(taus),max(taus)],[sum(spikfreqsAll[0])]*2,'k--',lw=0.4)
axinset.set_xlabel(r'$C_m$ (pF)',fontsize=5.5)
axinset.set_ylabel('AUC (Hz$\\times$pA)',fontsize=5.5)
print("min err: C_m = "+str(taus[argmin(errs)]))
print("min AUC err: C_m = "+str(taus[argmin([abs(sum(spikfreqsAll[0])-AUCs[i]) for i in range(0,len(AUCs))])]))
  
f.savefig('figcortfit.pdf')

#From analyzecort_meanMMN_0.5sec.py:
Nperpop = 40

gcortAMPAs = [2000.0, 3000.0, 4000.0, 5000.0, 6000.0]
noiseCoeffs = [1000.0,1250.0,1500.0,1750.0,2000.0,2250.0,2500.0] #3000.0,4000.0]

errs_all = []
for igcortAMPA in range(0,len(gcortAMPAs)):
  for inoiseCoeff in range(0,len(noiseCoeffs)):
    if (igcortAMPA+inoiseCoeff)%2==0:
      polygon = Polygon(array([[2.5*inoiseCoeff,2.5*inoiseCoeff+2.5,2.5*inoiseCoeff+2.5,2.5*inoiseCoeff],[30*igcortAMPA,30*igcortAMPA,30*igcortAMPA+30,30*igcortAMPA+30]]).T)
      p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
      p.set_facecolor('#EEEEEE')
      p.set_edgecolor(None)
      axarr[1].add_collection(p)

curves_all = []
Nspikes_target_all = []
Nspikes_nontarget_all = []
for igcortAMPA in range(0,len(gcortAMPAs)):
  axarr[1].text(-0.66,30*igcortAMPA+12,'{:.1f}'.format(gcortAMPAs[igcortAMPA]/1000),fontsize=5.5,rotation=90)

  axarr[1].plot([0,15],[30*igcortAMPA,30*igcortAMPA],'k-',lw=0.4)
  axarr[1].plot([0,15],[30*igcortAMPA+target_standard,30*igcortAMPA+target_standard],'b--',lw=0.4)
  axarr[1].plot([0,15],[30*igcortAMPA+target_deviant,30*igcortAMPA+target_deviant],'r--',color='#CC8800',lw=0.4)
  Nspikes_target_thisampa = []
  Nspikes_nontarget_thisampa = []
  errs_thisampa = []
  for inoiseCoeff in range(0,len(noiseCoeffs)):
    if igcortAMPA == 0:
      axarr[1].text(2.5*inoiseCoeff+1.0,30*len(gcortAMPAs),'{:.2f}'.format(noiseCoeffs[inoiseCoeff]/1000),fontsize=5.5)
    Nspikes_target_thisnoise = []
    Nspikes_nontarget_thisnoise = []
    for myseed in [1,2,3,4,5]:
      curves_thisseed = []
      myseedAdd = '' if myseed == 1 else '_seed'+str(myseed)
      Nspikes_target_thisseed = []
      Nspikes_nontarget_thisseed = []
      for imodel in range(0,16):
        curves_thismodel = []
        filename = 'MMNs_2pm_sep_noISDIDD_limtau_CTRL_gcortAMPA'+str(gcortAMPAs[igcortAMPA])+'_noise'+str(noiseCoeffs[inoiseCoeff])+'_model'+str(imodel)+myseedAdd+'.mat'
        if not exists(filename):
          print(filename+' does not exist')
          continue
        print('Loading '+filename)
        A = scipy.io.loadmat(filename)
        for q in ['standard', 'deviant', 'pacemaker', 'pacemaker2', 'output', 'standardBoost', 'deviantBoost','cortOutput']:

          try:
            shp = A[q].shape
            for iy in range(0,shp[0]):
              for ix in range(0,shp[1]):
                if A[q][iy,ix].shape[0] == 1 and A[q][iy,ix].shape[1] > 1:
                  A[q][iy,ix] = A[q][iy,ix][0]
          except:
            pass

        
        Nspikes_target_forFF_thismodel = []
        Nspikes_nontarget_forFF_thismodel = []
        for iMMN in range(0,4):
          A_iMMNtype_order = [1,0,2,3]
          spikes = A['cortOutput'][A_iMMNtype_order[iMMN],0]
          spikers = A['cortOutput'][A_iMMNtype_order[iMMN],1]
          #spikes = A['cortOutput'][iMMN,0]
          #spikers = A['cortOutput'][iMMN,1]
          if noiseCoeffs[inoiseCoeff] == targetNoise and gcortAMPAs[igcortAMPA] == targetAMPA: #Choose here the best param combination. Make sure that it is the same as ibestampa and ibestnoise!
            if myseed == 1 and imodel == 0:
              if iMMN == 0:
                axarr[2].plot(spikes,spikers,'k.',ms=0.5,mew=0.5,lw=0.5)
              it = its[i_target]
              axarr[3].plot([spikes[i]-it+iMMN*1400+0 for i in range(0,len(spikes)) if spikes[i] >= it and spikes[i] < it+500],[spikers[i] for i in range(0,len(spikes)) if spikes[i] >= it and spikes[i] < it+500],'k.',
                            ms=0.5,mew=0.5,lw=0.5)
              it = its[i_nontarget[2]]
              axarr[3].plot([spikes[i]-it+iMMN*1400+600 for i in range(0,len(spikes)) if spikes[i] >= it and spikes[i] < it+500],[spikers[i] for i in range(0,len(spikes)) if spikes[i] >= it and spikes[i] < it+500],'k.',
                            ms=0.5,mew=0.5,lw=0.5)
              
            thiscurve = zeros([4000])
            mysigma = 25 #25 ms std
            for ispike in range(0,len(spikes)):
              thiscurve = thiscurve + 1/mysigma/sqrt(2*pi)*exp(-0.5*((array(range(0,4000))-spikes[ispike])/mysigma)**2)
            curves_thismodel.append(thiscurve[:])
              
          if len(spikes) == 1:
            spikes = spikes[0]

          Nspikes_nontarget_this = 0
          Nspikes_target_forFF = []
          Nspikes_nontarget_forFF = []
          for it in its:
            Nspikes_nontarget_this = Nspikes_nontarget_this + len([1 for x in spikes if x >= it and x < it+500])
            #print(str([x for x in spikes if x >= it and x < it+500]))
            if it == its[i_target]:
              Nspikes_target_forFF.append(len([1 for x in spikes if x >= it and x < it+500]))
            else:
              Nspikes_nontarget_forFF.append(len([1 for x in spikes if x >= it and x < it+500]))
          
          Nspikes_target_forFF_thismodel.append(Nspikes_target_forFF[:])
          Nspikes_nontarget_forFF_thismodel.append(Nspikes_nontarget_forFF[:])
        Nspikes_target_thisseed.append(Nspikes_target_forFF_thismodel[:])
        Nspikes_nontarget_thisseed.append(Nspikes_nontarget_forFF_thismodel[:])
        if noiseCoeffs[inoiseCoeff] == targetNoise and gcortAMPAs[igcortAMPA] == targetAMPA:
          curves_thisseed.append(curves_thismodel[:])
      Nspikes_target_thisnoise.append(Nspikes_target_thisseed[:])
      Nspikes_nontarget_thisnoise.append(Nspikes_nontarget_thisseed[:])
      if noiseCoeffs[inoiseCoeff] == targetNoise and gcortAMPAs[igcortAMPA] == targetAMPA:
        curves_all.append(curves_thisseed[:])
    Nspikes_target_thisampa.append(Nspikes_target_thisnoise[:])
    Nspikes_nontarget_thisampa.append(Nspikes_nontarget_thisnoise[:])

    #axarr[1].bar(2.5*inoiseCoeff+0.5,mean([mean([mean([mean([Nspikes_nontarget_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for i in range(0,len(Nspikes_nontarget_thisnoise[iseed][imodel][iMMN]))]) for imodel in range(0,len(Nspikes_nontarget_thisnoise[iseed]))]) for iseed in range(0,len(Nspikes_nontarget_thisnoise))]) for iMMN in range(0,4)]),bottom=30*igcortAMPA,facecolor='#4444FF')
    #axarr[1].bar(2.5*inoiseCoeff+1.5,mean([mean([mean([mean([Nspikes_target_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for i in range(0,len(Nspikes_target_thisnoise[iseed][imodel][iMMN]))]) for imodel in range(0,len(Nspikes_target_thisnoise[iseed]))]) for iseed in range(0,len(Nspikes_target_thisnoise))]) for iMMN in range(0,4)]),bottom=30*igcortAMPA,facecolor='#DD9911')
    axarr[1].bar(2.5*inoiseCoeff+0.5,mean([mean([mean([mean([Nspikes_nontarget_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for imodel in range(0,len(Nspikes_nontarget_thisnoise[iseed]))]) for i in range(0,len(Nspikes_nontarget_thisnoise[iseed][0][iMMN]))]) for iseed in range(0,len(Nspikes_nontarget_thisnoise))]) for iMMN in range(0,4)]),bottom=30*igcortAMPA,facecolor='#4444FF')
    axarr[1].bar(2.5*inoiseCoeff+1.5,mean([mean([mean([mean([Nspikes_target_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for imodel in range(0,len(Nspikes_target_thisnoise[iseed]))]) for i in range(0,len(Nspikes_target_thisnoise[iseed][0][iMMN]))]) for iseed in range(0,len(Nspikes_target_thisnoise))]) for iMMN in range(0,4)]),bottom=30*igcortAMPA,facecolor='#DD9911')

    #errs_thisampa.append(abs(mean([mean([mean([mean([Nspikes_nontarget_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for i in range(0,len(Nspikes_nontarget_thisnoise[iseed][imodel][iMMN]))]) for imodel in range(0,len(Nspikes_nontarget_thisnoise[iseed]))]) for iseed in range(0,len(Nspikes_nontarget_thisnoise))]) for iMMN in range(0,4)])-target_standard)+abs(mean([mean([mean([mean([Nspikes_target_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for i in range(0,len(Nspikes_target_thisnoise[iseed][imodel][iMMN]))]) for imodel in range(0,len(Nspikes_target_thisnoise[iseed]))]) for iseed in range(0,len(Nspikes_target_thisnoise))]) for iMMN in range(0,4)])-target_deviant))
    errs_thisampa.append(abs(mean([mean([mean([mean([Nspikes_nontarget_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for imodel in range(0,len(Nspikes_nontarget_thisnoise[iseed]))]) for i in range(0,len(Nspikes_nontarget_thisnoise[iseed][0][iMMN]))]) for iseed in range(0,len(Nspikes_nontarget_thisnoise))]) for iMMN in range(0,4)])-target_standard)+abs(mean([mean([mean([mean([Nspikes_target_thisnoise[iseed][imodel][iMMN][i]/Nperpop/0.5 for imodel in range(0,len(Nspikes_target_thisnoise[iseed]))]) for i in range(0,len(Nspikes_target_thisnoise[iseed][0][iMMN]))]) for iseed in range(0,len(Nspikes_target_thisnoise))]) for iMMN in range(0,4)])-target_deviant))

  errs_all.append(errs_thisampa[:])
  Nspikes_target_all.append(Nspikes_target_thisampa[:])
  Nspikes_nontarget_all.append(Nspikes_nontarget_thisampa[:])

  
ibestampa = argmin([min(x) for x in errs_all])
ibestnoise = argmin(errs_all[ibestampa])
#highlight_ax = axs[ibestampa,ibestnoise]
axarr[1].plot([2.5*ibestnoise,2.5*ibestnoise,2.5*ibestnoise+2,2.5*ibestnoise+2,2.5*ibestnoise],[30*ibestampa,30*ibestampa+28,30*ibestampa+28,30*ibestampa,30*ibestampa],'r-')
axarr[1].set_xticks([])
axarr[1].set_yticks([])
axarr[1].set_xlim([0,15])
axarr[1].spines['bottom'].set_visible(False)
axarr[1].spines['left'].set_visible(False)

axarr[1].text(-1.9,90,'AMPA conductance to CO population (nS)',fontsize=6.5,rotation=90,va='center')
axarr[1].text(2.5*3,30*len(gcortAMPAs)+9,'Level of noise (nA)',fontsize=6.5,ha='center')

axnew = []
for iax in range(0,4):
  axnew.append(f.add_axes([0.24+0.21*iax,0.26,0.07,0.06]))
  axnew[iax].tick_params(axis='both', which='major', labelsize=4, direction='out', width=0.4, length=2)
  boxoff(axnew[iax])
  axnew[iax].set_yticks([])
  axnew[iax].set_ylim([0,3])
  axnew[iax].text(-1.5,2.5,r'$f_{\text{deviant}}$',fontsize=6.5,ha='right',va='center')
  axnew[iax].text(-1.5,1.5,r'$f_{\text{standard}}$',fontsize=6.5,ha='right',va='center')
  axnew[iax].text(-1.5,0.5,r'$f_{\text{dd}}$',fontsize=6.5,ha='right',va='center')

#Curves:    
for iMMNtype in range(0,4):
  print('MMN='+str(iMMNtype))
  it_target = its[i_target]
  curves = array([[mean([curves_all[iseed][imodel][iMMNtype][it] for imodel in range(0,len(curves_all[0]))]) for it in range(0,len(curves_all[0][imodel][iMMNtype]))] for iseed in range(0,len(curves_all))])
  for iseed in range(0,len(curves_all)):
    target_curve = curves[iseed][it_target:it_target+500]
    nontarget_curve = zeros([500])
    for iit in i_nontarget:
      it = its[iit]
      nontarget_curve = nontarget_curve +  curves[iseed][it:it+500]
    nontarget_curve = nontarget_curve/len(i_nontarget) #normalize by the number of stimuli (in [900,1400,1900,2900,3400])
    axarr[4].plot(range(1400*iMMNtype,1400*iMMNtype+500),target_curve,'-',lw=0.1,color='#AAAAAA')
    axarr[4].plot(range(1400*iMMNtype+600,1400*iMMNtype+1100),nontarget_curve,'-',lw=0.1,color='#AAAAAA')
  mean_curve = mean(curves,axis=0)
  target_curve = mean_curve[it_target:it_target+500]
  
  nontarget_curve = zeros([500])
  for iit in i_nontarget:
    it = its[iit]
    nontarget_curve = nontarget_curve +  mean_curve[it:it+500]
  nontarget_curve = nontarget_curve/len(i_nontarget) #normalize by the number of stimuli (in [900,1400,1900,2900,3400])                                                                                                                                           
        
  axarr[4].plot(range(1400*iMMNtype,1400*iMMNtype+500),target_curve,'k-',lw=0.6)
  axarr[4].plot(range(1400*iMMNtype+600,1400*iMMNtype+1100),nontarget_curve,'k-',lw=0.6)

  #mybar(axnew[iMMNtype],0.5,[[Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_target_all[ibestampa][ibestnoise]))] for imodel in range(0,16)])
  #mybar(axnew[iMMNtype],1.5,[[Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))] for imodel in range(0,16)])
  #mybar(axnew[iMMNtype],2.5,[[Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5-Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))] for imodel in range(0,16)])
  #myhorbarstd(axnew[iMMNtype],[mean([Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_target_all[ibestampa][ibestnoise]))]) for imodel in range(0,16)],2.5,facecolor='#DDDDFF')
  #myhorbarstd(axnew[iMMNtype],[mean([Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))]) for imodel in range(0,16)],1.5,facecolor='#EEEEEE')
  #myhorbarstd(axnew[iMMNtype],[mean([Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5-Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))]) for imodel in range(0,16)],0.5,facecolor='#FFFFFF')
  #pval = scipy.stats.ranksums([mean([Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_target_all[ibestampa][ibestnoise]))]) for imodel in range(0,16)],
  #                             [mean([Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))]) for imodel in range(0,16)])[1]
  myhorbarstd(axnew[iMMNtype],[mean([Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for imodel in range(0,16)]) for seed in range(0,len(Nspikes_target_all[ibestampa][ibestnoise]))],2.5,facecolor='#DDDDFF')
  myhorbarstd(axnew[iMMNtype],[mean([Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for imodel in range(0,16)]) for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))],1.5,facecolor='#EEEEEE')
  myhorbarstd(axnew[iMMNtype],[mean([Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5-Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for imodel in range(0,16)]) for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))],0.5,facecolor='#FFFFFF')
  pval = scipy.stats.ranksums([mean([Nspikes_target_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for imodel in range(0,16)]) for seed in range(0,len(Nspikes_target_all[ibestampa][ibestnoise]))],
                               [mean([Nspikes_nontarget_all[ibestampa][ibestnoise][iseed][imodel][iMMNtype][0]/Nperpop/0.5 for imodel in range(0,16)]) for iseed in range(0,len(Nspikes_nontarget_all[ibestampa][ibestnoise]))])[1]
  print('pval = '+str(pval))
  if pval < 0.05/4:
    axnew[iMMNtype].plot([26,27.5,27.5,26],[2.5,2.5,1.5,1.5],'k-',lw=0.4)
    axnew[iMMNtype].text(28.5,1.7,'*',fontsize=6.5)


axarr[3].plot([1050,1300],[35,35],'k-',lw=0.6)
axarr[3].text(1175,37,'250 ms',fontsize=6.5,ha='center',va='bottom')

axarr[2].set_xlabel('$t$ (ms)',fontsize=6.5)
#axarr[2].set_ylabel('Neuron ID')
axarr[2].set_yticks([])
axarr[2].set_ylim([0,40])
axarr[3].set_yticks([])
axarr[3].set_xticks([])
axarr[4].set_ylabel('$f$ (spikes/sec)',fontsize=6.5)

axarr[3].set_ylim([0,40])
axarr[4].set_xticks([])
axarr[4].set_ylim([0,5.5])

for iMMN in [0,1,2,3]:
  axarr[3].text(1400*iMMN+550,4.5+40,MMNtypes[iMMN],ha='center',va='bottom',fontsize=6.5,clip_on=False)
  axnew[iax].set_xlim([0,29])
  axnew[iax].set_xticks([0,10,20])

for iax in range(0,5):
  pos = axarr[iax].get_position()
  f.text(pos.x0 - 0.05, pos.y1 - 0.01, chr(ord('A')+iax), fontsize=11)

f.savefig('figcortfit.pdf')