559 lines
25 KiB
TeX
559 lines
25 KiB
TeX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\chapter{Spiketrain analysis}
|
|
\exercisechapter{Spiketrain analysis}
|
|
|
|
\enterm[action potential]{Action potentials} (\enterm{spikes}) are
|
|
the carriers of information in the nervous system. Thereby it is the
|
|
time at which the spikes are generated that is of importance for
|
|
information transmission. The waveform of the action potential is
|
|
largely stereotyped and therefore does not carry information.
|
|
|
|
The result of the pre-processing of electrophysiological recordings are
|
|
series of spike times, which are termed \enterm{spiketrains}. If
|
|
measurements are repeated we get 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-plots of
|
|
ten trials of data illustrating the times of the action
|
|
potentials. Each vertical dash illustrates the time at which an
|
|
action potential was observed. Each line displays the events of
|
|
one trial. Shown is a stationary point process (left, homogeneous
|
|
point process with a rate $\lambda=20$\;Hz, left) and an
|
|
non-stationary point process (right, perfect integrate-and-fire
|
|
neuron dirven by Ohrnstein-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 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}
|
|
|
|
A temporal \enterm{point process} is a stochastic process that
|
|
generates a sequence of events at times $\{t_i\}$, $t_i \in \reZ$. 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}. 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}).
|
|
|
|
\begin{exercise}{rasterplot.m}{}
|
|
Implement a function \code{rasterplot()} that displays the times of
|
|
action potentials within the first \code{tmax} seconds in a raster
|
|
plot. The spike times (in seconds) recorded in the individual trials
|
|
are stored as vectors of times within a \codeterm{cell-array}.
|
|
\end{exercise}
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Interval statistics}
|
|
|
|
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
|
|
interspike intervals are described using common measures for
|
|
describing the statistics of stochastic real-valued variables:
|
|
|
|
\begin{figure}[t]
|
|
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
|
|
\titlecaption{\label{isihexamplesfig}Interspike interval
|
|
histograms}{of the spike trains shown 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 spike times (in seconds) of each
|
|
trial are stored as vectors within a \codeterm{cell-array}.
|
|
\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]{Coefficient of variation}: $CV_{ISI} =
|
|
\frac{\sigma_{ISI}}{\mu_{ISI}}$.
|
|
\item \enterm[diffusion coefficient]{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
|
|
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}. 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 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}
|
|
Counting the number of events $n_i$ (counts) in time windows $i$ of 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 counts: $\mu_n = \langle n \rangle$.
|
|
\item Variance of counts: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
|
|
\item \determ{Fano Factor} (variance of counts divided by average count): $F = \frac{\sigma_n^2}{\mu_n}$.
|
|
\end{itemize}
|
|
Of particular interest is the average firing rate $r$ (spike count per
|
|
time interval , \determ{Feuerrate}) that is given in Hertz
|
|
\sindex[term]{firing rate!average 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.\pagebreak[4]
|
|
\end{exercise}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\section{Homogeneous Poisson process}
|
|
|
|
The Gaussian distribution is, because of the central limit theorem,
|
|
the standard distribution 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 event 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 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.}{The processes differ in their
|
|
rate (left: $\lambda=20$\,Hz, right: $\lambda=100$\,Hz). The red
|
|
lines indicate the corresponding exponential interval distribution
|
|
\eqnref{poissonintervals}.}
|
|
\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:
|
|
\begin{equation}
|
|
\label{poissoncounts}
|
|
P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!}
|
|
\end{equation}
|
|
(\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}Distribution of counts of a
|
|
Poisson spiketrain.}{The count statistics was generated for two
|
|
different windows of width $W=10$\,ms (left) and width $W=100$\,ms
|
|
(right). The red line illustrates the corresponding Poisson
|
|
distribution \eqnref{poissoncounts}.}
|
|
\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 showing 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{convratefig}
|
|
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}
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\printsolutions
|