scientificComputing/pointprocesses/lecture/pointprocesses.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Spiketrain analysis}

\selectlanguage{english}

\enterm[Actionspotentials]{Actionspotentials} (\enterm{spikes}) are
the carriers of information in the nervous system. Thereby it is
mainly the time at which the spikes are generated that is of
importance. The waveform of the action potential is largely
stereotyped and does not carry information.

The result of the processing of electrophysiological recordings are
series of spike times, which are then termed \enterm{spiketrains}. If
measurements are repeated we yield several \enterm{trials} of
spiketrains (\figref{rasterexamplesfig}).

Spiketrains are times of events, the action potentials. The analysis
of these leads into the realm of the so called \enterm[point
  process]{point processes}.

\begin{figure}[ht]
  \includegraphics[width=1\textwidth]{rasterexamples}
  \titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plot of
    ten realizations of a stationary point process (homogeneous point
    process with a rate $\lambda=20$\;Hz, left) and an inhomogeneous
    point process (perfect integrate-and-fire neuron dirven by
    Ohrnstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms,
    right). Each vertical dash illustrates the time at which the
    action potential was observed. Each line represents the event of
    each trial.}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Point processes}

A temporal \enterm{point process} is a stochastic process that
generates a sequence of events at times $\{t_i\}$, $t_i \in
\reZ$.

\begin{ibox}{Examples of point processes}
  Every point process is generated by a temporally continuously
  developing process. An event is generated whenever this process
  reaches a certain threshold. For example:
  \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 Evoked communication calls in crickets/frogs/birds: shaped by
    the dynamics of nervous system and the muscle appartus.
  \end{itemize}
\end{ibox}

\begin{figure}[t]
  \texpicture{pointprocessscetch}
  \titlecaption{\label{pointprocessscetchfig} Statistics of point
    processes.}{A point process is a sequence of instances in time
    $t_i$ that can be characterized through the inter-event-intervals
    $T_i=t_{i+1}-t_i$ and the number of events $n_i$. }
\end{figure}

In the neurosciences, the statistics of point processes is of
importance since the timing of the neuronal events (the action
potentials) is crucial for information transmission and can be treated
as such a process.

Point processes can be described using the intervals between
successive events $T_i=t_{i+1}-t_i$ and the number of observed events
within a certain time window $n_i$ (\figref{pointprocessscetchfig}).

The events originating from a point process can be illustrated in form
of a scatter- or 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}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\section{Intervalstatistics}

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 \enterm{interspike intervals}. The statistics of these
are described using the common measures.

\begin{figure}[t]
  \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
  \titlecaption{\label{isihexamplesfig}Interspike interval
    histogram}{of the spikes depicted in \figref{rasterexamplesfig}.}
\end{figure}

\begin{exercise}{isis.m}{}
  Implement a function \code{isis()} that calculates the interspike
  intervals from several spike trains. The function should return a
  single vector of intervals. The action potentials recorded in the
  individual trials are stored as vectors of spike times within a
  \codeterm{cell-array}. Spike times are given in seconds.
\end{exercise}

\subsection{First order interval statistics}
\begin{itemize}
\item Probability density $p(T)$ of the intervals $T$
  (\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \; dT
  = 1$.
\item Average interval: $\mu_{ISI} = \langle T \rangle =
  \frac{1}{n}\sum\limits_{i=1}^n T_i$.
\item Standard deviation of the interspike intervals: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
    \rangle)^2 \rangle}$\vspace{1ex}
\item \enterm{Coefficient of variation}: $CV_{ISI} =
  \frac{\sigma_{ISI}}{\mu_{ISI}}$.
\item \enterm{Diffusion coefficient}: $D_{ISI} =
  \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
\end{itemize}

