exercises for sta and reconstruction

pointprocesses/code/binnedRate.m

@@ -0,0 +1,19 @@
+function [time, rate] = binned_rate(spike_times, bin_width, dt, t_max)
+% Calculates the firing rate with the binning method. The hist funciton is
+% used to count the number of spikes in each bin. 
+% Arguments: 
+%           spike_times, vector containing the times of the spikes.
+%           bin_width, the width of the bins in seconds.
+%           dt, the temporal resolution.
+%           t_max, the tiral duration.
+%
+% Returns two vectors containing the time and the rate.
+
+time = 0:dt:t_max-dt;
+bins = 0:bin_width:t_max;
+rate = zeros(size(time));
+
+h = hist(spike_times, bins) ./ bin_width;
+for i = 2:length(bins)
+    rate(round(bins(i - 1) / dt) + 1:round(bins(i) / dt)) = h(i);
+end

pointprocesses/code/convolutionRate.m

@@ -0,0 +1,24 @@
+function [time, rate] = convolution_rate(spike_times, sigma, dt, t_max)
+% Calculates the firing rate with the convolution method. 
+% Arguments: 
+%           spike_times, a vector containing the spike times.
+%           sigma, the standard deviation of the Gaussian kernel in seconds.
+%           dt, the temporal resolution in seconds.
+%           t_max, the trial duration in seconds.
+%
+% Returns two vectors containing the time and the rate.
+
+time = 0:dt:t_max - dt;
+rate = zeros(size(time));
+spike_indices = round(spike_times / dt);
+rate(spike_indices) = 1;
+kernel = gauss_kernel(sigma, dt);
+
+rate = conv(rate, kernel, 'same');
+end
+
+
+function y = gauss_kernel(s, step)
+    x = -4 * s:step:4 * s;
+    y = exp(-0.5 .* (x ./ s) .^ 2) ./ sqrt(2 * pi) / s;
+end

pointprocesses/code/instantaneousRate.m

@@ -0,0 +1,22 @@
+function [time, rate] = instantaneous_rate(spike_times, dt, t_max)
+% Function calculates the firing rate as the inverse of the interspike
+% interval.
+% Arguments:
+%            spike_times, vector containing the times of the spikes.
+%            dt, the temporal resolutions of the recording.
+%            t_max, the duration of the trial.
+%
+% Returns the vector representing time and a vector containing the rate.
+
+time = 0:dt:t_max-dt;
+rate = zeros(size(time));
+
+isis = diff([0 spike_times]);
+inst_rate = 1 ./ isis;
+spike_indices = [1 round(spike_times ./ dt)];
+
+for i = 2:length(spike_indices)
+    rate(spike_indices(i - 1):spike_indices(i)) = inst_rate(i - 1);
+end
+    
+

pointprocesses/code/reconstructStimulus.m

@@ -0,0 +1,19 @@
+function s_est = reconstructStimulus(spike_times, sta, stim_duration, dt)
+% Function estimates the stimulus from the Spike-Triggered-Average
+% (sta).
+% Arguments:
+%         spike_times, a vector containing the spike times in seconds.
+%         sta, a vector containing the spike-triggered-average.
+%         stim_duration, the total duration of the stimulus.
+%         dt, the sampling interval given in seconds.
+%
+% Returns:
+%         the estimated stimulus.
+
+s_est = zeros(round(stim_duration / dt), 1);
+
+binary_spikes = zeros(size(s_est));
+binary_spikes(round(spike_times ./ dt)) = 1;
+
+s_est = conv(binary_spikes, sta, 'same');
+

pointprocesses/code/spikeTriggeredAverage.m

@@ -0,0 +1,32 @@
+function [sta, std_sta, valid_spikes] = spikeTriggeredAverage(stimulus, spike_times, count, sampling_rate)
+% Function estimates the Spike-Triggered-Average (sta).
+%
+% Arguments:
+%           stimulus, a vector containing stimulus intensities
+%           as a function of time.
+%           spike_times, a vector containing the spike times 
+%           in seconds.
+%           count, the number of datapoints that are taken around
+%           the spike times.
+%           sampling_rate, the sampling rate of the stimulus.
+%
+% Returns: 
+%           the sta, a vector containing the staandard deviation and 
+%           the number of spikes taken into account.
+
+snippets = zeros(numel(spike_times), 2*count);
+valid_spikes = 1;
+for i = 1:numel(spike_times)
+    t = spike_times(i);
+    index = round(t*sampling_rate);
+    if index <= count || (index + count) > length(stimulus)
+        continue
+    end
+    snippets(valid_spikes,:) = stimulus(index-count:index+count-1);
+    valid_spikes = valid_spikes + 1;
+end
+
+snippets(valid_spikes:end,:) = [];
+
+sta = mean(snippets, 1);
+std_sta = std(snippets,[],1);

