from pylab import *
import scipy.io
import pickle
from matplotlib.collections import PatchCollection
from os.path import exists
params_all = ['95.763,528.724,127.429,86.771,46.098,5.427,3.186,4.578,1393.875','106.712,760.671,100.983,8.947,47.145,3.772,1.699,2.329,2693.683','130.765,543.857,169.167,8.660,26.611,3.194,2.055,2.717,2809.427','99.210,706.796,145.360,5.214,20.302,4.551,2.016,3.098,642.686','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','83.635,130.840,78.502,19.798,62.214,5.000,2.017,2.581,664.426','125.670,130.620,352.533,3.608,17.885,3.812,2.647,3.117,2795.527','105.438,642.913,142.107,19.602,40.482,3.599,2.391,2.974,1386.488','114.411,522.838,77.007,2.422,57.527,3.983,1.366,1.767,2359.879','118.680,111.483,149.149,23.627,49.491,3.429,2.414,2.980,2102.036','120.485,234.492,89.890,62.929,64.518,3.716,2.159,2.792,3000.000','88.810,871.501,140.692,4.467,91.855,4.311,1.620,1.908,1708.201','118.267,784.214,90.639,86.903,90.008,3.459,2.299,2.764,2836.151','169.502,193.848,423.755,2.220,22.328,2.789,1.817,2.166,2105.334']
GIRKcond = 0.033 #nS
baclofen_std = 40.0
gabaMax = 50.0
somaticIscale = 4.5
tstop = 4000
T = 3000
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()
f,axs = subplots(2,9)
axarr = axs.reshape(prod(axs.shape),).tolist()
axnew = f.add_axes([0.09,0.73,0.09,0.25])
for iax in range(0,len(axarr)):
boxoff(axarr[iax])
for axis in ['top','bottom','left','right']:
axarr[iax].spines[axis].set_linewidth(0.3)
axarr[iax].tick_params(axis='x',direction='out',width=0.4,length=2)
axarr[iax].tick_params(axis='y',direction='out',width=0.4,length=2)
axarr[iax].tick_params(axis='both', which='major', labelsize=5)
axarr[0].set_position([0.01,0.44,0.18,0.55])
axarr[0].set_xticks([])
axarr[0].set_yticks([])
boxoff(axarr[0])
for axis in ['top','bottom','left','right']:
axnew.spines[axis].set_linewidth(0.3)
axnew.set_xticks([])
axnew.set_yticks([])
for iax in range(0,7):
axarr[1+iax].set_position([0.3+0.1*iax-0.05*(iax==0),0.75,0.09-0.04*(iax==0),0.21])
for iax in range(0,7):
axarr[8+iax].set_position([0.3+0.1*iax-0.05*(iax==0),0.44,0.09,0.21])
for iax in range(0,3):
axarr[15+iax].set_position([0.25+0.15*iax,0.15,0.09,0.21])
for iax in range(0,5):
axarr[3+iax].set_yticklabels([])
#axarr[0]: plot the morphology
if True:
#Pre-saved file with the data for plotting the morphologies
unpicklefile = open('morph_distfromapicalbranchcolor_radius45.sav', 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
segdata = unpickledlist[:]
for ipos in range(0,len(segdata)):
coord1 = segdata[ipos][0]
coord2 = segdata[ipos][1]
mystyle = segdata[ipos][2]
diam = segdata[ipos][3]
mycol = segdata[ipos][4]
if mystyle == 'bx':
axarr[0].plot(coord1, coord2,mystyle,color=mycol, linewidth=0.3, ms=3.3, mew=0.5)
axnew.plot(coord1, coord2,mystyle,color=mycol, linewidth=0.3, ms=3.3, mew=0.5)
else:
axarr[0].plot(coord1, coord2,mystyle,color=mycol, linewidth=diam)
axnew.plot(coord1, coord2,mystyle,color=mycol, linewidth=4*diam)
axarr[0].axis('equal')
axnew.axis('equal')
axnew.set_ylim([570,770])
axnew.set_xlim([-80,80])
zs = [5*i for i in range(0,11)]
for iz in range(0,len(zs)):
polygon = Polygon(array([[350+100*x for x in [1,2,2,1]],[-300+65*x for x in [iz+1,iz+1,iz+2,iz+2]]]).T)
mynum = int(zs[iz]/50.0*255.99)
if mynum < 16:
mychr = '0'+hex(mynum)[2]
else:
mychr = hex(mynum)[2:]
mycol = "#"+'00'+mychr+'00'
p = PatchCollection([polygon], cmap=matplotlib.cm.jet)
p.set_facecolor(mycol); p.set_edgecolor('none')
axarr[0].add_collection(p)
axarr[0].text(380,-200+iz*65,str(int(zs[iz]))+'$\mu$M',fontsize=4,ha='right',va='center')
axarr[0].set_xlim([-167,637])
ax0bbox = [axarr[0].get_xlim(),axarr[0].get_ylim()]
bbox = [axnew.get_xlim(),axnew.get_ylim()]
axarr[0].axis('auto')
axarr[0].set_xlim(ax0bbox[0])
axarr[0].set_ylim(ax0bbox[1])
axarr[0].plot([bbox[0][0],bbox[0][1],bbox[0][1],bbox[0][0],bbox[0][0]],[bbox[1][0],bbox[1][0],bbox[1][1],bbox[1][1],bbox[1][0]],'k-',lw=0.3,color='#888888')
axnewpos = axnew.get_position()
ax0pos = axarr[0].get_position()
axarr[0].plot([bbox[0][1],axarr[0].get_xlim()[0]+(axarr[0].get_xlim()[1]-axarr[0].get_xlim()[0])*(axnewpos.x1 - ax0pos.x0)/(ax0pos.x1 - ax0pos.x0)],
[bbox[1][0],axarr[0].get_ylim()[0]+(axarr[0].get_ylim()[1]-axarr[0].get_ylim()[0])*(axnewpos.y0 - ax0pos.y0)/(ax0pos.y1 - ax0pos.y0)],'k-',lw=0.3,color='#888888')
axarr[0].plot([bbox[0][0],axarr[0].get_xlim()[0]+(axarr[0].get_xlim()[1]-axarr[0].get_xlim()[0])*(axnewpos.x0 - ax0pos.x0)/(ax0pos.x1 - ax0pos.x0)],
[bbox[1][1],axarr[0].get_ylim()[0]+(axarr[0].get_ylim()[1]-axarr[0].get_ylim()[0])*(axnewpos.y1 - ax0pos.y0)/(ax0pos.y1 - ax0pos.y0)],'k-',lw=0.3,color='#888888')
axarr[0].set_xlim(ax0bbox[0])
axarr[0].set_ylim(ax0bbox[1])
f.savefig("fig3.pdf")
#axarr[1]: plot the error versus somatic I scaling factor:
filenames = ['AP_frq_Istep_110817a1.mat', 'AP_frq_Istep_110831a1.mat', 'AP_frq_Istep_110913a1.mat', 'AP_frq_Istep_110919a2.mat', 'AP_frq_Istep_111117a1.mat', 'AP_frq_Istep_111123a4.mat', 'AP_frq_Istep_111201a2.mat', 'AP_frq_Istep_111205a2.mat', 'AP_frq_Istep_111206a3.mat', 'AP_frq_Istep_111209a1.mat', 'AP_frq_Istep_111212a1.mat', 'AP_frq_Istep_111213a1.mat', 'AP_frq_Istep_120104a1.mat', 'AP_frq_Istep_120105a1.mat', 'AP_frq_Istep_120109a1.mat', 'AP_frq_Istep_120112a3.mat', 'AP_frq_Istep_120113a2.mat', 'AP_frq_Istep_120117a3.mat', 'AP_frq_Istep_120118a2.mat', 'AP_frq_Istep_120119a3.mat', 'AP_frq_Istep_120125a2.mat', 'AP_frq_Istep_120126b1.mat', 'AP_frq_Istep_120127a1.mat', 'AP_frq_Istep_120131a4.mat']
#This is how the Schulz data were retrieved:
#APfreqs_schulz = []
#APfreqs_schulz_baclofen = []
#for ifile in range(0,len(filenames)):
# A = scipy.io.loadmat('postsyn_data_schulz/'+filenames[ifile])
# print(str(A['AP_frq'].shape)+", "+str(A['x_soma'])+", "+str(A['y_dend']))
# myAPfreq = A['AP_frq'][:,:,0].reshape(A['AP_frq'].shape[0],A['AP_frq'].shape[1])
# myAPfreq_baclofen = A['AP_frq'][:,:,1].reshape(A['AP_frq'].shape[0],A['AP_frq'].shape[1])
# if myAPfreq.shape[1] < 6:
# myAPfreq = c_[myAPfreq, nan+zeros([5,6-myAPfreq.shape[1]])]
# myAPfreq_baclofen = c_[myAPfreq_baclofen, nan+zeros([5,6-myAPfreq_baclofen.shape[1]])]
# APfreqs_schulz.append(myAPfreq[:])
# APfreqs_schulz_baclofen.append(myAPfreq_baclofen[:])
#
#meanAPfreqs_schulz = nanmean(array(APfreqs_schulz),axis=0)
#stdAPfreqs_schulz = nanstd(array(APfreqs_schulz),axis=0)
#meanAPfreqs_schulz_baclofen = nanmean(array(APfreqs_schulz_baclofen),axis=0)
#stdAPfreqs_schulz_baclofen = nanstd(array(APfreqs_schulz_baclofen),axis=0)
meanAPfreqs_schulz = array([[ 0. , 0. , 0. , 2.375 , 6. , 9.57142857],
[ 0.83333333, 2.75 , 5.58333333, 13.54166667, 20.84090909, 23.28571429],
[ 9.5 , 14.41666667, 19.39583333, 24.08333333, 29.06818182, 32.71428571],
[21.89583333, 25.125 , 28.14583333, 31.64583333, 35.11363636, 37.85714286],
[30.8125 , 32.91666667, 35.14583333, 37.25 , 39.15909091, 41.14285714]])
stdAPfreqs_schulz = array([[0. , 0. , 0. , 4.63512405, 7.85667059, 7.10848875],
[2.82351239, 4.71919838, 6.96369554, 7.78341809, 7.5819558 , 6.27271383],
[7.71902412, 7.09998044, 6.39495108, 6.60755039, 6.27457977, 5.62429134],
[7.33425065, 6.5211994 , 6.11262213, 6.10409556, 6.05642375, 5.30306034],
[7.30412626, 6.49946579, 6.24412919, 5.94593699, 5.814141 , 3.90708409]])
meanAPfreqs_schulz_baclofen = array([[ 0. , 0. , 0. , 0. , 1.68181818, 0.28571429],
[ 0. , 0.79166667, 2.79166667, 4.70833333, 10.13636364, 10.42857143],
[ 5.04166667, 7.58333333, 12.41666667, 17.29166667, 21.59090909, 21.28571429],
[16.79166667, 19.66666667, 23.04166667, 26.91666667, 29.95454545, 29.57142857],
[26.29166667, 28.70833333, 31.29166667, 34.16666667, 36.5 , 35.57142857]])
stdAPfreqs_schulz_baclofen = array([[0. , 0. , 0. , 0. , 4.00438809, 0.69985421],
[0. , 2.54916925, 5.07427472, 7.24269268, 8.01303174, 5.97272713],
[7.34551545, 8.32624699, 8.36618525, 8.38390903, 7.43247842, 3.41066754],
[9.25103222, 9.04464237, 7.87125131, 7.39322137, 7.0419346 , 3.11022015],
[8.94650382, 8.27888462, 7.73778159, 7.16278965, 6.61094134, 3.20076521]])
Kv31factors = [1.0,2.0,2.5,3.0,4.0,5.0]
amps = [0.25*i for i in range(0,6)]
AUCs_control_schulz = sum(meanAPfreqs_schulz,axis=0)*0.25
somaticIscale=1.0
Errors = []
insides_all = []
for iKv31factor in range(0,len(Kv31factors)):
Kv31factor = Kv31factors[iKv31factor]
APfreqs_test = []
Ninsides = []
for st2ind in range(0,6):
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022_p0_gSKv3_1bar_SKv3_1x'+str(Kv31factor)+'_0.0_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+'{:.1f}'.format(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std45.0.sav'
if not exists(filename):
print(filename+' does not exist')
filename = 'fIbaclofen_constsurf_p0_gSKv3_1bar_SKv3_1x'+str(Kv31factor)+'_0.0_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+'{:.1f}'.format(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std50.0.sav'
print('Loading '+filename)
if not exists(filename):
print(filename+' does not exist')
continue
unpicklefile = open(filename, 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
spikes_all = unpickledlist[1]
APfreqs_test.append([sum([1 for y in x if y >= tstop-T])/T*1000 for x in spikes_all])
insides = [[APfreqs_test[st2ind][i] >= meanAPfreqs_schulz.T[st2ind][i] - stdAPfreqs_schulz.T[st2ind][i] and
APfreqs_test[st2ind][i] <= meanAPfreqs_schulz.T[st2ind][i] + stdAPfreqs_schulz.T[st2ind][i] for i in range(0,len(APfreqs_test[st2ind]))] for st2ind in range(0,6)]
insides_all.append(insides[:])
AUCs_control = sum(array(APfreqs_test).T,axis=0)*0.25
Error_tot =sum(abs( meanAPfreqs_schulz - array(APfreqs_test).T))
Errors.append(Error_tot)
if Kv31factor == 2.5:
APfreqs_chosen = APfreqs_test[:]
insides_all_noGABA = insides_all[:]
axarr[1].plot(Kv31factors,Errors,'k:o',lw=0.3,ms=2.0,mew=1.0)
axarr[1].plot([Kv31factors[i] for i in range(0,len(Kv31factors)) if Kv31factors[i]==3.0],[Errors[i] for i in range(0,len(Kv31factors)) if Kv31factors[i]==3.0],'k.',lw=0.5,mew=1.2,ms=1.2)
axarr[1].set_ylabel('AUC error (spikes/s)',fontsize=5)
axarr[1].set_xlabel('Kv3.1 conductance factor',fontsize=5)
axarr[1].set_xticks([1,3,5])
axarr[1].set_ylim([0,300])
axarr_right = axarr[1].twinx()
pos = axarr[1].get_position()
axarr_right.set_position(pos)
axarr_right.plot(Kv31factors, [30-sum(insides_all[i]) for i in range(0,len(Kv31factors))], 'r:o',lw=0.3,ms=2.0,mew=1.0,color='#FF0000')
axarr_right.set_ylim([0,25])
axarr_right.yaxis.label.set_color('#FF0000')
axarr_right.tick_params(axis='y', colors='#FF0000')
for axis in ['top','bottom','left','right']:
axarr_right.spines[axis].set_linewidth(0.3)
axarr_right.tick_params(axis='y',direction='out',width=0.4,length=2)
axarr_right.tick_params(axis='both', which='major', labelsize=5)
axarr_right.set_ylabel('Data points further than 1 SD',color='#FF0000',fontsize=5)
#axarr[2-7]: plot the f-I curves and how much within 1SD
axarr[2].set_ylabel('spiking freq. (spikes/s)',fontsize=5)
axarr[4].set_xlabel('somatic input strength (A.U.)',fontsize=5)
insides_noGABA_all = []
N_insides_noGABA = 0
for iapic in range(0,6):
axarr[2+iapic].set_title('Apic. '+str(0.25*iapic)+' nA',fontsize=5)
polygon1 = Polygon(array([[1.0*i for i in [0,1,2,3,4,4,3,2,1,0]],[meanAPfreqs_schulz.T[iapic][i]+j*stdAPfreqs_schulz.T[iapic][i] for i,j in zip([0,1,2,3,4,4,3,2,1,0],[-1,-1,-1,-1,-1,1,1,1,1,1])]]).T)
p1 = PatchCollection([polygon1], cmap=matplotlib.cm.jet)
p1.set_facecolor('#D0D0D0')
p1.set_alpha(0.5)
axarr[2+iapic].add_collection(p1)
axarr[2+iapic].plot([1.0*i for i in range(0,5)], meanAPfreqs_schulz.T[iapic],lw=0.5,color='#666666')
axarr[2+iapic].set_ylim([0,45])
axarr[2+iapic].plot([1.0*i for i in [0,1,2,3,4]], APfreqs_chosen[iapic], 'k-', lw=0.5)
insides_noGABA = [APfreqs_chosen[iapic][isoma] >= meanAPfreqs_schulz.T[iapic][isoma]-stdAPfreqs_schulz.T[iapic][isoma] and APfreqs_chosen[iapic][isoma] <= meanAPfreqs_schulz.T[iapic][isoma]+stdAPfreqs_schulz.T[iapic][isoma] for isoma in range(0,5)]
for isoma in range(0,5):
if insides_noGABA[isoma]:
axarr[2+iapic].plot(1.0*isoma, APfreqs_chosen[iapic][isoma], 'r*',ms=4,mew=0.5)
N_insides_noGABA = N_insides_noGABA+1
else:
axarr[2+iapic].plot(1.0*isoma, APfreqs_chosen[iapic][isoma], 'kx',ms=3.3,mew=0.5)
print("isoma="+str(isoma)+",iapic="+str(iapic)+": freq = "+str(APfreqs_chosen[iapic][isoma])+" vs "+str(meanAPfreqs_schulz.T[iapic][isoma]))
insides_noGABA_all.append(insides_noGABA[:])
#axarr[8]: plot number of correctly reproduced data points for each parameter set
for iax in [8,15,16,17]:
axarr[iax].set_ylim([0,30])
axarr[iax].set_xlim([0,16])
axarr[iax].plot([0,16],[sum(insides_noGABA_all)*0.66,sum(insides_noGABA_all)*0.66],'k--',lw=0.5)
axarr[iax].set_xlabel('Parameter set',fontsize=5)
axarr[15].set_ylabel('#accepted',fontsize=5)
axarr[8].set_xticklabels([])
axarr[15].set_xticklabels([])
somaticIscale = 1.0
rel_errs_all = []
abs_N_insides_all = []
insides_all = []
insides_baclofen_all = []
APfreqs_all = []
APfreqs_baclofen_all = []
for iparam in range(0,15):
par = '_p' +str(iparam)
rel_errs = []
abs_errs = []
abs_N_insides = []
APfreqs = []
APfreqs_baclofen = []
for st2ind in range(0,6):
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022'+par+'_gSKv3_1bar_SKv3_1x2.5_0.0_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+str(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std'+str(baclofen_std)+'.sav'
if exists(filename):
print('Loading '+filename)
else:
print(filename+' does not exist')
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022'+par+'_gSKv3_1bar_SKv3_1x2.5_0.0_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+str(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std45.0.sav'
if exists(filename):
print('Using '+filename+' instead')
else:
continue
unpicklefile = open(filename, 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
spikes_all = unpickledlist[1]
APfreqs.append([sum([1 for y in x if y >= tstop-T])/T*1000 for x in spikes_all])
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022'+par+'_gSKv3_1bar_SKv3_1x2.5_0.033_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+str(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std'+str(baclofen_std)+'.sav'
if exists(filename):
print('Loading '+filename)
else:
print(filename+' does not exist')
continue
unpicklefile = open(filename, 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
spikes_all = unpickledlist[1]
APfreqs_baclofen.append([sum([1 for y in x if y >= tstop-T])/T*1000 for x in spikes_all])
if len(APfreqs) < 6 or len(APfreqs_baclofen) < 6:
print('Not all files loaded for iparam='+str(iparam))
rel_errs.append(nan)
abs_errs.append(nan)
abs_N_insides.append(nan)
continue
insides = [[APfreqs[st2ind][i] >= meanAPfreqs_schulz.T[st2ind][i] - stdAPfreqs_schulz.T[st2ind][i] and
APfreqs[st2ind][i] <= meanAPfreqs_schulz.T[st2ind][i] + stdAPfreqs_schulz.T[st2ind][i] for i in range(0,len(APfreqs[st2ind]))] for st2ind in range(0,6)]
insides_baclofen = [[APfreqs_baclofen[st2ind][i] >= meanAPfreqs_schulz_baclofen.T[st2ind][i] - stdAPfreqs_schulz_baclofen.T[st2ind][i] and
APfreqs_baclofen[st2ind][i] <= meanAPfreqs_schulz_baclofen.T[st2ind][i] + stdAPfreqs_schulz_baclofen.T[st2ind][i] for i in range(0,len(APfreqs_baclofen[st2ind]))] for st2ind in range(0,6)]
insides_all.append(insides[:])
insides_baclofen_all.append(insides_baclofen[:])
AUCs_control = sum(array(APfreqs).T,axis=0)*0.25
AUCs_baclofen = sum(array(APfreqs_baclofen).T,axis=0)*0.25
N_insides = 0
N_insides_max = 0
for i_inside in range(0,len(insides)):
for i in range(0,len(insides[i_inside])):
if insides[i_inside][i]:
N_insides_max = N_insides_max + 1
if insides_baclofen[i_inside][i]:
N_insides = N_insides + 1
myfacecol = '#999999' if N_insides < sum(insides_noGABA_all)*0.666666 else '#4444FF'
axarr[8].bar(iparam+1,N_insides,facecolor=myfacecol)
APfreqs_all.append(APfreqs[:]) #This will be the same for all, saved here for debugging purpose
APfreqs_baclofen_all.append(APfreqs_baclofen[:])
f.savefig("fig3.pdf")
APfreqs_baclofen = APfreqs_baclofen_all[0]
axarr[9].set_ylabel('spiking freq. (spikes/s)',fontsize=5)
axarr[11].set_xlabel('somatic input strength (A.U.)',fontsize=5)
for iapic in range(0,6):
axarr[9+iapic].set_title('Apic. '+str(0.25*iapic)+' nA',fontsize=5)
polygon1 = Polygon(array([[1.0*i for i in [0,1,2,3,4,4,3,2,1,0]],[meanAPfreqs_schulz_baclofen.T[iapic][i]+j*stdAPfreqs_schulz_baclofen.T[iapic][i] for i,j in zip([0,1,2,3,4,4,3,2,1,0],[-1,-1,-1,-1,-1,1,1,1,1,1])]]).T)
p1 = PatchCollection([polygon1], cmap=matplotlib.cm.jet)
p1.set_facecolor('#D0D0D0')
p1.set_alpha(0.5)
axarr[9+iapic].add_collection(p1)
axarr[9+iapic].plot([1.0*i for i in range(0,5)], meanAPfreqs_schulz_baclofen.T[iapic],lw=0.5,color='#666666')
axarr[9+iapic].set_ylim([0,45])
axarr[9+iapic].plot([1.0*i for i in [0,1,2,3,4]], APfreqs_baclofen[iapic], 'k-', lw=0.5)
for isoma in range(0,5):
if APfreqs_baclofen[iapic][isoma] >= meanAPfreqs_schulz_baclofen.T[iapic][isoma]-stdAPfreqs_schulz_baclofen.T[iapic][isoma] and APfreqs_baclofen[iapic][isoma] <= meanAPfreqs_schulz_baclofen.T[iapic][isoma]+stdAPfreqs_schulz_baclofen.T[iapic][isoma]:
if insides[iapic][isoma]:
axarr[9+iapic].plot(1.0*isoma, APfreqs_baclofen[iapic][isoma], 'r*',ms=4,mew=0.5)
else:
axarr[9+iapic].plot(1.0*isoma, APfreqs_baclofen[iapic][isoma], 'b*',ms=4,mew=0.5)
else:
axarr[9+iapic].plot(1.0*isoma, APfreqs_baclofen[iapic][isoma], 'kx',ms=3.3,mew=0.5)
f.savefig("fig3.pdf")
#axarr[15-16]: plot number of correctly reproduced data points for each parameter set with alternative SD
altSDs = [20.0, 100.0, 200.0]
for iSD in range(0,len(altSDs)):
baclofen_std = altSDs[iSD]
for iparam in range(0,15):
par = '_p' +str(iparam)
rel_errs = []
abs_errs = []
abs_N_insides = []
APfreqs = []
APfreqs_baclofen = []
insides_all = []
insides_baclofen_all = []
for st2ind in range(0,6):
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022'+par+'_gSKv3_1bar_SKv3_1x2.5_0.0_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+str(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std'+str(baclofen_std)+'.sav'
if exists(filename):
print('Loading '+filename)
else:
print(filename+' does not exist')
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022'+par+'_gSKv3_1bar_SKv3_1x2.5_0.0_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+str(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std45.0.sav'
if exists(filename):
print('Using '+filename+' instead')
else:
continue
unpicklefile = open(filename, 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
spikes_all = unpickledlist[1]
APfreqs.append([sum([1 for y in x if y >= tstop-T])/T*1000 for x in spikes_all])
filename = 'fIbaclofen_constsurfGABAGIRK8.3e-05_0.00022'+par+'_gSKv3_1bar_SKv3_1x2.5_0.033_tstop'+str(tstop)+'_T'+str(T)+'_st2ind'+str(st2ind)+'_somIsc'+str(somaticIscale)+'_amp'+str(amps[st2ind])+'_max50.0_std'+str(baclofen_std)+'.sav'
if exists(filename):
print('Loading '+filename)
else:
print(filename+' does not exist')
continue
unpicklefile = open(filename, 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
spikes_all = unpickledlist[1]
APfreqs_baclofen.append([sum([1 for y in x if y >= tstop-T])/T*1000 for x in spikes_all])
if len(APfreqs) < 6 or len(APfreqs_baclofen) < 6:
print('Not all files loaded for iparam='+str(iparam))
rel_errs.append(nan)
abs_errs.append(nan)
abs_N_insides.append(nan)
continue
insides = [[APfreqs[st2ind][i] >= meanAPfreqs_schulz.T[st2ind][i] - stdAPfreqs_schulz.T[st2ind][i] and
APfreqs[st2ind][i] <= meanAPfreqs_schulz.T[st2ind][i] + stdAPfreqs_schulz.T[st2ind][i] for i in range(0,len(APfreqs[st2ind]))] for st2ind in range(0,6)]
insides_baclofen = [[APfreqs_baclofen[st2ind][i] >= meanAPfreqs_schulz_baclofen.T[st2ind][i] - stdAPfreqs_schulz_baclofen.T[st2ind][i] and
APfreqs_baclofen[st2ind][i] <= meanAPfreqs_schulz_baclofen.T[st2ind][i] + stdAPfreqs_schulz_baclofen.T[st2ind][i] for i in range(0,len(APfreqs_baclofen[st2ind]))] for st2ind in range(0,6)]
insides_all.append(insides[:])
insides_baclofen_all.append(insides_baclofen[:])
AUCs_control = sum(array(APfreqs).T,axis=0)*0.25
AUCs_baclofen = sum(array(APfreqs_baclofen).T,axis=0)*0.25
N_insides = 0
N_insides_max = 0
for i_inside in range(0,len(insides)):
for i in range(0,len(insides[i_inside])):
if insides[i_inside][i]:
N_insides_max = N_insides_max + 1
if insides_baclofen[i_inside][i]:
N_insides = N_insides + 1
myfacecol = '#999999' if N_insides < sum(insides_noGABA_all)*0.666666 else '#4444FF'
axarr[15+iSD].bar(iparam+1,N_insides,facecolor=myfacecol)
ipanelsLabeled = [0,1,2,8,9] #,15,16]
for ipanel in range(0,len(ipanelsLabeled)):
iax = ipanelsLabeled[ipanel]
pos = axarr[iax].get_position()
f.text(pos.x0 - 0.05+0.06*(iax==0), pos.y1 - 0.01 - 0.02*(iax==0), chr(ord('A')+ipanel), fontsize=10)
for iax in range(0,3):
axarr[8+6*(iax>0)+iax].set_position([0.25,0.65-0.07*(iax+1),0.09,0.07])
axarr[17].set_visible(False)
f.savefig("fig3.pdf")