diff --git a/header.tex b/header.tex index 2f2360d..73288b8 100644 --- a/header.tex +++ b/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]{} -% 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}}} diff --git a/pointprocesses/code/counthist.m b/pointprocesses/code/counthist.m index 67b4cae..2f5b746 100644 --- a/pointprocesses/code/counthist.m +++ b/pointprocesses/code/counthist.m @@ -1,16 +1,11 @@ -function [counts, bins] = counthist(spikes, w) -% Compute and plot histogram of spike counts. +function counthist(spikes, w) +% Plot histogram of spike counts. % -% [counts, bins] = counthist(spikes, w) +% counthist(spikes, w) % % Arguments: % spikes: a cell array of vectors of spike times in seconds -% w: observation window duration in seconds for computing the counts -% -% Returns: -% counts: the histogram of counts normalized to probabilities -% bins: the bin centers for the histogram - +% w: duration of window in seconds for computing the counts % collect spike counts: tmax = spikes{1}(end); n = []; @@ -21,12 +16,12 @@ function [counts, bins] = counthist(spikes, w) % for tk = 0:w:tmax-w % nn = sum((times >= tk) & (times < tk+w)); % %nn = length(find((times >= tk) & (times < tk+w))); -% n = [n nn]; +% n = [n, nn]; % end % alternative 2: use the hist() function to do that! tbins = 0.5*w:w:tmax-0.5*w; nn = hist(times, tbins); - n = [n nn]; + n = [n, nn]; end % histogram of spike counts: @@ -36,9 +31,7 @@ function [counts, bins] = counthist(spikes, w) counts = counts / sum(counts); % plot: - if nargout == 0 - bar(bins, counts); - xlabel('counts k'); - ylabel('P(k)'); - end + bar(bins, counts); + xlabel('counts k'); + ylabel('P(k)'); end diff --git a/pointprocesses/code/isihist.m b/pointprocesses/code/isihist.m index 00b02c3..a24989c 100644 --- a/pointprocesses/code/isihist.m +++ b/pointprocesses/code/isihist.m @@ -8,7 +8,7 @@ function [pdf, centers] = isihist(isis, binwidth) % binwidth: optional width in seconds to be used for the isi bins % % Returns: -% pdf: vector with probability density of interspike intervals in Hz +% pdf: vector with pdf of interspike intervals in Hertz % centers: vector with centers of interspikeintervalls in seconds if nargin < 2 diff --git a/pointprocesses/code/isis.m b/pointprocesses/code/isis.m index 65a346b..509dcc5 100644 --- a/pointprocesses/code/isis.m +++ b/pointprocesses/code/isis.m @@ -5,11 +5,9 @@ function isivec = isis(spikes) % % Arguments: % spikes: a cell array of vectors of spike times in seconds -% isivec: a column vector with all the interspike intervalls % % Returns: % isivec: a column vector with all the interspike intervalls - isivec = []; for k = 1:length(spikes) difftimes = diff(spikes{k}); diff --git a/pointprocesses/code/isiserialcorr.m b/pointprocesses/code/isiserialcorr.m index 21bb474..1b813fc 100644 --- a/pointprocesses/code/isiserialcorr.m +++ b/pointprocesses/code/isiserialcorr.m @@ -1,15 +1,15 @@ -function isicorr = isiserialcorr(isivec, maxlag) +function [isicorr, lags] = isiserialcorr(isivec, maxlag) % serial correlation of interspike intervals % % isicorr = isiserialcorr(isivec, maxlag) % % Arguments: % isivec: vector of interspike intervals in seconds -% maxlag: the maximum lag in seconds +% maxlag: the maximum lag % % Returns: % isicorr: vector with the serial correlations for lag 0 to maxlag - +% lags: vector with the lags corresponding to isicorr lags = 0:maxlag; isicorr = zeros(size(lags)); for k = 1:length(lags) @@ -21,14 +21,4 @@ function isicorr = isiserialcorr(isivec, maxlag) isicorr(k) = corr(isivec(1:end-lag), isivec(lag+1:end)); end end - - if nargout == 0 - % plot: - plot(lags, isicorr, '-b'); - hold on; - scatter(lags, isicorr, 50.0, 'b', 'filled'); - hold off; - xlabel('Lag k') - ylabel('\rho_k') - end end diff --git a/pointprocesses/code/plotisiserialcorr.m b/pointprocesses/code/plotisiserialcorr.m new file mode 100644 index 0000000..3f21418 --- /dev/null +++ b/pointprocesses/code/plotisiserialcorr.m @@ -0,0 +1,16 @@ +function isicorr = plotisiserialcorr(isivec, maxlag) +% plot serial correlation of interspike intervals +% +% plotisiserialcorr(isivec, maxlag) +% +% Arguments: +% isivec: vector of interspike intervals in seconds +% maxlag: the maximum lag + [isicorr, lags] = isiserialcorr(isivec, maxlag); + plot(lags, isicorr, '-b'); + hold on; + scatter(lags, isicorr, 20.0, 'b', 'filled'); + hold off; + xlabel('Lag k') + ylabel('\rho_k') +end diff --git a/pointprocesses/code/rasterplot.m b/pointprocesses/code/rasterplot.m index 7b28820..ccbfe60 100644 --- a/pointprocesses/code/rasterplot.m +++ b/pointprocesses/code/rasterplot.m @@ -2,21 +2,35 @@ function rasterplot(spikes, tmax) % Display a spike raster of the spike times given in spikes. % % rasterplot(spikes, tmax) +% +% Arguments: % spikes: a cell array of vectors of spike times in seconds -% tmax: plot spike raster upto tmax seconds - +% tmax: plot spike raster up to tmax seconds + in_msecs = tmax < 1.5 + spiketimes = []; + trials = []; ntrials = length(spikes); for k = 1:ntrials times = spikes{k}; times = times(times0$, 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 introduce +correlations between intervals. Such non-renewal dynamics is then +characterized by the non-zero serial correlations +(\figref{returnmapfig}). \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. + Implement a function \varcode{isiserialcorr()} that takes a vector + of interspike intervals as input argument and calculates the serial + correlations up to some maximum lag. +\end{exercise} + +\begin{exercise}{plotisiserialcorr.m}{} + Implement a function \varcode{plotisiserialcorr()} that takes a + vector of interspike intervals as input argument and generates a + plot of the serial correlations. \end{exercise} @@ -189,11 +348,11 @@ with itself and is always 1. \begin{figure}[t] \includegraphics{countexamples} - \titlecaption{\label{countstatsfig}Count statistics.}{The - distribution of the number of events $k$ (blue) counted within - windows of 20\,ms (left) or 200\,ms duration (right) of the - homogeneous Poisson spike train with a rate of 20\,Hz shown in - \figref{rasterexamplesfig}. For Poisson spike trains these + \titlecaption{\label{countstatsfig}Count statistics.}{Probability + distributions of counting $k$ events $k$ (blue) within windows of + 20\,ms (left) or 200\,ms duration (right) of a homogeneous Poisson + spike train with a rate of 20\,Hz + (\figref{rasterexamplesfig}). For Poisson spike trains these distributions are given by Poisson distributions (red).} \end{figure} @@ -204,7 +363,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,28 +375,44 @@ 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$ appropriately normalized to + probability distributions. 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$. \end{itemize} -Because spike counts are unitless and positive numbers, the +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 -$W$. The average spike count, for example, grows linearly with $W$ for -sufficiently large time windows: $\langle n \rangle = r W$, -\eqnref{firingrate}. Doubling the counting window doubles the spike -count. As does the spike-count variance (\figref{fanofig}). At smaller -time windows the statistics of the event counts might depend on the -particular duration of the counting window. There might be an optimal -time window for which the variance of the spike count is minimal. The -Fano factor plotted as a function of the time window illustrates such -properties of point processes (\figref{fanofig}). + +Note that all of these statistics depend in general on the chosen +window length $W$. The average spike count, for example, grows +linearly with $W$ for sufficiently large time windows: $\langle n +\rangle = r W$, \eqnref{firingrate}. Doubling the counting window +doubles the spike count. As does the spike-count variance +(\figref{fanofig}). At smaller time windows the statistics of the +event counts might depend on the particular duration of the counting +window. There might be an optimal time window for which the variance +of the spike count is minimal. The Fano factor plotted as a function +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 @@ -246,126 +421,33 @@ information encoded in the mean spike count is transmitted. \begin{figure}[t] \includegraphics{fanoexamples} - \titlecaption{\label{fanofig} - Count variance and Fano factor.}{Variance of event counts as a - function of mean counts obtained by varying the duration of the - count window (left). Dividing the count variance by the respective - mean results in the Fano factor that can be plotted as a function - of the count window (right). For Poisson spike trains the variance - 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 + \titlecaption{\label{fanofig} Fano factor.}{Counting events in time + windows of given duration and then dividing the variance of the + counts by their mean results in the Fano factor. Here, the Fano + factor is plotted as a function of the duration of the window used + to count events. For Poisson spike trains the variance always + equals the mean counts and consequently the Fano factor equals one + irrespective of the count window (left). A spike train with + positive correlations between interspike intervals (caused by an + Ornstein-Uhlenbeck process) has a minimum in the Fano factor, that is an analysis window for which the relative count variance is minimal somewhere close to the correlation time scale of the - interspike intervals (bottom).} + interspike intervals (right).} \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. + window. The function should take two input arguments: a cell-array + of vectors containing the spike times in seconds observed in a + number of trials, and the duration of the time window that is used + to evaluate the counts. \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} +\label{nonstationarysec} So far we have discussed stationary spike trains. The statistical properties of these did not change within the observation time (stationary point @@ -533,6 +615,7 @@ relevate time-scale. \end{exercise} \section{Spike-triggered Average} +\label{stasec} The graphical representation of the neuronal activity alone is not sufficient tot investigate the relation between the neuronal response diff --git a/pointprocesses/lecture/rasterexamples.py b/pointprocesses/lecture/rasterexamples.py index 07b8628..7a9f4f0 100644 --- a/pointprocesses/lecture/rasterexamples.py +++ b/pointprocesses/lecture/rasterexamples.py @@ -85,7 +85,7 @@ def plot_inhomogeneous_spikes(ax): if __name__ == "__main__": - fig, (ax1, ax2) = plt.subplots(1, 2, figsize=cm_size(figure_width, 0.5*figure_width)) + fig, (ax1, ax2) = plt.subplots(1, 2) fig.subplots_adjust(**adjust_fs(fig, left=4.0, right=1.0, top=1.2)) plot_homogeneous_spikes(ax1) plot_inhomogeneous_spikes(ax2) diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex index 13ca240..016a6a4 100644 --- a/scientificcomputing-script.tex +++ b/scientificcomputing-script.tex @@ -116,9 +116,10 @@ \includechapter{pointprocesses} % add a chapter on simulating point-processes/spike trains -% (Poisson spike trains, integrate-and-fire models) +% (Poisson spike trains, LMP models, integrate-and-fire models, full HH like models) % add chapter on information theory, mutual information, stimulus reconstruction, coherence +% move STA here! % add a chapter on Bayesian inference (the Neuroscience of it and a % bit of application for statistical problems).