562 lines
25 KiB
TeX
562 lines
25 KiB
TeX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
\chapter{Spiketrain analysis}
|
|
\exercisechapter{Spiketrain analysis}
|
|
|
|
\entermde[action potential]{Aktionspotential}{Action potentials}
|
|
(\enterm[spike|seealso{action potential}]{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. Analyzing
|
|
spike trains leads into the realm of the so called \entermde[point
|
|
process]{Punktprozess}{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 (homogeneous point
|
|
process with a rate $\lambda=20$\;Hz, left) and an non-stationary
|
|
point process (perfect integrate-and-fire neuron driven 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 apparatus.
|
|
\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 \entermde{Punktprozess}{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 \varcode{rasterplot()} that displays the times of
|
|
action potentials within the first \varcode{tmax} seconds in a raster
|
|
plot. The spike times (in seconds) recorded in the individual trials
|
|
are stored as vectors of times within a 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 \entermde{Interspikeintervalle}{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 \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}
|
|
|
|
%\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 \entermde[coefficient of variation]{Variationskoeffizient}{Coefficient of variation}:
|
|
$CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}$.
|
|
\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 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
|
|
\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.
|
|
\end{exercise}
|
|
|
|
\subsection{Interval correlations}
|
|
So called \entermde[return map]{return map}{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}
|
|
(\determ{Verz\"ogerung}) $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}
|
|
(\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 using the
|
|
\entermde[correlation!serial]{Korrelation!serielle}{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 \varcode{isiserialcorr()} that takes a vector of
|
|
interspike intervals as input argument and calculates the serial
|
|
correlation. The function should further plot the serial
|
|
correlation.
|
|
\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 \entermde{Fano Faktor}{Fano factor} (variance of counts divided by average count): $F = \frac{\sigma_n^2}{\mu_n}$.
|
|
\end{itemize}
|
|
Of particular interest is the \enterm[firing rate!average]{average
|
|
firing rate} $r$ (spike count per time interval, \determ{Feuerrate})
|
|
that is given in Hertz
|
|
\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 \varcode{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
|
|
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, because of the central limit theorem,
|
|
the standard distribution for continuous measures. The equivalent in
|
|
the realm of point processes is the
|
|
\entermde[distribution!Poisson]{Verteilung!Poisson-}{Poisson distribution}.
|
|
|
|
In a \entermde[Poisson process!homogeneous]{Poissonprozess!homogener}{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 \entermde[Poisson process!inhomogeneous]{Poissonprozess!inhomogener}{inhomogeneous Poisson
|
|
process}, however, the rate $\lambda$ depends on time: $\lambda =
|
|
\lambda(t)$.
|
|
|
|
\begin{exercise}{poissonspikes.m}{}
|
|
Implement a function \varcode{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} (\determ{Erneuerungsprozess}).
|
|
\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 \varcode{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 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 spiketrain
|
|
(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 spiketrain 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 spiketrain as in \figref{instratefig} (top). The
|
|
convolution of the spiketrain 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}
|
|
|
|
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
|