#! /usr/bin/python import numpy as np import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from matplotlib.mlab import griddata import glob import sys import os from joblib import Parallel, delayed import multiprocessing import pickle from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes from scipy.stats import linregress from scipy import stats import pandas as pd from sklearn.preprocessing import scale from sklearn.decomposition import PCA from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from IPython.display import display, HTML import seaborn as sns def save_obj(obj, name ): with open(name + '.pkl', 'wb') as f: pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL) def load_obj(name ): with open(name, 'rb') as f: return pickle.load(f) def pickle_the_pickle(file_string, outName): '''merges pickled data from multiple file into one large file The pickled files must be two "layered" dicts. The keys to the first layer will be lost in the merge ''' c = 1 new_dict = {} files = glob.glob(file_string) for f in files: res = pickle.load( open( f, "rb" ) ) for k in res: new_dict[c] = {} for k2 in res[k]: new_dict[c][k2] = res[k][k2] print k2, res[k][k2] c+=1 save_obj(new_dict, outName) def repickle_static(outName): ''' merges pickled data from static simulation in this case 'DRM*' and 'StatRandModInv*' ''' c = 1 new_dict = {} channels = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] I = [320, 330, 340] files = glob.glob('Beskow/save/StatRandModInv*') for f in files: res = pickle.load( open( f, "rb" ) ) for k1 in res: if isinstance(k1, int): new_dict[c] = {} factors = [int(round(i * 100)) for i in res[k1]['factors']] new_dict[c]['factors'] = factors new_dict[c]['channels'] = res['channels'] for i in I: new_dict[c][i] = res[k1][i] c+=1 print c files = glob.glob('Beskow/save/DRM*') for f in files: factors = [] # extract modulation factors-------- # extract ID from file name ID = f.split('-')[-1].split('.')[0] # loop over channels for n, channel in enumerate(channels): #extract factor from ID if n == len(channels)-1: mod_fact = int( ID.split(channels[n])[1] ) else: mod_fact = int( ID.split(channels[n])[1].split(channels[n+1])[0] ) factors.append(mod_fact) # ----- res = pickle.load( open( f, "rb" ) ) new_dict[c] = {} new_dict[c]['factors'] = factors new_dict[c]['channels'] = channels for i in I: new_dict[c][i] = res['vm'][i] #print i, res['vm'][i] c+=1 print c # repickle the combined files save_obj(new_dict, outName) def repickle_static2(inName, outName): ''' merges pickled data from static simulation ''' c = 1 new_dict = {} channels = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] I = [320, 330, 340] files = glob.glob('Beskow/save/'+inName) for f in files: res = pickle.load( open( f, "rb" ) ) for k1 in res: if isinstance(k1, int): new_dict[c] = {} factors = [int(round(i * 100)) for i in res[k1]['factors']] new_dict[c]['factors'] = factors new_dict[c]['channels'] = res['channels'] for i in I: new_dict[c][i] = res[k1][i] c+=1 print c # repickle the combined files save_obj(new_dict, outName) def adjust_spines(ax, spines, detache=False): for loc, spine in ax.spines.items(): if loc in spines: if detache: spine.set_position(('outward', 10)) # outward by 10 points spine.set_smart_bounds(True) else: spine.set_color('none') # don't draw spine # turn off ticks where there is no spine if 'left' in spines: ax.yaxis.set_ticks_position('left') else: # no yaxis ticks ax.yaxis.set_ticks([]) if 'bottom' in spines: ax.xaxis.set_ticks_position('bottom') else: # no xaxis ticks ax.xaxis.set_ticks([]) def pie_autopct(pct): #print pct return ('%1.0f' % pct) if pct > 20 else '' def load_file(f): x,y = np.loadtxt(f, unpack=True) dy = np.gradient(y) ''' dist = int( f.split('_')[2].split('u')[1] ) if dist < 115: c = 'b' elif dist < 140: c = 'r' else: c = 'k' #c = 'k' ''' #parts = f.split('_') #name = parts[2] + parts[3].split('.')[0] return [x, y, dy, f] def simple_plot(fString, color=None): ''' plots all files matching fString ''' files = glob.glob(fString) print fString print files num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) #fig,ax = plt.subplots(2,1, figsize=(8,6)) for x,y,dy,f in M: #ax[0].plot(x,y, label=f, lw=2, color=c) #ax[1].plot(x,dy, label=f, lw=2, color=c) if color: plt.plot(x,y, label=f, lw=2, color=color) else: plt.plot(x,y, label=f, lw=2, color='grey' ) print f if color: plt.plot([0,2000], [-70,-70], 'k--', lw=4) #plt.legend(loc='best') #[x1,y1] = np.loadtxt('../../results/plateau/Plotkin_plateau_D1.csv', unpack=True) #ax[0].plot((x1*1000-100), (y1*1000-3), 'g-.o', lw=3, ms=2, alpha=0.6) #fig.savefig('../../Dropbox/manuscript/Frontiers/Figures/plateau_dy.png', transparent=True) plt.show() def plot_modulation(fString, save=1): ''' plots all files matching fString ''' fig, ax = plt.subplots(3,1, figsize=(8,20) ) files = reversed(sorted(glob.glob(fString))) num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) #fig,ax = plt.subplots(2,1, figsize=(8,6)) color = ['r', 'g', 'b', 'orange', 'm' ] i = 0 for x,y,dy,f in M: path_file = os.path.split(f) channel = path_file[1].split('_')[-1].split('.')[0] lw = 2 if path_file[1] == 'vm_vm_315.out': ax[1].plot(x,y, label='control', ls='--', lw=4, color='k') ax[2].plot(x,y, label='control', ls='--', lw=4, color='k') elif path_file[1] == 'vm_modulation_315_.out': ax[1].plot(x,y, label='all mod', lw=4, color='k') ax[2].plot(x,y, label='all mod', lw=4, color='k') elif channel[0] == 'c': continue elif channel[0:2] == 'ka': ax[1].plot(x,y, label=channel, lw=2, color=color[i] ) i += 1 elif channel == 'KIR-kaf': ax[2].plot(x,y, label=channel, lw=5, color=color[i] ) i += 1 else: continue ax[2].plot(x,y, label=channel, lw=2, color=color[i] ) i += 1 ax[1].legend(loc=2) ax[2].legend(loc=2) fString = ''.join([path_file[0], '/[a,d,p]*.out']) files = glob.glob(fString) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) for i, m in enumerate(M): substrate = m[3].split('/')[-1].split('.')[0] ax[0].plot(m[0],m[1], color=color[i], lw=3, label=substrate ) ax[0].legend(loc=1) if save == 1: plt.savefig('../../../Dropbox/manuscript/Frontiers/Figures/modulation.png', transparent=True) if save == 2: plt.savefig('../../../Dropbox/manuscript/Frontiers/Figures/modulation_KIR-down_Kaf-zero_allMod.png', transparent=True) plt.show() def getSpikedata_x_y(x,y): ''' getSpikedata_x_y(x,y) -> return count Extracts and returns the number of spikes from spike trace data. # arguments x = time vector y = vm vector # returns count = number of spikes in trace (int) # extraction algorithm -threshold y and store list containing index for all points larger than 0 V -sorts out and counts the index that are the first one(s) crossing the threshold, i.e. the first index of each spike. This is done by looping over all index and check if the index is equal to the previous index + 1. If not it is the first index of a spike. If no point is above threshold in the trace the function returns 0. ''' count = 0 spikes = [] # pick out index for all points above zero potential for potential trace spikeData = [i for i,v in enumerate(y) if v > 0] # if no point above 0 if len(spikeData) == 0: return spikes else: # pick first point of each individaul transient (spike)... for j in range(0, len(spikeData)-1): if j==0: count += 1 spikes.append(x[spikeData[j]]) # ...by checking above stated criteria elif not spikeData[j] == spikeData[j-1]+1: count += 1 spikes.append(x[spikeData[j]]) return spikes def loadFile(f): x,y = np.loadtxt(f, unpack=True) amp = int(f.split('/')[-1].split('_')[2].split('.')[0]) #print f, amp return [getSpikedata_x_y(x,y), amp, [x,y] ] def plot_modulation_dynamical(): ''' plots instant firing f over time during dynamical modulation -loops over lists of wild card strings; for each string: -imports all files matching string using loadFile; for each: -extracting current injection amp and spike times using getSpikedata -saves all data in result dict -plots data from result dict ''' all_files_base = ['*_.out', \ '*_kaf.out', \ '*_kafKIR.out',\ '*_-*-.out', \ '*_kaf-*-.out', \ '*_kafKIR-*-.out' ] sim = 'modulation' if sim == 'directMod': path = 'Modulation/IncreasedNafKIR/Static/' elif sim == 'modulation': path = 'Modulation/IncreasedNafKIR/Dynamic/' all_files = [] for afb in all_files_base: # all mod files e.g. '*_.out' all_files.append( ''.join([path, 'vm_', sim, afb ]) ) # add control all_files.append( ''.join([path, 'vm_vm*_.out']) ) # global result vector res = {} vm = {} file_markers = [] counter = 1.0 for file_string in all_files: print file_string ID = file_string.split('_')[-1] files = glob.glob(file_string) #returns [getSpikedata_x_y(x,y), amp, [x,y] ] num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(loadFile)( f ) for f in files) least_spikes = 100 last_subThresh = 0 # update result vector if ID in res: ID = 'control' res[ ID ] = {} vm[ ID ] = {} for m in M: vm[ID][m[1]] = m[2] if len(m[0]) == 0: # make artificial result vector with zero freq spanning whole simulation range res[ID][m[1]] = {'isi': [0, 0], 't': [ 0, m[2][0][-1] ] } elif len(m[0]) == 1: # make artificial result vector assuming 2 Hz freq since only one spike ISI = [0.0, 0.0, counter, 0.0] counter += 1.0 freq = np.divide(1, ISI) res[ID][m[1]] = {'isi': freq, 't': [ 0, m[0][0], m[0][0], m[2][0][-1] ] } else: # extract each ISI and place in time vector (center of ISI) for i in range( 1, len(m[0]) ): if i == 1: # initiate local result vectors ISI = [0, 0] T = [0, m[0][0] ] isi = m[0][i] - m[0][i-1] t = m[0][i-1] + isi/2 ISI.append(isi) T.append(t) ISI.append(0) T.append(m[0][-1]) freq = np.divide(ISI, 1000) freq = np.divide(1, freq) res[ID][m[1]] = {'isi': freq, 't': T } if len(m[0]) < 3: print ID, m[1], m[0], res[ID][m[1]] colors = [[0,0,0], [1, 1, 0], [0.5, 0.5, 0.5], [1,0,0], [0,1,0], [0,0,1], [1,0,1] ] currents = np.arange(320,355,10) fig,ax = plt.subplots(len(currents)+1, 2, figsize=(25,18), sharex=True) legend = {} for i,ID in enumerate(res): for j, trace in enumerate(currents): print ID, i, trace if ID in legend: label = '' else: label = ID legend[ID] = 1 ax[len(currents)-j, 0].plot(vm[ID][trace][0], vm[ID][trace][1], color=colors[i], label=label ) if len(res[ID][trace]['t']) == 4: ax[len(currents)-j, 1].plot(res[ID][trace]['t'], res[ID][trace]['isi'], 'o', label=label, ms=20, color=colors[i], alpha=0.5) else: ax[len(currents)-j, 1].plot(res[ID][trace]['t'], res[ID][trace]['isi'], label=label, lw=4, color=colors[i]) ax[len(currents),0].legend(fontsize=20) # Hide the right and top spines for a in ax[:,1]: a.set_ylim([0,8]) a.spines['right'].set_visible(False) a.spines['top'].set_visible(False) a.spines['bottom'].set_visible(False) # Only show ticks on the left and bottom spines a.yaxis.set_ticks_position('left') a.xaxis.set_ticks([]) yticks = np.arange(1,8,3) a.set_yticks(yticks) a.set_yticklabels(yticks, fontsize=30) a.tick_params(width=2, length=4) ''' # set x-range ax.set_xlim([300, 500]) # set ticks yticks = np.arange(10,31,10) ax.set_yticks(yticks) ax.set_yticklabels(yticks, fontsize=30) xticks = np.arange(300,505,100) ax.set_xticks(xticks) ''' # size of frame for axis in ['left']: a.spines[axis].set_linewidth(4) xticks = [500, 1000, 2000, 3000] a.xaxis.set_ticks(xticks) a.set_xticklabels(xticks, fontsize=30) a.spines['bottom'].set_visible(True) #a.xaxis.set_ticks_position('bottom') a.spines['bottom'].set_linewidth(4) for a in ax[:,0]: a.set_ylim([-100,40]) #plot da and pka substrates = glob.glob('[!v]??[!r,!.]*ut') for sub in substrates: ax[0,1].plot( *np.loadtxt(sub,unpack=True), lw=5 ) ax[0,1].set_ylim([0,2500]) ax[0,1].set_axis_off() ax[0,0].set_axis_off() fig.savefig('../../../Dropbox/manuscript/Frontiers/Figures/modulation_FI_dynamic.png', bbox_inches='tight', transparent=True, pad_inches=0 ) #plt.figure() plt.show() def add2quant(quant, noSpikes, channels, mod_factors): ''' sorts into array ''' rheobase = 'more' if noSpikes: if noSpikes == 320: rheobase = 'equal' elif noSpikes > 320: rheobase = 'less' # get all channel index for i,chan in enumerate(channels): quant[rheobase][chan].append(mod_factors[i]) def format_boxplot(bp, median=False, lw=2, colors=None): '''adjust boxplot to stadardized format ''' if not colors: for box in bp['boxes']: box.set(color='k', linewidth=lw) for w in bp['whiskers']: w.set(color='k', linewidth=lw) for cap in bp['caps']: cap.set(color='k', linewidth=lw) else: for c,box in enumerate(bp['boxes']): box.set(color=colors[c], linewidth=lw) for c,w in enumerate(bp['whiskers']): c = c/2 w.set(color=colors[c], linewidth=lw) for c,cap in enumerate(bp['caps']): c=c/2 cap.set(color=colors[c], linewidth=lw) for c,flier in enumerate(bp["fliers"]): flier.set(marker='|', markersize=2, markerfacecolor=colors[c], markeredgecolor=colors[c]) if not median: for med in bp['medians']: med.set(color='w', linewidth=0) def sort_input( simulation, modulation, channels, res, i, j, exception=False ): flag = True not_in_range = '-' # check if all mod factors are within range for n,mod in enumerate(modulation): # get mod factor mod_fact = simulation['factors'][n] # check if factor is within range if (mod == '*') or ( np.abs(mod[0] - mod_fact) <= mod[1] ): if channels[n] in [ 'naf']: print channels[n], mod_fact else: not_in_range = channels[n] + str(mod_fact) + ' ' + str( np.abs(mod[0] - mod_fact) ) flag = False if exception and channels[n] in ['kaf', 'naf']: print channels[n], mod_fact flag = True elif channels[n] in [ 'naf']: print channels[n], mod_fact if flag: #print simulation['factors'] for n,mod in enumerate(modulation): # check if spikes and sort into dict if len(simulation['spikes']) > 0: res[j]['spikes'][channels[n]].append( int(mod_fact*100) ) res[j]['spikes'][i] = simulation['spikes'] else: res[j]['no_spikes'][channels[n]].append( int(mod_fact*100) ) res[j]['no_spikes']['runs'].append( i ) def plot_synMod_partial( sim='all-uniform_T1p', marker=None, modulation=[[0.7,0.1], \ [0.75,0.1], \ [0.8,0.05], \ [1.05,0.2], \ '*', '*', '*', \ [1.2, 0.2], \ [1.4, 0.2], \ [1.0, 0.4] ] ): ''' fraction of plot_synMod_dyn_pdc() analyse and plot time to spike from kinetic simulations using synaptic activation -How fast can DA trigger spikes? Mechanism? ''' # channels need to be in this order! ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] channels = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can', 'amp', 'nmd', 'gab' ] f = glob.glob('Beskow/save/'+sim+'.pkl') IN = load_obj(f[0]) i = 1 res = {} res[i] = { 'no_spikes':{ 'runs':[] }, 'spikes':{} } t_first_spike = [] for chan in channels: res[i]['no_spikes'][chan] = [] res[i]['spikes'][ chan] = [] # sort out simulations outside of range for j,k in enumerate(IN): sort_input( IN[k], modulation, channels, res, j, i) # how big proportion of traces spikes? N = [ len(res[i]['spikes']['kas']), len(res[i]['no_spikes']['kas']) ] plt.pie(N, colors=['r', 'grey'], \ shadow=False, \ autopct=pie_autopct ) plt.rcParams['font.size'] = 22 plt.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/highExcUniform_dynamic_PKA_proportion.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) # what's the time to first spike? for key in res[i]['spikes']: if key not in channels: t_first_spike.append( res[i]['spikes'][key][0] - 1000) # plot t_first_spike plt.figure() boxprops = dict(linewidth=2, color='grey') medianprops = dict(linestyle='-', linewidth=2, color='k') # over modulation setups bp1 = plt.boxplot(t_first_spike, \ vert=False, \ medianprops=medianprops, \ boxprops=boxprops) format_boxplot(bp1, median=True, lw=2, colors=['#bf5b17', '#bf5b17', '#bf5b17', '#bf5b17', '#bf5b17']) print 'substrates: ', bp1['medians'][0].get_ydata() plt.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/highExcUniform_dynamic_PKA_dtSpike_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) plt.xlim([0, 1000]) plt.show() def plot_synMod_dyn_pdc( sim='all_synMod', marker=None, modulation=[[0.7,0.1], \ [0.75,0.1], \ [0.8,0.05], \ [1.05,0.2], \ '*', '*', '*', \ [1.2, 0.2], \ [1.4, 0.2], \ [1.0, 0.4] ] ): ''' Produces figure 6 for the frontiers analyse and plot result from kinetic simulations using synaptic activation -How fast can DA trigger spikes? Mechanism? ''' # channels need to be in this order! ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] channels = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can', 'amp', 'nmd', 'gab' ] f2, a2 = plt.subplots(2,1, figsize=(4.0, 6), gridspec_kw = {'height_ratios':[2, 2]} ) f3, a3 = plt.subplots(1,1, figsize=(4.0, 3) ) res = {} t_first_spike = [] simulations = ['T1p-synModAll', \ 'cAMP_synModAll', \ 'Gbgolf_synModAll', \ 'D1RDAGolf_synModAll', \ 'T1p-glutModIonChan', \ 'T1p-synModGABA', \ 'T1p-glutMod', \ 'T1p-ionModOnly', \ 'T1p-Naf100_synModAll', \ 'T1p-Kaf100_synModAll', \ 'T1p-kafNaf100_synModAll', \ 'T1p-longer_synModAll' ] names = ['PKA', \ 'cAMP', \ r'G$_{\beta\gamma}$', \ 'D1R-G', \ 'No GABA', \ 'Synaptic', \ 'Excitatory', \ 'Intrinsic', \ 'no Naf', \ 'no Kaf', \ 'no Kaf,Naf', \ 'PKA 1.5 s' ] if not marker: f1, a1 = plt.subplots(4,3, figsize=(9.6,10) ) #fKaf, aKaf = plt.subplots(1,1, figsize=(3.2,2.4)) for i, sim in enumerate(simulations): f = glob.glob('Beskow/save/all_'+sim+'.pkl') print f IN = load_obj(f[0]) res[i] = { 'no_spikes':{ 'runs':[] }, 'spikes':{} } for chan in channels: res[i]['no_spikes'][chan] = [] res[i]['spikes'][ chan] = [] # sort out simulations outside of range for j,k in enumerate(IN): if i <= 7: sort_input( IN[k], modulation, channels, res, j, i ) else: sort_input( IN[k], modulation, channels, res, j, i, exception=True ) # how big proportion of traces spikes? N = [ len(res[i]['spikes']['kas']), len(res[i]['no_spikes']['kas']) ] if i <= 3: a1[i,2].pie(N, colors=['r', 'grey'], \ shadow=False, \ autopct=pie_autopct ) a1[i,2].set_title(names[i]) elif i <= 7: a1[i-4,0].pie(N, colors=['r', 'grey'], \ shadow=False, \ autopct=pie_autopct ) a1[i-4,0].set_title(names[i]) else: a1[i-8,1].pie(N, colors=['r', 'grey'], \ shadow=False, \ autopct=pie_autopct ) a1[i-8,1].set_title(names[i]) plt.rcParams['font.size'] = 22 #plt.axis('equal') # what's the time to first spike? t_first_spike.append( [] ) for key in res[i]['spikes']: if key not in channels: t_first_spike[i].append( res[i]['spikes'][key][0] - 1000) save_obj(t_first_spike, 'Beskow/save/fig_dynamic_t_first_spike') f1.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/dynamic_PKA_proportion_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) ''' fKaf.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/dynamic_PKA-noNafMod_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' else: t_first_spike = load_obj('Beskow/save/fig_dynamic_t_first_spike.pkl') colors = ['#bf5b17', '#4daf4a', '#984ea3', '#ff7f00'] # plot noise and example trace files = glob.glob('Beskow/save/spiking_*') num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) for x,y,dy,f in M: x = np.subtract(x,1000) a2[0].plot(x,y, label=f, lw=2, color='k') files = glob.glob('Beskow/save/test_*') num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) for x,y,dy,f in M: x = np.subtract(x,1000) a2[0].plot(x,y, label=f, lw=2, color=[0.5,0.5,0.5]) a2[0].plot([-250,1000], [-70,-70], 'k--', lw=3) # plot substrates files = ['Target1p', 'cAMP', 'Gbgolf', 'D1RDAGolf', ] for c,f in enumerate(files): x, y = np.loadtxt('Beskow/save/substrate_'+f+'.out', unpack=True) y = np.divide(y, max(y)/50) #y = np.subtract(y, 10) x = np.subtract(x, 1000) a2[0].plot(x,y, label=f, lw=3, color=colors[c]) # plot t_first_spike boxprops = dict(linewidth=2, color='grey') medianprops = dict(linestyle='-', linewidth=2, color='k') # over modulation setups bp1 = a3.boxplot([t_first_spike[i] for i in [0,4,5,6,7]], \ labels=[ names[i] for i in [0,4,5,6,7]], \ vert=False, \ medianprops=medianprops, \ boxprops=boxprops) format_boxplot(bp1, median=True, lw=2, colors=['#bf5b17', '#bf5b17', '#bf5b17', '#bf5b17', '#bf5b17']) print 'substrates: ', bp1['medians'][0].get_ydata() # over substrates (PKA, cAMP, Gbg, DIR)--reversed bp2 = a2[1].boxplot([t_first_spike[i] for i in [3,2,1,0]], \ labels=[ names[i] for i in [3,2,1,0]], \ vert=False, \ medianprops=medianprops, \ boxprops=boxprops) format_boxplot(bp2, median=True, lw=2, colors=colors[::-1]) for median in range(len(bp2['medians'])): print print bp2['medians'][median].get_xdata() for ax in a2[0:1]: ax.set_xlim([-250, 1000]) a2[1].set_xticks( [-250,0,500,1000] ) #a2[1].yaxis.tick_right() a2[1].spines['left'].set_visible(False) a2[1].spines['top'].set_visible(False) a2[1].set_yticks([]) #a2[1].yaxis.set_ticks_position('right') a2[1].xaxis.set_ticks_position('bottom') # size of frame for axis in ['right', 'bottom']: a2[1].spines[axis].set_linewidth(2) a3. spines[axis].set_linewidth(2) a2[0].set_ylim([-78, 55]) a2[0].axis('off') a3.set_xlim([750, 1000]) a3.set_xticks( [800,900,1000] ) a3.set_yticks([]) a3.spines['left'].set_visible(False) a3.spines['top'].set_visible(False) a3.set_yticks([]) a3.xaxis.set_ticks_position('bottom') # Show no spines #a2[0].tick_params(axis='both', bottom='off', top='off', left='off', right='off', labelleft='off') f2.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/dynamic_substrates_Example_dtSpike_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) f3.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/dynamic_PKA_dtSpike_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) plt.show() def plot_static_modulation_pdc(pickled_data, modulation=['*', '*', '*', '*', '*', '*', '*'], marker='test'): ''' remake of fig 3 for the Frontiers manuscript with additional "inverse-pyramidal distributed modulation factors merged with the old pyramidal distributed factors -> create uniform distribution??? ARGUMENTS: pickled_data: merged data set, e.g. 'Beskow/save/STATIC_grandTotal' modulation: restrict dataset to specific range of factors? marker: mark resulting figures, showing restrictions etc ''' # check excitability --------------------------------- # result structure quant = {'more':{}, 'equal':{}, 'less':{}} counter = 1 I = np.arange(320, 345, 10) # load pickle DATA = load_obj(''.join([pickled_data, '.pkl']) ) for n in DATA: # flags that checks that all modulation factors are within range and keeps track of rheobase current mod_within_range = True noSpikes = None mod_factors = DATA[n]['factors'] channels = DATA[n]['channels'] # if first iteration, setup of result structure if not 'can' in quant['more']: for channel in channels: for key in quant: quant[key][channel] = [] # check if all mod factors are within range ---- for m,mod in enumerate(modulation): mod_fact = mod_factors[m] if not mod == '*': # check if factor is outside of range and if so update flag if np.abs(mod[0] - mod_fact) > int(mod[1]): mod_within_range = False if mod_within_range: for amp in I: if len(DATA[n][amp]['spikes']) == 0: noSpikes = amp add2quant(quant, noSpikes, channels, mod_factors) # plot distribution of factors --------------------------------- f5, a5 = plt.subplots(len(channels), 1, figsize=(6,6), sharex=True) f6, a6 = plt.subplots(len(channels), 1, figsize=(6,6), sharex=True) bins = np.linspace(0, 200, 50) color = ['r', 'grey', 'b'] TOT = {} for i, chan in enumerate(channels): TOT[chan] = quant['less'][chan] + quant['equal'][chan] + quant['more'][chan] a6[i].hist( TOT[chan], bins=10, color='k' ) for c, cathegory in enumerate(['less', 'equal', 'more']): weights = np.ones_like(quant[cathegory][chan])/float(len(quant[cathegory][chan]))*100 histtype='step' a5[i].hist( quant[cathegory][chan], bins=10, color=color[2-c], alpha=0.2, histtype='stepfilled') a5[i].hist( quant[cathegory][chan], bins=10, color=color[2-c], histtype='step', lw=3) for a in [a5, a6]: a[i].set_yticks([]) #a[i].set_ylim([0,8]) a[i].set_xticks([0,100,200]) a[i].set_xticklabels([0,100,200], fontsize=24) a[i].spines['left'].set_visible(False) a[i].spines['top'].set_visible(False) a[i].spines['right'].set_visible(False) a[i].spines['bottom'].set_visible(True) a[i].xaxis.set_ticks_position('bottom') a[i].spines['bottom'].set_linewidth(2) ''' f5.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/extendedDATA_uniform_static_distributions_pdc_', marker, '.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' # procetage in each category plt.figure(figsize=(4,3)) N = [ len(quant['more']['naf']), len(quant['equal']['naf']), len(quant['less']['naf']) ] plt.pie(N, colors=color, autopct=pie_autopct, shadow=False) plt.rcParams['font.size'] = 40 plt.axis('equal') ''' plt.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/extendedDATA_uniform_static_pie_pdc_', marker, '.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' plt.show() def plot_modulation_dynamical_pdc(modulation=[[70,10], [75,10], [100,3], [102,22], '*', '*', '*'], sim='DRM', marker='test' ): ''' fig 2/3 frontiers - run from ipython using command: import plot_modulation as sp sp.plot_modulation_dynamical_pdc(modulation=[[70,10], [75,10], [100,3], [102,22], '*', '*', '*'], mark ...: er='all', sim='DRM') sp.plot_modulation_dynamical_pdc(modulation=[[70,10], [75,10], [80,3], [102,22], '*', '*', '*'], mark ...: er='kaf', sim='DRM') sp.plot_modulation_dynamical_pdc(modulation=[[70,10], [75,10], [80,3], [102,22], '*', '*', '*'], mark ...: er='all', sim='ADRM') to create figure for A, B and C respectively. The other figure can be created by manipulating the range of the modulation parameters. ADRM specifies the set without axon modulation DRM specifies the set with uniform modulation (including axon). ''' # channels need to be in this order! ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] channels = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] all_files = glob.glob( ''.join(['Beskow/save/',sim,'*']) ) res = {} quant = {'more':{}, 'equal':{}, 'less':{}} raster = {} mod_dist = {} for channel in channels: mod_dist[channel] = [] for key in quant: quant[key][channel] = [] counter = 1 I = np.arange(300, 355, 10) for file_string in all_files: mod_factors = [] noSpikes = None # flag that checks that all modulation factors are within range mod_within_range = True # extract ID from file name ID = file_string.split('-')[-1].split('.')[0] # check if all mod factors are within range for n,mod in enumerate(modulation): #extract factor from ID if n == len(channels)-1: mod_fact = int( ID.split(channels[n])[1] ) else: mod_fact = int( ID.split(channels[n])[1].split(channels[n+1])[0] ) mod_dist[channels[n]].append(mod_fact) mod_factors.append(mod_fact) if not mod == '*': # check if factor is outside of range and if so update flag if np.abs(mod[0] - mod_fact) > int(mod[1]): mod_within_range = False if mod_within_range: # load pickle RES = load_obj(file_string) res[ ID ] = {} raster[ID] = {} print ID for amp in I: spikes = RES['vm'][amp]['spikes'] raster[ID][amp] = RES['vm'][amp]['spikes'] if len(spikes) == 0: # make artificial result vector with zero freq spanning whole simulation range res[ID][amp] = {'isi': [0, 0], 't': [ 0, 2000 ] } noSpikes = amp elif len(spikes) == 1: # make artificial result vector assuming 2 Hz freq since only one spike ISI = [0.0, 0.0, counter, 0.0] counter += 1.0 freq = np.divide(1, ISI) res[ID][amp] = {'isi': freq, 't': [ 0, spikes[0], spikes[0], 2000 ] } else: # extract each ISI and place in time vector (center of ISI) for i in range( 1, len(spikes) ): if i == 1: # initiate local result vectors ISI = [0, 0] T = [0, spikes[0] ] isi = spikes[i] - spikes[i-1] t = spikes[i-1] + isi/2 ISI.append(isi) T.append(t) ISI.append(0) T.append(spikes[-1]) freq = np.divide(ISI, 1000) freq = np.divide(1, freq) res[ID][amp] = {'isi': freq, 't': T } add2quant(quant, noSpikes, channels, mod_factors) colors = [[0,0,0], [1, 1, 0], [0.5, 0.5, 0.5], [1,0,0], [0,1,0], [0,0,1], [1,0,1] ] currents = np.arange(300,355,10) '''fig,ax = plt.subplots(len(currents)+1, 1, figsize=(6,9), sharex=True) legend = {} for i,ID in enumerate(res): for j, trace in enumerate(currents): if ID in legend: label = '' else: label = ID legend[ID] = 1 #ax[len(currents)-j, 0].plot(vm[ID][trace][0], vm[ID][trace][1], label=label ) if len(res[ID][trace]['t']) == 4: ax[len(currents)-j].plot(res[ID][trace]['t'], res[ID][trace]['isi'], 'o', label=label, ms=20, alpha=0.5) else: ax[len(currents)-j].plot(res[ID][trace]['t'], res[ID][trace]['isi'], label=label, lw=4, ) #ax[len(currents),0].legend(fontsize=20) # Hide the right and top spines for a in ax: a.set_ylim([0,10]) a.spines['right'].set_visible(False) a.spines['top'].set_visible(False) a.spines['bottom'].set_visible(False) # Only show ticks on the left and bottom spines a.yaxis.set_ticks_position('left') a.xaxis.set_ticks([]) yticks = np.arange(1,8,3) a.set_yticks(yticks) a.set_yticklabels(yticks, fontsize=20) a.tick_params(width=2, length=4) # size of frame for axis in ['left']: a.spines[axis].set_linewidth(4) xticks = [500, 1000, 2000, 3000] a.xaxis.set_ticks(xticks) a.set_xticklabels(xticks, fontsize=30) a.spines['bottom'].set_visible(True) #a.xaxis.set_ticks_position('bottom') a.spines['bottom'].set_linewidth(4) #plot da and pka substrates = glob.glob('[!v]??[!r,!.]*ut') for sub in substrates: ax[0].plot( *np.loadtxt(sub,unpack=True), lw=5 ) ax[0].set_ylim([0,2500]) ax[0].set_axis_off() ax[0].set_axis_off() #ax[0].set_title('Instant spike Freq over time', fontsize=50) ax[-1].set_xlabel("Time (ms)", fontsize=20) ax[3].set_ylabel("Instant spike f (Hz)", fontsize=20) fig.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_instantF_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' #plt.figure() print 'fig 1 done' '''f2, a2 = plt.subplots(len(channels), 1, figsize=(6,6), sharex=True) for i,channel in enumerate(channels): adjust_spines(a2[i], ['bottom', 'left']) n, bins, patches = a2[i].hist(mod_dist[channel], 50, facecolor=[0.2,0.2,0.2], alpha=0.75) #a2[i].set_ylabel(channel, fontsize=20) a2[i].set_yticks([]) #a2[i].set_yticklabels([0,20], fontsize=18) #print n, channel #a2[0].set_title('Total sample distribution', fontsize=50) a2[-1].set_xlabel('Modulation (%)', fontsize=35) a2[-1].set_xticks(np.arange(0,210,50)) a2[-1].set_xticklabels(np.arange(0,210,50), fontsize=24) f2.savefig('../../../Dropbox/manuscript/Frontiers/Figures/static_chanDistAll_pdc.png', bbox_inches='tight', transparent=True, pad_inches=0 ) # only one distribution (example) plt.figure() adjust_spines(plt.gca(), ['bottom', 'left']) n, bins, patches = plt.hist(mod_dist['kas'], 50, facecolor=[0.2,0.2,0.2], alpha=0.75) plt.yticks([]) plt.xlabel('Modulation (%)', fontsize=35) plt.xticks(np.arange(0,210,50), fontsize=28) plt.yticks([0,1000,2000,3000], [0,1,2,3], fontsize=28) print n, np.sum(n) plt.savefig('../../../Dropbox/manuscript/Frontiers/Figures/static_chanDistKas_pdc.png', bbox_inches='tight', transparent=True, pad_inches=0 )''' print 'fig 2 done' f3, a3 = plt.subplots(1, 1, figsize=(5,4)) #axins = plt.axes([.60, .13, .3, .25]) medians = [[],[],[]] exc_list = ['more', 'equal', 'less'] y = [] for c,channel in enumerate(channels): y.append([]) for i,mod in enumerate(exc_list): medians[i].append(np.median(quant[mod][channel])) y[c] = y[c] +list(quant[mod][channel]) adjust_spines(a3, ['bottom', 'left']) bp = a3.boxplot(y, labels=channels, vert=False) #bpInl = axins.boxplot(y, labels=channels, vert=False) format_boxplot(bp) #format_boxplot(bpInl) color = ['r', 'grey', 'b'] for c,median in enumerate(medians): for m,chan in enumerate(median): y = [m+1-0.4, m+1+0.4] x = [chan, chan] a3.plot(x,y,color=color[c], lw=3) #axins.plot(x,y,color=color[c], lw=3) a3.set_xticks(np.arange(0,210,50)) a3.set_xticklabels(np.arange(0,210,50), fontsize=20) labels = a3.get_yticklabels() a3.set_yticklabels(labels, fontsize=20) a3.set_xlabel('Percentage of control', fontsize=20) #a3.set_title('Distribution of modulation factors', fontsize=40) #axins.set_ylim([0,4]) #axins.set_xlim([55,85]) #plt.setp(axins, xticks=[], yticks=[]) ''' f3.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_chanDist_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' print 'fig 3 done' f5, a5 = plt.subplots(len(channels), 1, figsize=(4,10), sharex=True) bins = np.linspace(0, 200, 30) for i, chan in enumerate(reversed(channels)): for c, cathegory in enumerate(['less', 'equal', 'more']): weights = np.ones_like(quant[cathegory][chan])/float(len(quant[cathegory][chan]))*100 a5[i].hist( quant[cathegory][chan], bins, weights=weights, color=color[2-c], alpha=0.6, label=cathegory ) a5[i].set_yticks([]) a5[i].set_ylim([0,8]) a5[i].set_xticks([0,100,200]) a5[i].set_xticklabels([0,100,200], fontsize=24) a5[i].spines['left'].set_linewidth(2) a5[i].spines['top'].set_visible(False) a5[i].spines['right'].set_visible(False) a5[i].spines['bottom'].set_visible(True) a5[i].xaxis.set_ticks_position('bottom') a5[i].spines['bottom'].set_linewidth(2) ''' f5.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_distributions_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' print 'fig 4 done' # plot correlation of kaf and naf channels # for the "more" excitable group #plt.figure( figsize=(4,3) ) fig = plt.figure( figsize=(4,3) ) ax = fig.add_subplot(1, 1, 1) ax.plot( quant['more']['naf'], quant['more']['kaf'], '.', color='k', alpha=0.2 ) ax.spines['left'].set_position('center') ax.spines['right'].set_color('none') ax.spines['bottom'].set_position('center') ax.spines['top'].set_color('none') ax.spines['left'].set_smart_bounds(True) ax.spines['bottom'].set_smart_bounds(True) ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') # plot naf,kaf pairs # calc and plot linear regression #fit = np.polyfit( quant['more']['naf'], quant['more']['kaf'], 1 ) #fit_fn = np.poly1d(fit) try: minv = min(quant['more']['naf']) maxv = max(quant['more']['naf']) x = [minv, maxv] #plt.plot(x, fit_fn(x), '-r') # recalc using scipy statistics = linregress(quant['more']['naf'], quant['more']['kaf']) print statistics y1 = statistics[0]*minv+statistics[1] y2 = statistics[0]*maxv+statistics[1] ax.plot( x, [y1,y2], '-r', lw=2 ) ax.set_xticks([0,200]) ax.set_xticklabels([0,200], fontsize=18) ax.set_yticks([0,200]) ax.set_yticklabels([0,200], fontsize=18) ''' fig.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_correlation_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' except: print 'could not save correlation figure--line ~1475' f6, a6 = plt.subplots(len(channels), len(channels), figsize=(10,10), sharex=True) for i1 in range(len(channels)): for i2 in range(len(channels)): a6[i1,i2].axis('off') bins = np.linspace(0, 200, 50) f7 = plt.figure( figsize=(4,3) ) ax = f7.add_subplot(1, 1, 1) ax.spines['left'].set_position('center') ax.spines['right'].set_color('none') ax.spines['bottom'].set_position('center') ax.spines['top'].set_color('none') ax.spines['left'].set_smart_bounds(True) ax.spines['bottom'].set_smart_bounds(True) ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') ax.set_xticks([0,200]) ax.set_xticklabels([0,200], fontsize=20) ax.set_yticks([0,200]) ax.set_yticklabels([0,200], fontsize=20) R = 0 try: for i1, chan1 in enumerate(channels): for i2 in range(i1,len(channels)): if i1==i2: for c, cathegory in enumerate(['less', 'equal', 'more']): weights = np.ones_like(quant[cathegory][chan1])/float(len(quant[cathegory][chan1]))*100 a6[i1,i2].hist( quant[cathegory][chan1], bins, weights=weights, color=color[2-c], edgecolor = "none", alpha=0.6 ) a6[i1,i2].set_ylim([0,8]) else: chan2 = channels[i2] # plot channel pairs using scipy--more excitable a6[i2,i1].plot( quant['more'][chan1], quant['more'][chan2], '.', color='k', alpha=0.02 ) minv = min(quant['more'][chan1]) maxv = max(quant['more'][chan1]) x = [minv, maxv] statistics = linregress(quant['more'][chan1], quant['more'][chan2]) R_square = statistics[2]**2 print 'more', chan1, chan2, R_square if R_square*1000 >= 1: if chan2 == 'kir': a6[i2,i1].text(100, 150, str.format('{0:.2f}', R_square), fontsize=12) else: a6[i2,i1].text(0, 150, str.format('{0:.2f}', R_square), fontsize=12) y1 = statistics[0]*minv+statistics[1] y2 = statistics[0]*maxv+statistics[1] a6[i2,i1].plot([0,200], [100,100], 'k') a6[i2,i1].plot([100,100], [0,200], 'k') a6[i2,i1].plot( x, [y1,y2], '-r', lw=3 ) # plot channel pairs using scipy--equal excitable a6[i1,i2].plot( quant['equal'][chan1], quant['equal'][chan2], '.', color='k', alpha=0.02 ) minv = min(quant['equal'][chan1]) maxv = max(quant['equal'][chan1]) x = [minv, maxv] statistics = linregress(quant['equal'][chan1], quant['equal'][chan2]) R_square = statistics[2]**2 print 'equal', chan1, chan2, R_square if R_square*1000 >= 1: if chan2 == 'kir': a6[i1,i2].text(100, 150, str.format('{0:.2f}', R_square), fontsize=12) else: a6[i1,i2].text(0, 150, str.format('{0:.2f}', R_square), fontsize=12) y1 = statistics[0]*minv+statistics[1] y2 = statistics[0]*maxv+statistics[1] a6[i1,i2].plot([0,200], [100,100], 'k') a6[i1,i2].plot([100,100], [0,200], 'k') a6[i1,i2].plot( x, [y1,y2], '-r', lw=3 ) # less excitable minv = min(quant['less'][chan1]) maxv = max(quant['less'][chan1]) x = [minv, maxv] statistics = linregress(quant['less'][chan1], quant['less'][chan2]) R_square = statistics[2]**2 y1 = statistics[0]*minv+statistics[1] y2 = statistics[0]*maxv+statistics[1] ax.plot([0,200], [100,100], 'k') ax.plot([100,100], [0,200], 'k') ax.plot( x, [y1,y2], '-r', lw=3 ) print 'less', chan1, chan2, R_square R = R + R_square ax.text(100, 150, str.format('{0:.2f}', R), fontsize=18) ''' f6.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_distAndScatter_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) f7.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_less_corr_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' print 'fig 4 done' except: print 'could not save correlation figure--line ~1599' ''' f4, a4 = plt.subplots(len(currents)-3, 1, figsize=(4,3)) for counter,ID in enumerate(raster): for c,I in enumerate(currents[2:-1]): i = len(currents) - (c+4) # create y vector y = np.ones( len(raster[ID][I]) ) + counter if len(raster[ID][I]) > 0: a4[i].plot(raster[ID][I],y, 'k|', mew=2, ms=3, alpha=1.0) ylimits = [-0.5, len(raster)+0.5] a4[i].set_ylim(ylimits) a4[i].set_yticks([]) a4[i].set_ylabel(str(I), fontsize=18) #a4[i].spines['right'].set_visible(False) #a4[i].spines['top'].set_visible(False) #a4[i].spines['bottom'].set_visible(False) a4[i].spines['left'].set_linewidth(2) a4[i].set_xlim([0,1000]) if not i == len(currents) - 1: a4[i].set_xticks([]) plt.tick_params( axis='x', # changes apply to the x-axis which='both', # both major and minor ticks are affected bottom='off', # ticks along the bottom edge are off top='off', # ticks along the top edge are off labelbottom='off') # labels along the bottom edge are off xticks = [0,500,1000] a4[-1].xaxis.set_ticks(xticks) a4[-1].set_xticklabels(xticks, fontsize=18) a4[-1].spines['bottom'].set_visible(True) a4[-1].xaxis.set_ticks_position('bottom') a4[-1].spines['bottom'].set_linewidth(2) a4[0].spines['top'].set_linewidth(1) #plt.gca().xaxis.set_minor_locator(plt.NullLocator()) #a4[0].set_title('Spike raster', fontsize=50) a4[-1].set_xlabel("Time (ms)", fontsize=20) files = glob.glob('Modulation/IncreasedNafKIR/Static/vm_vm_3[0-5]*') #returns [getSpikedata_x_y(x,y), amp, [x,y] ] num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(loadFile)( f ) for f in files) Max = len(raster) for c,I in enumerate(currents): if len(M[c][0]) > 0 and int(M[c][1]) < 350: indx = list(currents).index(M[c][1]) i = len(currents) - (indx+2) print indx, M[c][1], i, c y = np.ones( len(M[c][0]) ) * Max * 0.5 a4[i].plot(M[c][0], y, '|', color=[0.2,0.2,0.2], ms=Max*3, mew=4, alpha=0.5) f4.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_raster_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 ) ''' plt.figure(figsize=(4,3)) N = [ len(quant['more']['naf']), len(quant['equal']['naf']), len(quant['less']['naf']) ] plt.pie(N, colors=color, autopct=pie_autopct, shadow=False) plt.rcParams['font.size'] = 50 plt.axis('equal') ''' plt.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_pie_pdc.png']), bbox_inches='tight', transparent=True, pad_inches=0 )''' #plt.close('all') plt.show() def plot_vm(): ''' Reproduces the base plots (validation: figure 1 B and C) for the frontiers paper. ''' files = glob.glob('Results/FI/vm_vm_*') f1,a1 = plt.subplots(1,1) fig,ax = plt.subplots(1,1) res = {} num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(loadFile)( f ) for f in files) least_spikes = 100 last_subThresh = 0 inj = np.arange(-100,345,40) # sort traces in spiking vs non spiking for i,m in enumerate(M): if len(m[0]) > 0: res[m[1]] = (len(m[0]) -1) * ( 1000 / (m[0][-1]-m[0][0]) ) if len(m[0]) < least_spikes: least_spikes = len(m[0]) rheob_index = i else: if m[1] in inj: a1.plot(m[2][0], m[2][1], 'grey') if m[1] > last_subThresh: last_subThresh = m[1] # add last trace without spike res[last_subThresh] = 0 # plot first trace to spike a1.plot(M[rheob_index][2][0], M[rheob_index][2][1], 'k', lw=2) a1.plot([10,110], [-20, -20], 'k') a1.plot([10, 10], [-20, -10], 'k') a1.axis('off') #plt.savefig('../../../Dropbox/manuscript/Frontiers/Figures/example_traces.png', transparent=True) # FI curve --------------------------------------------------------------------------------------- for i,f in enumerate(glob.glob('Exp_data/FI/Planert2013-D1-FI-trace*') ): [x_i,y_i] = np.loadtxt(f, unpack=True) if i == 0: label = 'D1planert' else: label = '' ax.plot(x_i, y_i, 'brown', lw=5, alpha=1, label=label, clip_on=False) #plt.legend() AMP = [] inj = np.arange(last_subThresh,480,40) for I in inj: AMP.append(res[I]) ax.plot(inj, AMP, color='k', lw=7, label='neuron', clip_on=False) # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') # set x-range ax.set_xlim([150, 800]) # set ticks yticks = np.arange(10,31,10) ax.set_yticks(yticks) ax.set_yticklabels(yticks, fontsize=30) xticks = np.arange(150,760,150) ax.set_xticks(xticks) ax.set_xticklabels(xticks, fontsize=30) ax.tick_params(width=3, length=6) # size of frame for axis in ['bottom','left']: ax.spines[axis].set_linewidth(4) #fig.savefig('../../../Dropbox/manuscript/Frontiers/Figures/FI-planert.png', transparent=True) #plt.show() def plot_IF_modulation(): ''' function for comparing IF traces of different modulations ''' fig,ax = plt.subplots(1,1, figsize=(4,5) ) #cmap = plt.cm.get_cmap('Set1', 8) ''' all_files = ['vm_directMod*_.out', \ 'vm_directMod*_kaf.out', \ 'vm_directMod*_kafKIR.out',\ 'vm_directMod*_-*-.out', \ 'vm_directMod*_kaf-*-.out', \ 'vm_directMod*_kafKIR-*-.out', \ 'vm_vm*_.out' ]''' all_files_base = ['*_.out', \ '*_kaf.out', \ '*_kafKIR.out',\ '*_-*-.out', \ '*_kaf-*-.out', \ '*_kafKIR-*-.out' ] sim = 'directMod' if sim == 'directMod': path = 'Modulation/IncreasedNafKIR/Static/' elif sim == 'modulation': path = 'Modulation/IncreasedNafKIR/Dynamic/' all_files = [] for afb in all_files_base: # all mod files e.g. '*_.out' all_files.append( ''.join([path, 'vm_', sim, afb ]) ) # add control all_files.append( ''.join([path, 'vm_vm*_.out']) ) for file_string in all_files: print file_string files = glob.glob(file_string) res = {} num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(loadFile)( f ) for f in files) least_spikes = 100 last_subThresh = 0 for i,m in enumerate(M): if len(m[0]) == 1: res[m[1]] = 2 least_spikes = len(m[0]) rheob_index = i elif len(m[0]) > 1: # calc mean instant spike freq as (#spikes-1) / dt [1/s] # i.e. total number of ISIs / total time in seconds from first to last spike res[m[1]] = (len(m[0]) -1) * ( 1000 / (m[0][-1]-m[0][0]) ) if len(m[0]) < least_spikes: # update least_spikes and rheobase index least_spikes = len(m[0]) rheob_index = i else: # if no spikes in train... if m[1] > last_subThresh: # ...and current amp larger than for current last_subthreshold, update! print file_string, last_subThresh, m[1] last_subThresh = m[1] # add last trace without spike to results res[last_subThresh] = 0 AMP = [] INJ = [] inj = np.arange(last_subThresh,500,10) for I in inj: if I in res: AMP.append(res[I]) INJ.append(I) label = file_string.split('_')[-1] #print label, last_subThresh ax.plot(INJ, AMP, lw=4, label=label, clip_on=False) #ax.legend(fontsize=26) # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') # set x-range ax.set_xlim([300, 360]) # set ticks yticks = np.arange(3,10,3) ax.set_yticks(yticks) ax.set_yticklabels(yticks, fontsize=30) xticks = np.arange(320, 380, 40) ax.set_xticks(xticks) ax.set_xticklabels(xticks, fontsize=30) ax.tick_params(width=3, length=6) # size of frame for axis in ['bottom','left']: ax.spines[axis].set_linewidth(4) fig.savefig('../../../Dropbox/manuscript/Frontiers/Figures/modulation_FI.png', bbox_inches='tight', transparent=True, pad_inches=0 ) plt.show() def get_max(f): x,y = np.loadtxt(f, unpack=True) m = max(y[3300:-1]) path_file = os.path.split(f) dist = int(path_file[1].split('_')[1]) return [m, dist, x[3300:-1], y[3300:-1]] def plot_Ca(fString): ''' plots resulting Ca as distance dependent using all traces matching fString ''' fig, ax = plt.subplots(1,1, figsize=(6,8)) files = glob.glob(fString) N = len(files) gradient = [(1-1.0*x/(N-1), 0, 1.0*x/(N-1)) for x in range(N)] res = {} distances = np.arange(0,200, 10) num_cores = multiprocessing.cpu_count() #int(np.ceil(multiprocessing.cpu_count() / 2)) M = Parallel(n_jobs=num_cores)(delayed(get_max)( f ) for f in files) for m in M: for d in distances: if m[1] > d-5 and m[1] < d+5: if d not in res: res[d] = [] res[d].append(m[0]) break mean_amp = [] spread = [] x = [] for d in distances: if d in res: mean_amp.append( np.mean(res[d]) ) spread.append( np.std( res[d]) ) x.append(d) for m in M: if m[1] >= 40: ax.plot(m[1], np.divide(m[0], mean_amp[3]), '.', ms=20, color='k', alpha=0.2) #print x #print mean_amp mean_amp = np.divide(mean_amp, mean_amp[3]) #day et al 2008 [x1,y1] = np.loadtxt('Exp_data/bAP/bAP-DayEtAl2006-D1.csv', unpack=True) ax.plot(x1, y1, 'brown', lw=6) #plt.fill_between(x, np.subtract(mean_amp, spread), np.add(mean_amp, spread), color='r', alpha=0.5) ax.plot(x[3:], mean_amp[3:], lw=6, color='k') # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') # set ticks yticks = [0,1] ax.set_yticks(yticks) ax.set_yticklabels(yticks, fontsize=35) xticks = [40, 120, 200] ax.set_xticks(xticks) ax.set_xticklabels(xticks, fontsize=35) ax.tick_params(width=2, length=4) # size of frame for axis in ['bottom','left']: ax.spines[axis].set_linewidth(4) #plt.savefig('../../../Dropbox/manuscript/Frontiers/Figures/ca_bAP.png', transparent=True) f,a = plt.subplots(1,1, figsize=(2,4)) path = os.path.dirname(files[0]) path_file = ''.join([path, '/vm_ca_2000.out']) f = glob.glob( path_file ) lw=4 a.plot(*np.loadtxt(f[0], unpack=True), color='k', lw=lw ) # construct I curve base = -110 a.plot([0,100], [base, base], 'k', lw=lw) a.plot([100,100],[base, -90], 'k', lw=lw) a.plot([100,102], [-90, -90], 'k', lw=lw) a.plot([102,102],[-90, base], 'k', lw=lw) a.plot([102,155], [base, base], 'k', lw=lw) a.set_ylim([base-1, 50]) a.set_xlim([95, 120]) a.axis('off') #plt.savefig('../../../Dropbox/manuscript/Frontiers/Figures/ca_bAP_vm.png', transparent=True) #plt.show() #plt.plot(*np.loadtxt(f, unpack=True), color=gradient[i]) #plt.ylim([0.00003, 0.0001]) #plt.show() def plot_Ca_updated(fString, ax, i, ax2): ''' plots resulting Ca as distance dependent using all traces matching fString Reproduces Fig 5B and C for the frontiers paper, if combined with ipython command below. import matplotlib.pyplot as plt import plot_modulation as func fig, ax = plt.subplots(1,1, figsize=(6,8)) fig2, ax2 = plt.subplots(1,1, figsize=(6,8)) sim = ['Org', 'DA', 'Half', 'Block'] color = ['black', 'm', 'b', 'g'] for i,s in enumerate(sim): fs = 'Beskow/Ca/'+s+'/*' x,v,ax2 = func.plot_Ca_updated(fs, ax, i, ax2) if i == 0: lw = 6 ls = '-' else: lw = 6 ls = '--' ax.plot(x[3:], v[3:], lw=lw, ls=ls, color=color[i]) ax.savefig(' ../../../Dropbox/Manuscripts/Frontiers/Figures/Fig5/ca_bAP_updated.png', transparent=True) ax2.savefig('../../../Dropbox/Manuscripts/Frontiers/Figures/Fig5/ca_transient.png', transparent=True) ''' files = glob.glob(fString) N = len(files) gradient = [(1-1.0*x/(N-1), 0, 1.0*x/(N-1)) for x in range(N)] distances = np.arange(0,200, 10) color = ['k', 'm', 'b', 'g'] res = {} num_cores = multiprocessing.cpu_count() M = Parallel(n_jobs=num_cores)(delayed(get_max)( f ) for f in files) for m in M: for d in distances: if m[1] > d-5 and m[1] < d+5: if d not in res: res[d] = [] res[d].append(m[0]) break lw = 3 a = 0.3 if m[1] >= 170 and m[1] <= 200: # distally m2 = np.subtract(m[2],100) m3 = np.multiply(m[3],500) # transformation from mM to uM and div by 2 (1000/2) # to be average of two same size pools instead of summation, # (the traces are the sum of the two...). ax2[0,1].plot(m2,m3, lw=lw, color=color[i], alpha=a) if i <= 1: if i == 0: a = 0.7 ax2[1,1].plot(m2,m3, lw=lw, color=color[i], alpha=a) elif m[1] >= 30 and m[1] <= 50: # proximally m2 = np.subtract(m[2],100) m3 = np.multiply(m[3],500) # transformation from mM to uM and div by 2 (1000/2) # to be average of two same size pools instead of summation, # (the traces are the sum of the two...). ax2[0,0].plot(m2,m3, lw=lw, color=color[i], alpha=a) if i <= 1: if i == 0: a = 0.7 ax2[1,0].plot(m2,m3, lw=lw, color=color[i], alpha=a) mean_amp = [] spread = [] x = [] for d in distances: if d in res: mean_amp.append( np.mean(res[d]) ) spread.append( np.std( res[d]) ) x.append(d) if i == 0: for m in M: if m[1] >= 40: ax.plot(m[1], np.divide(m[0], mean_amp[3]), '.', ms=20, color='k', alpha=0.2) #day et al 2008 [x1,y1] = np.loadtxt('../../results/bAP/bAP-DayEtAl2006-D1.csv', unpack=True) ax.plot(x1, y1, 'brown', lw=6) # Hide the right and top spines ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines ax.yaxis.set_ticks_position('left') ax.xaxis.set_ticks_position('bottom') # set ticks yticks = [0,1] ax.set_yticks(yticks) ax.set_yticklabels(yticks, fontsize=35) xticks = [40, 120, 200] ax.set_xticks(xticks) ax.set_xticklabels(xticks, fontsize=35) ax.tick_params(width=2, length=4) ax.set_xlim([0,300]) ax.set_ylim([0,1.2]) # size of frame for axis in ['bottom','left']: ax.spines[axis].set_linewidth(4) elif i == 3: tag = fString.split('/')[2] for j in [0,1]: ax2[0,j].set_ylim([0,6.5]) ax2[1,j].set_ylim([0,1.3]) ax2[j,0].set_xlim([-5, 105]) ax2[j,1].set_xlim([-5, 105]) ax2[0,j].set_xticks([]) ax2[1,j].set_xticks([0,100]) ax2[1,j].set_xticklabels([0,100], fontsize=30) ax2[0,0].set_yticks([0,6]) ax2[0,0].set_yticklabels([0,6], fontsize=30) ax2[0,1].set_yticks([]) ax2[1,0].set_yticks([0,1.2]) ax2[1,0].set_yticklabels([0,1.2], fontsize=30) ax2[1,1].set_yticks([]) for a in ax2.reshape(-1): a.tick_params(width=4, length=6) # Hide the right and top spines a.spines['right'].set_visible(False) a.spines['top'].set_visible(False) # Only show ticks on the left and bottom spines a.yaxis.set_ticks_position('left') a.xaxis.set_ticks_position('bottom') # line width for axis in ['bottom','left']: a.spines[axis].set_linewidth(4) print mean_amp mean_amp = np.divide(mean_amp, mean_amp[3]) return x, mean_amp, ax2 def plot_HBP(): fig, axes = plt.subplots(2, 1, figsize=(12,24)) # sub fig 1 -------------------------------------------------------------------- files = glob.glob('[!v]*.out') for f in sorted(files): print f label=f.split('.')[0] if label == 'ach': label = 'ACh' ls = '--' c = 'k' a = 0.6 elif label == 'da': label = label.upper() ls = ':' c = 'k' a = 1 elif label == 'pka': label = label.upper() ls = '-' c = 'r' a = 1 # unpack, plot and set label axes[0].plot(*np.loadtxt(f,unpack=True), linestyle=ls, \ color=c, label=label, linewidth=4, alpha=1) # sub fig 2 ------------------------------------------------------------------- files = glob.glob('v*.out') for f in files: print f label=f.split('_')[1].split('.')[0] label = label[0].upper() + label[1:] # unpack, plot and set label axes[1].plot(*np.loadtxt(f,unpack=True), label=label, linewidth=4) # axes[0].legend(fontsize=26) axes[0].set_title('Substrates', fontsize=40) axes[0].set_ylabel('Concentrarion [nM]', fontsize=28) axes[0].set_xticks([]) axes[0].set_xlim([0, int(sys.argv[1])]) axes[0].set_yticks([0, 100, 200, 500, 1000, 1500]) axes[0].tick_params(labelsize=20) axes[1].plot([0,10000], [-50, -50], 'k--') axes[1].legend(fontsize=26, loc='best') axes[1].set_title('Somatic membrane potential', fontsize=40) axes[1].set_ylabel('Potential [mV]', fontsize=28) axes[1].set_xlabel('Time [ms]', fontsize=28) axes[1].set_xticks([0, 2000, 4000, 6000, 8000, 10000]) axes[1].set_xlim([0, int(sys.argv[1])]) axes[1].set_yticks([-90, -60, -30, 0, 30, 60]) axes[1].tick_params(labelsize=20) plt.show() #----------------------------------------------------------------------------------------- def hinton(matrix, max_weight=None, ax=None, col=None, row=None): """Draw Hinton diagram for visualizing a weight matrix. from: http://python-for-multivariate-analysis.readthedocs.io/a_little_book_of_python_for_multivariate_analysis.html """ ax = ax if ax is not None else plt.gca() if not max_weight: max_weight = 2**np.ceil(np.log(np.abs(matrix).max())/np.log(2)) ax.patch.set_facecolor('lightgray') ax.set_aspect('equal', 'box') ax.xaxis.set_major_locator(plt.NullLocator()) ax.yaxis.set_major_locator(plt.NullLocator()) if col and row: w = matrix[col][row] color = 'orange' if w > 0 else 'black' print w, color size = np.sqrt(np.abs(w)) rect = plt.Rectangle([0.5 - size / 2, 0.5 - size / 2], size, size, facecolor=color, edgecolor=color) ax.add_patch(rect) else: for (x, y), w in np.ndenumerate(matrix): color = 'orange' if w > 0 else 'black' size = np.sqrt(np.abs(w)) rect = plt.Rectangle([x - size / 2, y - size / 2], size, size, facecolor=color, edgecolor=color) ax.add_patch(rect) nticks = matrix.shape[0] #ax.xaxis.tick_top() ax.set_xticks(range(nticks)) ax.set_xticklabels(list(matrix.columns), rotation=45, fontsize=30) ax.set_yticks(range(nticks)) ax.set_yticklabels(matrix.columns, fontsize=30) ax.grid(False) ax.autoscale_view() ax.invert_yaxis() def rpredict(lda, X, y, out=False): ret = {"class": lda.predict(X), "posterior": pd.DataFrame(lda.predict_proba(X), columns=lda.classes_)} ret["x"] = pd.DataFrame(lda.fit_transform(X, y)) ret["x"].columns = ["LD"+str(i+1) for i in range(ret["x"].shape[1])] if out: print("class") print(ret["class"]) print() print("posterior") print(ret["posterior"]) print() print("x") print(ret["x"]) return ret def pretty_scalings(lda, X, out=False): ret = pd.DataFrame(lda.scalings_, index=X.columns, columns=["LD"+str(i+1) for i in range(lda.scalings_.shape[1])]) if out: print("Coefficients of linear discriminants:") display(ret) return ret def ldahist(data, g, sep=False): xmin = np.trunc(np.min(data)) - 1 xmax = np.trunc(np.max(data)) + 1 ncol = len(set(g)) binwidth = 0.5 print binwidth, (xmax + binwidth -xmin)/binwidth bins=np.arange(xmin, xmax + binwidth, binwidth) if sep: fig, axl = plt.subplots(ncol, 1, sharey=True, sharex=True) else: fig, axl = plt.subplots(1, 1, sharey=True, sharex=True) axl = [axl]*ncol for ax, (group, gdata) in zip(axl, data.groupby(g)): sns.distplot(gdata.values, bins, ax=ax, label="group "+str(group)) ax.set_xlim([xmin, xmax]) if sep: ax.set_xlabel("group"+str(group)) else: ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)) plt.tight_layout() def pca_scatter(pca, standardised_values, classifs): foo = pca.transform(standardised_values) bar = pd.DataFrame(zip(foo[:, 0], foo[:, 1], classifs), columns=["PC1", "PC2", "Class"]) sns.lmplot("PC1", "PC2", bar, hue="Class", fit_reg=False) def pca_summary(pca, standardised_data, out=True): names = ["PC"+str(i) for i in range(1, len(pca.explained_variance_ratio_)+1)] a = list(np.std(pca.transform(standardised_data), axis=0)) b = list(pca.explained_variance_ratio_) c = [np.sum(pca.explained_variance_ratio_[:i]) for i in range(1, len(pca.explained_variance_ratio_)+1)] columns = pd.MultiIndex.from_tuples([("sdev", "Standard deviation"), ("varprop", "Proportion of Variance"), ("cumprop", "Cumulative Proportion")]) summary = pd.DataFrame(zip(a, b, c), index=names, columns=columns) if out: print("Importance of components:") display(summary) return summary def run_PCA(df, channels, marker): ''' runs principal component analysis on data frame ''' colors = ['b', 'grey', 'r'] alpha = [0.004, 0.04, 0.005] #df.boxplot(by='group') f_corr, a_corr = plt.subplots(1,3) f_sc, a_sc = plt.subplots(2,3, sharey=True, sharex=True) f_sc.subplots_adjust(wspace=0.5) #hspace=0.5, f_sc3d, a_sc3d = plt.subplots(1,3, sharey=True) #f_sc3d.subplots_adjust(wspace=0.5) fig = plt.figure() gs = gridspec.GridSpec(2, 7, height_ratios=[4, 1], hspace = .6) ax1 = plt.subplot(gs[0, :]) axes = [] for i in range(7): axes.append( plt.subplot(gs[1, i]) ) ax1.plot([-1,27], [100, 100], '--k', lw=3) # loop over groups for g,group in enumerate([0.0, 1.0, 2.0]): # correlation between (Na and K+) channels df1 = df.loc[df['group'] == group, channels[0:4]] print df1 corrmat = df1.corr() #method='spearman') hinton(corrmat, ax=a_corr[g]) if g > 0: a_corr[g].set_yticks([]) kaf = df.loc[df['group'] == group, 'kaf'] naf = df.loc[df['group'] == group, 'naf'] kir = df.loc[df['group'] == group, 'kir'] kas = df.loc[df['group'] == group, 'kas'] chan_2b_paired = [naf, kir] position = [] mean = [[],[],[],[],[],[],[]] a_sc3d[g].plot([0,200], [100,100], 'k') a_sc3d[g].plot([100,100], [0,200], 'k') sc = a_sc3d[g].scatter(kaf, naf, c=kir, cmap='copper', s=kas, lw=0.5, edgecolor='g', alpha=1) fs = 14 if g == 2: # colorbar fc,ac = plt.subplots() cbar = fc.colorbar(sc, ticks=[0, 100, 200], orientation='horizontal') cbar.ax.set_xticklabels([0, 100, 200], fontsize=fs) ac.remove() a_sc3d[g].set_xticks([0,100,200]) a_sc3d[g].set_yticks([0,100,200]) a_sc3d[g].set_xlim([0,200]) a_sc3d[g].set_ylim([0,200]) a_sc3d[g].set(adjustable='box-forced', aspect='equal') a_sc3d[g].set_xticklabels([0,100,200], fontsize=fs) a_sc3d[g].set_yticklabels([0,100,200], fontsize=fs) # pick specific channel with value g for c,channel in enumerate(channels): if channel in ['naf','kir']: a = ['naf','kir'].index(channel) chan = chan_2b_paired[a] a_sc[a,g].plot([0,200], [100,100], 'k') a_sc[a,g].plot([100,100], [0,200], 'k') a_sc[a,g].plot(kaf, chan, 'o', \ ms=8, color=colors[g], alpha=alpha[g]) a_sc[a,g].set_xticks([0,100,200]) a_sc[a,g].set_yticks([0,100,200]) a_sc[a,g].set(adjustable='box-forced', aspect='equal') a_sc[a,g].set_xticklabels([0,100,200], fontsize=24) a_sc[a,g].set_yticklabels([0,100,200], fontsize=24) # add channel with specific group to list df1 = df.loc[df['group'] == group, channel] mean[c].append(df1) # append position position.append( g+c*4 ) # plot lw=3 boxprops = dict(linewidth=lw, color=colors[g]) medianprops = dict(linestyle='-', linewidth=lw, color=colors[g]) ax1.boxplot( mean, \ positions=position, \ boxprops=boxprops, \ medianprops=medianprops ) ax1.set_xlim([-1,27]) ax1.set_xticks(np.arange(1,27,4)) ax1.set_yticks([0,100,200]) ax1.set_xticklabels(channels, fontsize=24, rotation=45) ax1.set_yticklabels([0,100,200], fontsize=24) # calc plot correlation of channel with excitability corrmat = df.corr() print corrmat['group'] for i in range( len(axes) ): hinton(corrmat, ax=axes[i], col='group', row=i+1) # scatter plot kaf vs naf or kir fig.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_chan-exc-cor_pdc.png']), bbox_inches='tight', transparent=True) f_corr.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_chan-cor_pdc.png']), bbox_inches='tight', transparent=True) f_sc.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_scatter_pdc.png']), bbox_inches='tight') f_sc3d.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_scatter3d_pdc.png']), bbox_inches='tight') #uncomment to save colorbar fc.savefig(''.join(['../../../Dropbox/manuscript/Frontiers/Figures/static_', marker, '_scatter_colorbar.png']), bbox_inches='tight') plt.close('all') ''' # principal component analysis y = df.group # dependent variable data standardisedX = scale(df) standardisedX = pd.DataFrame(standardisedX, index=df.index, columns=df.columns) # principal component analysis (PCA) pca = PCA().fit(standardisedX) pca_scatter(pca, standardisedX, y) # Linear discriminant analysis (LDA) lda = LinearDiscriminantAnalysis().fit(df, y) #pretty_scalings_ = pretty_scalings(lda, df, out=True) lda_values = rpredict(lda, df, y, True) ldahist(lda_values["x"].LD1, y) sns.lmplot("LD1", "LD2", lda_values["x"].join(y), hue="group", fit_reg=False); plt.show() for group in [0.0, 1.0, 2.0]: df1 = df[df['group'] == group] df1 = df1.drop('group', 1) # standardize data (mean and std) standardisedX = scale(df1) standardisedX = pd.DataFrame(standardisedX, index=df1.index, columns=df1.columns) pca = PCA().fit(standardisedX) summary = pca_summary(pca, standardisedX) print summary.sdev**2 print pca.components_[0] #pca_scatter(pca, standardisedX, y)''' def read_and_reshape(pickled_string): ''' reads pickled data in the form created by repickle_static2, reshapes to array and stores in pandas DataFrame ''' files = glob.glob('Beskow/save/'+pickled_string+'*') for f in files: # result structure quant = {'more':{}, 'equal':{}, 'less':{}} counter = 1 I = np.arange(320, 345, 10) # load pickle DATA = load_obj(f ) # create array N = len(DATA) A = np.zeros((N,8)) for n in DATA: # flags that checks that all modulation factors are within range and keeps track of rheobase current noSpikes = None mod_factors = DATA[n]['factors'] channels = DATA[n]['channels'] # if first iteration, setup of result structure if not 'can' in quant['more']: for channel in channels: for key in quant: quant[key][channel] = [] for amp in I: if len(DATA[n][amp]['spikes']) == 0: noSpikes = amp add2quant(quant, noSpikes, channels, mod_factors) # counter c = 0 for i,group in enumerate(['less', 'equal', 'more']): # fill array # length of subset l = len(quant[group]['naf']) # create array with group A[c:c+l,0] = np.ones(l) * i for j,channel in enumerate(channels): A[c:c+l,j+1] = quant[group][channel] c += l df = pd.DataFrame(A, columns=['group']+channels) # analyse (pca) run_PCA(df, channels, pickled_string) def repickle_and_plot_uniform(): ''' repickles data from uniform distribution and plot what tag change (interval) save Tag ____________________________________________________________________________________________ full range StatUniFull STATuniform_Full +kaf StatUnikaf [0.97,1.03] -> [0.75,0.85] STATuniform_kaf kir relaxed StatUniKir [0.85,1.25] -> [0.00,1.25] STATuniform_Kir naf relaxed StatUniNaf [0.60,0.80] -> [0.60,1.00] STATuniform_Naf kas relaxed StatUniKas [0.65,0.85] -> [0.00,0.85] STATuniform_kas full-noAxon StatUniFull-noAxon- STATuniform_full-noAxon kaf-noAxon StatUnikaf-noAxon [0.97,1.03] -> [0.75,0.85] STATuniform_kaf-noAxon ''' # list of tags all_files = ['StatUniFull-na', \ 'StatUnikaf-na', \ 'StatUniKir', \ 'StatUniNaf', \ 'StatUniKas', \ 'StatUniFull-no', \ 'StatUnikaf-no' ] repickle_tag = ['STATuniform_Full', \ 'STATuniform_Kaf', \ 'STATuniform_Kir', \ 'STATuniform_Naf', \ 'STATuniform_Kas', \ 'STATuniform_Full-noAxon', \ 'STATuniform_Kaf-noAxon' ] for i,tag in enumerate(all_files): # repickle repickle_static2(tag+'*', 'Beskow/save/'+repickle_tag[i]) # plot plot_static_modulation_pdc('Beskow/save/'+repickle_tag[i], marker=repickle_tag[i]) def plot_substrates(): ''' plot substrate values (normalized to 1). copied from "plot_synMod_dyn_pdc" ''' f2, a2 = plt.subplots(1,1, figsize=(4.0, 3) ) #, gridspec_kw = {'height_ratios':[2, 2, 2]} ) colors = ['#bf5b17', '#4daf4a', '#984ea3', '#ff7f00'] # plot substrates files = ['Target1p', 'cAMP', 'Gbgolf', 'D1RDAGolf', ] peak_da = ['500', '1500'] #, '4500', '9500'] for i, peak in enumerate(peak_da): for c,f in enumerate(files): x, y = np.loadtxt('substrate_'+f+peak+'.out', unpack=True) y = np.divide(y, max(y)) #y = np.subtract(y, 10) x = np.subtract(x, 1000) a2.plot(x,y, label=f, lw=3, color=colors[c]) plt.show() def plot_fig6B(directory, SPIKES): ''' reporduces fig 6B, except for that the random traces are different copied from "plot_synMod_dyn_pdc" ''' f2, a2 = plt.subplots(2,1, figsize=(4.0, 6), gridspec_kw = {'height_ratios':[2, 2]} ) colors = ['#bf5b17', '#4daf4a', '#984ea3', '#ff7f00'] # plot noise and example trace num_cores = multiprocessing.cpu_count() files = glob.glob(''.join([directory, 'spiking_[0-9]_naf*.out'])) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) for x,y,dy,f in M: x = np.subtract(x,1000) a2[0].plot(x,y, label=f, lw=2, color='k') files = glob.glob(''.join([directory, 'spiking_[0-9]_control.out'])) M = Parallel(n_jobs=num_cores)(delayed(load_file)( f ) for f in files) for x,y,dy,f in M: x = np.subtract(x,1000) a2[0].plot(x,y, label=f, lw=2, color=[0.5,0.5,0.5]) a2[0].plot([-250,1000], [-70,-70], 'k--', lw=3) # plot substrates files = ['Target1p', 'cAMP', 'Gbgolf', 'D1RDAGolf', ] for c,f in enumerate(files): x, y = np.loadtxt(''.join([directory, 'substrate_', f, '.out']), unpack=True) y = np.divide(y, max(y)/50) x = np.subtract(x, 1000) a2[0].plot(x,y, label=f, lw=3, color=colors[c]) names = ['PKA', \ 'cAMP', \ r'G$_{\beta\gamma}$', \ 'D1R-G'] t_first_spike = [[],[],[],[]] for t,target in enumerate(files): for run in SPIKES[target]: if len( SPIKES[target][run] ) > 0: t_first_spike[t].append( SPIKES[target][run][0]-1000 ) print t_first_spike medianprops = dict(linestyle='-', linewidth=2, color='k') boxprops = dict(linewidth=2, color='grey') # over substrates (PKA, cAMP, Gbg, DIR)--reversed bp2 = a2[1].boxplot([t_first_spike[i] for i in [3,2,1,0]], \ labels=[ names[i] for i in [3,2,1,0]], \ vert=False, \ medianprops=medianprops, \ boxprops=boxprops) format_boxplot(bp2, median=True, lw=2, colors=colors[::-1]) for median in range(len(bp2['medians'])): print print bp2['medians'][median].get_xdata() for ax in a2[0:1]: ax.set_xlim([-250, 1000]) a2[1].set_xticks( [-250,0,500,1000] ) a2[1].spines['left'].set_visible(False) a2[1].spines['top'].set_visible(False) #a2[1].set_yticks([]) a2[1].set_ylabel('Time to spike', fontsize=20) a2[0].set_ylabel('Example Substrates') a2[1].xaxis.set_ticks_position('bottom') # size of frame for axis in ['right', 'bottom']: a2[1].spines[axis].set_linewidth(2) a2[0].set_ylim([-78, 55]) a2[0].axis('off') plt.show() def plot_fig4C(spike_dict): ''' "reproduce" the raster and pie chart from fig 4A (with smaller sample set). code base from functions "plot_static_modulation_pdc" and "plot_modulation_dynamical_pdc" ''' # raster plot------------------------------------------------------------------------- f, a = plt.subplots(3, 1, figsize=(4,3)) currents = np.arange(320,345,10) for c,I in enumerate(currents): i = len(currents) - (c+1) counter = 0 for ID in spike_dict: # create y vector and plot raster if ID == 0: y = np.ones( len(spike_dict[ID][I]) ) * (len(spike_dict)-1) * 0.5 if len(spike_dict[ID][I]) > 0: a[i].plot(spike_dict[ID][I], y, '|', color=[0.2,0.2,0.2], ms=len(spike_dict)*5, mew=4, alpha=0.5) else: y = np.ones( len(spike_dict[ID][I]) ) + counter - 1 if len(spike_dict[ID][I]) > 0: a[i].plot(spike_dict[ID][I], y, 'k|', mew=2, ms=3, alpha=1.0) counter += 1 ylimits = [-0.5, len(spike_dict)-0.5] a[i].set_ylim(ylimits) a[i].set_yticks([]) a[i].set_ylabel(str(I), fontsize=18) a[i].spines['left'].set_linewidth(2) a[i].set_xlim([0,1000]) if not i == len(currents) - 1: a[i].set_xticks([]) plt.tick_params( axis='x', # changes apply to the x-axis which='both', # both major and minor ticks are affected bottom='off', # ticks along the bottom edge are off top='off', # ticks along the top edge are off labelbottom='off') # labels along the bottom edge are off xticks = [0,500,1000] a[-1].xaxis.set_ticks(xticks) a[-1].set_xticklabels(xticks, fontsize=18) a[-1].spines['bottom'].set_visible(True) a[-1].xaxis.set_ticks_position('bottom') a[-1].spines['bottom'].set_linewidth(2) a[0].spines['top'].set_linewidth(1) a[-1].set_xlabel("Time (ms)", fontsize=20) # pie plot --------------------------------------------------------------------------- plt.figure(figsize=(4,3)) quant = {320:0, 330:0, 340:0} # quantify by finding the lowest current amplitude that make the cell spikes for n in range(1, len(spike_dict)): for i,I in enumerate(currents): if len(spike_dict[n][I]) > 0 or i == 2: quant[I] += 1 break N = [ quant[320], quant[330], quant[340] ] color = ['r', 'grey', 'b'] plt.pie(N, colors=color, autopct=pie_autopct, shadow=False) plt.rcParams['font.size'] = 50 plt.axis('equal') # channel distributions ----------------------------------------------------- channels = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ] f5, a5 = plt.subplots(len(channels), 1, figsize=(6,6), sharex=True) bins = np.linspace(0, 200, 50) TOT = {} quant = { 320: {}, 330:{}, 340:{} } for n in range(1, len(spike_dict)): for c, I in enumerate(currents): if 'kaf' not in quant[I]: for chan in channels: quant[I][chan] = [] if len(spike_dict[n][I]) > 0 or c == 2: for i, chan in enumerate(channels): quant[I][chan].append( spike_dict[n]['factors'][i]*100 ) break for i, chan in enumerate(channels): for c, cathegory in enumerate([320, 330, 340]): weights = np.ones_like(quant[cathegory][chan])/float(len(quant[cathegory][chan]))*100.0 a5[i].hist( quant[cathegory][chan], bins=bins, weights=weights, color=color[c], alpha=0.2, histtype='stepfilled') a5[i].hist( quant[cathegory][chan], bins=bins, weights=weights, color=color[c], histtype='step', lw=3) for a in [a5]: a[i].set_yticks([]) a[i].set_xticks([0,100,200]) a[i].set_xticklabels([0,100,200], fontsize=24) a[i].spines['left'].set_visible(False) a[i].spines['top'].set_visible(False) a[i].spines['right'].set_visible(False) a[i].spines['bottom'].set_visible(True) a[i].xaxis.set_ticks_position('bottom') a[i].spines['bottom'].set_linewidth(2) plt.show()