[pointprocesses] improved chapter and added rasterplot()

This commit is contained in:
Jan Benda 2018-11-26 19:05:57 +01:00
parent c72ea6ad37
commit acf7fc8c81
10 changed files with 107 additions and 86 deletions

View File

@ -9,7 +9,8 @@ pythonplots : $(PYPDFFILES)
$(PYPDFFILES) : %.pdf: %.py $(PYPDFFILES) : %.pdf: %.py
echo $$(which python) echo $$(which python)
python3 $< #python3 $<
python $<
cleanpythonplots : cleanpythonplots :
rm -f $(PYPDFFILES) rm -f $(PYPDFFILES)

View File

@ -1,13 +1,13 @@
subplot(1, 3, 1); subplot(1, 3, 1);
spikeraster(poissonspikes, 1.0); rasterplot(poissonspikes, 1.0);
title('Poisson'); title('Poisson');
subplot(1, 3, 2); subplot(1, 3, 2);
spikeraster(pifouspikes, 1.0); rasterplot(pifouspikes, 1.0);
title('PIF OU'); title('PIF OU');
subplot(1, 3, 3); subplot(1, 3, 3);
spikeraster(lifadaptspikes, 1.0); rasterplot(lifadaptspikes, 1.0);
title('LIF adapt'); title('LIF adapt');
savefigpdf(gcf, 'spikeraster.pdf', 15, 5); savefigpdf(gcf, 'spikeraster.pdf', 15, 5);

View File

@ -1,7 +1,7 @@
function spikeraster(spikes, tmax) function rasterplot(spikes, tmax)
% Display a spike raster of the spike times given in spikes. % Display a spike raster of the spike times given in spikes.
% %
% spikeraster(spikes, tmax) % rasterplot(spikes, tmax)
% spikes: a cell array of vectors of spike times in seconds % spikes: a cell array of vectors of spike times in seconds
% tmax: plot spike raster upto tmax seconds % tmax: plot spike raster upto tmax seconds

View File

@ -18,7 +18,7 @@
\section{TODO} \section{TODO}
\begin{itemize} \begin{itemize}
\item Add spikeraster function \item Explain difference stationary versus non-stationary point process
\item Multitrial firing rates \item Multitrial firing rates
\item Choice of bin width for PSTH, kernel width, also in relation sto \item Choice of bin width for PSTH, kernel width, also in relation sto
stimulus time scale stimulus time scale

View File

