In [60]:
#bigger plots
figsize(14,5)
In [61]:
import scipy.io
a1 = scipy.io.loadmat('data/A1_SCV_SPK_040528_cell1_001_sec333_end_nostruct.mat',squeeze_me=True)
pp1 = scipy.io.loadmat('data/PP1_Gamma_course_day2_nostruct.mat',squeeze_me=True)
pp1_select = 0 #0-2 select one of the 3 spiketrains from the data
If you detected nspike in your observation interval, how many inter-spike intervals do you have?
Compute the inter-spike intervals by simply subtracting from each spike time the previous spike time.
What is the average ISI length?
To estimate the distribution of ISIs construct a histogram of ISIs.Define a reasonable width of your histogram classes (bin width) and construct a vector of histogram edges (bin vector).
In [62]:
# Load the Data
spikes_rat_cortex = a1['TotalOriginalSpikeTrainMS']
#select one of the three spike trains (pp1_select) and convert to ms
spikes_gamma_renewal = pp1['Spikes'][pp1_select] * 1e3
# Compute the ISIs
ISIs_rat = diff(spikes_rat_cortex)
ISIs_gamma = diff(spikes_gamma_renewal)
In [63]:
print 'mean ISI cortical cell:', ISIs_rat.mean(),'±',ISIs_rat.std(), 'ms'
print 'mean ISI renewal gamma:', ISIs_gamma.mean(),'±',ISIs_gamma.std(), 'ms'
In [64]:
# Construct the bin vector
# min value: 0, max value: longest ISI (rat or gamma), number of bins: nbin
nbin = 20
bins = linspace(0, max(ISIs_rat.max(),ISIs_gamma.max()), nbin)
# Plot
subplot(1,2,1)
cr,binr,_ = hist(ISIs_rat, bins, normed=True)
xlabel('time[ms]')
title('ISI distribution cortical cell')
subplot(1,2,2)
cg,bing,_ = hist(ISIs_gamma, bins, normed=True)
title('ISI distribution gamma point process')
xlabel('time[ms]')
Out[64]:
The ISI CV is a measure of Interval Variability
In [65]:
CV_rat = ISIs_rat.std()/ISIs_rat.mean()
CV_gamma = ISIs_gamma.std()/ISIs_gamma.mean()
print 'ISI cortical cell:', CV_rat
print 'ISI renewal gamma:', CV_gamma
Can you think of biophysical mechanisms that explain why the output train of action potentials of a biological neuron cannot be memoryless as in the case of the Poisson process? Briefly describe them in your protocol.
In [66]:
# Refraction Time (relative & absolute)
The autocorrelogram measures for any spike the (average) probability of finding another spike in temporal distance dt. In other words, we ask for any given spike $s_i$ at time $t_i=0$, how likely is it to find a second spike $s_j$ at time $t_j \neq 0$.
Optional: Program your own function that performs a cross - correlation of discrete event times. Numerically you may follow three steps:
In [67]:
#shorten variable names (optional)
spikes_r = spikes_rat_cortex
spikes_g = spikes_gamma_renewal
In [68]:
#select the correlation time and resolution [in ms]
tcorr = 200 #ms
dtcorr = 20 #ms
#construct bin vector for histogram
nbin = tcorr/dtcorr
bins = linspace(1,tcorr,nbin+1)
#calculate pairwise spike-time differences
spikes_pdiff_r = subtract(spikes_r[:,newaxis],spikes_r[newaxis,:])
spikes_pdiff_g = subtract(spikes_g[:,newaxis],spikes_g[newaxis,:])
In [69]:
# Plot
subplot(1,2,1)
cr,bins,_ = hist(spikes_pdiff_r.flatten(), bins, normed=True)
xlim([0,tcorr])
xlabel('time[ms]')
title('Autocorrelation cortical cell')
subplot(1,2,2)
cg,bins,_ = hist(spikes_pdiff_g.flatten(), bins, normed=True)
title('Autocorrelation gamma point process')
xlim([0,tcorr])
xlabel('time[ms]')
Out[69]:
In [70]:
# Alternatively, use numpys correlate function on PSTH
#first, calculate the psth using a bin vector with desired resolution
bins_r = arange(0,spikes_rat_cortex.max(), dtcorr)
bins_g = arange(0,spikes_gamma_renewal.max(), dtcorr)
PSTH_r,bins_r = histogram(spikes_rat_cortex,bins_r)
PSTH_g,bins_g = histogram(spikes_gamma_renewal,bins_g)
# the autocorrelation at time t is the correlation of the signal with signal shifted by t
# roll(PSTH,i) shifts the time of the psth by i*dt
autoc_r = [correlate(PSTH_r,roll(PSTH_r,i)) for i in arange(1,nbin)]
autoc_g = [correlate(PSTH_g,roll(PSTH_g,i)) for i in arange(1,nbin)]
In [71]:
# Plot
subplot(1,2,1)
plot(linspace(dtcorr,tcorr,nbin-1),autoc_r)
xlim([0,tcorr])
xlabel('time[ms]')
title('Autocorrelation cortical cell')
subplot(1,2,2)
plot(linspace(dtcorr,tcorr,nbin-1),autoc_g,)
title('Autocorrelation gamma point process')
xlim([0,tcorr])
xlabel('time[ms]')
Out[71]:
Side Note: The naive way of calculating the PSTH autocorrelation would be:
In [72]:
ac_full = correlate(PSTH_r,PSTH_r,mode='full')
ac_full[ac_full==ac_full.max()]=0 #optional, remove the peak at 0
time = linspace(-len(PSTH_r)*dtcorr,len(PSTH_r)*dtcorr,len(ac_full))
plot(time,ac_full)
xlabel('correlation time [ms]')
Out[72]:
This is the full autocorrelation up to times of 100000ms (=10s). The triangular shape is actually an artefact due to finate data. However, we are only intersted in the biophysical time on the order of 100 ms. If we zoom in into this time scale we should see the same autocorrelation shape as above.
Thanks to Katharina and Diana for the code base, thanks to Christian for hinting towards ´np.subtract´ for pairwise differences