%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Analyse von Spiketrains} \determ{Aktionspotentiale} (\enterm{Spikes}) sind die Tr\"ager der Information in Nervensystemen. Dabei ist in erster Linie nur der Zeitpunkt des Auftretens eines Aktionspotentials von Bedeutung. Die genaue Form des Aktionspotentials spielt keine oder nur eine untergeordnete Rolle. Nach etwas Vorverarbeitung haben elektrophysiologische Messungen deshalb Listen von Spikezeitpunkten als Ergebniss --- sogenannte \enterm{Spiketrains}. Diese Messungen k\"onnen wiederholt werden und es ergeben sich mehrere \enterm{trials} von Spiketrains (\figref{rasterexamplesfig}). Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen --- und deren Analyse f\"allt daher in das Gebiet der Statistik von sogenannten \determ{Punktprozessen}. \begin{figure}[ht] \includegraphics[width=1\textwidth]{rasterexamples} \titlecaption{\label{rasterexamplesfig}Raster-Plot.}{Raster-Plot von jeweils 10 Realisierungen eines station\"arenen Punktprozesses (homogener Poisson Prozess mit Rate $\lambda=20$\;Hz, links) und eines nicht-station\"aren Punktprozesses (perfect integrate-and-fire Neuron getrieben mit Ohrnstein-Uhlenbeck Rauschen mit Zeitkonstante $\tau=100$\,ms, rechts).} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Punktprozesse} Ein zeitlicher Punktprozess (\enterm{point process}) ist ein stochastischer Prozess, der eine Abfolge von Ereignissen zu den Zeiten $\{t_i\}$, $t_i \in \reZ$, generiert. \begin{ibox}{Beispiele von Punktprozessen} Jeder Punktprozess wird durch einen sich in der Zeit kontinuierlich entwickelnden Prozess generiert. Wann immer dieser Prozess eine Schwelle \"uberschreitet wird ein Ereigniss des Punktprozesses erzeugt. Zum Beispiel: \begin{itemize} \item Aktionspotentiale/Herzschlag: wird durch die Dynamik des Membranpotentials eines Neurons/Herzzelle erzeugt. \item Erdbeben: wird durch die Dynamik des Druckes zwischen tektonischen Platten auf beiden Seiten einer geologischen Verwerfung erzeugt. \item Zeitpunkt eines Grillen/Frosch/Vogelgesangs: wird durch die Dynamik des Nervensystems und des Muskelapparates erzeugt. \end{itemize} \end{ibox} \begin{figure}[t] \texpicture{pointprocessscetch} \titlecaption{\label{pointprocessscetchfig} Statistik von Punktprozessesen.}{Ein Punktprozess ist eine Abfolge von Zeitpunkten $t_i$ die auch durch die Intervalle $T_i=t_{i+1}-t_i$ oder die Anzahl der Ereignisse $n_i$ beschrieben werden kann. } \end{figure} F\"ur die Neurowissenschaften ist die Statistik der Punktprozesse besonders wichtig, da die Zeitpunkte der Aktionspotentiale als zeitlicher Punktprozess betrachtet werden k\"onnen und entscheidend f\"ur die Informations\"ubertragung sind. Bei Punktprozessen k\"onnen wir die Zeitpunkte $t_i$ ihres Auftretens, die Intervalle zwischen diesen Zeitpunkten $T_i=t_{i+1}-t_i$, sowie die Anzahl der Ereignisse $n_i$ bis zu einer bestimmten Zeit betrachten (\figref{pointprocessscetchfig}). Zwei Punktprozesse mit verschiedenen Eigenschaften sind in \figref{rasterexamplesfig} als Rasterplot dargestellt, bei dem die Zeitpunkte der Ereignisse durch senkrechte Striche markiert werden. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Intervallstatistik} Die Intervalle $T_i=t_{i+1}-t_i$ zwischen aufeinanderfolgenden Ereignissen sind reelle, positive Zahlen. Bei Aktionspotentialen heisen die Intervalle auch \enterm{Interspikeintervalle}. Deren Statistik kann mit den \"ublichen Gr\"o{\ss}en beschrieben werden. \begin{figure}[t] \includegraphics[width=1\textwidth]{isihexamples}\hfill \titlecaption{\label{isihexamplesfig}Interspikeintervall Histogramme}{der in \figref{rasterexamplesfig} gezeigten Spikes.} \end{figure} \subsection{Intervallstatistik erster Ordnung} \begin{itemize} \item Wahrscheinlichkeitsdichte $p(T)$ der Intervalle $T$ (\figref{isihexamplesfig}). Normiert auf $\int_0^{\infty} p(T) \; dT = 1$. \item Mittleres Intervall: $\mu_{ISI} = \langle T \rangle = \frac{1}{n}\sum\limits_{i=1}^n T_i$. \item Varianz der Intervalle: $\sigma_{ISI}^2 = \langle (T - \langle T \rangle)^2 \rangle$\vspace{1ex} \item Variationskoeffizient (\enterm{coefficient of variation}): $CV_{ISI} = \frac{\sigma_{ISI}}{\mu_{ISI}}$. \item Diffusions Koeffizient: $D_{ISI} = \frac{\sigma_{ISI}^2}{2\mu_{ISI}^3}$. \end{itemize} \subsection{Korrelationen der Intervalle} In \enterm{return maps} werden die um das \enterm{lag} $k$ verz\"ogerten Intervalle $T_{i+k}$ gegen die Intervalle $T_i$ geplottet. Dies macht m\"ogliche Abh\"angigkeiten von aufeinanderfolgenden Intervallen sichtbar. \begin{figure}[t] \includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{serialcorrexamples} \titlecaption{\label{returnmapfig}Interspikeintervall return maps und serielle Korrelationen}{zwischen aufeinander folgenden Intervallen im Abstand des Lags $k$.} \end{figure} Solche Ab\"angigkeiten werden durch die serielle Korrelation der Intervalle quantifiziert. Das ist der Korrelationskoeffizient zwischen aufeinander folgenden Intervallen getrennt durch \enterm{lag} $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) \] \"Ublicherweise wird die Korrelation $\rho_k$ gegen den Lag $k$ aufgetragen (\figref{returnmapfig}). $\rho_0=1$ (Korrelation jedes Intervalls mit sich selber). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Z\"ahlstatistik} % \begin{figure}[t] % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} % \titlecaption{\label{countstatsfig}Count Statistik.}{} % \end{figure} Die Anzahl der Ereignisse $n_i$ in Zeifenstern $i$ der L\"ange $W$ ergeben ganzzahlige, positive Zufallsvariablen, die meist durch folgende Sch\"atzer charakterisiert werden: \begin{itemize} \item Histogramm der counts $n_i$. \item Mittlere Anzahl von Ereignissen: $\mu_N = \langle n \rangle$. \item Varianz der Anzahl: $\sigma_n^2 = \langle (n - \langle n \rangle)^2 \rangle$. \item Fano Faktor (Varianz geteilt durch Mittelwert): $F = \frac{\sigma_n^2}{\mu_n}$. \end{itemize} Insbesondere ist die mittlere Rate der Ereignisse $r$ (Spikes pro Zeit, \determ{Feuerrate}) gemessen 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} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Homogener Poisson Prozess} F\"ur kontinuierliche Me{\ss}gr\"o{\ss}en ist die Normalverteilung u.a. wegen dem Zentralen Grenzwertsatz die Standardverteilung. Eine \"ahnliche Rolle spielt bei Punktprozessen der \determ{Poisson Prozess}. Beim homogenen Poisson Prozess treten Ereignisse mit einer festen Rate $\lambda=\text{const.}$ auf und sind unabh\"angig von der Zeit $t$ und unabh\"angig von den Zeitpunkten fr\"uherer Ereignisse (\figref{hompoissonfig}). Die Wahrscheinlichkeit zu irgendeiner Zeit ein Ereigniss in einem kleinen Zeitfenster der Breite $\Delta t$ zu bekommen ist \[ P = \lambda \cdot \Delta t \; . \] Beim inhomogenen Poisson Prozess h\"angt die Rate $\lambda$ von der Zeit ab: $\lambda = \lambda(t)$. \begin{figure}[t] \includegraphics[width=1\textwidth]{poissonraster100hz} \titlecaption{\label{hompoissonfig}Rasterplot von Spikes eines homogenen Poisson Prozesses mit $\lambda=100$\,Hz.}{} \end{figure} \begin{figure}[t] \includegraphics[width=0.45\textwidth]{poissonisihexp20hz}\hfill \includegraphics[width=0.45\textwidth]{poissonisihexp100hz} \titlecaption{\label{hompoissonisihfig}Interspikeintervallverteilungen zweier Poissonprozesse.}{} \end{figure} Der homogne Poissonprozess hat folgende Eigenschaften: \begin{itemize} \item Die Intervalle $T$ sind exponentiell verteilt: $p(T) = \lambda e^{-\lambda T}$ (\figref{hompoissonisihfig}). \item Das mittlere Intervall ist $\mu_{ISI} = \frac{1}{\lambda}$ . \item Die Varianz der Intervalle ist $\sigma_{ISI}^2 = \frac{1}{\lambda^2}$ . \item Der Variationskoeffizient ist also immer $CV_{ISI} = 1$ . \item Die seriellen Korrelationen $\rho_k =0$ f\"ur $k>0$, da das Auftreten der Ereignisse unabh\"angig von der Vorgeschichte ist. Ein solcher Prozess wird auch \determ{Erneuerungsprozess} genannt (\enterm{renewal process}). \item Die Anzahl der Ereignisse $k$ innerhalb eines Fensters der L\"ange W ist Poissonverteilt: \[ P(k) = \frac{(\lambda W)^ke^{\lambda W}}{k!} \] (\figref{hompoissoncountfig}) \item Der Fano Faktor ist immer $F=1$ . \end{itemize} \begin{figure}[t] \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz10ms}\hfill \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms} \titlecaption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.}{} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Zeitabh\"angige Feuerraten} Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik sich innerhalb der Analysezeit nicht ver\"andert (station\"are Punktprozesse). Meistens jedoch \"andert sich die Statistik der Spiketrains eines Neurons mit der Zeit. Z.B. kann ein sensorisches Neuron auf einen Reiz hin mit einer erh\"ohten Feuerrate antworten (nichtstation\"arer Punktprozess). Wie die mittlere Anzahl der Spikes sich mit der Zeit ver\"andert, die \determ{Feuerrate} $r(t)$, ist die wichtigste Gr\"o{\ss}e bei nicht-station\"aren Spiketrains. Die Einheit der Feuerrate ist Hertz, also Anzahl Aktionspotentiale pro Sekunde. Es gibt verschiedene Methoden diese zu bestimmen. Drei solcher Methoden sind in Abbildung \ref{psthfig} dargestellt. Alle Methoden haben ihre Berechtigung und ihre Vor- und Nachteile. Im folgenden werden die drei Methoden aus Abbildung \ref{psthfig} n\"aher erl\"autert. \begin{figure}[tp] \includegraphics[width=\columnwidth]{firingrates} \titlecaption{Verschiedene Methoden die zeitabh\"angige Feuerrate zu bestimmen.}{\textbf{A)} Rasterplot einer einzelnen neuronalen Antwort. Jeder vertikale Strich notiert den Zeitpunkt eines Aktionspotentials. \textbf{B)} Feurerrate aus der instantanen Feuerrate bestimmt. \textbf{C)} klassisches PSTH mit der Binning Methode. \textbf{D)} Feuerrate durch Faltung mit einem Gauss Kern bestimmt.}\label{psthfig} \end{figure} \subsection{Instantane Feuerrate} \begin{figure}[tp] \includegraphics[width=\columnwidth]{isimethod} \titlecaption{Bestimmung des zeitabh\"angigen Feuerrate aus dem Interspike Intervall.}{\textbf{A)} Skizze eines Rasterplots einer einzelnen neuronalen Antwort. Jeder vertikale Strich notiert den Zeitpunkt eines Aktionspotentials. Die Pfeile zwischen aufeinanderfolgenden Aktionspotentialen illustrieren das Interspike Interval. \textbf{B)} Der Kehrwert des Interspike Intervalls ist die Feuerrate.}\label{instrate} \end{figure} Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist die sogenannte \determ{instantane Feuerrate}. Dabei wird die Feuerrate aus dem Kehrwert der Interspikeintervalle, der Zeit zwischen zwei aufeinander folgenden Aktionspotentialen (\figref{instrate} A), bestimmt. Die abgesch\"atzte Feuerrate (\figref{instrate} B) ist g\"ultig f\"ur das gesammte Interspikeintervall. Diese Methode hat den Vorteil, dass sie sehr einfach zu berechnen ist und keine Annahme \"uber eine relevante Zeitskala (der Kodierung oder des Auslesemechanismus der postsynaptischen Zelle) macht. $r(t)$ ist allerdings keine kontinuierliche Funktion, die Spr\"unge in der Feuerrate k\"onnen f\"ur manche Analysen nachteilig sein. Au{\ss}erdem wird die Feuerrate nie gleich Null, auch wenn lange keine Aktionspotentiale generiert wurden. \subsection{Peri-Stimulus-Zeit-Histogramm} W\"ahrend die Instantane Rate den Kehrwert der Zeit von einem bis zum n\"achsten Aktionspotential misst, sch\"atzt das sogenannte \determ{Peri-Stimulus-Zeit-Histogramm} (\enterm{peri stimulus time histogram}, PSTH) die Wahrscheinlichkeit ab, zu einem Zeitpunkt Aktionspotentiale anzutreffen. Es wird versucht die mittlere Rate \eqnref{firingrate} im Grenzwert kleiner Beobachtungszeiten abzusch\"atzen: \begin{equation} \label{psthrate} r(t) = \lim_{W \to 0} \frac{\langle n \rangle}{W} \; , \end{equation} wobei die Anzahl $n$ der Aktionspotentiale, die im Zeitintervall $(t,t+W)$ aufgetreten sind, \"uber trials gemittelt wird. Eine solche Rate enspricht der zeitabh\"angigen Rate $\lambda(t)$ des inhomogenen Poisson-Prozesses. Das PSTH \eqnref{psthrate} kann entweder \"uber die Binning-Methode oder durch Verfaltung mit einem Kern bestimmt werden. Beiden Methoden gemeinsam ist die Notwendigkeit der Wahl einer zus\"atzlichen Zeitskala, die der Beobachtungszeit $W$ in \eqnref{psthrate} entspricht. \subsubsection{Binning-Methode} Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in die jeweiligen Bins fallen, gez\"ahlt (\figref{binpsth} A). Um diese Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der Binweite normiert werden. Das ist fast so, wie beim Absch\"atzen einer Wahrscheinlichkeitsdichte. Es kann auch die \code{hist} Funktion zur Bestimmung des PSTHs verwendet werden. Die bestimmte Feuerrate gilt f\"ur das gesamte Bin (\figref{binpsth} B). Das so berechnete PSTH hat wiederum eine stufige Form, die von der Wahl der Binweite anh\"angt. $r(t)$ ist also keine stetige Funktion. Die Binweite bestimmt die zeitliche Aufl\"osung der Absch\"atzung. \"Anderungen in der Feuerrate, die innerhalb eines Bins vorkommen k\"onnen nicht aufgl\"ost werden. Mit der Wahl der Binweite wird somit eine Annahme \"uber die relevante Zeitskala des Spiketrains gemacht. \begin{figure}[tp] \includegraphics[width=\columnwidth]{binmethod} \titlecaption{Bestimmung des PSTH mit der Binning Methode.}{\textbf{A)} Skizze eines Rasterplots einer einzelnen neuronalen Antwort. Jeder vertikale Strich notiert den Zeitpunkt eines Aktionspotentials. Die roten gestrichelten Linien stellen die Grenzen der Bins dar und die Zahlen geben den Spike Count pro Bin an. \textbf{B)} Die Feuerrate erh\"alt man indem das Zeithistogramm mit der Binweite normiert.}\label{binpsth} \end{figure} \subsubsection{Faltungsmethode} Bei der Faltungsmethode werden die harten Kanten der Bins der Binning-Methode vermieden. Der Spiketrain wird mit einem Kern verfaltet, d.h. jedes Aktionspotential wird durch den Kern ersetzt. Zur Berechnung wird die Aktionspotentialfolge zun\"achst ``bin\"ar'' dargestellt. Dabei wird ein Spiketrain als (Zeit-)Vektor dargestellt, in welchem die Zeitpunkte der Aktionspotentiale als 1 notiert werden. Alle anderen Elemente des Vektors sind 0. Anschlie{\ss}end wir dieser bin\"are Spiketrain mit einem Gau{\ss}-Kern bestimmter Breite verfaltet: \[r(t) = \int_{-\infty}^{\infty} \omega(\tau) \, \rho(t-\tau) \, {\rm d}\tau \; , \] wobei $\omega(\tau)$ der Filterkern und $\rho(t)$ die bin\"are Antwort ist. Bildlich geprochen wird jede 1 in $\rho(t)$ durch den Filterkern ersetzt (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert wurde (Integral gleich Eins), ergibt sich die Feuerrate direkt aus der \"Uberlagerung der Kerne (Abb. \ref{convrate} B). \begin{figure}[tp] \includegraphics[width=\columnwidth]{convmethod} \titlecaption{Schematische Darstellung der Faltungsmethode.}{\textbf{A)} Rasterplot einer einzelnen neuronalen Antwort. Jeder vertikale Strich notiert den Zeitpunkt eines Aktionspotentials. In der Faltung werden die mit einer 1 notierten Aktionspotential durch den Faltungskern ersetzt. \textbf{B)} Bei korrekter Normierung des Kerns ergibt sich die Feuerrate direkt aus der \"Uberlagerung der Kerne.}\label{convrate} \end{figure} Die Faltungsmethode f\"uhrt, anders als die anderen Methoden, zu einer kontinuierlichen Funktion was insbesondere f\"ur spektrale Analysen von Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur Binweite, die zeitliche Aufl\"osung von $r(t)$. Die Breite des Kerns macht also auch wieder eine Annahme \"uber die relevante Zeitskala des Spiketrains. \section{Spike-triggered Average} Die graphischer Darstellung der Feuerrate allein reicht nicht aus um den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu analysieren. Eine Methode um mehr \"uber diesen Zusammenhang zu erfahren, ist der \enterm{Spike-triggered average} (STA). Der STA \begin{equation} STA(\tau) = \langle s(t - \tau) \rangle = \frac{1}{N} \sum_{i=1}^{N} s(t_i - \tau) \end{equation} der $N$ Aktionspotentiale zu den Zeiten $t_i$ in Anwort auf den Stimulus $s(t)$ ist der mittlere Stimulus, der zu einem Aktionspotential in der neuronalen Antwort f\"uhrt. Der STA l\"a{\ss}t sich relativ einfach berechnen, indem aus dem Stimulus f\"ur jeden beobachteten Spike ein entsprechender Abschnitt ausgeschnitten wird und diese dann gemittelt werde (\figref{stafig}). \begin{figure}[t] \includegraphics[width=0.5\columnwidth]{sta} \titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors.}{Der Rezeptor wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0 ist der Zeitpunkt des beobachteten Aktionspotentials.}\label{stafig} \end{figure} Aus dem STA k\"onnen verschiedene Informationen \"uber den Zusammenhang zwischen Stimulus und neuronaler Antwort gewonnen werden. Die Breite des STA repr\"asentiert die zeitliche Pr\"azision, mit der Stimulus und Antwort zusammenh\"angen und wie weit das Neuron zeitlich integriert. Die Amplitude des STA (gegeben in der gleichen Einheit wie der Stimulus) deutet auf die Empfindlichkeit des Neurons bez\"uglich des Stimulus hin. Eine hohe Amplitude besagt, dass es einen starken Stimulus ben\"otigt, um ein Aktionspotential hervorzurufen. Aus dem zeitlichen Versatz des STA kann die Zeit abgelesen werden, die das System braucht, um auf den Stimulus zu antworten. Der STA kann auch dazu benutzt werden, aus den Antworten der Zelle den Stimulus zu rekonstruieren (\figref{reverse_reconstruct_fig}). Bei der \determ{invertierten Rekonstruktion} wird die Zellantwort mit dem STA verfaltet. \begin{figure}[t] \includegraphics[width=\columnwidth]{reconstruction} \titlecaption{Rekonstruktion des Stimulus mittels STA.}{Die Zellantwort wird mit dem STA verfaltet, um eine Rekonstruktion des Stimulus zu erhalten.}\label{reverse_reconstruct_fig} \end{figure}