[pointprocesses] better incorporated Poisson spike train

header.tex

@@ -222,11 +222,26 @@
   {\vspace{-2ex}\lstset{#1}\noindent\minipage[t]{1\linewidth}}%
   {\endminipage}
 
-%%%%% english, german, code and file terms: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%% english and german terms: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \usepackage{ifthen}
 
+% \enterm[en-index]{term}: Typeset term and add it to the index of english terms
+%
+% \determ[de-index]{term}: Typeset term and add it to the index of german terms
+%
+% \entermde[en-index]{de-index}{term}: Typeset term and add it to the index of english terms
+%                                      and de-index to the index of german terms
+%
+% how to specificy an index:
+% \enterm{term}                                 - just put term into the index
+% \enterm[en-index]{term}                       - typeset term and put en-index into the index
+% \enterm[statistics!mean]{term}                - mean is a subentry of statistics
+% \enterm[statistics!average|see{statistics!mean}]{term} - cross reference to statistics mean
+% \enterm[statistics@\textbf{statistics}]{term} - put index at statistics but use 
+%                                               \textbf{statistics} for typesetting in the index
+
 % \enterm[english index entry]{<english term>}
-% typeset the term in italics and add it (or the optional argument) to
+% typeset the term in italics and add it (or rather the optional argument) to
 % the english index.
 \newcommand{\enterm}[2][]{\textit{#2}\ifthenelse{\equal{#1}{}}{\protect\sindex[enterm]{#2}}{\protect\sindex[enterm]{#1}}}
 

pointprocesses/lecture/isihexamples.py

@@ -93,6 +93,8 @@ def plot_hom_isih(ax):
     ax.set_ylim(0.0, 31.0)
     ax.set_xticks(np.arange(0.0, 151.0, 50.0))
     ax.set_yticks(np.arange(0.0, 31.0, 10.0))
+    tt = np.linspace(0.0, 0.15, 100)
+    ax.plot(1000.0*tt, rate*np.exp(-rate*tt), **lsB)
     plotisih(ax, isis(homspikes), 0.005)
 
     

pointprocesses/lecture/pointprocesses.tex

@@ -34,7 +34,7 @@ process]{Punktprozess}{point processes}.
     trial. Shown is a stationary point process (homogeneous point
     process with a rate $\lambda=20$\;Hz, left) and an non-stationary
     point process with a rate that varies in time (noisy perfect
-    integrate-and-fire neuron driven by Ohrnstein-Uhlenbeck noise with
+    integrate-and-fire neuron driven by Ornstein-Uhlenbeck noise with
     a time-constant $\tau=100$\,ms, right).}
 \end{figure}
 
@@ -87,15 +87,60 @@ certain time window $n_i$ (\figref{pointprocesssketchfig}).
   are stored as vectors of times within a cell array.
 \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{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}
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 \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
+referred to as \entermde[interspike
-  intervals}. The statistics of interspike intervals are described
+interval]{Interspikeintervall}{interspike intervals}, in short
-using common measures for describing the statistics of stochastic
+\entermde[ISI|see{interspike
-real-valued variables:
+  interval}]{ISI|see{Interspikeintervall}}{ISI}s. The statistics of
+interspike intervals are described using common measures for
+describing the statistics of real-valued variables:
 
 \begin{figure}[t]
   \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
@@ -110,19 +155,45 @@ real-valued variables:
   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
+  (\figref{isihexamplesfig}). Normalized to $\int_0^{\infty} p(T) \;
-  = 1$.
+  dT = 1$. Commonly referred to as the \enterm[interspike
-\item Average interval: $\mu_{ISI} = \langle T \rangle =
+  interval!histogram]{interspike interval histogram}. Its shape
-  \frac{1}{n}\sum\limits_{i=1}^n T_i$.
+  reveals many interesting aspects like locking or bursting that
-\item Standard deviation of the interspike intervals: $\sigma_{ISI} = \sqrt{\langle (T - \langle T
+  cannot be inferred from the mean or standard deviation. A particular
-    \rangle)^2 \rangle}$\vspace{1ex}
+  reference is the exponential distribution of intervals
-\item \entermde[coefficient of variation]{Variationskoeffizient}{Coefficient of variation}:
+  \begin{equation}
-  $CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}$.
+    \label{hompoissonexponential}
-\item \entermde[diffusion coefficient]{Diffusionskoeffizient}{Diffusion coefficient}: $D_{ISI} =
+    p_{exp}(T) = \lambda e^{-\lambda T}
-  \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
+  \end{equation}
+  of a homogeneous Poisson spike train with rate $\lambda$.
+\item Mean interval: $\mu_{ISI} = \langle T \rangle =
+  \frac{1}{n}\sum\limits_{i=1}^n T_i$. The average time it takes from
+  one event to the next. The inverse of the mean interval is identical
+  with the mean rate $\lambda$ (number of events per time, see below)
+  of the process.
+\item Standard deviation of intervals: $\sigma_{ISI} = \sqrt{\langle
+    (T - \langle T \rangle)^2 \rangle}$. Periodically spiking neurons
+  have little variability in their intervals, whereas many cortical
+  neurons cover a wide range with their intervals. The standard
+  deviation of homogeneous Poisson spike trains, $\sigma_{ISI} =
+  \frac{1}{\lambda}$, also equals the inverse rate.  Whether the
+  standard deviation of intervals is low or high, however, is better
+  quantified by the
+\item \entermde[coefficient of
+  variation]{Variationskoeffizient}{Coefficient of variation}, the
+  standard deviation of the ISIs relative to their mean:
+  \begin{equation}
+    \label{cvisi}
+    CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}
+  \end{equation}
+  Homogeneous Poisson spike trains have an CV of exactly one. The
+  lower the CV the more regularly firing a neuron is firing. CVs
+  larger than one are also possible in spike trains with small
+  intervals separated by really long ones.
+%\item \entermde[diffusion coefficient]{Diffusionskoeffizient}{Diffusion coefficient}: $D_{ISI} =
+%  \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$.
 \end{itemize}
 
 \begin{exercise}{isihist.m}{}
