'''
Jonathon Reed, Avrama Blackwell
Analysis of simulation of neuron morphologies
comparing original morphology with predicted radii
Using two types of current injection and 1 location of synaptic input
2021 March 12
'''
import matplotlib.pyplot as plt
plt.ion()
from matplotlib.ticker import StrMethodFormatter
import numpy as np
from scipy import optimize
import os
import glob
convert2milli=1000
def save_png(png, title): #default to save plots instead of plot to window
png.savefig(title, bbox_inches = 'tight')
print('File Created : ' + str(title))
png.close()
def exp2(t,tau1,tau2,A,B,C):
return A*np.exp(-t/tau1)+B*np.exp(-t/tau2)+C #0th tau --> single soma as RM*CM
#tau1 = 0th tau(1+(n*pi/elec_len)**2) where n = tau #
#tau --> short sims beginning at peak
#change in voltage --> long sims @ steady reach
def exp1(t,tau1,A,B):
return A*np.exp(-t/tau1)+B #0th tau --> single soma as RM*CM
def main_plot(data, var, colors,title = '', ax_titles = None, labels = None, save_as = None, where = 'best', size = None,xmax=None):
#plt.rc('font', size = 36) #default plot and font sizes
if save_as:
scale=2
else:
scale=1
plt.rc('font', size = scale*9)
fig=plt.figure(figsize = (scale*4,scale*2))
for key,traces in data.items():
if labels:
plt.plot(data[key][var[0]], data[key][var[1]], label = labels[key], linewidth = scale*1,color=colors[key])
else:
plt.plot(data[key][var[0]], data[key][var[1]], linewidth = scale*1,color=colors[key])
if labels:
plt.legend(loc = where)
if xmax:
plt.xlim([0,xmax])
plt.xlabel(ax_titles[0]) if ax_titles else plt.xlabel(var[0])
plt.ylabel(ax_titles[1]) if ax_titles else plt.ylabel(var[1])
#plt.gca().xaxis.set_major_formatter(StrMethodFormatter('{x:,.2f}'))
#plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.2f}'))
plt.tight_layout()
if save_as:
#save_png(plt,save_as + '.png')# if save_as else save_png(plt,title + '.png')
fig.savefig(save_as+'.tiff')
else:
plt.title(title)
plt.show()
return plt
def show_curvefit(func,sim_data,popt,rkey,start,comp='somaVm',title=None):
estimate = func(np.array(sim_data[rkey]['time'][start:]),*popt) #send in all curve_fit params for estimate
#show result of curve fit
plt.figure()
if title:
plt.title(title)
else:
plt.title(rkey)
plt.plot(sim_data[rkey]['time'][start:],sim_data[rkey][comp][start:],'o')
plt.plot(sim_data[rkey]['time'][start:], estimate)
#save_png(plt,'check_eq.png')
#plt.close()
_FUNCTIONS={
'exp1': exp1,
'exp2': exp2}
####### Uncomment the root file name (pattern) and synapse compartment of one of these #########
#only saving syntau and dVsyn in soma for one of the synapses
pattern = 'Purkinje-slice-ageP35-1.CNG_*' #matches quite well to pred+in. Looks very similar to original plot (shape, not values)
syn=['12_3'] #synaptic input ,,'1342_3',
#pattern = 'WT-0201MSN03.CNG_*' #Lai striatum - matches quite well, including synaptic input
#syn=['19_3']#,'35_3',
#pattern = 'WT-P270-20-14ak.CNG_*' #Lindroos striatum -
#syn=['35_3']#'194_3',
#pattern = 'ri06_mod.CNG_*' #Golding hippocampus - good fit, good synapse.
#syn=['31_4']#,'3400_3'] #
#pattern = 'Adol-20100419cell1.CNG_*' #Groen - hippocampus - good fit, good synapse.
#syn=['32_3']#,'112_3'] #response is ~0.1 mV
### to fit to single exponential, change to 'exp1'
fit='exp2'
param_bounds = ([0,0,-100,-100,-100],[1000,1000, 100,100,100])#msec,msec,mV,mV,mV
showfit=True
#parameters to control plotting, specific to the set of simulations, but same for all file names
inj_set={'short' : '0.001', 'long': '0.8'}
syn_set={'syn'+str(i):syn[i] for i in range(len(syn))}
stim_set=dict(syn_set,**inj_set)
labels={'org':'Original','pred':'Predict','pred_i':'Predict+In.','zdia2':'2.0 Diameter'}
colors={'org':'k','pred':'b','pred_i':'g','zdia2':'gray'}
ax_titles = ['Time (ms)','Vm (mV)']
#dictionary to store summary measures
measures={'tau':{},'tau2':{},'deltaV':{},'syntau':{},'syntau2':{},'dVsyn':{}}
fit_cov={}
for stimtype,stim in stim_set.items():
fnames=sorted(glob.glob(pattern+stim+'*.txt'))
sim_data = {}
for fname in fnames:
rkey=fname.split(stim)[0].split('CNG_')[-1]
header=open(fname).readline()
data_keys=header.split()[1:]
sim_data[rkey] = {dk:[] for dk in data_keys}
print('processing',fname,' .p file type=',rkey,' data columns',data_keys)
data=np.loadtxt(fname)
for i, dk in enumerate(data_keys):
sim_data[rkey][dk]=data[:,i]*convert2milli
if len(data_keys)>2:
if 'dVsyn'+'_'+data_keys[2] not in measures.keys():
for kk in ['dVsyn','syntau','syntau2']:
measures[kk+'_'+data_keys[2]]={}
if stimtype=='short': #determine time constants
start = np.argmax(sim_data[rkey]['somaVm']) #point at which is highest voltage
#print('exp start', start,sim_data[rkey]['time'][start])
if fit=='exp2':
#curve fit to double tau, needs parameter constraints
popt,pcov = optimize.curve_fit(exp2,sim_data[rkey]['time'][start:],sim_data[rkey]['somaVm'][start:],bounds=param_bounds,maxfev=5000)
measures['tau'][rkey]=max(popt[0],popt[1])
measures['tau2'][rkey]=min(popt[0],popt[1])
else:
popt,pcov = optimize.curve_fit(exp1,sim_data[rkey]['time'][start:],sim_data[rkey]['somaVm'][start:])
measures['tau'][rkey] = popt[0]
fit_cov[rkey]=pcov[0,0]
if showfit:
show_curvefit(_FUNCTIONS[fit],sim_data,popt,rkey,start)
elif stimtype=='long': #determine deltaV (proportional to input resistance)
delay = 50e-3
pt = np.max(np.where(np.array(sim_data[rkey]['time'])<delay))
baseV=sim_data[rkey]['somaVm'][pt]
peakpt=np.argmax(sim_data[rkey]['somaVm'])
peakV=sim_data[rkey]['somaVm'][peakpt]
measures['deltaV'][rkey] = peakV-baseV
elif stimtype.startswith('syn'): #determine peak time and peak Vm
delay=10e-3
for comp in list(sim_data[rkey].keys())[1:]:
peakpt=np.array(sim_data[rkey][comp]).argmax()
print(' analysis of',stimtype,comp, 'for', rkey, 'peakpt=',peakpt)
peakV=sim_data[rkey][comp][peakpt]
basept = np.max(np.where(np.array(sim_data[rkey]['time'])<delay))
baseV=sim_data[rkey][comp][basept]
if comp=='somaVm':
add_to_key=''
else:
add_to_key='_'+comp
if fit=='exp2':
popt,pcov = optimize.curve_fit(exp2,sim_data[rkey]['time'][peakpt:],sim_data[rkey][comp][peakpt:],bounds=param_bounds,maxfev=5000)
measures['syntau2'+add_to_key][rkey]=min(popt[0],popt[1])
measures['syntau'+add_to_key][rkey]=max(popt[0],popt[1])
else:
popt,pcov = optimize.curve_fit(exp1,sim_data[rkey]['time'][peakpt:],sim_data[rkey][comp][peakpt:])
measures['syntau'+add_to_key][rkey]=popt[0]
measures['dVsyn'+add_to_key][rkey]=peakV-baseV
if showfit:
show_curvefit(_FUNCTIONS[fit],sim_data,popt,rkey,peakpt,comp=comp,title=rkey+add_to_key)
#plot all morphologies for one type of stimulation
title=fname.split('.')[0]
if title == 'ri06_mod':
title = 'ri06'
print(title)
savename = title+'_'+stim_set[stimtype]
main_plot(sim_data, data_keys, colors, title=title+' '+stim_set[stimtype], ax_titles=ax_titles, labels = labels, where = 'upper right',save_as=savename)
if stimtype.startswith('syn'):
savename = title+'_'+'dend_'+stim_set[stimtype]
main_plot(sim_data, np.array(data_keys)[[0,2]], colors, title='dend '+stim_set[stimtype], ax_titles=ax_titles, labels = labels, where = 'upper right',xmax=200, save_as=savename)
#calculate the difference as ratio abs((orig-new)/orig)
diff={k:{} for k in measures.keys()}
for meastype,measure in measures.items():
if len(measure):
for key in measure.keys():
diff[meastype][key]=np.abs(measure[key]/measure['org']-1)
print ('*******', meastype,'******')
print([(k,round(m,5)) for k,m in measure.items()],'\ndiff = ',[(k,round(d,3)) for k,d in diff[meastype].items()])
if meastype=='tau':
print('fit', fit_cov[key])
##### Save measures to a file #####3
outfname=pattern.split('.')[0]+'.measures'
f=open(outfname,'w')
for key in diff:
f.write(key+'='+str(measures[key])+'\n')
f.write('diff='+str(diff[key])+'\n')
f.close()
####### to save figure as tif:
#import PIL
#fig.savefig('fname.tif')