%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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