\begin{exercise}{isihist.m}{}
  Implement a function \code{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 further returns the
  probability density as well as the centers of the bins.
\end{exercise}

\begin{exercise}{plotisihist.m}{}
  Implement a function that takes the return values of
  \code{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.
\end{exercise}

\subsection{Interval correlations}
So called \enterm{return maps} are used to illustrate
interdependencies between successive interspike intervals. The return
map plots the delayed interval $T_{i+k}$ against the interval
$T_i$. The parameter $k$ is called the \enterm{lag} $k$. Stationary
and non-stationary return maps are distinctly different
\figref{returnmapfig}.

\begin{figure}[t]
  \includegraphics[width=1\textwidth]{returnmapexamples}
  \includegraphics[width=1\textwidth]{serialcorrexamples}
  \titlecaption{\label{returnmapfig}Interspike interval analyses of a
    stationary and a non-stationary pointprocess.}{Upper plots show the
    return maps and the lower panels depict the serial correlation of
    successive intervals separated by the lag $k$.}
\end{figure}

Such dependencies can be further quantified using the \enterm{serial
  correlations} \figref{returnmapfig}. The serial correlation is the
correlation coefficient of the intervals $T_i$ and the intervals
delayed by the lag $T_{i+k}$:
\[ \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) \] The resulting correlation coefficient
$\rho_k$ is usually plotted against the lag $k$
\figref{returnmapfig}. $\rho_0=1$ is the correlation of each interval
with itself and is always 1.

\begin{exercise}{isiserialcorr.m}{}
  Implement a function \code{isiserialcorr()} that takes a vector of
  interspike intervals as input argument and calculates the serial
  correlation. The function should further plot the serial
  correlation. \pagebreak[4]
\end{exercise}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Count statistics}

% \begin{figure}[t]
%   \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill
%   \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
%   \titlecaption{\label{countstatsfig}Count Statistik.}{}
% \end{figure}
The number of events $n_i$ (counts) in a time window $i$ of the duration $W$
yields positive integer random numbers that are commonly quantified
using the following measures:
\begin{itemize}
\item Histogram of the counts $n_i$.
\item Average number of events: $\mu_N = \langle n \rangle$.
\item Variance of the counts: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
\item \determ{Fano Faktor} (The variance divided by the average): $F = \frac{\sigma_n^2}{\mu_n}$.
\end{itemize}
And in particular the average firing rate $r$ (spike count per time interval
, \determ{Feuerrate}) that is given in Hertz \sindex[term]{Feuerrate!mittlere Rate}
\begin{equation}
  \label{firingrate}
  r = \frac{\langle n \rangle}{W} \; .
\end{equation}

% \begin{figure}[t]
%   \begin{minipage}[t]{0.49\textwidth}
%     Poisson process $\lambda=100$\,Hz:\\
%     \includegraphics[width=1\textwidth]{poissonfano100hz}
%   \end{minipage}
%   \hfill
%   \begin{minipage}[t]{0.49\textwidth}
%     LIF $I=10$, $\tau_{adapt}=100$\,ms:\\
%     \includegraphics[width=1\textwidth]{lifadaptfano10-100ms}
%   \end{minipage}
%   \titlecaption{\label{fanofig}Fano factor.}{}
% \end{figure}

\begin{exercise}{counthist.m}{}
  Implement a function \code{counthist()} that calculates and plots
  the distribution of spike counts observed in a certain time
  window. The function should take two input arguments: (i) a
  \codeterm{cell-array} of vectors containing the spike times in
  seconds observed in a number of trials and (ii) the duration of the
  time window that is used to evaluate the counts.
\end{exercise}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Homogeneous Poisson process}

The Gaussian distribution is, due to the central limit theorem, the
standard for continuous measures. The equivalent in the realm of point
processes is the \enterm{Poisson distribution}.

In a \enterm[Poisson process!homogeneous]{homogeneous Poisson process}
the events occur at a fixed rate $\lambda=\text{const.}$ and are
independent of both the time $t$ and occurrence of previous events
(\figref{hompoissonfig}). The probability of observing an even within a
small time window of width $\Delta t$ is given by
\begin{equation}
  \label{hompoissonprob}
  P = \lambda \cdot \Delta t \; .
\end{equation}

In an \enterm[Poisson process!inhomogeneous]{inhomogeneous Poisson
  process}, however, the rate $\lambda$ depends on the time: $\lambda =
\lambda(t)$.

\begin{exercise}{poissonspikes.m}{}
  Implement a function \code{poissonspikes()} that uses a homogeneous
  Poisson process to generate events at a given rate for a certain
  duration and a number of trials. The rate should be given in Hertz
  and the duration of the trials is given in seconds. The function
  should return the event times in a cell-array. Each entry in this
  array represents the events observed in one trial. Apply
  \eqnref{hompoissonprob} to generate the event times.
\end{exercise}

\begin{figure}[t]
  \includegraphics[width=1\textwidth]{poissonraster100hz}
  \titlecaption{\label{hompoissonfig}Rasterplot of spikes of a
    homogeneous Poisson process with a rate $\lambda=100$\,Hz.}{}
\end{figure}

