import matplotlib
matplotlib.use('Agg')
import numpy
from pylab import *
import mytools
import pickle
import time
import sys
import random
import scipy.io
import mutation_stuff
MT = mutation_stuff.getMT()
geneNames = mutation_stuff.getgenenames()
defVals = mutation_stuff.getdefvals()
unpicklefile = open('scalings.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
theseCoeffsAll = unpickledlist[0]
theseMutValsAll = unpickledlist[2]
styles = ['b-','b-','b-','b-','b-','b-','b-','b-','b-','b-']
#cols = ['#666666','#012345','#aa00aa','#bbaa00','#ee6600','#ff0000', '#00aaaa','#772277','#00cc00']
cols = ['#444444','#012345','#aa00aa','#bbaa00','#ee6600','#ff0000', '#009999','#772277','#00cc00']
col_control = '#2222ff'
coeffCoeffs = [[0.25,0],[0.125,0],[0.5,0],[0.5,1.0/3],[0.5,2.0/3],[0.5,1.0],[-0.25,0],[-0.125,0],[-0.5,0]]
lw = 1.5
fs = 10
ispDef = 0 # Consider a local maximum above -50mV a spike only if after last spike the membrane potential came below -50mV
#ispDef = 1 # Consider every local maximum above -50mV a spike
##Is = [0.65+0.025*x for x in range(0,11)]
#Is = [0.7+0.0125*x for x in range(0,17)]
Is = [0.7+0.005*x for x in range(0,41)]
unpicklefile = open('ifcurves_0_0_0.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
spTimesThisMutVal = unpickledlist[1+ispDef]
nSpikes_control = [sum([1 for x in spTimesThisMutVal[5][j] if x >= 350]) for j in range(0,len(Is))]
unpicklefile = open('steadystate_0_0_0.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
times_control = unpickledlist[1][5]
Vsoma_control = unpickledlist[2][5]
unpicklefile = open('DCshortthreshs_0_0_0.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
threshIs_control = unpickledlist[0][5]
variants = [[0,0,0],[2,3,0],[1,2,13],[3,0,1],[5,0,0]]
close("all")
f,axarr = plt.subplots(2,len(variants))
iters = [0,2,5,6,8]
timesAll = []
VsomaAll = []
spikeFreqsAll = []
threshIsAll = []
for ivar in range(0,len(variants)):
igene = variants[ivar][0]
imut = variants[ivar][1]
nVals = len(MT[igene][imut])*[0]
thesemutvars = []
theseCoeffs = theseCoeffsAll[igene][imut]
for imutvar in range(0,len(MT[igene][imut])):
thesemutvars.append(MT[igene][imut][imutvar][0])
if type(MT[igene][imut][imutvar][1]) is int or type(MT[igene][imut][imutvar][1]) is float:
MT[igene][imut][imutvar][1] = [MT[igene][imut][imutvar][1]]
nVals[imutvar] = len(MT[igene][imut][imutvar][1])
cumprodnVals = cumprod(nVals)
allmutvars = cumprodnVals[len(MT[igene][imut])-1]*[thesemutvars]
allmutvals = []
for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
allmutvals.append([0]*len(thesemutvars))
for iallmutval in range(0,cumprodnVals[len(MT[igene][imut])-1]):
for imutvar in range(0,len(MT[igene][imut])):
if imutvar==0:
allmutvals[iallmutval][imutvar] = MT[igene][imut][imutvar][1][iallmutval%nVals[imutvar]]
else:
allmutvals[iallmutval][imutvar] = MT[igene][imut][imutvar][1][(iallmutval/cumprodnVals[imutvar-1])%nVals[imutvar]]
iallmutval = variants[ivar][2]
unpicklefile = open('ifcurves_'+str(igene)+'_'+str(imut)+'_'+str(iallmutval)+'.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
ISIsThisMutVal = unpickledlist[0]
spTimesThisMutVal = unpickledlist[1+ispDef]
unpicklefile = open('steadystate_'+str(igene)+'_'+str(imut)+'_'+str(iallmutval)+'.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
timesThisVal = unpickledlist[1]
VsomaThisVal = unpickledlist[2]
unpicklefile = open('DCshortthreshs_'+str(igene)+'_'+str(imut)+'_'+str(iallmutval)+'.sav', 'r')
unpickledlist = pickle.load(unpicklefile)
unpicklefile.close()
threshIsThisVal = unpickledlist[0]
mutText = ""
for imutvar in range(0,len(MT[igene][imut])):
if imutvar > 0 and imutvar%2==0:
mutText = mutText+"\n"
mutvars = allmutvars[iallmutval][imutvar]
mutvals = allmutvals[iallmutval][imutvar]
if type(mutvars) is str:
mutvars = [mutvars]
mutText = mutText + str(mutvars) + ": "
for kmutvar in range(0,len(mutvars)):
mutvar = mutvars[kmutvar]
if (mutvar.find('off') > -1 and mutvar.find('offc') < 0) or mutvar.find('eh') > -1:
if mutvals >= 0 and kmutvar==0:
mutText = mutText + "+" + str(mutvals) +" mV"
elif kmutvar==0:
mutText = mutText + str(mutvals) +" mV"
else:
if kmutvar==0:
mutText = mutText + "*" + str(mutvals)
if kmutvar < len(mutvars)-1:
mutText = mutText + ", "
print mutText
nSpsThis = []
spikeFreqsThisVal = []
for iiter in range(0,len(iters)):
iter = iters[iiter]
if iter == 5:
continue
nSpikes = [sum([1 for x in spTimesThisMutVal[iiter][j] if x >= 350]) for j in range(0,len(Is))]
spikeFreqsThisVal.append([x/15.65 for x in nSpikes])
axarr[0,ivar].plot(Is, [x/15.65 for x in nSpikes], styles[iter],color=cols[iter])
axarr[1,ivar].plot(timesThisVal[iiter], VsomaThisVal[iiter], styles[iter],color=cols[iter])
spikeFreqs_control = [x/15.65 for x in nSpikes_control]
axarr[0,ivar].plot(Is, [x/15.65 for x in nSpikes_control], styles[iter],color=col_control)
axarr[1,ivar].plot(times_control, Vsoma_control, styles[iter],color=col_control)
axarr[0,ivar].set_ylim([0,50])
axarr[1,ivar].set_xlim([190,600])
axarr[0,ivar].set_xlabel('I (nA)')
axarr[1,ivar].set_xlabel('t (ms)')
if ivar==0:
axarr[0,ivar].set_ylabel('F (Hz)')
axarr[0,ivar].set_ylabel('V_m (mV)')
#timesAll.append(timesThisVal[:])
#VsomaAll.append(VsomaThisVal[:])
timesAll.append([timesThisVal[i] for i in [0,1,3,4]])
VsomaAll.append([VsomaThisVal[i] for i in [0,1,3,4]])
threshIsAll.append([threshIsThisVal[i] for i in [0,1,3,4]])
spikeFreqsAll.append(spikeFreqsThisVal[:])
f.suptitle(mutText)
f.savefig("fig2.eps")
scipy.io.savemat('fig2_curves.mat', {'Is': Is, 'spikeFreqs': spikeFreqsAll, 'times': timesAll, 'Vsoma': VsomaAll, 'threshIs': threshIsAll,
'spikeFreqs_control': spikeFreqs_control, 'times_control': times_control, 'Vsoma_control': Vsoma_control, 'threshI_control': threshIs_control,
'variants': variants})
picklelist = [Is, spikeFreqsAll, timesAll, VsomaAll, threshIsAll, spikeFreqs_control, times_control, Vsoma_control, threshIs_control, variants]
file = open('fig2_curves.sav', 'w')
pickle.dump(picklelist,file)
file.close()