689 lines
32 KiB
TeX
689 lines
32 KiB
TeX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\chapter{Spike-train analysis}
|
|
\label{pointprocesseschapter}
|
|
\exercisechapter{Spike-train analysis}
|
|
|
|
\entermde[action potential]{Aktionspotential}{Action potentials}
|
|
(\enterm[spike|seealso{action potential}]{spikes}) carry information
|
|
within neural systems. More precisely, the times at which action
|
|
potentials are generated contain the information. The waveform of the
|
|
action potential is largely stereotyped and therefore conveys no
|
|
information. Analyzing the statistics of spike times and their
|
|
relation to sensory stimuli or motor actions is central to
|
|
neuroscientific research. With multi-electrode arrays it is nowadays
|
|
possible to record from hundreds or even thousands of neurons
|
|
simultaneously. The open challenge is how to analyze such huge data
|
|
sets in smart ways in order to gain insights into the way neural
|
|
systems work. Let's start with the basics in this chapter.
|
|
|
|
The result of the pre-processing of electrophysiological recordings
|
|
are series of spike times, which are termed \enterm[spike train]{spike
|
|
trains}. If measurements are repeated we get several \enterm{trials}
|
|
of spike trains (\figref{rasterexamplesfig}). Spike trains are lists
|
|
of times of events, the action potentials. Analyzing spike trains
|
|
leads into the realm of the statistics of so called \entermde[point
|
|
process]{Punktprozess}{point processes}.
|
|
|
|
\begin{figure}[bt]
|
|
\includegraphics[width=1\textwidth]{rasterexamples}
|
|
\titlecaption{\label{rasterexamplesfig}Raster plots of spike
|
|
trains.}{Raster plots of ten trials of data illustrating the times
|
|
of action potentials. Each vertical stroke illustrates the time at
|
|
which an action potential was observed. Each row displays the
|
|
events of one trial. Shown is a stationary point process
|
|
(homogeneous point process with a rate $\lambda=20$\;Hz, left) and
|
|
an non-stationary point process with a rate that varies in time
|
|
(noisy perfect integrate-and-fire neuron driven by
|
|
Ornstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms,
|
|
right).}
|
|
\end{figure}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Point processes}
|
|
|
|
\begin{ibox}{Examples of point processes}
|
|
Every point process is generated by a temporally continuously
|
|
developing process. An event is generated whenever this process
|
|
crosses some threshold. For example:\vspace{-1ex}
|
|
\begin{itemize}
|
|
\item Action potentials/heart beat: created by the dynamics of the
|
|
neuron/sinoatrial node.
|
|
\item Earthquake: defined by the dynamics of the pressure between
|
|
tectonical plates.
|
|
\item Communication calls in crickets/frogs/birds: shaped by
|
|
the dynamics of the nervous system and the muscle apparatus.
|
|
\end{itemize}
|
|
\end{ibox}
|
|
|
|
\begin{figure}[tb]
|
|
\includegraphics{pointprocesssketch}
|
|
\titlecaption{\label{pointprocesssketchfig} Statistics of point
|
|
processes.}{A temporal point process is a sequence of events in
|
|
time, $t_i$, that can be also characterized by the corresponding
|
|
inter-event intervals $T_i=t_{i+1}-t_i$ and event counts $n_i$,
|
|
i.e. the number of events that occurred so far.}
|
|
\end{figure}
|
|
|
|
\noindent
|
|
A temporal \entermde{Punktprozess}{point process} is a stochastic
|
|
process that generates a sequence of events at times $\{t_i\}$. In
|
|
the neurosciences, the statistics of point processes is of importance
|
|
since the timing of neuronal events (action potentials, post-synaptic
|
|
potentials, events in EEG or local-field recordings, etc.) is crucial
|
|
for information transmission and can be treated as such a process.
|
|
|
|
The events of a point process can be illustrated by means of a raster
|
|
plot in which each vertical line indicates the time of an event. The
|
|
event from two different point processes are shown in
|
|
\figref{rasterexamplesfig}. In addition to the event times, point
|
|
processes can be described using the intervals $T_i=t_{i+1}-t_i$
|
|
between successive events or the number of observed events $n_i$
|
|
within a certain time window (\figref{pointprocesssketchfig}).
|
|
|
|
In \enterm[point process!stationary]{stationary} point processes the
|
|
statistics does not change over time. In particular, the rate of the
|
|
process, the number of events per time, is constant. The homogeneous
|
|
point process shown in \figref{rasterexamplesfig} on the left is an
|
|
example of a stationary point process. Although locally within each
|
|
trial there are regions with many events per time and others with long
|
|
intervals between events, the average number of events within a small
|
|
time window over trials does not change in time. In the first sections
|
|
of this chapter we introduce various statistics for characterizing
|
|
stationary point processes.
|
|
|
|
On the other hand, the example shown in \figref{rasterexamplesfig} on
|
|
the right is a non-stationary point process. The rate of the events in
|
|
all trials first decreases and then increases again. This common
|
|
change in rate may have been caused by some sensory input to the
|
|
neuron. How to estimate temporal changes in event rates (firing rates)
|
|
is covered in section~\ref{nonstationarysec}.
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Homogeneous Poisson process}
|
|
|
|
Before we are able to start analyzing point processes, we need some
|
|
data. As before, because we are working with computers, we can easily
|
|
simulate them. To get us started we here briefly introduce the
|
|
homogeneous Poisson process. The Poisson process is to point processes
|
|
what the Gaussian distribution is to the statistics of real-valued
|
|
data. It is a standard against which everything is compared.
|
|
|
|
In a \entermde[Poisson
|
|
process!homogeneous]{Poissonprozess!homogener}{homogeneous Poisson
|
|
process} events at every given time, that is within every small time
|
|
window of width $\Delta t$ occur with the same constant
|
|
probability. This probability is independent of absolute time and
|
|
independent of any events occurring before. To observe an event right
|
|
after an event is as likely as to observe an event at some specific
|
|
time later on.
|
|
|
|
The simplest way to simulate events of a homogeneous Poisson process
|
|
is based on the exponential distribution \eqref{hompoissonexponential}
|
|
of event intervals. We randomly draw intervals from this distribution
|
|
and then sum them up to convert them to event times. The only
|
|
parameter of a homogeneous Poisson process is its rate. It defines how
|
|
many events per time are expected.
|
|
|
|
Here is a function that generates several trials of a homogeneous
|
|
Poisson process:
|
|
|
|
\lstinputlisting[caption={hompoissonspikes.m}]{hompoissonspikes.m}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Raster plot}
|
|
|
|
Let's generate some events with this function and display them in a
|
|
raster plot as in \figref{rasterexamplesfig}. For the raster plot we
|
|
need to draw for each event a line at the corresponding time of the
|
|
event and at the height of the corresponding trial. You can go with a
|
|
one for-loop through the trials and then with another for-loop through
|
|
the event times and every time plot a two-point line. This, however,
|
|
is very slow. The fastest way is to concatenate the coordinates of all
|
|
strokes into large vectors, separate the events by \varcode{nan}
|
|
entries, and pass this to a single call of \varcode{plot()}:
|
|
|
|
\pageinputlisting[caption={rasterplot.m}]{rasterplot.m}
|
|
|
|
Adapt this function to your needs and use it where ever possible to
|
|
illustrate your spike train data to your readers. They appreciate
|
|
seeing your raw data and being able to judge the data for themselves
|
|
before you go on analyzing firing rates, interspike-interval
|
|
correlations, etc.
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Interval statistics}
|
|
|
|
Let's start with the interval statistics of stationary point
|
|
processes. The intervals $T_i=t_{i+1}-t_i$ between successive events
|
|
are real positive numbers. In the context of action potentials they
|
|
are referred to as \entermde[interspike
|
|
interval]{Interspikeintervall}{interspike intervals}, in short
|
|
\entermde[ISI|see{interspike
|
|
interval}]{ISI|see{Interspikeintervall}}{ISI}s. For analyzing event
|
|
intervals we can use all the usual statistics that we know from
|
|
describing univariate data sets of real numbers:
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics[width=0.96\textwidth]{isihexamples}
|
|
\titlecaption{\label{isihexamplesfig}Interspike-interval
|
|
histograms}{of the spike trains shown in
|
|
\figref{rasterexamplesfig}. The intervals of the homogeneous
|
|
Poisson process (left) are exponentially distributed according to
|
|
\eqnref{hompoissonexponential} (red). Typically for many sensory
|
|
neurons that fire more regularly is an ISI histogram like the one
|
|
shown on the right. There is a preferred interval where the
|
|
distribution peaks. The distribution falls off quickly towards
|
|
smaller intervals, and very small intervals are absent, probably
|
|
because of refractoriness of the spike generator. The tail of the
|
|
distribution again approaches an exponential distribution as for
|
|
the Poisson process.}
|
|
\end{figure}
|
|
|
|
\begin{exercise}{isis.m}{}
|
|
Implement a function \varcode{isis()} that calculates the interspike
|
|
intervals from several spike trains. The function should return a
|
|
single vector of intervals. The spike times (in seconds) of each
|
|
trial are stored as vectors within a cell-array.
|
|
\end{exercise}
|
|
|
|
\begin{itemize}
|
|
\item Probability density $p(T)$ of the intervals $T$
|
|
(\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \;
|
|
dT = 1$. Commonly referred to as the \enterm[interspike
|
|
interval!histogram]{interspike interval histogram}. Its shape
|
|
reveals many interesting aspects like locking or bursting that
|
|
cannot be inferred from the mean or standard deviation. A particular
|
|
reference is the exponential distribution of intervals
|
|
\begin{equation}
|
|
\label{hompoissonexponential}
|
|
p_{exp}(T) = \lambda e^{-\lambda T}
|
|
\end{equation}
|
|
of a homogeneous Poisson spike train with rate $\lambda$.
|
|
\item Mean interval
|
|
\begin{equation}
|
|
\label{meanisi}
|
|
\mu_{ISI} = \langle T \rangle = \frac{1}{n}\sum_{i=1}^n T_i
|
|
\end{equation}
|
|
The average time it takes from one event to the next. For stationary
|
|
point processes the inverse of the mean interval is identical with
|
|
the mean rate $\lambda$ (number of events per time, see below) of
|
|
the process.
|
|
\item Standard deviation of intervals
|
|
\begin{equation}
|
|
\label{stdisi}
|
|
\sigma_{ISI} = \sqrt{\langle (T - \langle T \rangle)^2 \rangle}
|
|
\end{equation}
|
|
Periodically spiking neurons have little variability in their
|
|
intervals, whereas many cortical neurons cover a wide range with
|
|
their intervals. The standard deviation of homogeneous Poisson spike
|
|
trains also equals the inverse rate. Whether the standard deviation
|
|
of intervals is low or high, however, is better quantified by the
|
|
\item \entermde[coefficient of
|
|
variation]{Variationskoeffizient}{Coefficient of variation}, the
|
|
standard deviation of the ISIs relative to their mean:
|
|
\begin{equation}
|
|
\label{cvisi}
|
|
CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}
|
|
\end{equation}
|
|
Homogeneous Poisson spike trains have an $CV$ of exactly one. The
|
|
lower the $CV$ the more regularly a neuron is firing. $CV$s larger than
|
|
one are also possible in spike trains with small intervals separated
|
|
by really long ones.
|
|
%\item \entermde[diffusion coefficient]{Diffusionskoeffizient}{Diffusion coefficient}: $D_{ISI} =
|
|
% \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
|
|
\end{itemize}
|
|
|
|
\begin{exercise}{isihist.m}{}
|
|
Implement a function \varcode{isihist()} that calculates the
|
|
normalized interspike interval histogram. The function should take
|
|
two input arguments; (i) a vector of interspike intervals and (ii)
|
|
the width of the bins used for the histogram. It returns the
|
|
probability density as well as the centers of the bins.
|
|
\end{exercise}
|
|
|
|
\begin{exercise}{plotisihist.m}{}
|
|
Implement a function that takes the returned values of
|
|
\varcode{isihist()} as input arguments and then plots the data. The
|
|
plot should show the histogram with the x-axis scaled to
|
|
milliseconds and should be annotated with the average ISI, the
|
|
standard deviation, and the coefficient of variation of the ISIs
|
|
(\figref{isihexamplesfig}).
|
|
\end{exercise}
|
|
|
|
\subsection{Interval correlations}
|
|
Intervals are not just numbers without an order, like weights of
|
|
tigers. Intervals are temporally ordered and there could be temporal
|
|
structure in the sequence of intervals. For example, short intervals
|
|
could be followed by more longer ones, and vice versa. Such
|
|
dependencies in the sequence of intervals do not show up in the
|
|
interval histogram. We need additional measures to also quantify the
|
|
temporal structure of the sequence of intervals.
|
|
|
|
We can use the same techniques we know for visualizing and quantifying
|
|
correlations in bivariate data sets, i.e. scatter plots and
|
|
correlation coefficients. We form $(x,y)$ data pairs by taking the
|
|
series of intervals $T_i$ as $x$-data values and pairing them with the
|
|
$k$-th next intervals $T_{i+k}$ as $y$-data values. The parameter $k$
|
|
is called \enterm{lag} (\determ{Verz\"ogerung}). For lag one we pair
|
|
each interval with the next one. A \entermde{return map}{return map}
|
|
illustrates dependencies between successive intervals by simply
|
|
plotting the intervals $T_{i+k}$ against the intervals $T_i$ in a
|
|
scatter plot (\figref{returnmapfig}). For Poisson spike trains there
|
|
is no structure beyond the one expected from the exponential
|
|
interspike interval distribution, hinting at neighboring interspike
|
|
intervals being independent of each other. For the spike train based
|
|
on an Ornstein-Uhlenbeck process the return map is more clustered
|
|
along the diagonal, hinting at a positive correlation between
|
|
succeeding intervals. That is, short intervals are more likely to be
|
|
followed by short ones and long intervals more likely by long
|
|
ones. This temporal structure was already clearly visible in each
|
|
trial of the spike raster shown in \figref{rasterexamplesfig} on the
|
|
right.
|
|
|
|
\begin{figure}[tp]
|
|
\includegraphics[width=1\textwidth]{serialcorrexamples}
|
|
\titlecaption{\label{returnmapfig}Interspike-interval
|
|
correlations}{of the spike trains shown in
|
|
\figref{rasterexamplesfig}. Upper panels show the return maps and
|
|
lower panels the serial correlations of successive intervals
|
|
separated by lag $k$. All the interspike intervals of the
|
|
stationary spike trains are independent of each other --- this is
|
|
a so called \enterm{renewal process}
|
|
(\determ{Erneuerungsprozess}). In contrast, the ones of the
|
|
non-stationary spike trains show positive correlations that decay
|
|
for larger lags. The positive correlations in this example are
|
|
caused by a common stimulus that slowly increases and decreases
|
|
the mean firing rate of the spike trains.}
|
|
\end{figure}
|
|
|
|
Such dependencies can be further quantified by
|
|
\entermde[correlation!serial]{Korrelation!serielle}{serial
|
|
correlations}. These quantify the correlations between successive
|
|
intervals by Pearson's correlation coefficients between the intervals
|
|
$T_{i+k}$ and $T_i$ in dependence on lag $k$:
|
|
\begin{equation}
|
|
\label{serialcorrelation}
|
|
\rho_k = \frac{\langle (T_{i+k} - \langle T \rangle)(T_i - \langle T \rangle) \rangle}{\langle (T_i - \langle T \rangle)^2\rangle} = \frac{{\rm cov}(T_{i+k}, T_i)}{{\rm var}(T_i)}
|
|
= {\rm corr}(T_{i+k}, T_i)
|
|
\end{equation}
|
|
The serial correlations $\rho_k$ are usually plotted against the lag
|
|
$k$ for a small range of lags (\figref{returnmapfig}). $\rho_0=1$ is
|
|
the correlation of each interval with itself and always equals one.
|
|
|
|
If the serial correlations all equal zero, $\rho_k =0$ for $k>0$, then
|
|
the length of an interval is independent of all the previous
|
|
ones. Such a process is a \enterm{renewal process}
|
|
(\determ{Erneuerungsprozess}). Each event, each action potential,
|
|
erases the history. The occurrence of the next event is independent of
|
|
what happened before. To a first approximation an action potential
|
|
erases all information about the past from the membrane voltage and
|
|
thus spike trains may approximate renewal processes.
|
|
|
|
However, other variables like the intracellular calcium concentration
|
|
or the states of slowly switching ion channels may carry information
|
|
from one interspike interval to the next and thus introduce
|
|
correlations between intervals. Such non-renewal dynamics is then
|
|
characterized by the non-zero serial correlations
|
|
(\figref{returnmapfig}).
|
|
|
|
\begin{exercise}{isiserialcorr.m}{}
|
|
Implement a function \varcode{isiserialcorr()} that takes a vector
|
|
of interspike intervals as input argument and calculates the serial
|
|
correlations up to some maximum lag.
|
|
\end{exercise}
|
|
|
|
\begin{exercise}{plotisiserialcorr.m}{}
|
|
Implement a function \varcode{plotisiserialcorr()} that takes a
|
|
vector of interspike intervals as input argument and generates a
|
|
plot of the serial correlations.
|
|
\end{exercise}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Count statistics}
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics{countexamples}
|
|
\titlecaption{\label{countstatsfig}Count statistics.}{Probability
|
|
distributions of counting $k$ events $k$ (blue) within windows of
|
|
20\,ms (left) or 200\,ms duration (right) of a homogeneous Poisson
|
|
spike train with a rate of 20\,Hz
|
|
(\figref{rasterexamplesfig}). For Poisson spike trains these
|
|
distributions are given by Poisson distributions (red).}
|
|
\end{figure}
|
|
|
|
The most commonly used measure for characterizing spike trains is the
|
|
\enterm[firing rate!average]{average firing rate}. The firing rate $r$
|
|
is the average number of spikes counted within some time interval $W$
|
|
\begin{equation}
|
|
\label{firingrate}
|
|
r = \frac{\langle n \rangle}{W}
|
|
\end{equation}
|
|
and is measured in Hertz. The average of the spike counts is taken
|
|
over trials. For stationary spike trains (no change in statistics, in
|
|
particular the firing rate, over time), the firing rate based on the
|
|
spike count equals the inverse average interspike interval
|
|
$1/\mu_{ISI}$.
|
|
|
|
The firing rate based on an averaged spike counts is one example of
|
|
many statistics based on event counts. Stationary spike trains can be
|
|
split into many segments $i$, each of duration $W$, and the number of
|
|
events $n_i$ in each of the segments can be counted. The integer event
|
|
counts can be quantified in the usual ways:
|
|
\begin{itemize}
|
|
\item Histogram of the counts $n_i$ appropriately normalized to
|
|
probability distributions. For homogeneous Poisson spike trains with
|
|
rate $\lambda$ the resulting probability distributions follow a
|
|
Poisson distribution (\figref{countstatsfig}), where the probability
|
|
of counting $k$ events within a time window $W$ is given by
|
|
\begin{equation}
|
|
\label{poissondist}
|
|
P(k) = \frac{(\lambda W)^k e^{\lambda W}}{k!}
|
|
\end{equation}
|
|
\item Average number of counts: $\mu_n = \langle n \rangle$.
|
|
\item Variance of counts:
|
|
$\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
|
|
\end{itemize}
|
|
Because spike counts are unitless and positive numbers the
|
|
\begin{itemize}
|
|
\item \entermde{Fano Faktor}{Fano factor} (variance of counts divided
|
|
by average count)
|
|
\begin{equation}
|
|
\label{fano}
|
|
F = \frac{\sigma_n^2}{\mu_n}
|
|
\end{equation}
|
|
is a commonly used measure for quantifying the variability of event
|
|
counts relative to the mean number of events. In particular for
|
|
homogeneous Poisson processes the Fano factor equals one,
|
|
independently of the time window $W$.
|
|
\end{itemize}
|
|
|
|
Note that all of these statistics depend in general on the chosen
|
|
window length $W$. The average spike count, for example, grows
|
|
linearly with $W$ for sufficiently large time windows: $\langle n
|
|
\rangle = r W$, \eqnref{firingrate}. Doubling the counting window
|
|
doubles the spike count. As does the spike-count variance
|
|
(\figref{fanofig}). At smaller time windows the statistics of the
|
|
event counts might depend on the particular duration of the counting
|
|
window. There might be an optimal time window for which the variance
|
|
of the spike count is minimal. The Fano factor plotted as a function
|
|
of the time window illustrates such properties of point processes
|
|
(\figref{fanofig}).
|
|
|
|
This also has consequences for information transmission in neural
|
|
systems. The lower the variance in spike count relative to the
|
|
averaged count, the higher the signal-to-noise ratio at which
|
|
information encoded in the mean spike count is transmitted.
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics{fanoexamples}
|
|
\titlecaption{\label{fanofig} Fano factor.}{Counting events in time
|
|
windows of given duration and then dividing the variance of the
|
|
counts by their mean results in the Fano factor. Here, the Fano
|
|
factor is plotted as a function of the duration of the window used
|
|
to count events. For Poisson spike trains the variance always
|
|
equals the mean counts and consequently the Fano factor equals one
|
|
irrespective of the count window (left). A spike train with
|
|
positive correlations between interspike intervals (caused by an
|
|
Ornstein-Uhlenbeck process) has a minimum in the Fano factor, that
|
|
is an analysis window for which the relative count variance is
|
|
minimal somewhere close to the correlation time scale of the
|
|
interspike intervals (right).}
|
|
\end{figure}
|
|
|
|
\begin{exercise}{counthist.m}{}
|
|
Implement a function \varcode{counthist()} that calculates and plots
|
|
the distribution of spike counts observed in a certain time
|
|
window. The function should take two input arguments: a cell-array
|
|
of vectors containing the spike times in seconds observed in a
|
|
number of trials, and the duration of the time window that is used
|
|
to evaluate the counts.
|
|
\end{exercise}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Time-dependent firing rate}
|
|
\label{nonstationarysec}
|
|
|
|
So far we have discussed stationary spike trains. The statistical properties
|
|
of these did not change within the observation time (stationary point
|
|
processes). Most commonly, however, this is not the case. A sensory
|
|
neuron, for example, might respond to a stimulus by modulating its
|
|
firing rate (non-stationary point process).
|
|
|
|
How the firing rate $r(t)$ changes over time is the most important
|
|
measure, when analyzing non-stationary spike trains. The unit of the
|
|
firing rate is Hertz, i.e. the number of action potentials per
|
|
second. There are different ways to estimate the firing rate and three
|
|
of these are illustrated in \figref{psthfig}. All have their own
|
|
justifications, their pros- and cons.
|
|
|
|
\begin{figure}[tp]
|
|
\includegraphics[width=\columnwidth]{firingrates}
|
|
\titlecaption{Estimating time-dependent firing rates.}{Rasterplot
|
|
showing one trial of spiking activity of a neuron (top).
|
|
\emph{Instantaneous rate}, classical PSTH estimatd with the
|
|
\emph{binning} method and by \emph{convolution} of the spike train
|
|
with a Gaussian kernel (bottom).}\label{psthfig}
|
|
\end{figure}
|
|
|
|
|
|
\subsection{Instantaneous firing rate}
|
|
|
|
\begin{figure}[tp]
|
|
\includegraphics[width=\columnwidth]{isimethod}
|
|
\titlecaption{Instantaneous firing rate.}{The recorded spike train
|
|
(top). Arrows illustrate the interspike intervals and numbers
|
|
give the intervals in milliseconds. The inverse of the interspike
|
|
intervals is the \emph{instantaneous firing rate}
|
|
(bottom).}\label{instratefig}
|
|
\end{figure}
|
|
|
|
A very simple method for estimating the time-dependent firing rate is
|
|
the \entermde[firing rate!instantaneous]{Feuerrate!instantane}{instantaneous firing rate}.
|
|
The firing rate can be directly estimated as the inverse of the time
|
|
between successive spikes, the interspike-interval (\figref{instratefig}).
|
|
|
|
\begin{equation}
|
|
\label{instantaneousrateeqn}
|
|
r_i = \frac{1}{T_i} .
|
|
\end{equation}
|
|
|
|
The instantaneous rate $r_i$ is valid for the whole interspike
|
|
interval. The method has the advantage of being extremely easy to
|
|
compute and that it does not make any assumptions about the relevant
|
|
timescale (of the encoding in the neuron or the decoding of a
|
|
postsynaptic neuron). The resulting $r(t)$, however, is no continuous
|
|
function, the firing rate jumps from one level to the next. Since the
|
|
interspike interval between successive spikes is never infinitely
|
|
long, the firing rate never reaches zero despite that the neuron may
|
|
not fire an action potential for a long time.
|
|
|
|
\begin{exercise}{instantaneousRate.m}{}
|
|
Implement a function that computes the instantaneous firing
|
|
rate. Plot the firing rate as a function of time.
|
|
%\note{TODO: example data!!!}
|
|
\end{exercise}
|
|
|
|
|
|
\subsection{Peri-stimulus-time-histogram}
|
|
While the instantaneous firing rate is based on the interspike
|
|
intervals, the \enterm{peri stimulus time histogram} (PSTH) is based on
|
|
spike counts within observation windows of the duration $W$. It
|
|
estimates the probability of observing a spike within that observation
|
|
time. It tries to estimate the average rate in the limit of small
|
|
obersvation times:
|
|
\begin{equation}
|
|
\label{psthrate}
|
|
r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; ,
|
|
\end{equation}
|
|
where $\langle n \rangle$ is the across trial average number of action
|
|
potentials observed within the interval $(t, t+W)$. Such description
|
|
matches the time-dependent firing rate $\lambda(t)$ of an
|
|
inhomogeneous Poisson process.
|
|
|
|
The firing probability can be estimated using the \enterm[firing
|
|
rate!binning method]{binning method} or by using \enterm[firing
|
|
rate!kernel density estimation]{kernel density estimations}. Both
|
|
methods make an assumption about the relevant observation time-scale
|
|
($W$ in \eqnref{psthrate}).
|
|
|
|
\subsubsection{Binning-method}
|
|
|
|
\begin{figure}[tp]
|
|
\includegraphics[width=\columnwidth]{binmethod}
|
|
\titlecaption{Estimating the PSTH using the binning method.}{The
|
|
same spike train as shown in \figref{instratefig} (top). Vertical
|
|
gray lines indicate the borders between adjacent bins in which the
|
|
number of action potentials is counted (red numbers). The firing
|
|
rate is then the histogram normalized to the binwidth
|
|
(bottom).}\label{binpsthfig}
|
|
\end{figure}
|
|
|
|
The \entermde[firing rate!binning
|
|
method]{Feuerrate!Binningmethode}{binning method} separates the time
|
|
axis into regular bins of the bin width $W$ and counts for each bin
|
|
the number of observed action potentials (\figref{binpsthfig}
|
|
top). The resulting histogram is then normalized with the bin width
|
|
$W$ to yield the firing rate shown in the bottom trace of figure
|
|
\ref{binpsthfig}. The above sketched process is equivalent to
|
|
estimating the probability density. For computing a PSTH the
|
|
\code{hist()} function can be used.
|
|
|
|
The estimated firing rate is valid for the total duration of each
|
|
bin. This leads to the step-like plot shown in
|
|
\figref{binpsthfig}. $r(t)$ is thus not a contiunous function in
|
|
time. The binwidth defines the temporal resolution of the firing rate
|
|
estimation Changes that happen within a bin cannot be resolved. Thus
|
|
chosing a bin width implies an assumption about the relevant
|
|
time-scale.
|
|
|
|
\begin{exercise}{binnedRate.m}{}
|
|
Implement a function that estimates the firing rate using the
|
|
binning method. The method should take the spike-times as an
|
|
input argument and returns the firing rate. Plot the PSTH.
|
|
\end{exercise}
|
|
|
|
\subsubsection{Convolution method --- Kernel density estimation}
|
|
|
|
\begin{figure}[tp]
|
|
\includegraphics[width=\columnwidth]{convmethod}
|
|
\titlecaption{Estimating the firing rate using the convolution
|
|
method.}{The same spike train as in \figref{instratefig} (top). The
|
|
convolution of the spike train with a kernel replaces each spike
|
|
event with the kernel (red). A Gaussian kernel is used here, but
|
|
other kernels are also possible. If the kernel is properly
|
|
normalized the firing rate results directly form the superposition
|
|
of the kernels (bottom).}\label{convratefig}
|
|
\end{figure}
|
|
|
|
With the \entermde[firing rate!convolution
|
|
method]{Feuerrate!Faltungsmethode}{convolution method} we avoid the
|
|
sharp edges of the binning method. The spiketrain is convolved with a
|
|
\entermde{Faltungskern}{convolution kernel}. Technically speaking we
|
|
need to first create a binary representation of the spike train. This
|
|
binary representation is a series of zeros and ones in which the ones
|
|
denote the spike. Then this binary vector is convolved with a kernel
|
|
of a certain width:
|
|
\[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \]
|
|
where $\omega(\tau)$ represents the kernel and $\rho(t)$ the binary
|
|
representation of the response. The process of convolution can be
|
|
imagined as replacing each event of the spiketrain with the kernel
|
|
(figure \ref{convratefig} top). The superposition of the replaced
|
|
kernels is then the firing rate (if the kernel is correctly normalized
|
|
to an integral of one, figure \ref{convratefig}
|
|
bottom).
|
|
|
|
In contrast to the other methods the convolution methods leads to a
|
|
continuous function which is often desirable (in particular when
|
|
applying methods in the frequency domain). The choice of the kernel
|
|
width defines, similar to the bin width of the binning method, the
|
|
temporal resolution of the method and thus makes assumptions about the
|
|
relevate time-scale.
|
|
|
|
\begin{exercise}{convolutionRate.m}{}
|
|
Implement the function that estimates the firing rate using the
|
|
convolution method. The method takes the spiketrain, the temporal
|
|
resolution of the recording (as the stepsize $dt$, in seconds) and
|
|
the width of the kernel (the standard deviation $\sigma$ of the
|
|
Gaussian kernel, in seconds) as input arguments. It returns the
|
|
firing rate. Plot the result.
|
|
\end{exercise}
|
|
|
|
\section{Spike-triggered Average}
|
|
\label{stasec}
|
|
|
|
The graphical representation of the neuronal activity alone is not
|
|
sufficient tot investigate the relation between the neuronal response
|
|
and a stimulus. One method to do this is the \entermde{Spike-triggered
|
|
Average}{spike-triggered average}, \enterm[STA|see{spike-triggered
|
|
average}]{STA}. The STA
|
|
\begin{equation}
|
|
STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau)
|
|
\end{equation}
|
|
of $N$ action potentials observed at the times $t_i$ in response to
|
|
the stimulus $s(t)$ is the average stimulus that led to a spike in the
|
|
neuron. The STA can be easily extracted by cutting snippets out of
|
|
$s(t)$ that surround the times of the spikes. The resulting stimulus
|
|
snippets are then averaged (\figref{stafig}).
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics[width=\columnwidth]{sta}
|
|
\titlecaption{Spike-triggered average of a P-type electroreceptor
|
|
and the stimulus reconstruction.}{The neuron was driven by a
|
|
\enterm{white-noise} stimulus (blue, right). The STA (left) is the
|
|
average stimulus that surrounds the times of the recorded action
|
|
potentials (40\,ms before and 20\,ms after the spike). Using the
|
|
STA as a kernel for convolving the spiketrain we can reconstruct
|
|
the stimulus from the neuronal response. In this way we can get an
|
|
impression of the stimulus features that are linearly encoded in
|
|
the neuronal response (orange, right).}\label{stafig}
|
|
\end{figure}
|
|
|
|
From the STA we can extract several pieces of information about the
|
|
relation of stimulus and response. The width of the STA represents the
|
|
temporal precision with which the neuron encodes the stimulus waveform
|
|
and in which temporal window the neuron integrates the (sensory)
|
|
input. The amplitude of the STA tells something about the sensitivity
|
|
of the neuron. The STA is given in the same units as the stimulus and
|
|
a small amplitude indicates that the neuron needs only a small
|
|
stimulus amplitude to create a spike, a large amplitude on the
|
|
contrary suggests the opposite. The temporal delay between the STA and
|
|
the time of the spike is a consequence of the time the system (neuron)
|
|
needs to process the stimulus.
|
|
|
|
We can further use the STA to do a \enterm{reverse reconstruction} and
|
|
estimate the stimulus from the neuronal response (\figref{stafig},
|
|
right). For this, the spiketrain is convolved with the STA as a
|
|
kernel.
|
|
|
|
\begin{exercise}{spikeTriggeredAverage.m}{}
|
|
Implement a function that calculates the STA. Use the dataset
|
|
\file{sta\_data.mat}. The function expects the spike train, the
|
|
stimulus and the temporal resolution of the recording as input
|
|
arguments and should return the following information:
|
|
\vspace{-1ex}
|
|
\begin{itemize}
|
|
\setlength{\itemsep}{0ex}
|
|
\item the spike-triggered average.
|
|
\item the standard deviation of the STA across the individual snippets.
|
|
\item The number of action potentials used to estimate the STA.
|
|
\end{itemize}
|
|
|
|
\end{exercise}
|
|
|
|
\begin{exercise}{reconstructStimulus.m}{}
|
|
Do the reverse reconstruction using the STA and the spike times. The
|
|
function should return the estimated stimulus in a vector that has
|
|
the same size as the original stimulus contained in file
|
|
\file{sta\_data.mat}.
|
|
\end{exercise}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\printsolutions
|