diff --git a/Makefile b/Makefile index dd62125..40ba6d3 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ BASENAME=scientificcomputing-script SUBDIRS=designpattern statistics bootstrap regression likelihood pointprocesses -#programming spike_trains +#programming SUBTEXS=$(foreach subd, $(SUBDIRS), $(subd)/lecture/$(subd).tex) pdf : chapters $(BASENAME).pdf diff --git a/header.tex b/header.tex index b88c763..0dab30c 100644 --- a/header.tex +++ b/header.tex @@ -170,14 +170,15 @@ \newcommand{\koZ}{\mathds{C}} %%%%% code/matlab commands: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\usepackage{textcomp} \newcommand{\code}[1]{\setlength{\fboxsep}{0.5ex}\colorbox{blue!10}{\texttt{#1}}} -\newcommand{\matlab}{MATLAB} +\newcommand{\matlab}{MATLAB$^{\copyright}$} \newcommand{\matlabfun}[1]{(\tr{\matlab{}-function}{\matlab-Funktion} \setlength{\fboxsep}{0.5ex}\colorbox{blue!10}{\texttt{#1}})} %%%%% exercises environment: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % usage: % -% \begin{exercise}[somecode.m] +% \begin{exercise}{somecode.m}{someoutput.out} % Write a function that computes the median. % \end{exercise} % @@ -187,6 +188,8 @@ % Write a function that computes the median. % Listing 1: somecode.m % the listing +% Output: +% content of someoutput.out % % Innerhalb der exercise Umgebung ist enumerate umdefiniert, um (a), (b), (c), .. zu erzeugen. \usepackage{ifthen} @@ -194,8 +197,7 @@ \newcounter{maxexercise} \setcounter{maxexercise}{10000} % show listings up to exercise maxexercise \newcounter{exercise}[chapter] -\setcounter{exercise}{0} -\newcommand{\theexercise}{\arabic{exercise}} +\renewcommand{\theexercise}{\arabic{exercise}} \newcommand{\codepath}{} \newenvironment{exercise}[2]% {\newcommand{\exercisesource}{#1}% diff --git a/plotting/lecture/plotting-chapter.tex b/plotting/lecture/plotting-chapter.tex index 323d374..2f260f3 100644 --- a/plotting/lecture/plotting-chapter.tex +++ b/plotting/lecture/plotting-chapter.tex @@ -11,7 +11,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} -\include{psth_sta} +\include{plotting} \end{document} diff --git a/pointprocesses/lecture/firingrates.py b/pointprocesses/lecture/firingrates.py new file mode 100644 index 0000000..e302ac8 --- /dev/null +++ b/pointprocesses/lecture/firingrates.py @@ -0,0 +1,268 @@ +import numpy as np +import matplotlib.pyplot as plt +import scipy.io as spio +import scipy.stats as spst +import scipy as sp +from IPython import embed + + +def set_axis_fontsize(axis, label_size, tick_label_size=None, legend_size=None): + """ + Sets axis, tick label and legend font sizes to the desired size. + + :param axis: the axes object + :param label_size: the size of the axis label + :param tick_label_size: the size of the tick labels. If None, lable_size is used + :param legend_size: the size of the font used in the legend.If None, the label_size is used. + + """ + if not tick_label_size: + tick_label_size = label_size + if not legend_size: + legend_size = label_size + axis.xaxis.get_label().set_fontsize(label_size) + axis.yaxis.get_label().set_fontsize(label_size) + for tick in axis.xaxis.get_major_ticks() + axis.yaxis.get_major_ticks(): + tick.label.set_fontsize(tick_label_size) + + l = axis.get_legend() + if l: + for t in l.get_texts(): + t.set_fontsize(legend_size) + + +def get_instantaneous_rate(times, max_t=30., dt=1e-4): + time = np.arange(0., max_t, dt) + indices = np.asarray(times / dt, dtype=int) + intervals = np.diff(np.hstack(([0], times))) + inst_rate = np.zeros(time.shape) + + for i, index in enumerate(indices[1:]): + inst_rate[indices[i-1]:indices[i]] = 1/intervals[i] + return time, inst_rate + + +def plot_isi_rate(spike_times, max_t=30, dt=1e-4): + times = np.squeeze(spike_times[0][0])[:50000] + time, rate = get_instantaneous_rate(times, max_t=50000*dt) + + rates = np.zeros((len(rate), len(spike_times))) + for i in range(len(spike_times)): + _, rates[:, i] = get_instantaneous_rate(np.squeeze(spike_times[i][0])[:50000], + max_t=50000*dt) + avg_rate = np.mean(rates, axis=1) + rate_std = np.std(rates, axis=1) + + fig = plt.figure() + ax1 = fig.add_subplot(311) + ax2 = fig.add_subplot(312) + ax3 = fig.add_subplot(313) + + ax1.vlines(times[times < (50000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("skpikes", fontsize=12) + set_axis_fontsize(ax1, 12) + ax1.set_xlim([0, 5]) + ax2.plot(time, rate, label="instantaneous rate, trial 1") + ax2.set_ylabel("firing rate [Hz]", fontsize=12) + ax2.legend(fontsize=12) + set_axis_fontsize(ax2, 12) + + ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue", + alpha=0.5, label="standard deviation") + ax3.plot(time, avg_rate, label="average rate") + ax3.set_xlabel("times [s]", fontsize=12) + ax3.set_ylabel("firing rate [Hz]", fontsize=12) + ax3.legend(fontsize=12) + ax3.set_ylim([0, 450]) + set_axis_fontsize(ax3, 12) + + fig.set_size_inches(15, 10) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95) + fig.set_facecolor("white") + fig.savefig("figures/instantaneous_rate.png") + plt.close() + + +def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4): + time = np.arange(0., max_t, dt) + bins = np.arange(0., max_t, bin_width) + bin_indices = bins / dt + hist, _ = sp.histogram(times, bins) + rate = np.zeros(time.shape) + + for i, b in enumerate(bin_indices[1:]): + rate[bin_indices[i-1]:b] = hist[i-1]/bin_width + return time, rate + + +def plot_bin_rate(spike_times, bin_width, max_t=30, dt=1e-4): + times = np.squeeze(spike_times[0][0]) + time, rate = get_binned_rate(times) + rates = np.zeros((len(rate), len(spike_times))) + for i in range(len(spike_times)): + _, rates[:, i] = get_binned_rate(np.squeeze(spike_times[i][0])) + avg_rate = np.mean(rates, axis=1) + rate_std = np.std(rates, axis=1) + + fig = plt.figure() + ax1 = fig.add_subplot(311) + ax2 = fig.add_subplot(312) + ax3 = fig.add_subplot(313) + + ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("skpikes", fontsize=12) + ax1.set_xlim([0, 5]) + set_axis_fontsize(ax1, 12) + + ax2.plot(time, rate, label="binned rate, trial 1") + ax2.set_ylabel("firing rate [Hz]", fontsize=12) + ax2.legend(fontsize=12) + ax2.set_xlim([0, 5]) + set_axis_fontsize(ax2, 12) + + ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue", + alpha=0.5, label="standard deviation") + ax3.plot(time, avg_rate, label="average rate") + ax3.set_xlabel("times [s]", fontsize=12) + ax3.set_ylabel("firing rate [Hz]", fontsize=12) + ax3.legend(fontsize=12) + ax3.set_xlim([0, 5]) + ax3.set_ylim([0, 450]) + set_axis_fontsize(ax3, 12) + + fig.set_size_inches(15, 10) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95) + fig.set_facecolor("white") + fig.savefig("figures/binned_rate.png") + plt.close() + + +def get_convolved_rate(times, sigma, max_t=30., dt=1.e-4): + time = np.arange(0., max_t, dt) + kernel = spst.norm.pdf(np.arange(-8*sigma, 8*sigma, dt),loc=0,scale=sigma) + indices = np.asarray(times/dt, dtype=int) + rate = np.zeros(time.shape) + rate[indices] = 1.; + conv_rate = np.convolve(rate, kernel, mode="same") + return time, conv_rate + + +def plot_conv_rate(spike_times, sigma=0.05, max_t=30, dt=1e-4): + times = np.squeeze(spike_times[0][0]) + time, rate = get_convolved_rate(times, sigma) + + rates = np.zeros((len(rate), len(spike_times))) + for i in range(len(spike_times)): + _, rates[:, i] = get_convolved_rate(np.squeeze(spike_times[i][0]), sigma) + avg_rate = np.mean(rates, axis=1) + rate_std = np.std(rates, axis=1) + + fig = plt.figure() + ax1 = fig.add_subplot(311) + ax2 = fig.add_subplot(312) + ax3 = fig.add_subplot(313) + + ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("skpikes", fontsize=12) + ax1.set_xlim([0, 5]) + set_axis_fontsize(ax1, 12) + + ax2.plot(time, rate, label="convolved rate, trial 1") + ax2.set_ylabel("firing rate [Hz]", fontsize=12) + ax2.legend(fontsize=12) + ax2.set_xlim([0, 5]) + set_axis_fontsize(ax2, 12) + + ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue", + alpha=0.5, label="standard deviation") + ax3.plot(time, avg_rate, label="average rate") + ax3.set_xlabel("times [s]", fontsize=12) + ax3.set_ylabel("firing rate [Hz]", fontsize=12) + ax3.legend(fontsize=12) + ax3.set_xlim([0, 5]) + ax3.set_ylim([0, 450]) + set_axis_fontsize(ax3, 12) + + fig.set_size_inches(15, 10) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95) + fig.set_facecolor("white") + fig.savefig("figures/convolved_rate.png") + plt.close() + + +def plot_comparison(spike_times, bin_width, sigma, max_t=30., dt=1e-4): + times = np.squeeze(spike_times[0][0]) + time, conv_rate = get_convolved_rate(times, bin_width/np.sqrt(12.)) + time, inst_rate = get_instantaneous_rate(times) + time, binn_rate = get_binned_rate(times, bin_width) + + fig = plt.figure() + ax1 = fig.add_subplot(411) + ax2 = fig.add_subplot(412) + ax3 = fig.add_subplot(413) + ax4 = fig.add_subplot(414) + ax1.spines["right"].set_visible(False) + ax1.spines["top"].set_visible(False) + ax1.yaxis.set_ticks_position('left') + ax1.xaxis.set_ticks_position('bottom') + ax2.spines["right"].set_visible(False) + ax2.spines["top"].set_visible(False) + ax2.yaxis.set_ticks_position('left') + ax2.xaxis.set_ticks_position('bottom') + ax3.spines["right"].set_visible(False) + ax3.spines["top"].set_visible(False) + ax3.yaxis.set_ticks_position('left') + ax3.xaxis.set_ticks_position('bottom') + ax4.spines["right"].set_visible(False) + ax4.spines["top"].set_visible(False) + ax4.yaxis.set_ticks_position('left') + ax4.xaxis.set_ticks_position('bottom') + + + + ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("spikes", fontsize=10) + ax1.set_xlim([1.5, 3.5]) + ax1.set_ylim([0, 1]) + ax1.set_yticks([0, 1]) + set_axis_fontsize(ax1, 10) + ax1.set_xticklabels([]) + + ax2.plot(time, inst_rate, lw=1.5, label="instantaneous rate") + ax2.set_ylabel("firing rate [Hz]", fontsize=10) + ax2.legend(fontsize=10) + ax2.set_xlim([1.5, 3.5]) + ax2.set_ylim([0, 300]) + set_axis_fontsize(ax2, 10) + ax2.set_xticklabels([]) + + ax3.plot(time, binn_rate, lw=1.5, label="binned rate") + ax3.set_ylabel("firing rate [Hz]", fontsize=10) + ax3.legend(fontsize=10) + ax3.set_xlim([1.5, 3.5]) + ax3.set_ylim([0, 300]) + set_axis_fontsize(ax3, 10) + ax3.set_xticklabels([]) + + ax4.plot(time, conv_rate, lw=1.5, label="convolved rate") + 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]) + ax4.set_ylim([0, 300]) + set_axis_fontsize(ax4, 10) + + fig.set_size_inches(7.5, 5) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95, ) + fig.set_facecolor("white") + fig.savefig("firingrates.pdf") + plt.close() + + +if __name__ == "__main__": + spike_times = spio.loadmat('lifoustim.mat')["spikes"] + # plot_isi_rate(spike_times) + # plot_bin_rate(spike_times, 0.05) + # plot_conv_rate(spike_times, 0.025) + plot_comparison(spike_times, 0.05, 0.025) + diff --git a/pointprocesses/lecture/isimethod.py b/pointprocesses/lecture/isimethod.py new file mode 100644 index 0000000..71c16e3 --- /dev/null +++ b/pointprocesses/lecture/isimethod.py @@ -0,0 +1,158 @@ +import matplotlib.pyplot as plt +import numpy as np +from IPython import embed + +def set_rc(): + plt.rcParams['xtick.labelsize'] = 8 + plt.rcParams['ytick.labelsize'] = 8 + plt.rcParams['xtick.major.size'] = 5 + plt.rcParams['xtick.minor.size'] = 5 + plt.rcParams['xtick.major.width'] = 2 + plt.rcParams['xtick.minor.width'] = 2 + plt.rcParams['ytick.major.size'] = 5 + plt.rcParams['ytick.minor.size'] = 5 + plt.rcParams['ytick.major.width'] = 2 + plt.rcParams['ytick.minor.width'] = 2 + plt.rcParams['xtick.direction'] = "out" + plt.rcParams['ytick.direction'] = "out" + + +def create_spikes(isi=0.08, duration=0.5): + times = np.arange(0., duration, isi) + times += np.random.randn(len(times)) * (isi / 2.5) + times = np.delete(times, np.nonzero(times < 0)) + times = np.delete(times, np.nonzero(times > duration)) + times = np.sort(times) + return times + + +def gaussian(sigma, dt): + x = np.arange(-4*sigma, 4*sigma, dt) + y = np.exp(-0.5 * (x / sigma)**2)/np.sqrt(2*np.pi)/sigma; + return x, y + + +def setup_axis(spikes_ax, rate_ax): + spikes_ax.spines["right"].set_visible(False) + spikes_ax.spines["top"].set_visible(False) + spikes_ax.yaxis.set_ticks_position('left') + spikes_ax.xaxis.set_ticks_position('bottom') + spikes_ax.set_yticks([0, 1.0]) + spikes_ax.set_ylim([0, 1.05]) + spikes_ax.set_ylabel("spikes", fontsize=9) + spikes_ax.text(-0.125, 1.2, "A", transform=spikes_ax.transAxes, size=10) + + rate_ax.spines["right"].set_visible(False) + rate_ax.spines["top"].set_visible(False) + rate_ax.yaxis.set_ticks_position('left') + rate_ax.xaxis.set_ticks_position('bottom') + rate_ax.set_xlabel('time[s]', fontsize=9) + rate_ax.set_ylabel('firing rate [Hz]', fontsize=9) + rate_ax.text(-0.125, 1.15, "B", transform=rate_ax.transAxes, size=10) + + +def plot_bin_method(): + dt = 1e-5 + duration = 0.5 + + spike_times = create_spikes(0.018, duration) + t = np.arange(0., duration, dt) + + bins = np.arange(0, 0.55, 0.05) + count, _ = np.histogram(spike_times, bins) + + plt.xkcd() + set_rc() + fig = plt.figure() + fig.set_size_inches(5., 2.5) + fig.set_facecolor('white') + spikes = plt.subplot2grid((7,1), (0,0), rowspan=3, colspan=1) + rate_ax = plt.subplot2grid((7,1), (3,0), rowspan=4, colspan=1) + setup_axis(spikes, rate_ax) + spikes.set_ylim([0., 1.25]) + + spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.25) + spikes.vlines(np.hstack((0,bins)), 0., 1.25, color="red", lw=1.5, linestyles='--') + for i,c in enumerate(count): + spikes.text(bins[i] + bins[1]/2, 1.05, str(c), fontdict={'color':'red'}) + spikes.set_xlim([0, duration]) + + rate = count / 0.05 + rate_ax.step(bins, np.hstack((rate, rate[-1])), where='post') + rate_ax.set_xlim([0., duration]) + rate_ax.set_ylim([0., 100.]) + rate_ax.set_yticks(np.arange(0,105,25)) + fig.tight_layout() + fig.savefig("binmethod.pdf") + plt.close() + + +def plot_conv_method(): + dt = 1e-5 + duration = 0.5 + spike_times = create_spikes(0.05, duration) + kernel_time, kernel = gaussian(0.02, dt) + + t = np.arange(0., duration, dt) + rate = np.zeros(t.shape) + rate[np.asarray(np.round(spike_times/dt), dtype=int)] = 1 + rate = np.convolve(rate, kernel, mode='same') + rate = np.roll(rate, -1) + + plt.xkcd() + set_rc() + fig = plt.figure() + fig.set_size_inches(5., 2.5) + fig.set_facecolor('white') + spikes = plt.subplot2grid((7,1), (0,0), rowspan=3, colspan=1) + rate_ax = plt.subplot2grid((7,1), (3,0), rowspan=4, colspan=1) + setup_axis(spikes, rate_ax) + + spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.5, zorder=2) + for i in spike_times: + spikes.plot(kernel_time + i, kernel/np.max(kernel), color="orange", lw=0.75, zorder=1) + spikes.set_xlim([0, duration]) + + rate_ax.plot(t, rate, color="darkblue", lw=1, zorder=2) + rate_ax.fill_between(t, rate, np.zeros(len(rate)), color="red", alpha=0.5) + rate_ax.set_xlim([0, duration]) + rate_ax.set_ylim([0, 50]) + rate_ax.set_yticks(np.arange(0,75,25)) + fig.tight_layout() + fig.savefig("convmethod.pdf") + + +def plot_isi_method(): + spike_times = create_spikes(0.09, 0.5) + + plt.xkcd() + set_rc() + fig = plt.figure() + fig.set_size_inches(5., 2.5) + fig.set_facecolor('white') + + spikes = plt.subplot2grid((7,1), (0,0), rowspan=3, colspan=1) + rate = plt.subplot2grid((7,1), (3,0), rowspan=4, colspan=1) + setup_axis(spikes, rate) + + spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.25) + spike_times = np.hstack((0, spike_times)) + for i in range(1, len(spike_times)): + t_start = spike_times[i-1] + t = spike_times[i] + spikes.annotate(s='', xy=(t_start, 0.5), xytext=(t,0.5), arrowprops=dict(arrowstyle='<->'), color='red') + + i_rate = 1./np.diff(spike_times) + + rate.step(spike_times, np.hstack((i_rate, i_rate[-1])),color="darkblue", lw=1.25, where="post") + rate.set_ylim([0, 75]) + rate.set_yticks(np.arange(0,100,25)) + + fig.tight_layout() + fig.savefig("isimethod.pdf") + + +if __name__ == '__main__': + plot_isi_method() + plot_conv_method() + plot_bin_method() diff --git a/pointprocesses/lecture/lifoustim.mat b/pointprocesses/lecture/lifoustim.mat new file mode 100644 index 0000000..b5fb8bd Binary files /dev/null and b/pointprocesses/lecture/lifoustim.mat differ diff --git a/pointprocesses/lecture/p-unit_spike_times.mat b/pointprocesses/lecture/p-unit_spike_times.mat new file mode 100644 index 0000000..75f53e6 Binary files /dev/null and b/pointprocesses/lecture/p-unit_spike_times.mat differ diff --git a/pointprocesses/lecture/p-unit_stimulus.mat b/pointprocesses/lecture/p-unit_stimulus.mat new file mode 100644 index 0000000..b6d0b25 Binary files /dev/null and b/pointprocesses/lecture/p-unit_stimulus.mat differ diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex index de9c1b6..79fe905 100644 --- a/pointprocesses/lecture/pointprocesses.tex +++ b/pointprocesses/lecture/pointprocesses.tex @@ -1,6 +1,35 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\chapter{\tr{Point processes}{Punktprozesse}} +\chapter{Analyse von Spiketrains} + +\begin{figure}[t] + \includegraphics[width=1\textwidth]{rasterexamples} + \caption{\label{rasterexamplesfig}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} + +Aktionspotentiale (``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 +spielt keine oder nur eine untergeordnete Rolle. + +Nach etwas Vorverarbeitung haben elektrophysiologischer Messungen +deshalb Listen von Spikezeitpunkten als Ergebniss - sogenannte +``Spiketrains''. Diese Messungen k\"onnen wiederholt werden und es +ergeben sich mehrere ``trials'' von Spiketrains +(\figref{rasterexamples}). + +Spiketrains sind Zeitpunkte von Ereignissen --- den Aktionspotentialen +--- und deren Analyse f\"allt daher in das Gebiet der Statistik von +sogenannten Punktprozessen. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Punktprozesse} \begin{figure}[t] \texpicture{pointprocessscetchB} @@ -28,16 +57,6 @@ erzeugt. Zum Beispiel: Dynamik des Nervensystems und des Muskelapparates erzeugt. \end{itemize} -\begin{figure}[t] - \includegraphics[width=1\textwidth]{rasterexamples} - \caption{\label{rasterexamplesfig}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} - 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 @@ -116,9 +135,9 @@ Intervalls mit sich selber). % \caption{\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: +Die Anzahl der Ereignisse (Spikes) $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$. @@ -146,7 +165,7 @@ Insbesondere ist die mittlere Rate der Ereignisse $r$ (``Spikes pro Zeit'', Feue \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 spilet bei Punktprozessen der ``Poisson Prozess''. +\"ahnliche Rolle spielt bei Punktprozessen der ``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 @@ -192,3 +211,176 @@ Der homogne Poissonprozess hat folgende Eigenschaften: \includegraphics[width=0.48\textwidth]{poissoncounthistdist100hz100ms} \caption{\label{hompoissoncountfig}Z\"ahlstatistik von Poisson Spikes.} \end{figure} + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Zeitabh\"angiger Feuerraten} + +Bisher haben wir station\"are Spiketrains betrachtet, deren Statistik +innerhlab der Analysezeit nicht ver\"andert. 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. Im folgenden Stellen wir drei Methoden +vor, mit denen eine sich zeitlich ver\"andernde Rate des Punktprozesses +abgesch\"atzt werden kann. + +Eine klassische Darstellung zeitabh\"angiger neuronaler Aktivit\"at +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. Es 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}[t] + \includegraphics[width=\columnwidth]{firingrates} + \caption{\textbf{Verschiedene Methoden die zeitabh\"angige Feuerrate + zu bestimmen. 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} + + +\paragraph{Instantane Feuerrate} +Ein sehr einfacher Weg, die zeitabh\"angige Feuerrate zu bestimmen ist +die sogenannte \textit{instantane Feuerrate}. Dabei wird die Feuerrate +aus dem Kehrwert des \textit{Interspike Intervalls}, der Zeit zwischen +zwei aufeinander folgenden Aktionspotentialen (Abbildung \ref{instrate} +A), bestimmt. Die abgesch\"atzte Feuerrate (Abbildung \ref{instrate} B) +ist g\"ultig f\"ur das gesammte Interspike Intervall +\ref{instrate}. 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. Des Weiteren ist die Feuerrate nie null, auch wenn lange keine +Aktionspotentiale generiert wurden. + +\begin{figure}[!htb] + \includegraphics[width=\columnwidth]{isimethod} + \caption{\textbf{Bestimmung des zeitabh\"angigen Feuerrate aus dem + Interspike Interval. 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} + + +\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 (Abbildung \ref{binpsth} A). Um +diese Z\"ahlungen in die Feuerrate umzurechnen muss noch mit der +Binweite normiert werden. Die bestimmte Feuerrate gilt f\"ur das +gesamte Bin (Abbildung \ref{binpsth} B). \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\"osung der Darstellung. \"Anderungen in der Feuerrate, die +innerhalb eines Bins vorkommen k\"onnen nicht aufgl\"ost werden. Die +Wahl der Binweite stellt somit eine Annahme \"uber die relevante +Zeitskala der Verarbeitung dar. Auch hier ist $r(t)$ keine +koninuierliche Funktion. + +\begin{figure}[h!] + \includegraphics[width=\columnwidth]{binmethod} + \caption{\textbf{Bestimmung des PSTH mit der Binning Methode. 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} + + +\paragraph{Faltungsmethode} +Bei der Faltungsmethode geht wird etwas anders vorgegangen um die +Feuerrate zeitaufgel\"ost zu berechnen. Die Aktionspotentialfolge wird +zun\"achst ``bin\"ar'' dargestellt. Das heisst, dass eine Antwort als +(Zeit-)Vektor dargestellt wird, 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 Gauss Kern 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 (Abbildung \ref{convrate} A). Wenn der Kern richtig normiert +wurde (Integral 1), ergibt sich die Feuerrate direkt aus der +\"Uberlagerung der Kerne (Abb. \ref{convrate} B). 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 also eine Annahme \"uber die +relevante Zeitskala. + +\begin{figure}[h!] + \includegraphics[width=\columnwidth]{convmethod} + \caption{\textbf{Schematische Darstellung der Faltungsmethode. 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} + + + +\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 mehr \"uber diesen Zusammenhang 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]{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]{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} diff --git a/pointprocesses/lecture/sta.py b/pointprocesses/lecture/sta.py new file mode 100644 index 0000000..fb0c19f --- /dev/null +++ b/pointprocesses/lecture/sta.py @@ -0,0 +1,88 @@ +import numpy as np +import matplotlib.pyplot as plt +import scipy.io as scio +from IPython import embed + + +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", 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 + + +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='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') + + 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('reconstruction.pdf') + plt.close() + + +if __name__ == "__main__": + punit_data = scio.loadmat('p-unit_spike_times.mat') + 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) diff --git a/regression/lecture/Makefile b/regression/lecture/Makefile index 71741c9..9a3574d 100644 --- a/regression/lecture/Makefile +++ b/regression/lecture/Makefile @@ -1,22 +1,29 @@ BASENAME=regression + PYFILES=$(wildcard *.py) PYPDFFILES=$(PYFILES:.py=.pdf) -pdf : $(BASENAME).pdf $(PYPDFFILES) +all : pdf + +# script: +pdf : $(BASENAME)-chapter.pdf -$(BASENAME).pdf : $(BASENAME)-chapter.tex +$(BASENAME)-chapter.pdf : $(BASENAME)-chapter.tex $(BASENAME).tex $(PYPDFFILES) pdflatex -interaction=scrollmode $< | tee /dev/stderr | fgrep -q "Rerun to get cross-references right" && pdflatex -interaction=scrollmode $< || true $(PYPDFFILES) : %.pdf : %.py python $< clean : - rm -f *~ $(BASENAME).aux $(BASENAME).log $(BASENAME).out $(BASENAME).toc $(BASENAME).nav $(BASENAME).snm $(BASENAME).vrb + rm -f *~ + rm -f $(BASENAME).aux $(BASENAME).log + rm -f $(BASENAME)-chapter.aux $(BASENAME)-chapter.log $(BASENAME)-chapter.out + rm -f $(PYPDFFILES) $(GPTTEXFILES) cleanall : clean - rm -f $(BASENAME).pdf + rm -f $(BASENAME)-chapter.pdf -watch : +watchpdf : while true; do ! make -q pdf && make pdf; sleep 0.5; done diff --git a/scientificcomputing-script.tex b/scientificcomputing-script.tex index aa66a08..6c2350a 100644 --- a/scientificcomputing-script.tex +++ b/scientificcomputing-script.tex @@ -19,10 +19,13 @@ \lstset{inputpath=programming/code} \include{programming/lectures/programming} -\chapter{Graphische Darstellung} -\graphicspath{{designpattern/lecture}{designpattern/lecture/figures}} -\lstset{inputpath=designpattern/code/} +\graphicspath{{plotting/lecture/}{plotting/lecture/images/}} +\lstset{inputpath=plotting/code/} +\include{plotting/lecture/plotting} + +\graphicspath{{designpattern/lecture/}{designpattern/lecture/figures/}} +\lstset{inputpath=designpattern/code} \include{designpattern/lecture/designpattern} \chapter{Cheat-Sheet} @@ -52,8 +55,4 @@ \renewcommand{\texinputpath}{pointprocesses/lecture/} \include{pointprocesses/lecture/pointprocesses} -\graphicspath{{spike_trains/lecture/images/}} -\lstset{inputpath=spike_trains/code/} -\include{spike_trains/lecture/psth_sta} - \end{document} diff --git a/spike_trains/lecture/psth_sta.tex b/spike_trains/lecture/psth_sta.tex index 96e5d71..c97b924 100644 --- a/spike_trains/lecture/psth_sta.tex +++ b/spike_trains/lecture/psth_sta.tex @@ -20,6 +20,8 @@ werden. Dar\"uber hinaus wird versucht die Beziehung zwischen der zeitabh\"angigen neuronalen Antwort und dem zugrundeliegenden Stimulus zu analysieren. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% The following moved to the pointprocesses chapter! \section{Darstellung der zeitabh\"angigen Feuerrate} Eine klassische Darstellung zeitabh\"angiger neuronaler Aktivit\"at