diff --git a/spike_trains/code/firing_rates.py b/spike_trains/code/firing_rates.py index 0f02e87..097250e 100644 --- a/spike_trains/code/firing_rates.py +++ b/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]) diff --git a/spike_trains/code/sta.py b/spike_trains/code/sta.py new file mode 100644 index 0000000..5b14083 --- /dev/null +++ b/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) diff --git a/spike_trains/lecture/images/reconstruction.pdf b/spike_trains/lecture/images/reconstruction.pdf new file mode 100644 index 0000000..7f5f663 Binary files /dev/null and b/spike_trains/lecture/images/reconstruction.pdf differ diff --git a/spike_trains/lecture/images/sta.pdf b/spike_trains/lecture/images/sta.pdf index 625ec9e..db21b29 100644 Binary files a/spike_trains/lecture/images/sta.pdf and b/spike_trains/lecture/images/sta.pdf differ diff --git a/spike_trains/lecture/psth_sta.tex b/spike_trains/lecture/psth_sta.tex index 57ea7a3..73d8fe1 100644 --- a/spike_trains/lecture/psth_sta.tex +++ b/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}