################################################################################
# Import, Set up MPI Variables, Load Necessary Files
################################################################################
from __future__ import division
import os
from os.path import join
import sys
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
np.seterr(divide='ignore', invalid='ignore')
import scipy
import scipy.fftpack
from scipy import signal as ss
from scipy import stats as st
import bootstrapped.bootstrap as bs
import bootstrapped.stats_functions as bs_stats
import itertools
import pandas as pd
font = {'family' : 'normal',
'weight' : 'normal',
'size' : 24}
matplotlib.rc('font', **font)
matplotlib.rc('legend',**{'fontsize':16})
controlm = ['y','o']
TuningType = 'both'
colors = ['dimgray', 'crimson', 'green', 'darkorange']
N_cells = 1000.
N_HL5PN = int(0.70*N_cells)
N_HL5PNso = int(N_HL5PN*0.5)
dt = 0.025
tstop = 7000
transient = 2000
stimtime = 6500 # Absolute time during stimulation that is set for stimtime
stimbegin = 3 # Time (ms) after stimtime where stimulation begins (i.e. the system delay)
stimend = 103 # Stimbegin + duration of stimulus
startslice = stimtime+stimbegin
endslice = stimtime+stimend
totalstim = endslice-startslice
tvec = np.arange(startslice,endslice+1,dt,dtype=np.float64)
synnums = [85,95]
N_seeds = 80
rndinds = np.linspace(1,N_seeds,N_seeds, dtype=int)
def cohen_d(y,x):
nx = len(x)
ny = len(y)
dof = nx + ny - 2
return (np.mean(x) - np.mean(y)) / np.sqrt(((nx-1)*np.std(x, ddof=1) ** 2 + (ny-1)*np.std(y, ddof=1) ** 2) / dof)
def barplot_annotate_brackets(num1, num2, data, center, height, yerr=None, dh=.05, barh=.05, fs=None, maxasterix=None):
"""
Annotate barplot with p-values.
:param num1: number of left bar to put bracket over
:param num2: number of right bar to put bracket over
:param data: string to write or number for generating asterixes
:param center: centers of all bars (like plt.bar() input)
:param height: heights of all bars (like plt.bar() input)
:param yerr: yerrs of all bars (like plt.bar() input)
:param dh: height offset over bar / bar + yerr in axes coordinates (0 to 1)
:param barh: bar height in axes coordinates (0 to 1)
:param fs: font size
:param maxasterix: maximum number of asterixes to write (for very small p-values)
"""
if type(data) is str:
text = data
else:
# * is p < 0.05
# ** is p < 0.005
# *** is p < 0.0005
# etc.
text = ''
p = .05
while data < p:
text += '*'
p /= 10.
if maxasterix and len(text) == maxasterix:
break
if len(text) == 0:
text = 'n. s.'
lx, ly = center[num1], height[num1]
rx, ry = center[num2], height[num2]
if yerr:
ly += yerr[num1]
ry += yerr[num2]
ax_y0, ax_y1 = plt.gca().get_ylim()
dh *= (ax_y1 - ax_y0)
barh *= (ax_y1 - ax_y0)
y = max(ly, ry) + dh
barx = [lx, lx, rx, rx]
bary = [y, y+barh/2, y+barh/2, y]
mid = ((lx+rx)/2, y+barh*(3/4))
plt.plot(barx, bary, c='black')
kwargs = dict(ha='center', va='bottom')
if fs is not None:
kwargs['fontsize'] = fs
plt.text(*mid, text, **kwargs)
def gaussian(x, mu, sig):
return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
Ratesy_allseeds1 = [[] for _ in range(0,len(rndinds))]
Rateso_allseeds1 = [[] for _ in range(0,len(rndinds))]
Ratesy_allseeds2 = [[] for _ in range(0,len(rndinds))]
Rateso_allseeds2 = [[] for _ in range(0,len(rndinds))]
SNRy_allseeds1 = [[] for _ in range(0,len(rndinds))]
SNRo_allseeds1 = [[] for _ in range(0,len(rndinds))]
SNRy_allseeds2 = [[] for _ in range(0,len(rndinds))]
SNRo_allseeds2 = [[] for _ in range(0,len(rndinds))]
BaseRatey_allseeds1 = [[] for _ in range(0,len(rndinds))]
BaseRateo_allseeds1 = [[] for _ in range(0,len(rndinds))]
BaseRatey_allseeds2 = [[] for _ in range(0,len(rndinds))]
BaseRateo_allseeds2 = [[] for _ in range(0,len(rndinds))]
for control in controlm:
if control == 'y': colc = 'k'
if control == 'o': colc = 'r'
path1 = control + '_' + str(synnums[0])+ '/HL5netLFPy/Circuit_output/'
path2 = control + '_' + str(synnums[1])+ '/HL5netLFPy/Circuit_output/'
for idx in range(0,len(rndinds)):
print('Seed #'+str(idx))
temp_sn1 = np.load(path1 + 'SPIKES_CircuitSeed1234StimSeed'+str(rndinds[idx])+'.npy',allow_pickle=True)
temp_sn2 = np.load(path2 + 'SPIKES_CircuitSeed1234StimSeed'+str(rndinds[idx])+'.npy',allow_pickle=True)
# Only analyze PNs here
SPIKES1 = [x for _,x in sorted(zip(temp_sn1.item()['gids'][0],temp_sn1.item()['times'][0]))]
SPIKES2 = [x for _,x in sorted(zip(temp_sn2.item()['gids'][0],temp_sn2.item()['times'][0]))]
for j in range(N_HL5PN):
tv_base = np.around(np.array(SPIKES1[j][(SPIKES1[j]>transient) & (SPIKES1[j]<stimtime)]),3)
rate_base = (len(tv_base)*1000)/(stimtime-transient) # convert to 1/seconds
tv = np.around(np.array(SPIKES1[j][(SPIKES1[j]>startslice) & (SPIKES1[j]<endslice)]),3)
rate = (len(tv)*1000)/(endslice-startslice) # convert to 1/seconds
SNR = rate/rate_base if (rate_base > 0) else 0
if control == 'y': Ratesy_allseeds1[idx].append(rate)
if control == 'o': Rateso_allseeds1[idx].append(rate)
if control == 'y': BaseRatey_allseeds1[idx].append(rate_base)
if control == 'o': BaseRateo_allseeds1[idx].append(rate_base)
if (control == 'y') & (rate_base > 0): SNRy_allseeds1[idx].append(SNR)
if (control == 'o') & (rate_base > 0): SNRo_allseeds1[idx].append(SNR)
tv_base = np.around(np.array(SPIKES2[j][(SPIKES2[j]>transient) & (SPIKES2[j]<stimtime)]),3)
rate_base = (len(tv_base)*1000)/(stimtime-transient) # convert to 1/seconds
tv = np.around(np.array(SPIKES2[j][(SPIKES2[j]>startslice) & (SPIKES2[j]<endslice)]),3)
rate = (len(tv)*1000)/(endslice-startslice) # convert to 1/seconds
SNR = rate/rate_base if (rate_base > 0) else 0
if control == 'y': Ratesy_allseeds2[idx].append(rate)
if control == 'o': Rateso_allseeds2[idx].append(rate)
if control == 'y': BaseRatey_allseeds2[idx].append(rate_base)
if control == 'o': BaseRateo_allseeds2[idx].append(rate_base)
if (control == 'y') & (rate_base > 0): SNRy_allseeds2[idx].append(SNR)
if (control == 'o') & (rate_base > 0): SNRo_allseeds2[idx].append(SNR)
# Initialize dataframe containing stats
df = pd.DataFrame(columns=["Cells", "Metric", "Mean Young", "SD Young", "Mean Old", "SD Old", "t-stat", "p-value", "Cohen's d"])
# All mean rates
MeanRatey = []
MeanRateo = []
MeanBaseRatey = []
MeanBaseRateo = []
PercentActivey = []
PercentActiveo = []
RateCVy = []
RateCVo = []
SNRy = []
SNRo = []
for i in range(0,len(rndinds)):
rvy = np.mean(Ratesy_allseeds1[i])
rvo = np.mean(Rateso_allseeds1[i])
MeanRatey.append(rvy)
MeanRateo.append(rvo)
# rvy = np.mean(Ratesy_allseeds2[i])
# rvo = np.mean(Rateso_allseeds2[i])
# MeanRatey.append(rvy)
# MeanRateo.append(rvo)
rvy = np.mean(BaseRatey_allseeds1[i])
rvo = np.mean(BaseRateo_allseeds1[i])
MeanBaseRatey.append(rvy)
MeanBaseRateo.append(rvo)
# rvy = np.mean(BaseRatey_allseeds2[i])
# rvo = np.mean(BaseRateo_allseeds2[i])
# MeanBaseRatey.append(rvy)
# MeanBaseRateo.append(rvo)
rvy = np.mean(SNRy_allseeds1[i])
rvo = np.mean(SNRo_allseeds1[i])
SNRy.append(rvy)
SNRo.append(rvo)
# rvy = np.mean(SNRy_allseeds2[i])
# rvo = np.mean(SNRo_allseeds2[i])
# SNRy.append(rvy)
# SNRo.append(rvo)
rvy = np.array(Ratesy_allseeds1[i])
rvo = np.array(Rateso_allseeds1[i])
PercentActivey.append((len(rvy[rvy>0])/N_HL5PN)*100)
PercentActiveo.append((len(rvo[rvo>0])/N_HL5PN)*100)
# rvy = np.array(Ratesy_allseeds2[i])
# rvo = np.array(Rateso_allseeds2[i])
# PercentActivey.append((len(rvy[rvy>0])/N_HL5PN)*100)
# PercentActiveo.append((len(rvo[rvo>0])/N_HL5PN)*100)
rvy = np.std(Ratesy_allseeds1[i])/np.mean(Ratesy_allseeds1[i])
rvo = np.std(Rateso_allseeds1[i])/np.mean(Rateso_allseeds1[i])
RateCVy.append(rvy)
RateCVo.append(rvo)
# rvy = np.std(Ratesy_allseeds2[i])/np.mean(Ratesy_allseeds2[i])
# rvo = np.std(Rateso_allseeds2[i])/np.mean(Rateso_allseeds2[i])
# RateCVy.append(rvy)
# RateCVo.append(rvo)
# Use this code for bootstrapping instead
bsMeans = False
if bsMeans:
x = bs.bootstrap(np.array(MeanRatey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey = x.value
lowerratey = x.lower_bound
upperratey = x.upper_bound
x = bs.bootstrap(np.array(MeanRateo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo = x.value
lowerrateo = x.lower_bound
upperrateo = x.upper_bound
x = bs.bootstrap(np.array(MeanBaseRatey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanbaseratey = x.value
lowerbaseratey = x.lower_bound
upperbaseratey = x.upper_bound
x = bs.bootstrap(np.array(MeanBaseRateo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanbaserateo = x.value
lowerbaserateo = x.lower_bound
upperbaserateo = x.upper_bound
x = bs.bootstrap(np.array(SNRy), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanSNRy = x.value
lowerSNRy = x.lower_bound
upperSNRy = x.upper_bound
x = bs.bootstrap(np.array(SNRo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanSNRo = x.value
lowerSNRo = x.lower_bound
upperSNRo = x.upper_bound
x = bs.bootstrap(np.array(PercentActivey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanactivey = x.value
loweractivey = x.lower_bound
upperactivey = x.upper_bound
x = bs.bootstrap(np.array(PercentActiveo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanactiveo = x.value
loweractiveo = x.lower_bound
upperactiveo = x.upper_bound
x = bs.bootstrap(np.array(RateCVy), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanCVy = x.value
lowerCVy = x.lower_bound
upperCVy = x.upper_bound
x = bs.bootstrap(np.array(RateCVo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanCVo = x.value
lowerCVo = x.lower_bound
upperCVo = x.upper_bound
else:
meanratey = np.mean(MeanRatey)
lowerratey = np.std(MeanRatey)
upperratey = np.std(MeanRatey)
meanrateo = np.mean(MeanRateo)
lowerrateo = np.std(MeanRateo)
upperrateo = np.std(MeanRateo)
meanbaseratey = np.mean(MeanBaseRatey)
lowerbaseratey = np.std(MeanBaseRatey)
upperbaseratey = np.std(MeanBaseRatey)
meanbaserateo = np.mean(MeanBaseRateo)
lowerbaserateo = np.std(MeanBaseRateo)
upperbaserateo = np.std(MeanBaseRateo)
meanSNRy = np.mean(SNRy)
lowerSNRy = np.std(SNRy)
upperSNRy = np.std(SNRy)
meanSNRo = np.mean(SNRo)
lowerSNRo = np.std(SNRo)
upperSNRo = np.std(SNRo)
meanactivey = np.mean(PercentActivey)
loweractivey = np.std(PercentActivey)
upperactivey = np.std(PercentActivey)
meanactiveo = np.mean(PercentActiveo)
loweractiveo = np.std(PercentActiveo)
upperactiveo = np.std(PercentActiveo)
meanCVy = np.mean(RateCVy)
lowerCVy = np.std(RateCVy)
upperCVy = np.std(RateCVy)
meanCVo = np.mean(RateCVo)
lowerCVo = np.std(RateCVo)
upperCVo = np.std(RateCVo)
df = pd.DataFrame(columns=["Cells", "Metric", "Mean Young", "SD Young", "Mean Old", "SD Old", "t-stat", "p-value", "Cohen's d"])
lay0 = abs(meanbaseratey-lowerbaseratey)
uay0 = abs(meanbaseratey-upperbaseratey)
lao0 = abs(meanbaserateo-lowerbaserateo)
uao0 = abs(meanbaserateo-upperbaserateo)
lay = abs(meanratey-lowerratey)
uay = abs(meanratey-upperratey)
lao = abs(meanrateo-lowerrateo)
uao = abs(meanrateo-upperrateo)
tstat_S0, pval_S0 = st.ttest_rel(MeanBaseRatey,MeanRatey)
tstat_S1, pval_S1 = st.ttest_rel(MeanBaseRateo,MeanRateo)
tstat_S2, pval_S2 = st.ttest_rel(MeanBaseRatey,MeanBaseRateo)
tstat_S3, pval_S3 = st.ttest_rel(MeanRatey,MeanRateo)
cd0 = cohen_d(MeanBaseRatey,MeanRatey)
cd1 = cohen_d(MeanBaseRateo,MeanRateo)
cd2 = cohen_d(MeanBaseRatey,MeanBaseRateo)
cd3 = cohen_d(MeanRatey,MeanRateo)
df2 = pd.DataFrame(columns=["Cells",
"Mean Base Young",
"SD Base Young",
"Mean Response Young",
"SD Response Young",
"Mean Base Old",
"SD Base Old",
"Mean Response Old",
"SD Response Old",
"t-stat (base-resp Y vs Y)",
"p-value (base-resp Y vs Y)",
"Cohen's d (base-resp Y vs Y)",
"t-stat (base-resp O vs O)",
"p-value (base-resp O vs O)",
"Cohen's d (base-resp O vs O)",
"t-stat (base-base Y vs O)",
"p-value (base-base Y vs O)",
"Cohen's d (base-base Y vs O)",
"t-stat (resp-resp Y vs O)",
"p-value (resp-resp Y vs O)",
"Cohen's d (resp-resp Y vs O)"])
df2 = df2.append({"Cells" : 'All',
"Mean Base Young" : meanbaseratey,
"SD Base Young" : [lay0,lay0] if bsMeans else lowerbaseratey,
"Mean Response Young" : meanratey,
"SD Response Young" : [lay,lay] if bsMeans else lowerratey,
"Mean Base Old" : meanbaserateo,
"SD Base Old" : [lao0,lao0] if bsMeans else lowerbaserateo,
"Mean Response Old" : meanrateo,
"SD Response Old" : [lao,lao] if bsMeans else lowerrateo,
"t-stat (base-resp Y vs Y)" : tstat_S0,
"p-value (base-resp Y vs Y)" : pval_S0,
"Cohen's d (base-resp Y vs Y)" : cd0,
"t-stat (base-resp O vs O)" : tstat_S1,
"p-value (base-resp O vs O)" : pval_S1,
"Cohen's d (base-resp O vs O)" : cd1,
"t-stat (base-base Y vs O)" : tstat_S2,
"p-value (base-base Y vs O)" : pval_S2,
"Cohen's d (base-base Y vs O)" : cd2,
"t-stat (resp-resp Y vs O)" : tstat_S3,
"p-value (resp-resp Y vs O)" : pval_S3,
"Cohen's d (resp-resp Y vs O)" : cd3},
ignore_index = True)
df2.to_csv('figs_tuning/stats_BaseVsResponseRates_'+str(totalstim)+'ms.csv')
x = [0.1,0.9,2.1,2.9]
fig = plt.figure(figsize=(7,5))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanbaseratey,meanbaserateo,meanratey,meanrateo], yerr=[[lay0,lao0],[uay0,uao0],[lay,lao],[uay,uao]] if bsMeans else [lowerbaseratey,lowerbaserateo,lowerratey,lowerrateo], color=['darkgray','lightcoral','darkgray','lightcoral'], width=0.8,tick_label=['Young','Old','Young','Old'],linewidth=1,ecolor='black',error_kw={'elinewidth':4,'markeredgewidth':4},capsize=12,edgecolor='k')
# if pval_S0 < 0.001:
# barplot_annotate_brackets(0, 1, '***', [x[i] for i in [0, 2]], [meanratey+2,meanratey+2], yerr=[uay0,uay] if bsMeans else [lowerbaseratey,lowerratey])
# elif pval_S0 < 0.01:
# barplot_annotate_brackets(0, 1, '**', [x[i] for i in [0, 2]], [meanratey+2,meanratey+2], yerr=[uay0,uay] if bsMeans else [lowerbaseratey,lowerratey])
# elif pval_S0 < 0.05:
# barplot_annotate_brackets(0, 1, '*', [x[i] for i in [0, 2]], [meanratey+2,meanratey+2], yerr=[uay0,uay] if bsMeans else [lowerbaseratey,lowerratey])
# if pval_S1 < 0.001:
# barplot_annotate_brackets(0, 1, '***', [x[i] for i in [1, 3]], [meanratey+3.5,meanratey+3.5], yerr=[uao0,uao] if bsMeans else [lowerbaserateo,lowerrateo])
# elif pval_S1 < 0.01:
# barplot_annotate_brackets(0, 1, '**', [x[i] for i in [1, 3]], [meanratey+3.5,meanratey+3.5], yerr=[uao0,uao] if bsMeans else [lowerbaserateo,lowerrateo])
# elif pval_S1 < 0.05:
# barplot_annotate_brackets(0, 1, '*', [x[i] for i in [1, 3]], [meanratey+3.5,meanratey+3.5], yerr=[uao0,uao] if bsMeans else [lowerbaserateo,lowerrateo])
if pval_S2 < 0.001:
barplot_annotate_brackets(0, 1, '***', [x[i] for i in [0, 1]], [meanbaseratey,meanbaserateo], yerr=[uay0,uao0] if bsMeans else [lowerbaseratey,lowerbaserateo])
elif pval_S2 < 0.01:
barplot_annotate_brackets(0, 1, '**', [x[i] for i in [0, 1]], [meanbaseratey,meanbaserateo], yerr=[uay0,uao0] if bsMeans else [lowerbaseratey,lowerbaserateo])
elif pval_S2 < 0.05:
barplot_annotate_brackets(0, 1, '*', [x[i] for i in [0, 1]], [meanbaseratey,meanbaserateo], yerr=[uay0,uao0] if bsMeans else [lowerbaseratey,lowerbaserateo])
if pval_S3 < 0.001:
barplot_annotate_brackets(0, 1, '***', [x[i] for i in [2, 3]], [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S3 < 0.01:
barplot_annotate_brackets(0, 1, '**', [x[i] for i in [2, 3]], [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S3 < 0.05:
barplot_annotate_brackets(0, 1, '*', [x[i] for i in [2, 3]], [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
ax1.set_ylabel('Mean Pyr Rate (Hz)')
ax1.set_xlim(-0.6,3.6)
ax1.set_xticks([0.5,2.5])
ax1.set_xticklabels(['Pre-stimulus','Post-stimulus'])
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/MeanBaseVsResponseRates_'+str(totalstim)+'ms_AllCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanactivey-loweractivey)
uay = abs(meanactivey-upperactivey)
lao = abs(meanactiveo-loweractiveo)
uao = abs(meanactiveo-upperactiveo)
tstat_S, pval_S = st.ttest_rel(PercentActivey,PercentActiveo)
cd = cohen_d(PercentActivey,PercentActiveo)
df = df.append({"Cells" : 'All',
"Metric" : '% Active',
"Mean Young" : meanactivey,
"SD Young" : [lay,lay] if bsMeans else loweractivey,
"Mean Old" : meanactiveo,
"SD Old" : [lao,lao] if bsMeans else loweractiveo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanactivey,meanactiveo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [loweractivey,loweractiveo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('% Silent p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
ax1.set_ylabel('Pyr % Active During Response')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/PercentActive_'+str(totalstim)+'ms_AllCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanratey-lowerratey)
uay = abs(meanratey-upperratey)
lao = abs(meanrateo-lowerrateo)
uao = abs(meanrateo-upperrateo)
tstat_S, pval_S = st.ttest_rel(MeanRatey,MeanRateo)
cd = cohen_d(MeanRatey,MeanRateo)
df = df.append({"Cells" : 'All',
"Metric" : 'Response Rate',
"Mean Young" : meanratey,
"SD Young" : [lay,lay] if bsMeans else lowerratey,
"Mean Old" : meanrateo,
"SD Old" : [lao,lao] if bsMeans else lowerrateo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanratey,meanrateo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerratey,lowerrateo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('Mean Rate p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
ax1.set_ylabel('Mean Pyr Response Rate (Hz)')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/MeanRate_'+str(totalstim)+'ms_AllCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanSNRy-lowerSNRy)
uay = abs(meanSNRy-upperSNRy)
lao = abs(meanSNRo-lowerSNRo)
uao = abs(meanSNRo-upperSNRo)
tstat_S, pval_S = st.ttest_rel(SNRy,SNRo)
cd = cohen_d(SNRy,SNRo)
df = df.append({"Cells" : 'All',
"Metric" : 'SNR',
"Mean Young" : meanSNRy,
"SD Young" : [lay,lay] if bsMeans else lowerSNRy,
"Mean Old" : meanSNRo,
"SD Old" : [lao,lao] if bsMeans else lowerSNRo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [-0.6,1.6]
fig = plt.figure(figsize=(4,5))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanSNRy,meanSNRo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerSNRy,lowerSNRo], color=['darkgray','lightcoral'], width=2.2,tick_label=['Younger','Older'],linewidth=1,ecolor='black',error_kw={'elinewidth':4,'markeredgewidth':4},capsize=12,edgecolor='k')
print('SNR p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', [x[i] for i in [0, 1]], [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', [x[i] for i in [0, 1]], [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', [x[i] for i in [0, 1]], [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
ax1.set_ylabel('Pyr SNR')
ax1.set_xlim(-2,3)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/SNR_'+str(totalstim)+'ms_AllCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanCVy-lowerCVy)
uay = abs(meanCVy-upperCVy)
lao = abs(meanCVo-lowerCVo)
uao = abs(meanCVo-upperCVo)
tstat_S, pval_S = st.ttest_rel(RateCVy,RateCVo)
cd = cohen_d(RateCVy,RateCVo)
df = df.append({"Cells" : 'All',
"Metric" : 'Response Rate CV',
"Mean Young" : meanCVy,
"SD Young" : [lay,lay] if bsMeans else lowerCVy,
"Mean Old" : meanCVo,
"SD Old" : [lao,lao] if bsMeans else lowerCVo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanCVy,meanCVo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerCVy,lowerCVo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('CV p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
ax1.set_ylabel('Response Rate CV')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/RateCV_'+str(totalstim)+'ms_AllCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
# Stimulated mean rate plots
MeanRatey = []
MeanRateo = []
PercentActivey = []
PercentActiveo = []
RateCVy = []
RateCVo = []
SNRy = []
SNRo = []
for i in range(0,len(rndinds)):
rvy = np.mean(Ratesy_allseeds1[i][:N_HL5PNso])
rvo = np.mean(Rateso_allseeds1[i][:N_HL5PNso])
MeanRatey.append(rvy)
MeanRateo.append(rvo)
rvy = np.mean(Ratesy_allseeds2[i][:N_HL5PNso])
rvo = np.mean(Rateso_allseeds2[i][:N_HL5PNso])
MeanRatey.append(rvy)
MeanRateo.append(rvo)
rvy = np.mean(SNRy_allseeds1[i][:N_HL5PNso])
rvo = np.mean(SNRo_allseeds1[i][:N_HL5PNso])
SNRy.append(rvy)
SNRo.append(rvo)
rvy = np.mean(SNRy_allseeds2[i][:N_HL5PNso])
rvo = np.mean(SNRo_allseeds2[i][:N_HL5PNso])
SNRy.append(rvy)
SNRo.append(rvo)
rvy = np.array(Ratesy_allseeds1[i][:N_HL5PNso])
rvo = np.array(Rateso_allseeds1[i][:N_HL5PNso])
PercentActivey.append((len(rvy[rvy>0])/N_HL5PNso)*100)
PercentActiveo.append((len(rvo[rvo>0])/N_HL5PNso)*100)
rvy = np.array(Ratesy_allseeds2[i][:N_HL5PNso])
rvo = np.array(Rateso_allseeds2[i][:N_HL5PNso])
PercentActivey.append((len(rvy[rvy>0])/N_HL5PNso)*100)
PercentActiveo.append((len(rvo[rvo>0])/N_HL5PNso)*100)
rvy = np.std(Ratesy_allseeds1[i][:N_HL5PNso])/np.mean(Ratesy_allseeds1[i][:N_HL5PNso])
rvo = np.std(Rateso_allseeds1[i][:N_HL5PNso])/np.mean(Rateso_allseeds1[i][:N_HL5PNso])
RateCVy.append(rvy)
RateCVo.append(rvo)
rvy = np.std(Ratesy_allseeds2[i][:N_HL5PNso])/np.mean(Ratesy_allseeds2[i][:N_HL5PNso])
rvo = np.std(Rateso_allseeds2[i][:N_HL5PNso])/np.mean(Rateso_allseeds2[i][:N_HL5PNso])
RateCVy.append(rvy)
RateCVo.append(rvo)
# Use this code for bootstrapping instead
if bsMeans:
x = bs.bootstrap(np.array(MeanRatey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey = x.value
lowerratey = x.lower_bound
upperratey = x.upper_bound
x = bs.bootstrap(np.array(MeanRateo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo = x.value
lowerrateo = x.lower_bound
upperrateo = x.upper_bound
x = bs.bootstrap(np.array(SNRy), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanSNRy = x.value
lowerSNRy = x.lower_bound
upperSNRy = x.upper_bound
x = bs.bootstrap(np.array(SNRo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanSNRo = x.value
lowerSNRo = x.lower_bound
upperSNRo = x.upper_bound
x = bs.bootstrap(np.array(PercentActivey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanactivey = x.value
loweractivey = x.lower_bound
upperactivey = x.upper_bound
x = bs.bootstrap(np.array(PercentActiveo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanactiveo = x.value
loweractiveo = x.lower_bound
upperactiveo = x.upper_bound
x = bs.bootstrap(np.array(RateCVy), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanCVy = x.value
lowerCVy = x.lower_bound
upperCVy = x.upper_bound
x = bs.bootstrap(np.array(RateCVo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanCVo = x.value
lowerCVo = x.lower_bound
upperCVo = x.upper_bound
else:
meanratey = np.mean(MeanRatey)
lowerratey = np.std(MeanRatey)
upperratey = np.std(MeanRatey)
meanrateo = np.mean(MeanRateo)
lowerrateo = np.std(MeanRateo)
upperrateo = np.std(MeanRateo)
meanSNRy = np.mean(SNRy)
lowerSNRy = np.std(SNRy)
upperSNRy = np.std(SNRy)
meanSNRo = np.mean(SNRo)
lowerSNRo = np.std(SNRo)
upperSNRo = np.std(SNRo)
meanactivey = np.mean(PercentActivey)
loweractivey = np.std(PercentActivey)
upperactivey = np.std(PercentActivey)
meanactiveo = np.mean(PercentActiveo)
loweractiveo = np.std(PercentActiveo)
upperactiveo = np.std(PercentActiveo)
meanCVy = np.mean(RateCVy)
lowerCVy = np.std(RateCVy)
upperCVy = np.std(RateCVy)
meanCVo = np.mean(RateCVo)
lowerCVo = np.std(RateCVo)
upperCVo = np.std(RateCVo)
lay = abs(meanactivey-loweractivey)
uay = abs(meanactivey-upperactivey)
lao = abs(meanactiveo-loweractiveo)
uao = abs(meanactiveo-upperactiveo)
tstat_S, pval_S = st.ttest_rel(PercentActivey,PercentActiveo)
cd = cohen_d(PercentActivey,PercentActiveo)
df = df.append({"Cells" : 'Stimulated',
"Metric" : '% Active',
"Mean Young" : meanactivey,
"SD Young" : [lay,lay] if bsMeans else loweractivey,
"Mean Old" : meanactiveo,
"SD Old" : [lao,lao] if bsMeans else loweractiveo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanactivey,meanactiveo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [loweractivey,loweractiveo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('% Silent p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
ax1.set_ylabel('Stimulated Pyr % Active During Response')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/PercentActive_'+str(totalstim)+'ms_StimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanratey-lowerratey)
uay = abs(meanratey-upperratey)
lao = abs(meanrateo-lowerrateo)
uao = abs(meanrateo-upperrateo)
tstat_S, pval_S = st.ttest_rel(MeanRatey,MeanRateo)
cd = cohen_d(MeanRatey,MeanRateo)
df = df.append({"Cells" : 'Stimulated',
"Metric" : 'Response Rate',
"Mean Young" : meanratey,
"SD Young" : [lay,lay] if bsMeans else lowerratey,
"Mean Old" : meanrateo,
"SD Old" : [lao,lao] if bsMeans else lowerrateo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanratey,meanrateo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerratey,lowerrateo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('Mean Rate p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
ax1.set_ylabel('Mean Stimulated Pyr Response Rate (Hz)')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/MeanRate_'+str(totalstim)+'ms_StimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanSNRy-lowerSNRy)
uay = abs(meanSNRy-upperSNRy)
lao = abs(meanSNRo-lowerSNRo)
uao = abs(meanSNRo-upperSNRo)
tstat_S, pval_S = st.ttest_rel(SNRy,SNRo)
cd = cohen_d(SNRy,SNRo)
df = df.append({"Cells" : 'Stimulated',
"Metric" : 'SNR',
"Mean Young" : meanSNRy,
"SD Young" : [lay,lay] if bsMeans else lowerSNRy,
"Mean Old" : meanSNRo,
"SD Old" : [lao,lao] if bsMeans else lowerSNRo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanSNRy,meanSNRo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerSNRy,lowerSNRo], color=['dimgray','r'], width=0.8,tick_label=['Younger','Older'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('SNR p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
ax1.set_ylabel('Stimulated Pyr SNR')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/SNR_'+str(totalstim)+'ms_StimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanCVy-lowerCVy)
uay = abs(meanCVy-upperCVy)
lao = abs(meanCVo-lowerCVo)
uao = abs(meanCVo-upperCVo)
tstat_S, pval_S = st.ttest_rel(RateCVy,RateCVo)
cd = cohen_d(RateCVy,RateCVo)
df = df.append({"Cells" : 'Stimulated',
"Metric" : 'Response Rate CV',
"Mean Young" : meanCVy,
"SD Young" : [lay,lay] if bsMeans else lowerCVy,
"Mean Old" : meanCVo,
"SD Old" : [lao,lao] if bsMeans else lowerCVo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanCVy,meanCVo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerCVy,lowerCVo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('CV p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
ax1.set_ylabel('Stimulated Pyr Response Rate CV')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/RateCV_'+str(totalstim)+'ms_StimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
# Untimulated mean rate plots
MeanRatey = []
MeanRateo = []
PercentActivey = []
PercentActiveo = []
RateCVy = []
RateCVo = []
SNRy = []
SNRo = []
for i in range(0,len(rndinds)):
rvy = np.mean(Ratesy_allseeds1[i][N_HL5PNso:])
rvo = np.mean(Rateso_allseeds1[i][N_HL5PNso:])
MeanRatey.append(rvy)
MeanRateo.append(rvo)
rvy = np.mean(Ratesy_allseeds2[i][N_HL5PNso:])
rvo = np.mean(Rateso_allseeds2[i][N_HL5PNso:])
MeanRatey.append(rvy)
MeanRateo.append(rvo)
rvy = np.mean(SNRy_allseeds1[i][N_HL5PNso:])
rvo = np.mean(SNRo_allseeds1[i][N_HL5PNso:])
SNRy.append(rvy)
SNRo.append(rvo)
rvy = np.mean(SNRy_allseeds2[i][N_HL5PNso:])
rvo = np.mean(SNRo_allseeds2[i][N_HL5PNso:])
SNRy.append(rvy)
SNRo.append(rvo)
rvy = np.array(Ratesy_allseeds1[i][N_HL5PNso:])
rvo = np.array(Rateso_allseeds1[i][N_HL5PNso:])
PercentActivey.append((len(rvy[rvy>0])/(N_HL5PN-N_HL5PNso))*100)
PercentActiveo.append((len(rvo[rvo>0])/(N_HL5PN-N_HL5PNso))*100)
rvy = np.array(Ratesy_allseeds2[i][N_HL5PNso:])
rvo = np.array(Rateso_allseeds2[i][N_HL5PNso:])
PercentActivey.append((len(rvy[rvy>0])/(N_HL5PN-N_HL5PNso))*100)
PercentActiveo.append((len(rvo[rvo>0])/(N_HL5PN-N_HL5PNso))*100)
rvy = np.std(Ratesy_allseeds1[i][N_HL5PNso:])/np.mean(Ratesy_allseeds1[i][N_HL5PNso:])
rvo = np.std(Rateso_allseeds1[i][N_HL5PNso:])/np.mean(Rateso_allseeds1[i][N_HL5PNso:])
RateCVy.append(rvy)
RateCVo.append(rvo)
rvy = np.std(Ratesy_allseeds2[i][N_HL5PNso:])/np.mean(Ratesy_allseeds2[i][N_HL5PNso:])
rvo = np.std(Rateso_allseeds2[i][N_HL5PNso:])/np.mean(Rateso_allseeds2[i][N_HL5PNso:])
RateCVy.append(rvy)
RateCVo.append(rvo)
# Use this code for bootstrapping instead
if bsMeans:
x = bs.bootstrap(np.array(MeanRatey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey = x.value
lowerratey = x.lower_bound
upperratey = x.upper_bound
x = bs.bootstrap(np.array(MeanRateo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo = x.value
lowerrateo = x.lower_bound
upperrateo = x.upper_bound
x = bs.bootstrap(np.array(SNRy), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanSNRy = x.value
lowerSNRy = x.lower_bound
upperSNRy = x.upper_bound
x = bs.bootstrap(np.array(SNRo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanSNRo = x.value
lowerSNRo = x.lower_bound
upperSNRo = x.upper_bound
x = bs.bootstrap(np.array(PercentActivey), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanactivey = x.value
loweractivey = x.lower_bound
upperactivey = x.upper_bound
x = bs.bootstrap(np.array(PercentActiveo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanactiveo = x.value
loweractiveo = x.lower_bound
upperactiveo = x.upper_bound
x = bs.bootstrap(np.array(RateCVy), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanCVy = x.value
lowerCVy = x.lower_bound
upperCVy = x.upper_bound
x = bs.bootstrap(np.array(RateCVo), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanCVo = x.value
lowerCVo = x.lower_bound
upperCVo = x.upper_bound
else:
meanratey = np.mean(MeanRatey)
lowerratey = np.std(MeanRatey)
upperratey = np.std(MeanRatey)
meanrateo = np.mean(MeanRateo)
lowerrateo = np.std(MeanRateo)
upperrateo = np.std(MeanRateo)
meanSNRy = np.mean(SNRy)
lowerSNRy = np.std(SNRy)
upperSNRy = np.std(SNRy)
meanSNRo = np.mean(SNRo)
lowerSNRo = np.std(SNRo)
upperSNRo = np.std(SNRo)
meanactivey = np.mean(PercentActivey)
loweractivey = np.std(PercentActivey)
upperactivey = np.std(PercentActivey)
meanactiveo = np.mean(PercentActiveo)
loweractiveo = np.std(PercentActiveo)
upperactiveo = np.std(PercentActiveo)
meanCVy = np.mean(RateCVy)
lowerCVy = np.std(RateCVy)
upperCVy = np.std(RateCVy)
meanCVo = np.mean(RateCVo)
lowerCVo = np.std(RateCVo)
upperCVo = np.std(RateCVo)
lay = abs(meanactivey-loweractivey)
uay = abs(meanactivey-upperactivey)
lao = abs(meanactiveo-loweractiveo)
uao = abs(meanactiveo-upperactiveo)
tstat_S, pval_S = st.ttest_rel(PercentActivey,PercentActiveo)
cd = cohen_d(PercentActivey,PercentActiveo)
df = df.append({"Cells" : 'Recurrent',
"Metric" : '% Active',
"Mean Young" : meanactivey,
"SD Young" : [lay,lay] if bsMeans else loweractivey,
"Mean Old" : meanactiveo,
"SD Old" : [lao,lao] if bsMeans else loweractiveo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanactivey,meanactiveo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [loweractivey,loweractiveo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('% Silent p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanactivey,meanactiveo], yerr=[uay,uao]if bsMeans else [loweractivey,loweractiveo])
ax1.set_ylabel('Recurrent Pyr % Active During Response')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/PercentActive_'+str(totalstim)+'ms_UnstimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanratey-lowerratey)
uay = abs(meanratey-upperratey)
lao = abs(meanrateo-lowerrateo)
uao = abs(meanrateo-upperrateo)
tstat_S, pval_S = st.ttest_rel(MeanRatey,MeanRateo)
cd = cohen_d(MeanRatey,MeanRateo)
df = df.append({"Cells" : 'Recurrent',
"Metric" : 'Response Rate',
"Mean Young" : meanratey,
"SD Young" : [lay,lay] if bsMeans else lowerratey,
"Mean Old" : meanrateo,
"SD Old" : [lao,lao] if bsMeans else lowerrateo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanratey,meanrateo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerratey,lowerrateo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('Mean Rate p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanratey,meanrateo], yerr=[uay,uao] if bsMeans else [lowerratey,lowerrateo])
ax1.set_ylabel('Mean Recurrent Pyr Response Rate (Hz)')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/MeanRate_'+str(totalstim)+'ms_UnstimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanSNRy-lowerSNRy)
uay = abs(meanSNRy-upperSNRy)
lao = abs(meanSNRo-lowerSNRo)
uao = abs(meanSNRo-upperSNRo)
tstat_S, pval_S = st.ttest_rel(SNRy,SNRo)
cd = cohen_d(SNRy,SNRo)
df = df.append({"Cells" : 'Recurrent',
"Metric" : 'SNR',
"Mean Young" : meanSNRy,
"SD Young" : [lay,lay] if bsMeans else lowerSNRy,
"Mean Old" : meanSNRo,
"SD Old" : [lao,lao] if bsMeans else lowerSNRo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanSNRy,meanSNRo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerSNRy,lowerSNRo], color=['dimgray','r'], width=0.8,tick_label=['Younger','Older'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('SNR p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanSNRy,meanSNRo], yerr=[uay,uao] if bsMeans else [lowerSNRy,lowerSNRo])
ax1.set_ylabel('Recurrent Pyr SNR')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/SNR_'+str(totalstim)+'ms_UnstimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
lay = abs(meanCVy-lowerCVy)
uay = abs(meanCVy-upperCVy)
lao = abs(meanCVo-lowerCVo)
uao = abs(meanCVo-upperCVo)
tstat_S, pval_S = st.ttest_rel(RateCVy,RateCVo)
cd = cohen_d(RateCVy,RateCVo)
df = df.append({"Cells" : 'Recurrent',
"Metric" : 'Response Rate CV',
"Mean Young" : meanCVy,
"SD Young" : [lay,lay] if bsMeans else lowerCVy,
"Mean Old" : meanCVo,
"SD Old" : [lao,lao] if bsMeans else lowerCVo,
"t-stat" : tstat_S,
"p-value" : pval_S,
"Cohen's d" : cd},
ignore_index = True)
x = [0,1]
fig = plt.figure(figsize=(6,9))
ax1 = fig.add_subplot(111)
ax1.bar(x, [meanCVy,meanCVo], yerr=[[lay,lao],[uay,uao]] if bsMeans else [lowerCVy,lowerCVo], color=['dimgray','r'], width=0.8,tick_label=['Young','Old'],linewidth=5,ecolor='black',error_kw={'elinewidth':3,'markeredgewidth':3},capsize=12,edgecolor='k')
print('CV p-value='+str(pval_S))
if pval_S < 0.001:
barplot_annotate_brackets(0, 1, '***', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
elif pval_S < 0.01:
barplot_annotate_brackets(0, 1, '**', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
elif pval_S < 0.05:
barplot_annotate_brackets(0, 1, '*', x, [meanCVy,meanCVo], yerr=[uay,uao] if bsMeans else [lowerCVy,lowerCVo])
ax1.set_ylabel('Unstimulated Pyr Response Rate CV')
ax1.set_xlim(-0.6,1.6)
ax1.spines["right"].set_visible(False)
ax1.spines["top"].set_visible(False)
fig.tight_layout()
fig.savefig('figs_tuning/RateCV_'+str(totalstim)+'ms_UnstimulatedCellsDuringStim_MergedOrientations.png',bbox_inches='tight', dpi=300, transparent=True)
plt.close(fig)
# Save stats dataframe to csv
df.to_csv('figs_tuning/stats_'+str(totalstim)+'ms.csv')
# Stimulated mean rates
MeanRatey1 = []
MeanRateo1 = []
MeanRatey2 = []
MeanRateo2 = []
for i in range(0,len(rndinds)):
rvy = np.mean(Ratesy_allseeds1[i][:N_HL5PNso])
rvo = np.mean(Rateso_allseeds1[i][:N_HL5PNso])
MeanRatey1.append(rvy)
MeanRateo1.append(rvo)
rvy = np.mean(Ratesy_allseeds2[i][:N_HL5PNso])
rvo = np.mean(Rateso_allseeds2[i][:N_HL5PNso])
MeanRatey2.append(rvy)
MeanRateo2.append(rvo)
x = bs.bootstrap(np.array(MeanRatey1), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey1 = x.value
x = bs.bootstrap(np.array(MeanRateo1), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo1 = x.value
x = bs.bootstrap(np.array(MeanRatey2), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey2 = x.value
x = bs.bootstrap(np.array(MeanRateo2), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo2 = x.value
mrys1 = meanratey1
mros1 = meanrateo1
mrys2 = meanratey2
mros2 = meanrateo2
# Unstimulated mean rates
MeanRatey1 = []
MeanRateo1 = []
MeanRatey2 = []
MeanRateo2 = []
for i in range(0,len(rndinds)):
rvy = np.mean(Ratesy_allseeds1[i][N_HL5PNso:])
rvo = np.mean(Rateso_allseeds1[i][N_HL5PNso:])
MeanRatey1.append(rvy)
MeanRateo1.append(rvo)
rvy = np.mean(Ratesy_allseeds2[i][N_HL5PNso:])
rvo = np.mean(Rateso_allseeds2[i][N_HL5PNso:])
MeanRatey2.append(rvy)
MeanRateo2.append(rvo)
x = bs.bootstrap(np.array(MeanRatey1), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey1 = x.value
x = bs.bootstrap(np.array(MeanRateo1), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo1 = x.value
x = bs.bootstrap(np.array(MeanRatey2), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanratey2 = x.value
x = bs.bootstrap(np.array(MeanRateo2), stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
meanrateo2 = x.value
mryu1 = meanratey1
mrou1 = meanrateo1
mryu2 = meanratey2
mrou2 = meanrateo2
# Plot tuning curves
CI_means_PNy1 = []
CI_lower_PNy1 = []
CI_upper_PNy1 = []
CI_means_PNo1 = []
CI_lower_PNo1 = []
CI_upper_PNo1 = []
CI_means_PNy2 = []
CI_lower_PNy2 = []
CI_upper_PNy2 = []
CI_means_PNo2 = []
CI_lower_PNo2 = []
CI_upper_PNo2 = []
bsMeans = False
if bsMeans:
for l in range(0,N_HL5PN):
print('Starting Bootstrapping')
x = bs.bootstrap(np.transpose(Ratesy_allseeds1)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNy1.append(x.value)
CI_lower_PNy1.append(x.lower_bound)
CI_upper_PNy1.append(x.upper_bound)
x = bs.bootstrap(np.transpose(Rateso_allseeds1)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNo1.append(x.value)
CI_lower_PNo1.append(x.lower_bound)
CI_upper_PNo1.append(x.upper_bound)
x = bs.bootstrap(np.transpose(Ratesy_allseeds2)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNy2.append(x.value)
CI_lower_PNy2.append(x.lower_bound)
CI_upper_PNy2.append(x.upper_bound)
x = bs.bootstrap(np.transpose(Rateso_allseeds2)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNo2.append(x.value)
CI_lower_PNo2.append(x.lower_bound)
CI_upper_PNo2.append(x.upper_bound)
else:
for l in range(0,N_HL5PN):
CI_means_PNy1.append(np.mean(np.transpose(Ratesy_allseeds1)[l]))
CI_lower_PNy1.append(np.mean(np.transpose(Ratesy_allseeds1)[l])-np.std(np.transpose(Ratesy_allseeds1)[l]))
CI_upper_PNy1.append(np.mean(np.transpose(Ratesy_allseeds1)[l])+np.std(np.transpose(Ratesy_allseeds1)[l]))
x = bs.bootstrap(np.transpose(Rateso_allseeds1)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNo1.append(np.mean(np.transpose(Rateso_allseeds1)[l]))
CI_lower_PNo1.append(np.mean(np.transpose(Rateso_allseeds1)[l])-np.std(np.transpose(Rateso_allseeds1)[l]))
CI_upper_PNo1.append(np.mean(np.transpose(Rateso_allseeds1)[l])+np.std(np.transpose(Rateso_allseeds1)[l]))
x = bs.bootstrap(np.transpose(Ratesy_allseeds2)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNy2.append(np.mean(np.transpose(Ratesy_allseeds2)[l]))
CI_lower_PNy2.append(np.mean(np.transpose(Ratesy_allseeds2)[l])-np.std(np.transpose(Ratesy_allseeds2)[l]))
CI_upper_PNy2.append(np.mean(np.transpose(Ratesy_allseeds2)[l])+np.std(np.transpose(Ratesy_allseeds2)[l]))
x = bs.bootstrap(np.transpose(Rateso_allseeds2)[l], stat_func=bs_stats.mean, alpha=0.05, num_iterations=500)
CI_means_PNo2.append(np.mean(np.transpose(Rateso_allseeds2)[l]))
CI_lower_PNo2.append(np.mean(np.transpose(Rateso_allseeds2)[l])-np.std(np.transpose(Rateso_allseeds2)[l]))
CI_upper_PNo2.append(np.mean(np.transpose(Rateso_allseeds2)[l])+np.std(np.transpose(Rateso_allseeds2)[l]))
gscalery1 = np.percentile(CI_means_PNy1,90)-mryu1
gscalero1 = np.percentile(CI_means_PNo1,90)-mrou1
mu_ind1 = int(N_HL5PNso*synnums[0]/180) # find index with peak input depending on orientation
gscalery2 = np.percentile(CI_means_PNy2,90)-mryu2
gscalero2 = np.percentile(CI_means_PNo2,90)-mrou2
mu_ind2 = int(N_HL5PNso*synnums[1]/180) # find index with peak input depending on orientation
sigma = int(N_HL5PNso*42/180) # Set sigma to 42 degrees (Bensmaia et al 2008) and scale by number of PN neurons
x_values = np.linspace(0, N_HL5PN-1, N_HL5PN)
gfuncy1 = gaussian(x_values, mu_ind1, sigma)*gscalery1+mryu1 # Gaussian vector (y=[0,1]) of length = cellnum
gfuncy1[x_values>N_HL5PNso]=mryu1
gfunco1 = gaussian(x_values, mu_ind1, sigma)*gscalero1+mrou1 # Gaussian vector (y=[0,1]) of length = cellnum
gfunco1[x_values>N_HL5PNso]=mrou1
gfuncy2 = gaussian(x_values, mu_ind2, sigma)*gscalery2+mryu2 # Gaussian vector (y=[0,1]) of length = cellnum
gfuncy2[x_values>N_HL5PNso]=mryu2
gfunco2 = gaussian(x_values, mu_ind2, sigma)*gscalero2+mrou2 # Gaussian vector (y=[0,1]) of length = cellnum
gfunco2[x_values>N_HL5PNso]=mrou2
# Plot Heatmap of Tuning for each Cell
x_values = np.linspace(0, N_HL5PN-1, N_HL5PN)
y_values = np.linspace(0, 180, N_HL5PNso)
sigma = int(N_HL5PNso*42/180) # Set sigma to 42 degrees (Bensmaia et al 2008) and scale by number of PN neurons
all_gfs = []
for y in y_values:
yin = int(N_HL5PNso*y/180)
gf = gaussian(x_values,yin,sigma)*45
gf[x_values>=N_HL5PNso]=0
all_gfs.append(gf)
fig, axarr = plt.subplots(nrows=1,ncols=1,figsize=(11,3.2))
im = axarr.imshow(all_gfs,origin='lower',extent = [0 , N_HL5PN-1, 0 , 180], aspect='auto')
axarr.set_ylabel('Angle ($^\circ$)')
axarr.set_xlabel('Neuron Index')
axarr.set_xlim(0,N_HL5PN)
axarr.set_xticks([0,100,200,300,400,500,600,700])
cbar = fig.colorbar(im)
cbar.set_label(r'$N_{inputs}$')
fig.tight_layout()
fig.savefig('figs_tuning/tuning_schematic.png',dpi=300,transparent=True)
xcells = np.arange(0,N_HL5PN,1)
fig, axarr = plt.subplots(nrows=2,ncols=1,figsize=(9,7),sharex=True,sharey=True)
fig.text(0.02, 0.5, 'Pyr Response Rate (Hz)', ha='center', va='center', rotation='vertical')
# Add background shading to highlight feature subsets
N_PNs = 700
N_PNs_stim = int(N_PNs/2)
NFeats = 100
SubsetFeatures1 = [0,N_PNs] # [Logical, Cell Start Index, Cell End Index] - Only triggers if logical = True
SubsetFeatures2 = [0,N_PNs_stim]
NFeats = 200
SubsetFeatures3 = [int(N_PNs_stim/2)-int(NFeats/2),int(N_PNs_stim/2)+int(NFeats/2)]
NFeats = 100
SubsetFeatures4 = [int(N_PNs_stim/2)-int(NFeats/2),int(N_PNs_stim/2)+int(NFeats/2)]
ymax = np.max([np.max(CI_upper_PNy1),np.max(CI_upper_PNy2),np.max(CI_upper_PNo1),np.max(CI_upper_PNo2)])
lb = [0,0]
ub = [ymax,ymax]
# axarr[0].fill_between(SubsetFeatures1,lb,ub,facecolor='k',alpha=0.1)
# axarr[0].fill_between(SubsetFeatures2,lb,ub,facecolor='k',alpha=0.2)
# axarr[0].fill_between(SubsetFeatures3,lb,ub,facecolor='k',alpha=0.3)
# axarr[0].fill_between(SubsetFeatures4,lb,ub,facecolor='k',alpha=0.4)
# axarr[1].fill_between(SubsetFeatures1,lb,ub,facecolor='k',alpha=0.1)
# axarr[1].fill_between(SubsetFeatures2,lb,ub,facecolor='k',alpha=0.2)
# axarr[1].fill_between(SubsetFeatures3,lb,ub,facecolor='k',alpha=0.3)
# axarr[1].fill_between(SubsetFeatures4,lb,ub,facecolor='k',alpha=0.4)
axarr[0].plot(xcells,CI_means_PNy1,color='k',alpha=1,label=str(synnums[0])+r'$\degree$', linewidth=0.5)
axarr[0].fill_between(xcells,CI_lower_PNy1,CI_upper_PNy1,facecolor='k',alpha=0.5)
axarr[0].plot(xcells,CI_means_PNy2,color='tab:brown',alpha=1,label=str(synnums[1])+r'$\degree$', linewidth=0.5)
axarr[0].fill_between(xcells,CI_lower_PNy2,CI_upper_PNy2,facecolor='tab:brown',alpha=0.3)
axarr[0].legend()
axarr[1].plot(xcells,CI_means_PNo1,color='r',alpha=1,label=str(synnums[0])+r'$\degree$', linewidth=0.5)
axarr[1].fill_between(xcells,CI_lower_PNo1,CI_upper_PNo1,facecolor='r',alpha=0.4)
axarr[1].plot(xcells,CI_means_PNo2,color='tab:orange',alpha=1,label=str(synnums[1])+r'$\degree$', linewidth=0.5)
axarr[1].fill_between(xcells,CI_lower_PNo2,CI_upper_PNo2,facecolor='tab:orange',alpha=0.3)
axarr[1].legend()
# axarr[0].plot(x_values,gfuncy1,color='k',alpha=1,lw=3,ls=':')
# axarr[0].plot(x_values,gfuncy2,color='tab:brown',alpha=1,lw=3,ls=':')
# axarr[1].plot(x_values,gfunco1,color='r',alpha=1,lw=3,ls=':')
# axarr[1].plot(x_values,gfunco2,color='tab:orange',alpha=1,lw=3,ls=':')
axarr[0].set_xlim(0,N_HL5PN)
axarr[0].set_ylim(0,ymax)
ylims=axarr[0].get_ylim()
axarr[0].plot(np.array([N_HL5PNso,N_HL5PNso]),ylims,color='dimgray',alpha=1,lw=2,ls='dashed')
axarr[1].plot(np.array([N_HL5PNso,N_HL5PNso]),ylims,color='dimgray',alpha=1,lw=2,ls='dashed')
axarr[0].set_ylim(0,ymax)
# axarr[1].set_ylabel('Spike Rate (Hz)')
secax = axarr[0].secondary_xaxis('top')
xpos = np.arange(0,N_HL5PNso+1,N_HL5PNso/8)
rot = 180/8
secax.set_xticks(xpos)
secax.set_xticklabels([u'\u2014']*len(xpos),fontsize=35,verticalalignment='bottom')
for i, t in enumerate(secax.get_xticklabels()):
t.set_rotation(i*rot)
axarr[1].set_xlim(0,N_HL5PN)
axarr[1].set_xlabel('Neuron Index')
fig.tight_layout()
fig.savefig('figs_tuning/tuningcurve_'+str(totalstim)+'ms_allcells_Compare.png',dpi=300,transparent=True)
xcells = np.arange(0,N_HL5PN,1)
fig, axarr = plt.subplots(nrows=1,ncols=1,figsize=(10,5),sharex=True,sharey=True)
# Add background shading to highlight feature subsets
ymax = np.max(CI_upper_PNy1)
lb = [0,0]
ub = [ymax,ymax]
axarr.scatter(xcells,CI_means_PNy1,c='k',alpha=0.25,label=str(synnums[0])+r'$\degree$', linewidths=1, edgecolors='none')
# axarr.plot(x_values,gfuncy1,color='k',alpha=1,lw=3,ls=':')
# axarr.fill_between(xcells,CI_lower_PNy1,CI_upper_PNy1,facecolor='k',alpha=0.5)
axarr.set_xlim(0,N_HL5PN)
axarr.set_ylim(0,ymax)
ylims=axarr.get_ylim()
axarr.plot(np.array([N_HL5PNso,N_HL5PNso]),ylims,color='dimgray',alpha=1,lw=2,ls='dashed')
axarr.set_ylim(0,ymax)
axarr.set_ylabel('Pyr Response Rate (Hz)')
# axins = inset_axes(axarr, width="40%", height="40%",loc=1)
# axins.plot(x_values,gfuncy1,color='k',alpha=1,lw=3,ls=':')
# axins.fill_between(xcells[0:N_HL5PNso],CI_lower_PNy1[0:N_HL5PNso],CI_upper_PNy1[0:N_HL5PNso],facecolor='k',alpha=0.5)
# axins.set_ylim(0,ymax)
# secax = axarr.secondary_xaxis('top')
xpos = np.arange(0,N_HL5PNso+1,N_HL5PNso/8)
rot = 180/8
# secax.set_xticks(xpos)
# secax.set_xticklabels([u'\u2014']*len(xpos),fontsize=35,verticalalignment='bottom')
# for i, t in enumerate(secax.get_xticklabels()):
# t.set_rotation(i*rot)
# axins.set_xticks(xpos)
# axins.set_xticklabels([u'\u2014']*len(xpos),fontsize=35)
# for i, t in enumerate(axins.get_xticklabels()):
# t.set_rotation(i*rot)
# axins.set_xlim(5,N_HL5PNso-5)
axarr.set_xlim(0,N_HL5PN)
axarr.set_xlabel('Neuron Index')
fig.tight_layout()
fig.savefig('figs_tuning/tuningcurve_scatter_YoungOnlyA1_'+str(totalstim)+'ms_allcells_Compare.png',dpi=300,transparent=True)
# Cross-correlations
xpos = np.arange(-N_HL5PN+1,N_HL5PN,1)
AllMeans = [CI_means_PNy1,CI_means_PNy2,CI_means_PNo1,CI_means_PNo2]
labels = ['YA1','YA2','OA1','OA2']
fig, axarr = plt.subplots(nrows=len(AllMeans),ncols=len(AllMeans),figsize=(14,14),sharex=True,sharey=True)
i = 0
for a in AllMeans:
j = 0
for b in AllMeans:
# Normalize
a = (a - np.mean(a)) / (np.std(a) * len(a))
b = (b - np.mean(b)) / (np.std(b))
c = np.correlate(a, b, 'full')
axarr[i][j].plot(xpos,c,color='k')
axarr[i][j].set_xlim(xpos[0],xpos[-1])
if i == 0 : axarr[i][j].set_title(labels[j],fontsize=25)
if j == 0 : axarr[i][j].set_ylabel(labels[i],fontsize=25)
j += 1
i += 1
fig.tight_layout()
fig.savefig('figs_tuning/tcurvexcorrs_'+str(totalstim)+'ms_allcells_Compare.png',dpi=300,transparent=True)
axarr[0][0].set_xlim(-100,100)
fig.savefig('figs_tuning/tcurvexcorrs_'+str(totalstim)+'ms_allcells_Compare_Zoomed.png',dpi=300,transparent=True)