In [1]:
import scipy.io
In [2]:
# please edit the data path if necessary
# squeeze_me=True removes unnessecary dimensions
data = scipy.io.loadmat('data/A1_data_set_040528_boucsein_nostruct.mat',
squeeze_me=True)
V = data['V']
dt = data['SampleIntervalSeconds']
print V.shape, dt
ONE
Plot only a short piece (e.g. 2s) of the complete voltage trace.
Set the range of the yaxis to (-75mV, 20mV). Compare to Fig. 1A in [1] which shows the same neuron.
How long is the complete voltage trace in seconds?
In [3]:
plot(V[0:2/dt])
ylim(-75,20)
T = len(V)*dt
print T, 'seconds for the complete trace'
TWO
How many action potentials did you detect? What is the average firing rate (i.e. action potentials per second) of this neuron?
In [4]:
threshold = - 30
V_above_threshold = V>threshold
# we need the first element of threshold crossing for each spike
# convert the bool array to int (0 and 1), diff is the elementwise difference
V_above_th_diff = diff(V_above_threshold.astype(int)) == 1
# 1 for spike onset and -1 for spike offset; select only onsets
spike_bins = V_above_th_diff == 1
N_spikes = spike_bins.sum()
print 'number of spikes:', N_spikes.sum(), '; firing rate:', N_spikes/T, 'Hz'
plot(spike_bins[0:2/dt])
Out[4]:
THREE - Spike Triggered Average (STA)
Now determine the average waveform of an action potential. To this end you need to loop across all AP onset times ti andobtain a cut-out of about +/- 10 ms around the onset time and average across all the cutouts.
In [5]:
cutout = int(10e-3/dt)
window = arange(-cutout,+cutout)
# create an array with the indices of all cut-outs around spikes onset
# first, repeat the cut-out window N_spike times
spike_select_idx = tile(window,(N_spikes,1)).T
# then, add a spike index to every window
spike_select_idx = spike_select_idx + find(spike_bins)
# finally in case the first or the last spike are to close to boundary,
# clip the indices to stay within limits
spike_select_idx = spike_select_idx.clip(0,len(V)-1)
print window.shape, N_spikes, spike_select_idx.shape
In [6]:
# now we can use the index array to get the voltage for all the cut-outs
V_STA = V[spike_select_idx].mean(1)
plot(V_STA)
Out[6]:
In [7]:
# Optional: Add to the plot the average waveform ± 1 standard deviation across all APs.
V_STA_std = V[spike_select_idx].std(1)
t = window * dt * 1e3
plot(t, V_STA, 'k')
fill_between(x = t,
y1 = V_STA + V_STA_std,
y2 = V_STA - V_STA_std,
facecolor='blue', alpha=0.4)
xlabel('time [ms]')
ylabel('V_mem [mV]')
Out[7]:
FOUR
What is the mean and variance of the membrane potential distribution when excluding the action potential cutouts?
In [8]:
# get the exlusive of all cut-out windows and the indices of V
no_ap_idx = setxor1d(spike_select_idx.flatten(),
indices(V.shape))
V_no_ap = V[no_ap_idx]
plot(V_no_ap[0:2/dt])
#make a histogram plot of the membrane potential distribution
figure(); hist(V_no_ap)
print V_no_ap.mean(),'±' ,V_no_ap.std()
The standard integrate-and-fire (I&F) model of a neuron assumes a fixed voltage threshold for initiating a spike. Can you determine a fixed voltage threshold that determines the onset of each spike? Repeat exercise 3 but now with a longer cut-out before the spike detection of 500ms.
Plot the mean action potential in red and in addition plot two functions as mean ± standard deviation.
In [9]:
# please edit the data path if necessary
data = scipy.io.loadmat('data/spontaneous_recording_nostruct.mat',
squeeze_me=True)
I = data['I']
FI = data['FI']
dt = data['TimeResolutionS']
print I.shape, FI.shape, dt
TWO
plot the raw signal (I) and the band-pass filtered signal (FI) in 2 different colors on top of each other (graph 1). Use an appropriate time scale (s or ms) on the x-axis and choose a time window that allows to observe the structure of the data (e.g. 1s length).
THREE
use a threshold to detect PSCs in the band-pass filtered data. The threshold should be chosen such that only obvious synaptic events cross the threshold, while noise fluctuations do not reach this threshold. Add the horizontal threshold line to the previous graph 1.
Assign a variable “onset_idx” to the vector indices at which threshold crossing occurred. How many PSCs did you detect? How many did your neighbors detect? What is the reason for differing numbers?
In [10]:
window = slice(0,1/dt)
threshold = 4
In [11]:
plot(I[window])
plot(FI[window], 'k', linewidth=2)
axhline(y=threshold, linestyle='--', color='k')
Out[11]:
In [12]:
# detecting the trhreshold indices explained in part I.2
I_above_threshold = I>threshold
I_above_th_diff = diff(I_above_threshold.astype(int)) == 1
EPSC_bins = I_above_th_diff == 1
onset_idx = find(EPSC_bins)
How many PSCs did you detect? How many did your neighbors detect? What is the reason for differing numbers?
In [13]:
N_EPSC = len(onset_idx)
print N_EPSC
FOUR
Make data cutouts of stretches of the raw signal around the onsets of the detected PSCs. Choose a reasonable time interval before and after the onset such that the time course of the PSCs is well captured. Assign the variable “PSCs” to the matrix of cutouts.
In [14]:
# couting out data around an event is exokained in Part I.3
cutout = int(1e-3/dt)
window = arange(-cutout,+cutout)
EPSC_select_idx = tile(window,(N_EPSC,1)).T
EPSC_select_idx += onset_idx
EPSC_select_idx = EPSC_select_idx.clip(0,len(I)-1)
print window.shape, N_EPSC, EPSC_select_idx.shape
FIVE
Plot the average PSC and two additional lines indicating the standard deviation in positive and negative direction with respect to the average PSC.
What is the peak amplitude of the avg. PSC?
In [15]:
I_EPSC = I[EPSC_select_idx]
I_EPSC_mean = I_EPSC.mean(1)
I_EPSC_std = I_EPSC.std(1)
In [16]:
# Plot all PSCs in one axis. On top plot the average PSC with a stronger line.
# Well, plotting over 100k EPSCs take forever. plot only the first 10 maybe
print I_EPSC.shape
plot(I_EPSC[:,:10], color='grey')
plot(I_EPSC_mean, linewidth=2)
Out[16]:
In [17]:
t = window * dt * 1e3
plot(t, I_EPSC_mean, 'k')
fill_between(x = t,
y1 = I_EPSC_mean + I_EPSC_std,
y2 = I_EPSC_mean - I_EPSC_std,
facecolor='blue', alpha=0.4)
xlabel('time [ms]')
ylabel('I_mem [?]')
print 'The peak amplitude of average EPSC is:',I_EPSC_mean.max()
SIX
Statistics: Determine the peak amplitude of each single PSC and produce a histogram (i.e. show the distribution) of PSC peak amplitudes (Figure 8). Is this distribution symmetric or skewed?
Produce a histogram on the logarithmic axis, or, equivalently, produce a histogram of the logarithmic values of peak amplitude. Distribution of synaptic strength should be well approximated by the so-called log-normal distribution (i.e. normally distributed on the logarithmic axis).
In [18]:
peaks_EPSC = I_EPSC.max(0)
p,bins,_ = hist(peaks_EPSC)
title('Distribution of peak EPSC amplitudes')
Out[18]:
In [19]:
p_log,bins,_ = hist(log(peaks_EPSC))
title('Distribution of peak EPSC amplitudes (log scale)')
Out[19]: