Hippocampus CA1 Interneuron Specific 3 (IS3) in vivo-like virtual NN simulations (Luo et al 2020)

 Download zip file 
Help downloading and running models
Accession:265523
"Disinhibition is a widespread circuit mechanism for information selection and transfer. In the hippocampus, disinhibition of principal cells is provided by the interneuron-specific interneurons that express the vasoactive intestinal polypeptide (VIP-IS) and innervate selectively inhibitory interneurons. By combining optophysiological experiments with computational models, we determined the impact of synaptic inputs onto the network state-dependent recruitment of VIP-IS cells. We found that VIP-IS cells fire spikes in response to both the Schaffer collateral and the temporoammonic pathway activation. Moreover, by integrating their intrinsic and synaptic properties into computational models, we predicted recruitment of these cells between the rising phase and peak of theta oscillation and during ripples. Two-photon Ca2+-imaging in awake mice supported in part the theoretical predictions, revealing a significant speed modulation of VIP-IS cells and their preferential albeit delayed recruitment during theta-run epochs, with estimated firing at the rising phase and peak of the theta cycle. However, it also uncovered that VIP-IS cells are not activated during ripples. Thus, given the preferential theta-modulated firing of VIP-IS cells in awake hippocampus, we postulate that these cells may be important for information gating during spatial navigation and memory encoding."
Reference:
1 . Luo X, Guet-McCreight A, Villette V, Francavilla R, Marino B, Chamberland S, Skinner FK, Topolnik L (2020) Synaptic Mechanisms Underlying the Network State-Dependent Recruitment of VIP-Expressing Interneurons in the CA1 Hippocampus. Cereb Cortex [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Synapse; Dendrite; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampal CA1 CR/VIP cell;
Channel(s): I Na,t; I Na,p; I A;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s): Glutamate; Gaba;
Simulation Environment: NEURON;
Model Concept(s): Spatial Navigation; Oscillations; Activity Patterns;
Implementer(s): Guet-McCreight, Alexandre [alexandre.guet.mccreight at mail.utoronto.ca];
Search NeuronDB for information about:  I Na,p; I Na,t; I A; Gaba; Glutamate;
/
LuoEtAl2020Code
Theta
SDprox2
NPYfiles
PLOTfiles
IKa.mod *
ingauss.mod *
Ksoma.mod *
Nap.mod *
Nasoma.mod *
vecevent.mod *
init.py *
IS3_M2_Case9StarRevised.hoc *
model_decaytimeinhvec.dat *
model_decaytimevec.dat *
model_dendsectionvec.dat *
model_distvec.dat *
model_minweightinhvec.dat *
model_minweightvec.dat *
model_risetimeinhvec.dat *
model_risetimevec.dat *
PlotResults.py
ranstream.hoc *
SynParamSearch.hoc
vecevent.hoc *
vecevent.ses *
                            
### Test Script for a file found in the SDprox2 results from initial Parallel Simulations
from __future__ import division
import numpy
import matplotlib
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D
import scipy
from scipy import signal
from scipy import stats

Case = 'SDprox2_E_COM_I_COM'
EXCSLM = 1
EXCSR = 1

Inh90SR = 1
Inh180SR = 1
Inh270SR = 1
Inh90SLM = 1
Inh180SLM = 1
Inh270SLM = 1

EXC = 1
INH = 1

numExcThetaSyns = 27 # i.e. 3 connections so 9 synapses per connection
numInhThetaSyns = 8 # i.e. 2 connections so 4 synapses per connection
SaveExample = 1
synspikesrandseedvec1 = numpy.array([10,1,89189,76511,23884976]) # Note that I use the first index in future simulations
synspikesrandseedvec2 = numpy.array([15,77716,267826,47246,28])
synlocrandseedvec1 = numpy.array([5,889829,294,76161,299835])
synlocrandseedvec2 = numpy.array([2,26556284,42,189,9817])

# HC Treshold Measurement Values
tstop = 10 # seconds
font_size = 13

Examples = numpy.load('NPYfiles/' + Case + '_ExampleHCModelParams.npy')
ExampleStrings = numpy.load('NPYfiles/' + Case + '_ExampleHCModelStrings.npy')
thetaSynMultiplier = numpy.array([0, 1, 2, 3])

tstop = h.tstop/1000
dt = h.dt