@ -2,17 +2,15 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Spiketrain analysis} \chapter{Spiketrain analysis}
\selectlanguage{english} \enterm[Action potentials]{action potentials} (\enterm{spikes}) are
the carriers of information in the nervous system. Thereby it is the
\enterm[Actionspotentials]{Actionspotentials} (\enterm{spikes}) are time at which the spikes are generated that is of importance for
the carriers of information in the nervous system. Thereby it is information transmission. The waveform of the action potential is
mainly the time at which the spikes are generated that is of largely stereotyped and therefore does not carry information.
importance. The waveform of the action potential is largely
stereotyped and does not carry information. The result of the pre-processing of electrophysiological recordings are
series of spike times, which are termed \enterm{spiketrains}. If
The result of the processing of electrophysiological recordings are measurements are repeated we get several \enterm{trials} of
series of spike times, which are then termed \enterm{spiketrains}. If
measurements are repeated we yield several \enterm{trials} of
spiketrains (\figref{rasterexamplesfig}). spiketrains (\figref{rasterexamplesfig}).
Spiketrains are times of events, the action potentials. The analysis Spiketrains are times of events, the action potentials. The analysis
@ -21,34 +19,32 @@ of these leads into the realm of the so called \enterm[point
\begin{figure}[ht] \begin{figure}[ht]
\includegraphics[width=1\textwidth]{rasterexamples} \includegraphics[width=1\textwidth]{rasterexamples}
\titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plot of \titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plots of
ten realizations of a stationary point process (homogeneous point ten trials of data illustrating the times of the action
process with a rate $\lambda=20$\;Hz, left) and an inhomogeneous potentials. Each vertical dash illustrates the time at which an
point process (perfect integrate-and-fire neuron dirven by action potential was observed. Each line displays the events of
Ohrnstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms, one trial. Shown is a stationary point process (left, homogeneous
right). Each vertical dash illustrates the time at which the point process with a rate $\lambda=20$\;Hz, left) and an
action potential was observed. Each line represents the event of non-stationary point process (right, perfect integrate-and-fire
each trial.} neuron dirven by Ohrnstein-Uhlenbeck noise with a time-constant
$\tau=100$\,ms, right).}
\end{figure} \end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Point processes} \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} \begin{ibox}{Examples of point processes}
Every point process is generated by a temporally continuously Every point process is generated by a temporally continuously
developing process. An event is generated whenever this process developing process. An event is generated whenever this process
reaches a certain threshold. For example: crosses some threshold. For example:\vspace{-1ex}
\begin{itemize} \begin{itemize}
\item Action potentials/heart beat: created by the dynamics of the \item Action potentials/heart beat: created by the dynamics of the
neuron/sinoatrial node neuron/sinoatrial node
\item Earthquake: defined by the dynamics of the pressure between \item Earthquake: defined by the dynamics of the pressure between
tectonical plates. tectonical plates.
\item Evoked communication calls in crickets/frogs/birds: shaped by \item Communication calls in crickets/frogs/birds: shaped by
the dynamics of nervous system and the muscle appartus. the dynamics of the nervous system and the muscle appartus.
\end{itemize} \end{itemize}
\end{ibox} \end{ibox}
@ -60,43 +56,51 @@ generates a sequence of events at times $\{t_i\}$, $t_i \in
$T_i=t_{i+1}-t_i$ and the number of events $n_i$. } $T_i=t_{i+1}-t_i$ and the number of events $n_i$. }
\end{figure} \end{figure}
In the neurosciences, the statistics of point processes is of A temporal \enterm{point process} is a stochastic process that
importance since the timing of the neuronal events (the action generates a sequence of events at times $\{t_i\}$, $t_i \in \reZ$. In
potentials) is crucial for information transmission and can be treated the neurosciences, the statistics of point processes is of importance
as such a process. since the timing of neuronal events (action potentials, post-synaptic
potentials, events in EEG or local-field recordings, etc.) is crucial
Point processes can be described using the intervals between for information transmission and can be treated as such a process.
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 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
The events originating from a point process can be illustrated in form event from two different point processes are shown in
of a scatter- or raster plot in which each vertical line indicates the \figref{rasterexamplesfig}. Point processes can be described using
time of an event. The event from two different point processes are the intervals between successive events $T_i=t_{i+1}-t_i$ and the
shown in \figref{rasterexamplesfig}. 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} \section{Interval statistics}
The intervals $T_i=t_{i+1}-t_i$ between successive events are real The intervals $T_i=t_{i+1}-t_i$ between successive events are real
positive numbers. In the context of action potentials they are positive numbers. In the context of action potentials they are
referred to as \enterm{interspike intervals}. The statistics of these referred to as \enterm{interspike intervals}. The statistics of
are described using the common measures. interspike intervals are described using common measures for
describing the statistics of stochastic real-valued variables:
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex} \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
\titlecaption{\label{isihexamplesfig}Interspike interval \titlecaption{\label{isihexamplesfig}Interspike interval
histogram}{of the spikes depicted in \figref{rasterexamplesfig}.} histograms}{of the spike trains shown in \figref{rasterexamplesfig}.}
\end{figure} \end{figure}
\begin{exercise}{isis.m}{} \begin{exercise}{isis.m}{}
Implement a function \code{isis()} that calculates the interspike Implement a function \code{isis()} that calculates the interspike
intervals from several spike trains. The function should return a intervals from several spike trains. The function should return a
single vector of intervals. The action potentials recorded in the single vector of intervals. The spike times (in seconds) of each
individual trials are stored as vectors of spike times within a trial are stored as vectors within a \codeterm{cell-array}.
\codeterm{cell-array}. Spike times are given in seconds.
\end{exercise} \end{exercise}
\subsection{First order interval statistics} %\subsection{First order interval statistics}
\begin{itemize} \begin{itemize}
\item Probability density $p(T)$ of the intervals $T$ \item Probability density $p(T)$ of the intervals $T$
(\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \; dT (\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \; dT
@ -138,10 +142,17 @@ and non-stationary return maps are distinctly different
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{returnmapexamples}
\includegraphics[width=1\textwidth]{serialcorrexamples} \includegraphics[width=1\textwidth]{serialcorrexamples}
\titlecaption{\label{returnmapfig}Interspike interval analyses of a \titlecaption{\label{returnmapfig}Interspike interval
stationary and a non-stationary pointprocess.}{Upper plots show the correlations}{of the spike trains shown in
return maps and the lower panels depict the serial correlation of \figref{rasterexamplesfig}. Upper panels show the return maps and
successive intervals separated by the lag $k$.} 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} \end{figure}
Such dependencies can be further quantified using the \enterm{serial Such dependencies can be further quantified using the \enterm{serial
@ -170,17 +181,18 @@ with itself and is always 1.
% \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms}
% \titlecaption{\label{countstatsfig}Count Statistik.}{} % \titlecaption{\label{countstatsfig}Count Statistik.}{}
% \end{figure} % \end{figure}
The number of events $n_i$ (counts) in a time window $i$ of the duration $W$ Counting the number of events $n_i$ (counts) in time windows $i$ of duration $W$
yields positive integer random numbers that are commonly quantified yields positive integer random numbers that are commonly quantified
using the following measures: using the following measures:
\begin{itemize} \begin{itemize}
\item Histogram of the counts $n_i$. \item Histogram of the counts $n_i$.
\item Average number of events: $\mu_N = \langle n \rangle$. \item Average number of counts: $\mu_n = \langle n \rangle$.
\item Variance of the counts: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$. \item Variance of 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}$. \item \determ{Fano Factor} (variance of counts divided by average count): $F = \frac{\sigma_n^2}{\mu_n}$.
\end{itemize} \end{itemize}
And in particular the average firing rate $r$ (spike count per time interval Of particular interest is the average firing rate $r$ (spike count per
, \determ{Feuerrate}) that is given in Hertz \sindex[term]{Feuerrate!mittlere Rate} time interval , \determ{Feuerrate}) that is given in Hertz
\sindex[term]{firing rate!average rate}
\begin{equation} \begin{equation}
\label{firingrate} \label{firingrate}
r = \frac{\langle n \rangle}{W} \; . r = \frac{\langle n \rangle}{W} \; .
@ -204,30 +216,30 @@ And in particular the average firing rate $r$ (spike count per time interval
the distribution of spike counts observed in a certain time the distribution of spike counts observed in a certain time
window. The function should take two input arguments: (i) a window. The function should take two input arguments: (i) a
\codeterm{cell-array} of vectors containing the spike times in \codeterm{cell-array} of vectors containing the spike times in
seconds observed in a number of trials and (ii) the duration of the seconds observed in a number of trials, and (ii) the duration of the
time window that is used to evaluate the counts. time window that is used to evaluate the counts.\pagebreak[4]
\end{exercise} \end{exercise}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Homogeneous Poisson process} \section{Homogeneous Poisson process}
The Gaussian distribution is, due to the central limit theorem, the The Gaussian distribution is, because of the central limit theorem,
standard for continuous measures. The equivalent in the realm of point the standard distribution for continuous measures. The equivalent in
processes is the \enterm{Poisson distribution}. the realm of point processes is the \enterm{Poisson distribution}.
In a \enterm[Poisson process!homogeneous]{homogeneous Poisson process} In a \enterm[Poisson process!homogeneous]{homogeneous Poisson process}
the events occur at a fixed rate $\lambda=\text{const.}$ and are the events occur at a fixed rate $\lambda=\text{const}$ and are
independent of both the time $t$ and occurrence of previous events independent of both the time $t$ and occurrence of previous events
(\figref{hompoissonfig}). The probability of observing an even within a (\figref{hompoissonfig}). The probability of observing an event within
small time window of width $\Delta t$ is given by a small time window of width $\Delta t$ is given by
\begin{equation} \begin{equation}
\label{hompoissonprob} \label{hompoissonprob}
P = \lambda \cdot \Delta t \; . P = \lambda \cdot \Delta t \; .
\end{equation} \end{equation}
In an \enterm[Poisson process!inhomogeneous]{inhomogeneous Poisson In an \enterm[Poisson process!inhomogeneous]{inhomogeneous Poisson
process}, however, the rate $\lambda$ depends on the time: $\lambda = process}, however, the rate $\lambda$ depends on time: $\lambda =
\lambda(t)$. \lambda(t)$.
\begin{exercise}{poissonspikes.m}{} \begin{exercise}{poissonspikes.m}{}
@ -249,7 +261,11 @@ In an \enterm[Poisson process!inhomogeneous]{inhomogeneous Poisson
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill
\includegraphics[width=0.45\textwidth]{poissonisihexp100hz} \includegraphics[width=0.45\textwidth]{poissonisihexp100hz}
\titlecaption{\label{hompoissonisihfig}Distribution of interspike intervals of two Poisson processes.}{} \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} \end{figure}
The homogeneous Poisson process has the following properties: The homogeneous Poisson process has the following properties:
@ -267,7 +283,10 @@ The homogeneous Poisson process has the following properties:
process is also called a \enterm{renewal process}. process is also called a \enterm{renewal process}.
\item The number of events $k$ within a temporal window of duration \item The number of events $k$ within a temporal window of duration
$W$ is Poisson distributed: $W$ is Poisson distributed:
\[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \] \begin{equation}
\label{poissoncounts}
P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!}
\end{equation}
(\figref{hompoissoncountfig}) (\figref{hompoissoncountfig})
\item The Fano Faktor is always $F=1$ . \item The Fano Faktor is always $F=1$ .
\end{itemize} \end{itemize}
@ -286,8 +305,11 @@ The homogeneous Poisson process has the following properties:
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill
\includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms} \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms}
\titlecaption{\label{hompoissoncountfig}Count statistics of Poisson \titlecaption{\label{hompoissoncountfig}Distribution of counts of a
spiketrains.}{} 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} \end{figure}
@ -312,7 +334,7 @@ closely.
\begin{figure}[tp] \begin{figure}[tp]
\includegraphics[width=\columnwidth]{firingrates} \includegraphics[width=\columnwidth]{firingrates}
\titlecaption{Estimating the time-dependent firing \titlecaption{Estimating the time-dependent firing
rate.}{\textbf{A)} Rasterplot depicting the spiking activity of a rate.}{\textbf{A)} Rasterplot showing the spiking activity of a
neuron. \textbf{B)} Firing rate calculated from the neuron. \textbf{B)} Firing rate calculated from the
\emph{instantaneous rate}. \textbf{C)} classical PSTH with the \emph{instantaneous rate}. \textbf{C)} classical PSTH with the
\emph{binning} method. \textbf{D)} Firing rate estimated by \emph{binning} method. \textbf{D)} Firing rate estimated by
@ -529,5 +551,3 @@ kernel.
the same size as the original stimulus contained in file the same size as the original stimulus contained in file
\file{sta\_data.mat}. \file{sta\_data.mat}.
\end{exercise} \end{exercise}
\selectlanguage{english}