@@ -142,12 +213,33 @@ real-valued variables:
 \end{exercise}
 
 \subsection{Interval correlations}
-So called \entermde[return map]{return map}{return maps} are used to
+Intervals are not just numbers without an order, like weights of
-illustrate interdependencies between successive interspike
+tigers.  Intervals are temporally ordered and there could be temporal
-intervals. The return map plots the delayed interval $T_{i+k}$ against
+structure in the sequence of intervals. For example, short intervals
-the interval $T_i$. The parameter $k$ is called the \enterm{lag}
+could be followed by more longer ones, and vice versa. Such
-(\determ{Verz\"ogerung}) $k$. Stationary and non-stationary return
+dependencies in the sequence of intervals do not show up in the
-maps are distinctly different \figref{returnmapfig}.
+interval histogram.  We need additional measures to also quantify the
+temporal structure of the sequence of intervals.
+
+We can use the same techniques we know for visualizing and quantifying
+correlations in bivariate data sets, i.e. scatter plots and
+correlation coefficients. We form $(x,y)$ data pairs by taking the
+series of intervals $T_i$ as $x$-data values and pairing them with the
+$k$-th next intervals $T_{i+k}$ as $y$-data values. The parameter $k$
+is called \enterm{lag} (\determ{Verz\"ogerung}). For lag one we pair
+each interval with the next one. A \entermde[return map]{return
+  map}{Return map} illustrates dependencies between successive
+intervals by simply plotting the intervals $T_{i+k}$ against the
+intervals $T_i$ in a scatter plot (\figref{returnmapfig}). For Poisson
+spike trains there is no structure beyond the one expected from the
+exponential interspike interval distribution, hinting at neighboring
+interspike intervals being independent of each other. For the spike
+train based on an Ornstein-Uhlenbeck process the return map is more
+clustered along the diagonal, hinting at a positive correlation
+between succeeding intervals. That is, short intervals are more likely
+to be followed by short ones and long intervals more likely by long
+ones. This temporal structure was already clearly visible in the spike
+raster shown in \figref{rasterexamplesfig}.
 
 \begin{figure}[tp]
   \includegraphics[width=1\textwidth]{serialcorrexamples}
@@ -165,16 +257,34 @@ maps are distinctly different \figref{returnmapfig}.
     the mean firing rate of the spike trains.}
 \end{figure}
 
-Such dependencies can be further quantified using the
+Such dependencies can be further quantified by
 \entermde[correlation!serial]{Korrelation!serielle}{serial
-  correlations} \figref{returnmapfig}. The serial correlation is the
+  correlations}. These are the correlation coefficients between the
-correlation coefficient of the intervals $T_i$ and the intervals
+intervals $T_{i+k}$ and $T_i$ in dependence on lag $k$:
-delayed by the lag $T_{i+k}$:
+\begin{equation}
-\[ \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)} 
+  \label{serialcorrelation}
-= {\rm corr}(T_{i+k}, T_i) \] The resulting correlation coefficient
+  \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)} 
-$\rho_k$ is usually plotted against the lag $k$
+= {\rm corr}(T_{i+k}, T_i)
-\figref{returnmapfig}. $\rho_0=1$ is the correlation of each interval
+\end{equation}
-with itself and is always 1.
+The serial correlations $\rho_k$ are usually plotted against the lag
+$k$ for a range small range of lags
+(\figref{returnmapfig}). $\rho_0=1$ is the correlation of each
+interval with itself and always equals one.
+
+If the serial correlations all equal zero, $\rho_k =0$ for $k>0$, then
+the length of an interval is independent of all the previous
+ones. Such a process is a \enterm{renewal process}
+(\determ{Erneuerungsprozess}). Each event, each action potential,
+erases the history. The occurrence of the next event is independent of
+what happened before. To a first approximation an action potential
+erases all information about the past from the membrane voltage and
+thus spike trains may approximate renewal processes.
+
+However, other variables like the intracellular calcium concentration
+or the states of slowly switching ion channels may carry information
+from one interspike interval to the next and thus introducing
+correlations. Such non-renewal dynamics can then be described by the
+non-zero serial correlations (\figref{returnmapfig}).
 
 \begin{exercise}{isiserialcorr.m}{}
   Implement a function \varcode{isiserialcorr()} that takes a vector of