for randix in range(0,5):
	x = 0
	randseedsynspikes1 = synspikesrandseedvec1[randix]
	randseedsynspikes2 = synspikesrandseedvec2[randix]
	randseedsynloc1 = synlocrandseedvec1[randix]
	randseedsynloc2 = synlocrandseedvec2[randix]
	if Examples[0][x] == 0:
		print(ExampleStrings[x].decode("utf-8") + ' is empty')
		continue
	for y in range(0,4):
		print('Simulating... ' + str(ExampleStrings[x].decode("utf-8")) + ' #' + str(y+1))
		ExampleString = ExampleStrings[x].decode("utf-8")
		HCNumber = x
		# Run Simulation of Example
		h.randomize_syns(randseedsynloc1,randseedsynloc2)
		h.f(Examples[0][x],Examples[1][x],Examples[2][x],Examples[3][x],SaveExample,randseedsynspikes1,randseedsynspikes2,1,INH*numInhThetaSyns*thetaSynMultiplier[y],EXC*numExcThetaSyns*thetaSynMultiplier[y],EXCSLM,EXCSR,Inh90SR,Inh180SR,Inh270SR,Inh90SLM,Inh180SLM,Inh270SLM) # i.e. same random seeds when comparing runs
		
		HC_Trace = numpy.fromfile("%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s" % ('model_',str(Examples[0][x]),'_NumInh_',str(Examples[1][x]),'_NumExc_',str(Examples[2][x]),'_InhSpikes_',str(Examples[3][x]),'_ExcSRSpikes_',str(Examples[3][x]),'_ExcSLMSpikes_',str(9),'_NumExcCommon_',str(4),'_NumInhCommon_X',str(thetaSynMultiplier[y]),'_ThetaMultiplier.dat'),dtype=float)
		voltvec = HC_Trace[1:len(HC_Trace)]
		
		if y == 0:
			HC_Trace_Baseline = HC_Trace[10000:100001]
			HC_SpikeTimes_Baseline = numpy.zeros((len(HC_Trace),), dtype=numpy.float)
			for i in range(0,len(h.apctimes)): HC_SpikeTimes_Baseline[int(h.apctimes.x[i]/dt)] = h.apctimes.x[i]
			HC_SpikeTimes_Baseline = HC_SpikeTimes_Baseline[10000:100001]
			HC_SpikeTimes_Baseline2 = numpy.array(h.apctimes,dtype=numpy.float)
		elif y == 1:
			HC_Trace_Rhythm = HC_Trace[10000:100001]
			HC_SpikeTimes_Rhythm = numpy.zeros((len(HC_Trace),), dtype=numpy.float)
			for i in range(0,len(h.apctimes)): HC_SpikeTimes_Rhythm[int(h.apctimes.x[i]/dt)] = h.apctimes.x[i]
			HC_SpikeTimes_Rhythm = HC_SpikeTimes_Rhythm[10000:100001]
			HC_SpikeTimes_Rhythm2 = numpy.array(h.apctimes,dtype=numpy.float)
		elif y == 2:
			HC_Trace_X2Rhythm = HC_Trace[10000:100001]
			HC_SpikeTimes_X2Rhythm = numpy.zeros((len(HC_Trace),), dtype=numpy.float)
			for i in range(0,len(h.apctimes)): HC_SpikeTimes_X2Rhythm[int(h.apctimes.x[i]/dt)] = h.apctimes.x[i]
			HC_SpikeTimes_X2Rhythm = HC_SpikeTimes_X2Rhythm[10000:100001]
			HC_SpikeTimes_X2Rhythm2 = numpy.array(h.apctimes,dtype=numpy.float)
		elif y == 3:
			HC_Trace_X3Rhythm = HC_Trace[10000:100001]
			HC_SpikeTimes_X3Rhythm = numpy.zeros((len(HC_Trace),), dtype=numpy.float)
			for i in range(0,len(h.apctimes)): HC_SpikeTimes_X3Rhythm[int(h.apctimes.x[i]/dt)] = h.apctimes.x[i]
			HC_SpikeTimes_X3Rhythm = HC_SpikeTimes_X3Rhythm[10000:100001]
			HC_SpikeTimes_X3Rhythm2 = numpy.array(h.apctimes,dtype=numpy.float)
		
		timevec = numpy.arange(0,10000.1,0.1)
		if y == 0:
			f, axarr = matplotlib.pyplot.subplots(4, sharex=True)
			axarr[0].plot(timevec,voltvec,'b')
			axarr[0].set_title('Baseline High-Conductance',fontsize = font_size-2)
			# axarr[0].set_title("%s%s%s%s%s%s%s%s%s" % ('InhSyns: ',str(Examples[0][x]),', ExcSyns: ',str(Examples[1][x]),',\nInhSpikes: ',str(Examples[2][x]/tstop),'Hz, ExcSpikes: ',str(Examples[3][x]/tstop),'Hz'))
			axarr[0].set_ylim(-85,30)
			axarr[0].set_xlim(8000,10000)
			axarr[0].spines['right'].set_visible(False)
			axarr[0].spines['top'].set_visible(False)
			axarr[0].spines['bottom'].set_visible(False)
			for tic in axarr[0].xaxis.get_major_ticks():
				tic.tick1line.set_visible = tic.tick2line.set_visible = False
		if y == 1:
			axarr[1].plot(timevec,voltvec,'r')
			axarr[1].set_title('Theta Input X1',fontsize=font_size-2)
			axarr[1].set_ylim(-85,30)
			axarr[1].set_ylabel('Voltage (mV)     ',fontsize=font_size-2)
			axarr[1].set_xlim(8000,10000)
			axarr[1].spines['right'].set_visible(False)
			axarr[1].spines['top'].set_visible(False)
			axarr[1].spines['bottom'].set_visible(False)
			for tic in axarr[1].xaxis.get_major_ticks():
				tic.tick1line.set_visible = tic.tick2line.set_visible = False
		if y == 2:
			axarr[2].plot(timevec,voltvec,'g')
			axarr[2].set_title('Theta Input X2',fontsize=font_size-2)
			axarr[2].set_ylim(-85,30)
			axarr[2].set_xlim(8000,10000)
			axarr[2].spines['right'].set_visible(False)
			axarr[2].spines['top'].set_visible(False)
			axarr[2].spines['bottom'].set_visible(False)
			for tic in axarr[2].xaxis.get_major_ticks():
				tic.tick1line.set_visible = tic.tick2line.set_visible = False
		if y == 3:
			axarr[3].plot(timevec,voltvec,'c')
			axarr[3].set_title('Theta Input X3',fontsize=font_size-2)
			axarr[3].set_ylim(-85,30)
			axarr[3].set_xlim(8000,10000)
			axarr[3].set_xlabel('Time (ms)',fontsize=font_size-2)
			axarr[3].spines['right'].set_visible(False)
			axarr[3].spines['top'].set_visible(False)
			
			pyplot.savefig('PLOTfiles/' + Case + '_Trace_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
			pyplot.savefig('PLOTfiles/' + Case + '_Trace_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
			pyplot.gcf().clear()
			pyplot.cla()
			pyplot.clf()
			pyplot.close()
	f, axarr = matplotlib.pyplot.subplots(2, sharex=False)
	f1, Pxx_den1 = signal.welch(HC_SpikeTimes_Baseline, 1/(dt/1000), nperseg=25000)
	f2, Pxx_den2 = signal.welch(HC_SpikeTimes_Rhythm, 1/(dt/1000), nperseg=25000)
	f3, Pxx_den3 = signal.welch(HC_SpikeTimes_X2Rhythm, 1/(dt/1000), nperseg=25000)
	f4, Pxx_den4 = signal.welch(HC_SpikeTimes_X3Rhythm, 1/(dt/1000), nperseg=25000)
	numpy.save('NPYfiles/' + Case + '_SpikeTimesBaseline_' + str(randix) + '_randix_' + str(HCNumber) + '_HCNumber.npy',HC_SpikeTimes_Baseline2)
	numpy.save('NPYfiles/' + Case + '_SpikeTimesRhythm_' + str(randix) + '_randix_' + str(HCNumber) + '_HCNumber.npy',HC_SpikeTimes_Rhythm2)
	numpy.save('NPYfiles/' + Case + '_SpikeTimesX2Rhythm_' + str(randix) + '_randix_' + str(HCNumber) + '_HCNumber.npy',HC_SpikeTimes_X2Rhythm2)
	numpy.save('NPYfiles/' + Case + '_SpikeTimesX3Rhythm_' + str(randix) + '_randix_' + str(HCNumber) + '_HCNumber.npy',HC_SpikeTimes_X3Rhythm2)
	axarr[0].loglog(f1, Pxx_den1,'b')
	axarr[0].loglog(f2, Pxx_den2,'r')
	axarr[0].loglog(f3, Pxx_den3,'g')
	axarr[0].loglog(f4, Pxx_den4,'c')
	axarr[0].hlines(numpy.amax(Pxx_den4),4,12,'k',linestyles='solid')
	axarr[0].text(4.05,numpy.amax(Pxx_den4)+30,'Theta (4-12Hz)')
	axarr[0].axvline(numpy.array([4]),ymin=0,ymax=0.95,color='k',linestyle='solid')
	axarr[0].axvline(numpy.array([12]),ymin=0,ymax=0.95,color='k',linestyle='solid')
	axarr[0].axvline(numpy.array([8]),ymin=0,ymax=0.95,color='k',linestyle='dashed')
	axarr[0].set_xlim(0,100)
	axarr[0].set_xlabel('frequency (Hz)')
	axarr[0].set_ylabel(r'$PSD (Spikes^2 / Hz)$')
	axarr[0].spines['right'].set_visible(False)
	axarr[0].spines['top'].set_visible(False)
	
	ind = numpy.arange(8)
	width = 0.4
	Area1 = numpy.trapz(Pxx_den1[(f1>4) & (f1<12)],x=f1[(f1>4) & (f1<12)])
	Area2 = numpy.trapz(Pxx_den2[(f2>4) & (f2<12)],x=f2[(f2>4) & (f2<12)])
	Area3 = numpy.trapz(Pxx_den3[(f3>4) & (f3<12)],x=f3[(f3>4) & (f3<12)])
	Area4 = numpy.trapz(Pxx_den4[(f4>4) & (f4<12)],x=f4[(f4>4) & (f4<12)])
	e8Hz1 = Pxx_den1[f1==8]
	e8Hz2 = Pxx_den2[f2==8]
	e8Hz3 = Pxx_den3[f3==8]
	e8Hz4 = Pxx_den4[f4==8]
	axarr[1].bar(ind+width, [Area1, Area2, Area3, Area4, e8Hz1[0], e8Hz2[0], e8Hz3[0], e8Hz4[0]], width, color='k')
	axarr[1].set_xticks(ind+width)
	axarr[1].set_xticklabels(('Base\n(4-12Hz)','ThetaX1\n(4-12Hz)', 'ThetaX2\n(4-12Hz)', 'ThetaX3\n(4-12Hz)', 'Base\n(8Hz)','ThetaX1\n(8Hz)','ThetaX2\n(8Hz)','ThetaX3\n(8Hz)'),fontsize=font_size-3, fontweight='bold', rotation=45)
	axarr[1].set_ylabel('PSD Magnitude')
	axarr[1].set_xlim(0,8+width)
	axarr[1].spines['right'].set_visible(False)
	axarr[1].spines['top'].set_visible(False)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PSD_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PSD_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	# Instantaneous Frequency Analyses (Note this is converting the ISIs into seconds and then instantaneous frequencies)
	IF_Baseline = numpy.concatenate([numpy.array([0],dtype=numpy.float),1000/numpy.diff(HC_SpikeTimes_Baseline2)])
	IF_ThetaX1 = numpy.concatenate([numpy.array([0],dtype=numpy.float),1000/numpy.diff(HC_SpikeTimes_Rhythm2)])
	IF_ThetaX2 = numpy.concatenate([numpy.array([0],dtype=numpy.float),1000/numpy.diff(HC_SpikeTimes_X2Rhythm2)])
	IF_ThetaX3 = numpy.concatenate([numpy.array([0],dtype=numpy.float),1000/numpy.diff(HC_SpikeTimes_X3Rhythm2)])
	
	heights1,bins1 = numpy.histogram(IF_Baseline,bins=100)
	heights2,bins2 = numpy.histogram(IF_ThetaX1,bins=100)
	heights3,bins3 = numpy.histogram(IF_ThetaX2,bins=100)
	heights4,bins4 = numpy.histogram(IF_ThetaX3,bins=100)
	
	# Normalize
	heights1 = heights1/float(sum(heights1))
	heights2 = heights2/float(sum(heights2))
	heights3 = heights3/float(sum(heights3))
	heights4 = heights4/float(sum(heights4))
	bin1=bins1[:-1]+numpy.diff(bins1)/2.
	bin2=bins2[:-1]+numpy.diff(bins2)/2.
	bin3=bins3[:-1]+numpy.diff(bins3)/2.
	bin4=bins4[:-1]+numpy.diff(bins4)/2.
	binMids1 = bin1[~numpy.isnan(bin1)]
	binMids2 = bin2[~numpy.isnan(bin2)]
	binMids3 = bin3[~numpy.isnan(bin3)]
	binMids4 = bin4[~numpy.isnan(bin4)]
	heights1 = heights1[~numpy.isnan(bin1)]
	heights2 = heights2[~numpy.isnan(bin2)]
	heights3 = heights3[~numpy.isnan(bin3)]
	heights4 = heights4[~numpy.isnan(bin4)]
	
	f, axarr = matplotlib.pyplot.subplots(2, sharex=False)
	axarr[0].semilogx(binMids1,heights1,'b')
	axarr[0].semilogx(binMids2,heights2,'r')
	axarr[0].semilogx(binMids3,heights3,'g')
	axarr[0].semilogx(binMids4,heights4,'c')
	axarr[0].vlines(numpy.array([4,12]),0,numpy.amax(heights4)+0.01,'k',linestyles='solid')
	axarr[0].vlines(numpy.array([8]),0,numpy.amax(heights4)+0.01,'k',linestyles='dashed')
	axarr[0].hlines(numpy.amax(heights4)+0.01,4,12,'k',linestyles='solid')
	axarr[0].text(4.3,numpy.amax(heights4)+0.02,'Theta (4-12Hz)')
	axarr[0].set_xlabel('Frequency (Hz)')
	axarr[0].set_ylabel('Probability')
	axarr[0].set_xlim(0,130)
	axarr[0].set_ylim(0,0.19)
	axarr[0].spines['right'].set_visible(False)
	axarr[0].spines['top'].set_visible(False)
	
	ind = numpy.arange(8)
	width = 0.4
	Area11 = numpy.trapz(heights1[(binMids1>4) & (binMids1<8)],x=binMids1[(binMids1>4) & (binMids1<8)])
	Area22 = numpy.trapz(heights2[(binMids2>4) & (binMids2<8)],x=binMids2[(binMids2>4) & (binMids2<8)])
	Area33 = numpy.trapz(heights3[(binMids3>4) & (binMids3<8)],x=binMids3[(binMids3>4) & (binMids3<8)])
	Area44 = numpy.trapz(heights4[(binMids4>4) & (binMids4<8)],x=binMids4[(binMids4>4) & (binMids4<8)])
	e8Hz11 = numpy.trapz(heights1[(binMids1>8) & (binMids1<12)],x=binMids1[(binMids1>8) & (binMids1<12)]) #IF NOT THIS TRY NUMPY.AMAX BETWEEN 7 AND 9
	e8Hz22 = numpy.trapz(heights2[(binMids2>8) & (binMids2<12)],x=binMids2[(binMids2>8) & (binMids2<12)])
	e8Hz33 = numpy.trapz(heights3[(binMids3>8) & (binMids3<12)],x=binMids3[(binMids3>8) & (binMids3<12)])
	e8Hz44 = numpy.trapz(heights4[(binMids4>8) & (binMids4<12)],x=binMids4[(binMids4>8) & (binMids4<12)])
	axarr[1].bar(ind+width, [Area11, Area22, Area33, Area44, e8Hz11, e8Hz22, e8Hz33, e8Hz44], width, color='k')
	axarr[1].set_xticks(ind+width)
	axarr[1].set_xticklabels(('Base\n(4-8Hz)','ThetaX1\n(4-8Hz)', 'ThetaX2\n(4-8Hz)', 'ThetaX3\n(4-8Hz)', 'Base\n(8-12Hz)','ThetaX1\n(8-12Hz)','ThetaX2\n(8-12Hz)','ThetaX3\n(8-12Hz)'),fontsize=font_size-3, fontweight='bold', rotation=45)
	axarr[1].set_ylabel('Probability')
	axarr[1].set_xlim(0,8+width)
	axarr[1].spines['right'].set_visible(False)
	axarr[1].spines['top'].set_visible(False)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_IFDistribution_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_IFDistribution_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	# Polar Distribution Plots
	STX0 = HC_SpikeTimes_Baseline2/125
	STX1 = HC_SpikeTimes_Rhythm2/125
	STX2 = HC_SpikeTimes_X2Rhythm2/125
	STX3 = HC_SpikeTimes_X3Rhythm2/125
	Rads_Baseline = 2*numpy.pi*(STX0-STX0.astype(int))
	Rads_1XRhythm = 2*numpy.pi*(STX1-STX1.astype(int))
	Rads_2XRhythm = 2*numpy.pi*(STX2-STX2.astype(int))
	Rads_3XRhythm = 2*numpy.pi*(STX3-STX3.astype(int))
	
	Rads_Baseline_Sorted = sorted(Rads_Baseline)
	Rads_1XRhythm_Sorted = sorted(Rads_1XRhythm)
	Rads_2XRhythm_Sorted = sorted(Rads_2XRhythm)
	Rads_3XRhythm_Sorted = sorted(Rads_3XRhythm)
	
	IF_Baseline_Sorted = [x for i, x in sorted(zip(Rads_Baseline,IF_Baseline))]
	IF_ThetaX1_Sorted = [x for i, x in sorted(zip(Rads_1XRhythm,IF_ThetaX1))]
	IF_ThetaX2_Sorted = [x for i, x in sorted(zip(Rads_2XRhythm,IF_ThetaX2))]
	IF_ThetaX3_Sorted = [x for i, x in sorted(zip(Rads_3XRhythm,IF_ThetaX3))]
	
	bin_size = 25
	range_if = (0,100)
	range_rads = (0,2*numpy.pi)
	
	heights11,be11 = numpy.histogram(Rads_Baseline,bins=bin_size,range=range_rads)
	bins11 = be11[:-1]+numpy.diff(be11)/2.
	PrefPhase_Baseline = bins11[heights11 == numpy.amax(heights11)]
	heights22,be22 = numpy.histogram(Rads_1XRhythm,bins=bin_size,range=range_rads)
	bins22 = be22[:-1]+numpy.diff(be22)/2.
	PrefPhase_ThetaX1 = bins22[heights22 == numpy.amax(heights22)]
	heights33,be33 = numpy.histogram(Rads_2XRhythm,bins=bin_size,range=range_rads)
	bins33 = be33[:-1]+numpy.diff(be33)/2.
	PrefPhase_ThetaX2 = bins33[heights33 == numpy.amax(heights33)]
	heights44,be44 = numpy.histogram(Rads_3XRhythm,bins=bin_size,range=range_rads)
	bins44 = be44[:-1]+numpy.diff(be44)/2.
	PrefPhase_ThetaX3 = bins44[heights44 == numpy.amax(heights44)]
	
	Rads_Baseline_MeanBins,be11,bn11 = scipy.stats.binned_statistic(Rads_Baseline_Sorted,Rads_Baseline_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	_BMeans = Rads_Baseline_MeanBins
	Rads_Baseline_MeanBins = Rads_Baseline_MeanBins[~numpy.isnan(Rads_Baseline_MeanBins)]
	Rads_Baseline_MeanBins = numpy.append(Rads_Baseline_MeanBins,Rads_Baseline_MeanBins[0])
	Rads_1XRhythm_MeanBins,be21,bn21 = scipy.stats.binned_statistic(Rads_1XRhythm_Sorted,Rads_1XRhythm_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	_1Means = Rads_1XRhythm_MeanBins
	Rads_1XRhythm_MeanBins = Rads_1XRhythm_MeanBins[~numpy.isnan(Rads_1XRhythm_MeanBins)]
	Rads_1XRhythm_MeanBins = numpy.append(Rads_1XRhythm_MeanBins,Rads_1XRhythm_MeanBins[0])
	Rads_2XRhythm_MeanBins,be31,bn31 = scipy.stats.binned_statistic(Rads_2XRhythm_Sorted,Rads_2XRhythm_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	_2Means = Rads_2XRhythm_MeanBins
	Rads_2XRhythm_MeanBins = Rads_2XRhythm_MeanBins[~numpy.isnan(Rads_2XRhythm_MeanBins)]
	Rads_2XRhythm_MeanBins = numpy.append(Rads_2XRhythm_MeanBins,Rads_2XRhythm_MeanBins[0])
	Rads_3XRhythm_MeanBins,be41,bn41 = scipy.stats.binned_statistic(Rads_3XRhythm_Sorted,Rads_3XRhythm_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	_3Means = Rads_3XRhythm_MeanBins
	Rads_3XRhythm_MeanBins = Rads_3XRhythm_MeanBins[~numpy.isnan(Rads_3XRhythm_MeanBins)]
	Rads_3XRhythm_MeanBins = numpy.append(Rads_3XRhythm_MeanBins,Rads_3XRhythm_MeanBins[0])
	
	IF_Baseline_MeanBins,be12,bn12 = scipy.stats.binned_statistic(Rads_Baseline_Sorted,IF_Baseline_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	IF_Baseline_MeanBins = IF_Baseline_MeanBins[~numpy.isnan(_BMeans)]
	IF_Baseline_MeanBins = numpy.append(IF_Baseline_MeanBins,IF_Baseline_MeanBins[0])
	IF_1XRhythm_MeanBins,be22,bn22 = scipy.stats.binned_statistic(Rads_1XRhythm_Sorted,IF_ThetaX1_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	IF_1XRhythm_MeanBins = IF_1XRhythm_MeanBins[~numpy.isnan(_1Means)]
	IF_1XRhythm_MeanBins = numpy.append(IF_1XRhythm_MeanBins,IF_1XRhythm_MeanBins[0])
	IF_2XRhythm_MeanBins,be32,bn32 = scipy.stats.binned_statistic(Rads_2XRhythm_Sorted,IF_ThetaX2_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	IF_2XRhythm_MeanBins = IF_2XRhythm_MeanBins[~numpy.isnan(_2Means)]
	IF_2XRhythm_MeanBins = numpy.append(IF_2XRhythm_MeanBins,IF_2XRhythm_MeanBins[0])
	IF_3XRhythm_MeanBins,be42,bn42 = scipy.stats.binned_statistic(Rads_3XRhythm_Sorted,IF_ThetaX3_Sorted,statistic='mean',bins=bin_size,range=range_rads)
	IF_3XRhythm_MeanBins = IF_3XRhythm_MeanBins[~numpy.isnan(_3Means)]
	IF_3XRhythm_MeanBins = numpy.append(IF_3XRhythm_MeanBins,IF_3XRhythm_MeanBins[0])
	
	IF_Baseline_MinsBins,be111,bn111 = scipy.stats.binned_statistic(Rads_Baseline_Sorted,IF_Baseline_Sorted,statistic='min',bins=bin_size,range=range_rads)
	IF_Baseline_MinsBins = IF_Baseline_MinsBins[~numpy.isnan(_BMeans)]
	IF_Baseline_MinsBins = numpy.append(IF_Baseline_MinsBins,IF_Baseline_MinsBins[0])
	IF_1XRhythm_MinsBins,be211,bn211 = scipy.stats.binned_statistic(Rads_1XRhythm_Sorted,IF_ThetaX1_Sorted,statistic='min',bins=bin_size,range=range_rads)
	IF_1XRhythm_MinsBins = IF_1XRhythm_MinsBins[~numpy.isnan(_1Means)]
	IF_1XRhythm_MinsBins = numpy.append(IF_1XRhythm_MinsBins,IF_1XRhythm_MinsBins[0])
	IF_2XRhythm_MinsBins,be311,bn311 = scipy.stats.binned_statistic(Rads_2XRhythm_Sorted,IF_ThetaX2_Sorted,statistic='min',bins=bin_size,range=range_rads)
	IF_2XRhythm_MinsBins = IF_2XRhythm_MinsBins[~numpy.isnan(_2Means)]
	IF_2XRhythm_MinsBins = numpy.append(IF_2XRhythm_MinsBins,IF_2XRhythm_MinsBins[0])
	IF_3XRhythm_MinsBins,be411,bn411 = scipy.stats.binned_statistic(Rads_3XRhythm_Sorted,IF_ThetaX3_Sorted,statistic='min',bins=bin_size,range=range_rads)
	IF_3XRhythm_MinsBins = IF_3XRhythm_MinsBins[~numpy.isnan(_3Means)]
	IF_3XRhythm_MinsBins = numpy.append(IF_3XRhythm_MinsBins,IF_3XRhythm_MinsBins[0])
	
	IF_Baseline_MaxBins,be111,bn111 = scipy.stats.binned_statistic(Rads_Baseline_Sorted,IF_Baseline_Sorted,statistic='max',bins=bin_size,range=range_rads)
	IF_Baseline_MaxBins = IF_Baseline_MaxBins[~numpy.isnan(_BMeans)]
	IF_Baseline_MaxBins = numpy.append(IF_Baseline_MaxBins,IF_Baseline_MaxBins[0])
	IF_1XRhythm_MaxBins,be211,bn211 = scipy.stats.binned_statistic(Rads_1XRhythm_Sorted,IF_ThetaX1_Sorted,statistic='max',bins=bin_size,range=range_rads)
	IF_1XRhythm_MaxBins = IF_1XRhythm_MaxBins[~numpy.isnan(_1Means)]
	IF_1XRhythm_MaxBins = numpy.append(IF_1XRhythm_MaxBins,IF_1XRhythm_MaxBins[0])
	IF_2XRhythm_MaxBins,be311,bn311 = scipy.stats.binned_statistic(Rads_2XRhythm_Sorted,IF_ThetaX2_Sorted,statistic='max',bins=bin_size,range=range_rads)
	IF_2XRhythm_MaxBins = IF_2XRhythm_MaxBins[~numpy.isnan(_2Means)]
	IF_2XRhythm_MaxBins = numpy.append(IF_2XRhythm_MaxBins,IF_2XRhythm_MaxBins[0])
	IF_3XRhythm_MaxBins,be411,bn411 = scipy.stats.binned_statistic(Rads_3XRhythm_Sorted,IF_ThetaX3_Sorted,statistic='max',bins=bin_size,range=range_rads)
	IF_3XRhythm_MaxBins = IF_3XRhythm_MaxBins[~numpy.isnan(_3Means)]
	IF_3XRhythm_MaxBins = numpy.append(IF_3XRhythm_MaxBins,IF_3XRhythm_MaxBins[0])
	
	IF_Baseline_StdBins,be111,bn111 = scipy.stats.binned_statistic(Rads_Baseline_Sorted,IF_Baseline_Sorted,statistic='std',bins=bin_size,range=range_rads)
	IF_Baseline_StdBins = IF_Baseline_StdBins[~numpy.isnan(_BMeans)]
	IF_Baseline_StdBins = numpy.append(IF_Baseline_StdBins,IF_Baseline_StdBins[0])
	IF_1XRhythm_StdBins,be211,bn211 = scipy.stats.binned_statistic(Rads_1XRhythm_Sorted,IF_ThetaX1_Sorted,statistic='std',bins=bin_size,range=range_rads)
	IF_1XRhythm_StdBins = IF_1XRhythm_StdBins[~numpy.isnan(_1Means)]
	IF_1XRhythm_StdBins = numpy.append(IF_1XRhythm_StdBins,IF_1XRhythm_StdBins[0])
	IF_2XRhythm_StdBins,be311,bn311 = scipy.stats.binned_statistic(Rads_2XRhythm_Sorted,IF_ThetaX2_Sorted,statistic='std',bins=bin_size,range=range_rads)
	IF_2XRhythm_StdBins = IF_2XRhythm_StdBins[~numpy.isnan(_2Means)]
	IF_2XRhythm_StdBins = numpy.append(IF_2XRhythm_StdBins,IF_2XRhythm_StdBins[0])
	IF_3XRhythm_StdBins,be411,bn411 = scipy.stats.binned_statistic(Rads_3XRhythm_Sorted,IF_ThetaX3_Sorted,statistic='std',bins=bin_size,range=range_rads)
	IF_3XRhythm_StdBins = IF_3XRhythm_StdBins[~numpy.isnan(_3Means)]
	IF_3XRhythm_StdBins = numpy.append(IF_3XRhythm_StdBins,IF_3XRhythm_StdBins[0])
	
	axarr = pyplot.subplot(111, projection='polar')
	axarr.plot(Rads_Baseline_Sorted,IF_Baseline_Sorted,'b',label='Baseline')
	axarr.plot(Rads_1XRhythm_Sorted,IF_ThetaX1_Sorted,'r',label='ThetaX1')
	axarr.plot(Rads_2XRhythm_Sorted,IF_ThetaX2_Sorted,'g',label='ThetaX2')
	axarr.plot(Rads_3XRhythm_Sorted,IF_ThetaX3_Sorted,'c',label='ThetaX3')
	position = 292.5
	axarr._r_label_position._t = (position, 0)
	axarr._r_label_position.invalidate()
	axarr.set_xticklabels(['0$^\circ$', '45$^\circ$', '90$^\circ$\nPeak', '135$^\circ$', '180$^\circ$', '225$^\circ$', '270$^\circ$\nTrough', '315$^\circ$'])
	# axarr.legend(loc=1)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlot_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlot_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	axarr = pyplot.subplot(111, projection='polar')
	axarr.plot(Rads_Baseline_MeanBins,IF_Baseline_MeanBins,'b',label='Baseline')
	# axarr.fill_between(Rads_Baseline_MeanBins,IF_Baseline_MinsBins,IF_Baseline_MaxBins,alpha=0.5, edgecolor='b', facecolor='b')
	axarr.fill_between(Rads_Baseline_MeanBins,numpy.clip(IF_Baseline_MeanBins-IF_Baseline_StdBins,0,1000),numpy.clip(IF_Baseline_MeanBins+IF_Baseline_StdBins,0,1000),alpha=0.5, edgecolor='b', facecolor='b')
	axarr.plot(Rads_1XRhythm_MeanBins,IF_1XRhythm_MeanBins,'r',label='ThetaX1')
	# axarr.fill_between(Rads_1XRhythm_MeanBins,IF_1XRhythm_MinsBins,IF_1XRhythm_MaxBins,alpha=0.5, edgecolor='r', facecolor='r')
	axarr.fill_between(Rads_1XRhythm_MeanBins,numpy.clip(IF_1XRhythm_MeanBins-IF_1XRhythm_StdBins,0,1000),numpy.clip(IF_1XRhythm_MeanBins+IF_1XRhythm_StdBins,0,1000),alpha=0.5, edgecolor='r', facecolor='r')
	axarr.plot(Rads_2XRhythm_MeanBins,IF_2XRhythm_MeanBins,'g',label='ThetaX2')
	# axarr.fill_between(Rads_2XRhythm_MeanBins,IF_2XRhythm_MinsBins,IF_2XRhythm_MaxBins,alpha=0.5, edgecolor='g', facecolor='g')
	axarr.fill_between(Rads_2XRhythm_MeanBins,numpy.clip(IF_2XRhythm_MeanBins-IF_2XRhythm_StdBins,0,1000),numpy.clip(IF_2XRhythm_MeanBins+IF_2XRhythm_StdBins,0,1000),alpha=0.5, edgecolor='g', facecolor='g')
	axarr.plot(Rads_3XRhythm_MeanBins,IF_3XRhythm_MeanBins,'c',label='ThetaX3')
	# axarr.fill_between(Rads_3XRhythm_MeanBins,IF_3XRhythm_MinsBins,IF_3XRhythm_MaxBins,alpha=0.5, edgecolor='c', facecolor='c')
	axarr.fill_between(Rads_3XRhythm_MeanBins,numpy.clip(IF_3XRhythm_MeanBins-IF_3XRhythm_StdBins,0,1000),numpy.clip(IF_3XRhythm_MeanBins+IF_3XRhythm_StdBins,0,1000),alpha=0.5, edgecolor='c', facecolor='c')
	position = 292.5
	axarr._r_label_position._t = (position, 0)
	axarr._r_label_position.invalidate()
	axarr.set_xticklabels(['0$^\circ$', '45$^\circ$', '90$^\circ$\nPeak', '135$^\circ$', '180$^\circ$', '225$^\circ$', '270$^\circ$\nTrough', '315$^\circ$'])
	# axarr.legend(loc=1)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	axarr = pyplot.subplot(111, projection='polar')
	axarr.plot(Rads_Baseline_MeanBins,IF_Baseline_MeanBins,'b',label='Baseline')
	# axarr.fill_between(Rads_Baseline_MeanBins,IF_Baseline_MinsBins,IF_Baseline_MaxBins,alpha=0.5, edgecolor='b', facecolor='b')
	axarr.fill_between(Rads_Baseline_MeanBins,numpy.clip(IF_Baseline_MeanBins-IF_Baseline_StdBins,0,1000),numpy.clip(IF_Baseline_MeanBins+IF_Baseline_StdBins,0,1000),alpha=0.5, edgecolor='b', facecolor='b')
	position = 292.5
	bottom = numpy.amax(numpy.clip(IF_Baseline_MeanBins+IF_Baseline_StdBins,0,1000))+5
	width = range_rads[1]/(bin_size+1)
	axarr.bar(_BMeans, heights11, width=width, bottom=bottom,color='b')
	axarr._r_label_position._t = (position, 0)
	axarr._r_label_position.invalidate()
	axarr.set_xticklabels(['0$^\circ$', '45$^\circ$', '90$^\circ$\nPeak', '135$^\circ$', '180$^\circ$', '225$^\circ$', '270$^\circ$\nTrough', '315$^\circ$'])
	# axarr.legend(loc=1)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinnedBaseline_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinnedBaseline_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	axarr = pyplot.subplot(111, projection='polar')
	axarr.plot(Rads_1XRhythm_MeanBins,IF_1XRhythm_MeanBins,'r',label='ThetaX1')
	# axarr.fill_between(Rads_1XRhythm_MeanBins,IF_1XRhythm_MinsBins,IF_1XRhythm_MaxBins,alpha=0.5, edgecolor='r', facecolor='r')
	axarr.fill_between(Rads_1XRhythm_MeanBins,numpy.clip(IF_1XRhythm_MeanBins-IF_1XRhythm_StdBins,0,1000),numpy.clip(IF_1XRhythm_MeanBins+IF_1XRhythm_StdBins,0,1000),alpha=0.5, edgecolor='r', facecolor='r')
	position = 292.5
	bottom = numpy.amax(numpy.clip(IF_1XRhythm_MeanBins+IF_1XRhythm_StdBins,0,1000))+5
	width = range_rads[1]/(bin_size+1)
	axarr.bar(_1Means, heights22, width=width, bottom=bottom,color='r')
	axarr._r_label_position._t = (position, 0)
	axarr._r_label_position.invalidate()
	axarr.set_xticklabels(['0$^\circ$', '45$^\circ$', '90$^\circ$\nPeak', '135$^\circ$', '180$^\circ$', '225$^\circ$', '270$^\circ$\nTrough', '315$^\circ$'])
	# axarr.legend(loc=1)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned1XTheta_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned1XTheta_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	axarr = pyplot.subplot(111, projection='polar')
	axarr.plot(Rads_2XRhythm_MeanBins,IF_2XRhythm_MeanBins,'g',label='ThetaX2')
	# axarr.fill_between(Rads_2XRhythm_MeanBins,IF_2XRhythm_MinsBins,IF_2XRhythm_MaxBins,alpha=0.5, edgecolor='g', facecolor='g')
	axarr.fill_between(Rads_2XRhythm_MeanBins,numpy.clip(IF_2XRhythm_MeanBins-IF_2XRhythm_StdBins,0,1000),numpy.clip(IF_2XRhythm_MeanBins+IF_2XRhythm_StdBins,0,1000),alpha=0.5, edgecolor='g', facecolor='g')
	position = 292.5
	bottom = numpy.amax(numpy.clip(IF_2XRhythm_MeanBins+IF_2XRhythm_StdBins,0,1000))+5
	width = range_rads[1]/(bin_size+1)
	axarr.bar(_2Means, heights33, width=width, bottom=bottom,color='g')
	axarr._r_label_position._t = (position, 0)
	axarr._r_label_position.invalidate()
	axarr.set_xticklabels(['0$^\circ$', '45$^\circ$', '90$^\circ$\nPeak', '135$^\circ$', '180$^\circ$', '225$^\circ$', '270$^\circ$\nTrough', '315$^\circ$'])
	# axarr.legend(loc=1)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned2XTheta_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned2XTheta_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()
	
	axarr = pyplot.subplot(111, projection='polar')
	axarr.plot(Rads_3XRhythm_MeanBins,IF_3XRhythm_MeanBins,'c',label='ThetaX3')
	# axarr.fill_between(Rads_3XRhythm_MeanBins,IF_3XRhythm_MinsBins,IF_3XRhythm_MaxBins,alpha=0.5, edgecolor='c', facecolor='c')
	axarr.fill_between(Rads_3XRhythm_MeanBins,numpy.clip(IF_3XRhythm_MeanBins-IF_3XRhythm_StdBins,0,1000),numpy.clip(IF_3XRhythm_MeanBins+IF_3XRhythm_StdBins,0,1000),alpha=0.5, edgecolor='c', facecolor='c')
	position = 292.5
	bottom = numpy.amax(numpy.clip(IF_3XRhythm_MeanBins+IF_3XRhythm_StdBins,0,1000))+5
	width = range_rads[1]/(bin_size+1)
	axarr.bar(_3Means, heights44, width=width, bottom=bottom,color='c')
	axarr._r_label_position._t = (position, 0)
	axarr._r_label_position.invalidate()
	axarr.set_xticklabels(['0$^\circ$', '45$^\circ$', '90$^\circ$\nPeak', '135$^\circ$', '180$^\circ$', '225$^\circ$', '270$^\circ$\nTrough', '315$^\circ$'])
	# axarr.legend(loc=1)
	pyplot.tight_layout()
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned3XTheta_' + str(randix) + '_randix_' + ExampleString + '.pdf', bbox_inches='tight')
	pyplot.savefig('PLOTfiles/' + Case + '_PolarPlotBinned3XTheta_' + str(randix) + '_randix_' + ExampleString + '.png', bbox_inches='tight')
	pyplot.gcf().clear()
	pyplot.cla()
	pyplot.clf()
	pyplot.close()

numpy.save('NPYfiles/' + Case + '_Areas.npy',numpy.array([Area1,Area2,Area3,Area4],dtype=numpy.float))
numpy.save('NPYfiles/' + Case + '_e8Hz.npy',numpy.array([e8Hz1,e8Hz2,e8Hz3,e8Hz4],dtype=numpy.float))

# LNI_LNE_LIS_LES
# HNI_LNE_LIS_LES
# LNI_HNE_LIS_LES
# HNI_HNE_LIS_LES
# LNI_LNE_HIS_LES
# HNI_LNE_HIS_LES
# LNI_HNE_HIS_LES
# HNI_HNE_HIS_LES
# LNI_LNE_LIS_HES
# HNI_LNE_LIS_HES
# LNI_HNE_LIS_HES
# HNI_HNE_LIS_HES
# LNI_LNE_HIS_HES
# HNI_LNE_HIS_HES
# LNI_HNE_HIS_HES
# HNI_HNE_HIS_HES

Loading data, please wait...