View File

@ -1,4 +1,4 @@
set term epslatex size 11.4cm, 6.5cm set term epslatex size 11.4cm, 6.4cm
set out 'pointprocessscetch.tex' set out 'pointprocessscetch.tex'
set border 0 set border 0

View File

@ -1,7 +1,7 @@
%!PS-Adobe-2.0 EPSF-2.0 %!PS-Adobe-2.0 EPSF-2.0
%%Title: pointprocessscetchA.tex %%Title: pointprocessscetchA.tex
%%Creator: gnuplot 4.6 patchlevel 4 %%Creator: gnuplot 4.6 patchlevel 4
%%CreationDate: Tue Nov 28 10:02:49 2017 %%CreationDate: Mon Nov 26 17:57:49 2018
%%DocumentFonts: %%DocumentFonts:
%%BoundingBox: 50 50 373 135 %%BoundingBox: 50 50 373 135
%%EndComments %%EndComments
@ -433,7 +433,7 @@ SDict begin [
/Author (benda) /Author (benda)
% /Producer (gnuplot) % /Producer (gnuplot)
% /Keywords () % /Keywords ()
/CreationDate (Tue Nov 28 10:02:49 2017) /CreationDate (Mon Nov 26 17:57:49 2018)
/DOCINFO pdfmark /DOCINFO pdfmark
end end
} ifelse } ifelse

View File

@ -1,7 +1,7 @@
%!PS-Adobe-2.0 EPSF-2.0 %!PS-Adobe-2.0 EPSF-2.0
%%Title: pointprocessscetchB.tex %%Title: pointprocessscetchB.tex
%%Creator: gnuplot 4.6 patchlevel 4 %%Creator: gnuplot 4.6 patchlevel 4
%%CreationDate: Tue Nov 28 10:02:49 2017 %%CreationDate: Mon Nov 26 17:57:49 2018
%%DocumentFonts: %%DocumentFonts:
%%BoundingBox: 50 50 373 237 %%BoundingBox: 50 50 373 237
%%EndComments %%EndComments
@ -433,7 +433,7 @@ SDict begin [
/Author (benda) /Author (benda)
% /Producer (gnuplot) % /Producer (gnuplot)
% /Keywords () % /Keywords ()
/CreationDate (Tue Nov 28 10:02:49 2017) /CreationDate (Mon Nov 26 17:57:49 2018)
/DOCINFO pdfmark /DOCINFO pdfmark
end end
} ifelse } ifelse