spiketrain chapter work

spike_trains/code/firing_rates.py

@@ -228,7 +228,7 @@ def plot_comparison(spike_times, bin_width, sigma, max_t=30., dt=1e-4):
     ax3.set_xticklabels([])
 
     ax4.plot(time, conv_rate, label="convolved rate")
-    ax4.set_xlabel("times [s]", fontsize=10)
+    ax4.set_xlabel("time [s]", fontsize=10)
     ax4.set_ylabel("firing rate [Hz]", fontsize=10)
     ax4.legend(fontsize=10)
     ax4.set_xlim([1.5, 3.5])

spike_trains/code/sta.py

@@ -0,0 +1,81 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.io as scio
+import seaborn as sb
+from IPython import embed
+sb.set_context("paper")
+
+
+def plot_sta(times, stim, dt, t_min=-0.1, t_max=.1):
+    count = 0
+    sta = np.zeros((abs(t_min) + abs(t_max))/dt)
+    time = np.arange(t_min, t_max, dt)
+    if len(stim.shape) > 1 and stim.shape[1] > 1:
+        stim = stim[:,1]
+    for i in range(len(times[0])):
+        times = np.squeeze(spike_times[0][i])
+        for t in times:
+            if (int((t + t_min)/dt) < 0) or ((t + t_max)/dt > len(stim)):
+                continue;
+            
+            min_index = int(np.round((t+t_min)/dt))
+            max_index = int(np.round((t+t_max)/dt))
+            snippet = np.squeeze(stim[ min_index : max_index])            
+            sta += snippet
+            count += 1
+    sta /= count
+    
+    fig = plt.figure()
+    fig.set_size_inches(5, 5)
+    fig.subplots_adjust(left=0.15, bottom=0.125, top=0.95, right=0.95, )
+    fig.set_facecolor("white")
+    
+    ax = fig.add_subplot(111)
+    ax.plot(time, sta, color="darkblue")
+    ax.set_xlabel("time [s]")
+    ax.set_ylabel("stimulus")
+    ax.xaxis.grid('off')
+    ylim = ax.get_ylim()
+    xlim = ax.get_xlim()
+    ax.plot(list(xlim), [0., 0.], zorder=1, color='darkgray', ls='--')
+    ax.plot([0., 0.], list(ylim), zorder=1, color='darkgray', ls='--')
+    ax.set_xlim(list(xlim))
+    ax.set_ylim(list(ylim))
+    fig.savefig("../lecture/images/sta.pdf")
+    plt.close()
+    return sta
+
+
+def reconstruct_stimulus(spike_times, sta, stimulus, t_max=30., dt=1e-4):
+    s_est = np.zeros((spike_times.shape[1], len(stimulus)))
+    for i in range(10):
+        times = np.squeeze(spike_times[0][i])
+        indices = np.asarray((np.round(times/dt)), dtype=int)
+        y = np.zeros(len(stimulus))
+        y[indices] = 1
+        s_est[i, :] = np.convolve(y, sta, mode='same')
+
+    plt.plot(np.arange(0, t_max, dt), stimulus[:,1], label='stimulus', color='darkblue', lw=2.)
+    plt.plot(np.arange(0, t_max, dt), np.mean(s_est, axis=0), label='reconstruction', color='silver', lw=1.5)
+    plt.xlabel('time[s]')
+    plt.ylabel('stimulus')
+    plt.xlim([0.0, 0.25])
+    plt.ylim([-1., 1.])
+    plt.legend()
+    plt.plot([0.0, 0.25], [0., 0.], color="darkgray", lw=1, ls='--', zorder=1)
+
+    fig = plt.gcf()
+    fig.set_size_inches(7.5, 5)
+    fig.subplots_adjust(left=0.15, bottom=0.125, top=0.95, right=0.95, )
+    fig.set_facecolor("white")
+    fig.savefig('../lecture/images/reconstruction.pdf')
+    plt.close()
+
+
+if __name__ == "__main__":
+    punit_data =  scio.loadmat('../../programming/exercises/p-unit_spike_times.mat')
+    punit_stim = scio.loadmat('../../programming/exercises/p-unit_stimulus.mat')
+    spike_times = punit_data["spike_times"]
+    stimulus = punit_stim["stimulus"]
+    sta = plot_sta(spike_times, stimulus, 5e-5, -0.05, 0.05)
+    reconstruct_stimulus(spike_times, sta, stimulus, 10, 5e-5)

spike_trains/lecture/psth_sta.tex