@@ -204,7 +314,7 @@ is the average number of spikes counted within some time interval $W$
   \label{firingrate}
   r = \frac{\langle n \rangle}{W}
 \end{equation}
-and is neasured in Hertz. The average of the spike counts is taken
+and is measured in Hertz. The average of the spike counts is taken
 over trials. For stationary spike trains (no change in statistics, in
 particular the firing rate, over time), the firing rate based on the
 spike count equals the inverse average interspike interval
@@ -216,7 +326,14 @@ split into many segments $i$, each of duration $W$, and the number of
 events $n_i$ in each of the segments can be counted. The integer event
 counts can be quantified in the usual ways:
 \begin{itemize}
-\item Histogram of the counts $n_i$ (\figref{countstatsfig}).
+\item Histogram of the counts $n_i$. For homogeneous Poisson spike
+  trains with rate $\lambda$ the resulting probability distributions
+  follow a Poisson distribution (\figref{countstatsfig}), where the
+  probability of counting $k$ events within a time window $W$ is given by
+  \begin{equation}
+    \label{poissondist}
+    P(k) = \frac{(\lambda W)^k e^{\lambda W}}{k!}
+  \end{equation}
 \item Average number of counts: $\mu_n = \langle n \rangle$.
 \item Variance of counts:
   $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$.
@@ -224,20 +341,29 @@ counts can be quantified in the usual ways:
 Because spike counts are unitless and positive numbers, the
 \begin{itemize}
 \item \entermde{Fano Faktor}{Fano factor} (variance of counts divided
-  by average count): $F = \frac{\sigma_n^2}{\mu_n}$.
+  by average count)
+  \begin{equation}
+    \label{fano}
+    F = \frac{\sigma_n^2}{\mu_n}
+  \end{equation}
+  is a commonly used measure for quantifying the variability of event
+  counts relative to the mean number of events. In particular for
+  homogeneous Poisson processes the Fano factor equals one,
+  independently of the time window $W$.
 \end{itemize}
 is an additional measure quantifying event counts.
 
-Note that all of these statistics depend on the chosen window length
+Note that all of these statistics depend in general on the chosen
-$W$. The average spike count, for example, grows linearly with $W$ for
+window length $W$. The average spike count, for example, grows
-sufficiently large time windows: $\langle n \rangle = r W$,
+linearly with $W$ for sufficiently large time windows: $\langle n
-\eqnref{firingrate}. Doubling the counting window doubles the spike
+\rangle = r W$, \eqnref{firingrate}. Doubling the counting window
-count. As does the spike-count variance (\figref{fanofig}). At smaller
+doubles the spike count. As does the spike-count variance
-time windows the statistics of the event counts might depend on the
+(\figref{fanofig}). At smaller time windows the statistics of the
-particular duration of the counting window. There might be an optimal
+event counts might depend on the particular duration of the counting
-time window for which the variance of the spike count is minimal. The
+window. There might be an optimal time window for which the variance
-Fano factor plotted as a function of the time window illustrates such
+of the spike count is minimal. The Fano factor plotted as a function
-properties of point processes (\figref{fanofig}).
+of the time window illustrates such properties of point processes
+(\figref{fanofig}).
 
 This also has consequences for information transmission in neural
 systems. The lower the variance in spike count relative to the
@@ -255,9 +381,9 @@ information encoded in the mean spike count is transmitted.
     always equals the mean counts and consequently the Fano factor
     equals one irrespective of the count window (top). A spike train
     with positive correlations between interspike intervals (caused by
-    Ohrnstein-Uhlenbeck noise) has a minimum in the Fano factor, that
+    an Ornstein-Uhlenbeck process) has a minimum in the Fano factor,
-    is an analysis window for which the relative count variance is
+    that is an analysis window for which the relative count variance
-    minimal somewhere close to the correlation time scale of the
+    is minimal somewhere close to the correlation time scale of the
     interspike intervals (bottom).}
 \end{figure}
 
@@ -271,99 +397,6 @@ information encoded in the mean spike count is transmitted.
 \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 spike train.}{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}