pointprocesses/lecture/pointprocesses.tex

@@ -301,6 +301,11 @@ 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.
 
+\begin{exercise}{instantaneousRate.m}{}
+  Implementiere die Absch\"atzung der Feuerrate auf Basis der
+  instantanen Feuerrate. Plotte die Feuerrate als Funktion der Zeit.
+\end{exercise}
+
 \subsubsection{Binning-Methode}
 Bei der Binning-Methode wird die Zeitachse in gleichm\"aßige
 Abschnitte (Bins) eingeteilt und die Anzahl Aktionspotentiale, die in
@@ -330,6 +335,10 @@ gemacht.
     Zeithistogramm mit der Binweite normiert.}\label{binpsth}
 \end{figure}
 
+\begin{exercise}{binnedRate.m}{}
+  Implementiere die Absch\"atzung der Feuerrate mit der ``binning''
+  Methode. Plotte das PSTH.
+\end{exercise}
 
 \subsubsection{Faltungsmethode}
 Bei der Faltungsmethode werden die harten Kanten der Bins der
@@ -365,6 +374,11 @@ 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.
 
+\begin{exercise}{convolutionRate.m}{}
+  Verwende die Faltungsmethode um die Feuerrate zu bestimmen. Plotte
+  das Ergebnis.
+\end{exercise}
+
 \section{Spike-triggered Average}
 Die graphischer Darstellung der Feuerrate allein reicht nicht aus um
 den Zusammenhang zwischen neuronaler Antwort und einem Stimulus zu
@@ -382,11 +396,15 @@ 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}
+  \includegraphics[width=\columnwidth]{sta}
+  \titlecaption{Spike-triggered Average eines P-Typ Elektrorezeptors
+    und Stimulusrekonstruktion.}{\textbf{A)} Der STA: der Rezeptor
+    wurde mit einem ``white-noise'' Stimulus getrieben. Zeitpunkt 0
+    ist der Zeitpunkt des beobachteten Aktionspotentials. Die Kurve
+    ergibt sich aus dem gemittelten Stimulus je 50\,ms vor und nach
+    einem Aktionspotential. \textbf{B)} Stimulusrekonstruktion mittels
+    STA. Die Zellantwort wird mit dem STA gefaltet um eine
+    Rekonstruktion des Stimulus zu erhalten.}\label{stafig}
 \end{figure}
 
 Aus dem STA k\"onnen verschiedene Informationen \"uber den
@@ -402,13 +420,26 @@ 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
+Stimulus zu rekonstruieren (\figref{stafig} B). 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}
+\begin{exercise}{spikeTriggeredAverage.m}{}
+  Implementiere eine Funktion, die den STA ermittelt. Verwende dazu
+  den Datensatz \codeterm{sta\_data.mat}. Die Funktion sollte folgende
+  R\"uckgabewerte haben:
+  \begin{itemize}
+  \item den Spike-Triggered-Average.
+  \item die Standardabweichung der individuellen STAs.
+  \item die Anzahl Aktionspotentiale, die dem STA zugrunde liegen.
+  \end{itemize}
+\end{exercise}
+
+\begin{exercise}{reconstructStimulus.m}{}
+  Rekonstruiere den Stimulus mithilfe des STA und der Spike
+  Zeiten. Die Funktion soll Vektor als R\"uckgabewert haben, der
+  genauso gro{\ss} ist wie der Originalstimulus aus der Datei
+  \codeterm{sta\_data.mat}.
+\end{exercise}
+
+

pointprocesses/lecture/sta.py

@@ -22,31 +22,7 @@ def plot_sta(times, stim, dt, t_min=-0.1, t_max=.1):
             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", lw=1)
-    ax.set_xlabel("time [s]")
-    ax.set_ylabel("stimulus")
-    ax.xaxis.grid('off')
-    ax.spines["right"].set_visible(False)
-    ax.spines["top"].set_visible(False)
-    ax.yaxis.set_ticks_position('left')
-    ax.xaxis.set_ticks_position('bottom')
-
-    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("sta.pdf")
-    plt.close()
-    return sta
+    return time, sta
 
 
 def reconstruct_stimulus(spike_times, sta, stimulus, t_max=30., dt=1e-4):