@@ -23,7 +23,7 @@ Stimulus zu analysieren.
 \section{Darstellung der zeitabh\"angigen Feuerrate}
 
 Eine klassische Darstellung zeitabh\"angiger neuronaler Aktivit\"at
-ist das sog. Peri Stimulus Zeithistogramm (peri stimulus time
+ist das sogenannte Peri Stimulus Zeithistogramm (peri stimulus time
 histogram, PSTH). Es wird der zeitliche Verlauf der Feuerrate $r(t)$
 dargestellt. Die Einheit der Feuerrate ist Hertz, das heisst, die
 Anzahl Aktionspotentiale pro Sekunde. Dabei gibt es verschiedene
@@ -32,16 +32,102 @@ Methoden diese zu bestimmen. Drei solcher Methoden sind in Abbildung
 
 \begin{figure}
   \includegraphics[width=\columnwidth]{images/psth_comparison}
-  \caption{}\label{psthfig}
+  \caption{\textbf{Verschiedene Methoden das PSTH zu bestimmen. A)}
+    Rasterplot einer einzelnen neuronalen Antwort. Jeder vertikale
+    Strich notiert den Zeitpunkt eines Aktionspotentials. \textbf{B)}
+    PSTH aus der instantanen Feuerrate bestimmt. \textbf{C)} PSTH mit
+    der Binning Methode. \textbf{D)} Feuerrate durch Faltung mit einem
+    Gauss Kern bestimmt.}\label{psthfig}
 \end{figure}
 
+\paragraph{Instantane Feuerrate}
+Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist
+die sogenannte \textit{instantane Feuerrate}. Dabie wird die Feuerrate
+aus dem Kehrwert des \textit{Interspike Intervalls}, der Zeit zwischen
+zwei aufeinanderfolgender Aktionspotentiale, bestimmt. Die
+abgesch\"atzte Feuerrate ist g\"ultig f\"ur das gesammte Interspike
+Intervall. Sie ist sehr einfach zu berechnen und hat den Vorteil keine
+Annahme \"uber eine relevante Zeitskala (der codierung oder des
+Auslesemechanismus der postsynaptischen Zelle) zu machen. $r(t)$ ist
+keine kontinuierliche Funktion, die Spr\"unge in der Feuerrate können f\"ur
+manche Analysen nachteilig sein.
 
+\paragraph{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. Um diese Z\"ahlungen in die Feuerrate
+umzurechnen muss noch mit der Binweite normiert werden. \textbf{Tipp:}
+Um die Anzahl Spikes pro Bin zu berechnen kann die \code{hist}
+Funktion benutzt werden. Das so berechnete PSTH hat wiederum eine
+stufige Form, die von der Wahl der Binweite anh\"angt. Die Binweite
+bestimmt die zeitliche Auflösung mit der Darstellung. \"Anderungen in
+der Feuerrate, die innerhalb eines Bins vorkommen koennen nicht
+aufglöst werden. Die Wahl der Binweite stellt somit eine Annahme \"uber
+die relevante Zeitskala der Verarbeitung dar. Auch hier ist $r(t)$
+keine koninuierliche Funktion.
 
-\section{Spike triggered Average}
 
-Der Spike triggered average (STA) ist der mittlere Stimulus, der zu
-einem Aktionspotential in der neuronalen Antwort f\"uhrt.
+\paragraph{Faltungsmethode}
+Bei der Faltungsmethode geht man anders vor. Die Aktionspotentialfolge
+wird ``bin\"ar'' dargestellt. Eine Antwort wird als Vektor dargestellt,
+in dem 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 Gausskern bestimmter Breite gefaltet. 
+ 
+\[r(t) = \int_{-\infty}^{\infty}d\tau \omega(\tau)\rho(t-\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.  Die Faltungsmethode f\"uhrt, anders als die anderen Methoden,
+zu einer kontinuierlichen Funktion was f\"ur spektrale Analysen von
+Vorteil sein kann. Die Wahl der Kernbreite bestimmt, \"ahnlich zur
+Binweite, die zeitliche Aufl\"osung von $r(t)$. Man macht eine Annahme
+\"uber die relevante Zeitskala.
+
+
+\section{Spike triggered Average}
+Die graphischer Darstellung der Feuerrate reicht nicht aus um den
+Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu
+analysieren.  Eine Methode mehr \"uber diesen zu erfahren ist der
+Spike triggered average (STA). Der STA ist der mittlere Stimulus, der
+zu einem Aktionspotential in der neuronalen Antwort f\"uhrt.
 
 \begin{equation}
   STA(\tau) = \frac{1}{\langle n \rangle} \left\langle  \displaystyle\sum_{i=1}^{n}{s(t_i - \tau)} \right\rangle
 \end{equation}
+
+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 mittelt (Abbildung
+\ref{stafig}). 
+
+\begin{figure}
+  \includegraphics[width=0.5\columnwidth]{images/sta}
+  \caption{\textbf{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 (Abbildung
+\ref{reverse_reconstruct_fig}). Bei der \textit{invertierten
+  Rekonstruktion} wird die Zellantwort mit dem STA gefaltet.
+
+\begin{figure}
+  \includegraphics[width=\columnwidth]{images/reconstruction}
+  \caption{\textbf{Rekonstruktion des Stimulus mittels STA.} Die
+    Zellantwort wird mit dem STA gefaltet um eine Rekonstruktion des
+    Stimulus zu erhalten.}\label{reverse_reconstruct_fig}
+\end{figure}