\begin{figure}[t]
  \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
  \includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
  \titlecaption{\label{hompoissonisihfig}Distribution of interspike intervals of two Poisson processes.}{}
\end{figure}

The homogeneous Poisson process has the following properties:
\begin{itemize}
\item Intervals $T$ are exponentially distributed (\figref{hompoissonisihfig}):
  \begin{equation}
    \label{poissonintervals}
    p(T) = \lambda e^{-\lambda T} \; .
  \end{equation}
\item The average interval is $\mu_{ISI} = \frac{1}{\lambda}$ .
\item The variance of the intervals is $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ .
\item Thus, the coefficient of variation is always $CV_{ISI} = 1$ .
\item The serial correlation is $\rho_k =0$ for $k>0$, since the
  occurrence of an event is independent of all previous events. Such a
  process is also called a \enterm{renewal process}.
\item The number of events $k$ within a temporal window of duration
  $W$ is Poisson distributed:
\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \]
(\figref{hompoissoncountfig})
\item The Fano Faktor is always $F=1$ .
\end{itemize}

\begin{exercise}{hompoissonspikes.m}{}
  Implement a function \code{hompoissonspikes()} that uses a
  homogeneous Poisson process to generate spike events at a given rate
  for a certain duration and a number of trials. The rate should be
  given in Hertz and the duration of the trials is given in
  seconds. The function should return the event times in a
  cell-array. Each entry in this array represents the events observed
  in one trial. Apply \eqnref{poissonintervals} to generate the event
  times.
\end{exercise}

\begin{figure}[t]
  \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
  \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
  \titlecaption{\label{hompoissoncountfig}Count statistics of Poisson
    spiketrains.}{}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Time-dependent firing rate}

So far we have discussed stationary spiketrains. 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 methods are illustrated in \figref{psthfig}. All of
these have their own justifications and pros- and cons. In the
following we will discuss these methods more
closely.

\begin{figure}[tp]
  \includegraphics[width=\columnwidth]{firingrates}
  \titlecaption{Estimating the time-dependent firing
    rate.}{\textbf{A)} Rasterplot depicting the spiking activity of a
    neuron. \textbf{B)} Firing rate calculated from the
    \emph{instantaneous rate}. \textbf{C)} classical PSTH with the
    \emph{binning} method. \textbf{D)} Firing rate estimated by
    \emph{convolution} of the activity with a Gaussian
    kernel.}\label{psthfig}
\end{figure}


\subsection{Instantaneous firing rate}

\begin{figure}[tp]
  \includegraphics[width=\columnwidth]{isimethod}
  \titlecaption{Instantaneous firing rate.}{Sketch of the recorded
    spiketrain (top).  Arrows illustrate the interspike intervals and
    number give the intervals in milliseconds. The inverse of the
    interspike interval is the \emph{instantaneous firing
      rate}.}\label{instratefig}
\end{figure}

A very simple method for estimating the time-dependent firing rate is
the \enterm[firing rate!instantaneous]{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 \emph{instantaneous firing rate} uses the interspike
interval, the \enterm{peri stimulus time histogram} (PSTH) uses the
spike count within observation windows of the duration $W$. It
estimates the probability of observing a spike within that observation
time. It tries to estimat the average rate in the limit of small
obersvation times \eqnref{psthrate}:
\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 \emph{binning
  method} or by using \emph{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 spiketrain as shown in \figref{instratefig} is used
    here. 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.}\label{binpsthfig}
\end{figure}

The 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. It is
possible to estimate the PSTH using the \code{hist()} method
\sindex[term]{Feuerrate!Binningmethode}

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.

\pagebreak[4]
\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 spiketrain as in \figref{instratefig}. The
    convolution of the spiketrain with a convolution kernel basically
    replaces each spike event with the kernel (top). 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.}\label{convratefig}
\end{figure}

With the convolution method we avoid the sharp edges of the binning
method. The spiketrain is convolved with a \enterm{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 kerel is correctly normalized
to an integral of one, figure \ref{convreffig}
bottom). \sindex[term]{Feuerrate!Faltungsmethode}

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.

\pagebreak[4]
\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}

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 \enterm[STA|see
  spike-triggered average]{spike-triggered average}. 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 electroreceptors
    and the stimulus reconstruction.}{The neuron was driven by a
    ``white-noise'' stimulus (blue curve, right panel). 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 convolution 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 encoded in the neuronal response (right
    panel).}\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}

\selectlanguage{english}