@@ -57,25 +33,57 @@ def reconstruct_stimulus(spike_times, sta, stimulus, t_max=30., dt=1e-4):
         y = np.zeros(len(stimulus))
         y[indices] = 1
         s_est[i, :] = np.convolve(y, sta, mode='same')
+    time = np.arange(0, t_max, dt)
+    return time, np.mean(s_est, axis=0) 
 
-    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='gray', 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.5, ls='--', zorder=1)
-    plt.gca().spines["right"].set_visible(False)
-    plt.gca().spines["top"].set_visible(False)
-    plt.gca().yaxis.set_ticks_position('left')
-    plt.gca().xaxis.set_ticks_position('bottom')
+
+def plot_results(sta_time, st_average, stim_time, s_est, stimulus, duration, dt):
+    sta_ax = plt.subplot2grid((1, 3), (0, 0), rowspan=1, colspan=1)
+    stim_ax = plt.subplot2grid((1, 3), (0, 1), rowspan=1, colspan=2)
     
     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_size_inches(15, 5)
+    fig.subplots_adjust(left=0.075, bottom=0.12, top=0.92, right=0.975)
     fig.set_facecolor("white")
-    fig.savefig('reconstruction.pdf')
+
+    sta_ax.plot(sta_time * 1000, st_average, color="dodgerblue", lw=2.)
+    sta_ax.set_xlabel("time [ms]", fontsize=12)
+    sta_ax.set_ylabel("stimulus", fontsize=12)
+    sta_ax.set_xlim([-50, 50])
+    #  sta_ax.xaxis.grid('off')
+    sta_ax.spines["right"].set_visible(False)
+    sta_ax.spines["top"].set_visible(False)
+    sta_ax.yaxis.set_ticks_position('left')
+    sta_ax.xaxis.set_ticks_position('bottom')
+    sta_ax.spines["bottom"].set_linewidth(2.0)
+    sta_ax.spines["left"].set_linewidth(2.0)
+    sta_ax.tick_params(direction="out", width=2.0)
+
+    ylim = sta_ax.get_ylim()
+    xlim = sta_ax.get_xlim()
+    sta_ax.plot(list(xlim), [0., 0.], zorder=1, color='darkgray', ls='--', lw=0.75)
+    sta_ax.plot([0., 0.], list(ylim), zorder=1, color='darkgray', ls='--', lw=0.75)
+    sta_ax.set_xlim(list(xlim))
+    sta_ax.set_ylim(list(ylim))
+    sta_ax.text(-0.225, 1.05, "A", transform=sta_ax.transAxes, size=14)
+
+    stim_ax.plot(stim_time * 1000, stimulus[:,1], label='stimulus', color='dodgerblue', lw=2.)
+    stim_ax.plot(stim_time * 1000, s_est, label='reconstruction', color='red', lw=2)
+    stim_ax.set_xlabel('time[ms]', fontsize=12)
+    stim_ax.set_xlim([0.0, 250])
+    stim_ax.set_ylim([-1., 1.])
+    stim_ax.legend()
+    stim_ax.plot([0.0, 250], [0., 0.], color="darkgray", lw=0.75, ls='--', zorder=1)
+    stim_ax.spines["right"].set_visible(False)
+    stim_ax.spines["top"].set_visible(False)
+    stim_ax.yaxis.set_ticks_position('left')
+    stim_ax.xaxis.set_ticks_position('bottom')
+    stim_ax.spines["bottom"].set_linewidth(2.0)
+    stim_ax.spines["left"].set_linewidth(2.0)
+    stim_ax.tick_params(direction="out", width=2.0)
+    stim_ax.text(-0.075, 1.05, "B", transform=stim_ax.transAxes, size=14)
+    
+    fig.savefig("sta.pdf")
     plt.close()
 
 
@@ -84,5 +92,6 @@ if __name__ == "__main__":
     punit_stim = scio.loadmat('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)
+    sta_time, sta = plot_sta(spike_times, stimulus, 5e-5, -0.05, 0.05)
+    stim_time, s_est = reconstruct_stimulus(spike_times, sta, stimulus, 10, 5e-5)
+    plot_results(sta_time, sta, stim_time, s_est, stimulus, 10, 5e-5)