diff --git a/programming/exercises/STA/p-unit_spike_times.mat b/programming/exercises/STA/p-unit_spike_times.mat new file mode 100644 index 0000000..75f53e6 Binary files /dev/null and b/programming/exercises/STA/p-unit_spike_times.mat differ diff --git a/programming/exercises/STA/p-unit_stimulus.mat b/programming/exercises/STA/p-unit_stimulus.mat new file mode 100644 index 0000000..b6d0b25 Binary files /dev/null and b/programming/exercises/STA/p-unit_stimulus.mat differ diff --git a/programming/exercises/STA/sta.m b/programming/exercises/STA/sta.m new file mode 100644 index 0000000..0b4e3cb --- /dev/null +++ b/programming/exercises/STA/sta.m @@ -0,0 +1,18 @@ +function [st_avg, std_sta, valid_spikes]= sta(stimulus, spike_times, count, sampling_rate) + +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(end-(end-valid_spikes):end,:) = []; + +st_avg = mean(snippets, 1); +std_sta = std(snippets,[],1); \ No newline at end of file diff --git a/programming/exercises/STA/sta_script.m b/programming/exercises/STA/sta_script.m new file mode 100644 index 0000000..432f971 --- /dev/null +++ b/programming/exercises/STA/sta_script.m @@ -0,0 +1,67 @@ +clear +clc + +%% load data +load p-unit_spike_times.mat +load p-unit_stimulus.mat +sample_rate = 20000; + +%% calculate STA +all_times = []; +for i = 1:length(spike_times) + all_times = cat(1, all_times, spike_times{i}); +end + +[st_average, sta_sd, num] = sta(stimulus(:,2), all_times, 1000, sample_rate); +fig = figure(); +set(fig, 'PaperUnits', 'centimeters'); +set(fig, 'PaperSize', [11.7 9.0]); +set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]); +set(fig,'Color', 'white') +plot(((1:length(st_average))-1000)/sample_rate, st_average) +xlabel('time [s]') +ylabel('stimulus') + +%% reverse reconstruction + +% make binary representation of the spike times +binary_spikes = zeros(size(stimulus, 1), length(spike_times)); +estimated_stims = zeros(size(binary_spikes)); +for i = 1:length(spike_times) + binary_spikes(round(spike_times{i}*sample_rate), i) = 1; + estimated_stims(:,i) = conv(binary_spikes(:,i), st_average, 'same'); +end + +fig = figure(); +set(fig, 'PaperUnits', 'centimeters'); +set(fig, 'PaperSize', [11.7 9.0]); +set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]); +set(fig,'Color', 'white') +plot(stimulus(:,1), stimulus(:,2), 'displayname','original') +hold on +plot(stimulus(:,1), mean(estimated_stims,2), 'r', 'displayname', 'reconstruction') +xlabel('time [s]') +ylabel('stimulus') +legend show + +%% calculate STC + +% we need to downsample the data otherwise the covariance matrixs gets too +% large 20Khz to 1kHz + +% downsampled_binary = zeros(size(stimulus, 1)/20, length(spike_times)); +downsampled_stim = zeros(size(stimulus,1)/20,1); + +% for i = 1:length(spike_times) +% indices = round(spike_times{i}.*1000); +% indices(indices < 1) = []; +% downsampled_binary(indices, i) = 1; +% end +for i = 1:length(downsampled_stim) + start_index = (i-1) * 20 + 1; + downsampled_stim(i) = mean(stimulus(start_index:start_index+19,2)); +end + +[st_average, ~, ~] = sta(downsampled_stim, all_times, 50, 1000); + + diff --git a/programming/lectures/sta_stc.tex b/programming/lectures/sta_stc.tex new file mode 100644 index 0000000..e789e20 --- /dev/null +++ b/programming/lectures/sta_stc.tex @@ -0,0 +1,380 @@ +\documentclass{beamer} +\usepackage{xcolor} +\usepackage{listings} +\usepackage{pgf} +%\usepackage{pgf,pgfarrows,pgfnodes,pgfautomata,pgfheaps,pgfshade} +%\usepackage{multimedia} + +\usepackage[english]{babel} +\usepackage{movie15} +\usepackage[latin1]{inputenc} +\usepackage{times} +\usepackage{amsmath} +\usepackage{bm} +\usepackage[T1]{fontenc} +\usepackage[scaled=.90]{helvet} +\usepackage{scalefnt} +\usepackage{tikz} +\usepackage{ textcomp } +\usepackage{soul} +\usepackage{hyperref} +\definecolor{lightblue}{rgb}{.7,.7,1.} +\definecolor{mygreen}{rgb}{0,1.,0} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\mode +{ + \usetheme{Singapore} + \setbeamercovered{opaque} + \usecolortheme{tuebingen} + \setbeamertemplate{navigation symbols}{} + \usefonttheme{default} + \useoutertheme{infolines} + % \useoutertheme{miniframes} +} + +\AtBeginSection[] +{ + \begin{frame} + \begin{center} + \Huge \insertsectionhead + \end{center} + % \frametitle{\insertsectionhead} + % \tableofcontents[currentsection,hideothersubsections] + \end{frame} +} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 + +\setbeamertemplate{blocks}[rounded][shadow=true] + +\title[]{Introduction to Scientific Computing -- \\ +Cross-Correlation, Spike--Triggered--Average and Reverse Reconstruction} +\author[]{Jan Grewe\\Abteilung f\"ur Neuroethologie\\ + Universit\"at T\"ubingen} + +\institute[Introduction to Scientific Computing]{} + \date{2014/10/13 - 2014/11/07} + %\logo{\pgfuseimage{../../resources/UT_BM_Rot_RGB.pdf}} + +\subject{Einf\"uhrung in wissenschaftliche Datenverarbeitung} +\vspace{1em} +\titlegraphic{ + \includegraphics[width=0.5\linewidth]{../../resources/UT_WBMW_Rot_RGB} +} +%%%%%%%%%% configuration for code +\lstset{ + basicstyle=\ttfamily, + numbers=left, + showstringspaces=false, + language=Matlab, + commentstyle=\itshape\color{darkgray}, + keywordstyle=\color{blue}, + stringstyle=\color{green}, + backgroundcolor=\color{blue!10}, + breaklines=true, + breakautoindent=true, + columns=flexible, + frame=single, + captionpos=b, + xleftmargin=1em, + xrightmargin=1em, + aboveskip=10pt + } +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\newcommand{\mycite}[1]{ +\begin{flushright} +\tiny \color{black!80} #1 +\end{flushright} +} + +\input{../../latex/environments.tex} +\makeatother + +\begin{document} + +\begin{frame}[plain] + \frametitle{} + \vspace{-1cm} + \titlepage % erzeugt Titelseite +\end{frame} + + +\begin{frame}[plain] + \frametitle{Rekapitulation} + \begin{enumerate} + \item PSTH\pause + \end{enumerate} +\end{frame} + + +\begin{frame} + \frametitle{Introduction to scientific computing} + \frametitle{Menue} + \begin{enumerate} + \item Cross-correlation + \item Estimation of the Spike--Triggered--Average -- STA + \item Reverse reconstruction + \end{enumerate} +\end{frame} + + +\begin{frame}[plain] + \huge{1. Recapitulation: PSTH} +\end{frame} + + +\begin{frame} + \frametitle{Recapitulation} + \only<1>{ + \framesubtitle{Displaying the neuronal response over time - Rasterplot} + \begin{figure} + \centering + \includegraphics[width=0.375\columnwidth]{images/rasterplot} + \end{figure} + } + \only<2>{ + \framesubtitle{Displaying the neuronal response over time - PSTH} + \begin{figure} + \centering + \includegraphics[width=0.5\columnwidth]{images/conv} + \end{figure} + } +\end{frame} + + +\begin{frame}[plain] + \huge{2. Relating stimulus and response} +\end{frame} + + +\begin{frame} + \frametitle{Relating stimulus and response} + \framesubtitle{How can we relate the response to the stimulus?} + \begin{figure} + \centering + \includegraphics[height=0.9\textheight]{images/conv_stim} + \end{figure} +\end{frame} + + +\begin{frame} + \frametitle{Relating stimulus and response} + \framesubtitle{Cross--correlation} + The cross - correlation for (time) discrete data is defined as: + \begin{equation} + Q_{rs}(\tau) = \frac{1}{2\cdot N-|\tau|} \displaystyle\sum^{T}_{\tau=-T}{r(t) \cdot s(t-\tau)} \cdot \Delta \tau + \end{equation} + \pause + Don't worry, this is already implemented in MATLAB! +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Cross--correlation} + \footnotesize + \begin{lstlisting} +[c,lags] = xcorr(stimulus, response, 10000, 'unbiased'); +fig = figure(); +set(fig, 'PaperUnits', 'centimeters'); +set(fig, 'PaperSize', [11.7 9.0]); +set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]); +set(fig,'Color', 'white') +plot(lags/sample_rate, c) +xlabel('lag [s]') +ylabel('correlation') + \end{lstlisting} +\end{frame} + + +\begin{frame} + \frametitle{Relating stimulus and response} + \framesubtitle{Cross--correlation} + \begin{figure} + \includegraphics[width=0.5\linewidth]{images/correlation} + \end{figure} +\end{frame} + + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Cross--correlation - Exercises} + \begin{enumerate} + \item calculate the cross-correlation between two vectors of random + numbers. + \item Calculate the cross-correlation between one of these vectors + and itself (auto-correlation). + \item Calculate the cross-correlation between one vector and a + time-shifted version of itself (use \verb+circshift+ to do this). + \end{enumerate} + \textbf{Note:} Select max\_lag to be less than 10\% of the length of + your vectors! + \begin{enumerate} + \item Create the cross correlation of the p-unit data and stimulus. + \item \textbf{Note:} you have to convert the spike\_times to a PSTH! + \item Find out the position of the correlation peak. + \item What does this tell you? + \end{enumerate} +\end{frame} + + +\begin{frame}[plain] + \huge{2. Spike--Triggered--Average --- STA} +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Spike--Triggered--Average --- STA} + + The \textbf{STA} is the average stimulus that led to a spike in the + neuron. + + \begin{equation} + STA(\tau) = \left\langle \frac{1}{n} \displaystyle\sum^{i=1}_{n}{s(t_i - \tau)} \right\rangle + \end{equation} + + Sadly, we need to implement this ourselves.\newline\newline\pause + \vspace{1em} + \noindent + \textbf{Algorithm:} + \begin{enumerate} + \item For each spike a snippet is taken from the respective stimulus. + \item The snippets are averaged. + \end{enumerate} +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Spike--Triggered--Average -- STA} + \vspace{-1em} + \begin{figure} + \centering + \includegraphics[width=0.65\columnwidth]{images/sta} + \end{figure} + \pause + \vspace{-0.5em} + What does the \textbf{STA} tell us? + \begin{enumerate} + \item Is there a relation between stimulus and response? \pause + \item Is there a lag between them and how large is it? \pause + \item How far in the past does a neuron encode? \pause + \item Can it see into the future? + \end{enumerate} +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{STA -- Exercises} + \textbf{Exercise:} + \begin{enumerate} + \item Write a function \verb+sta(x, y, count, sample_rate)+ that + takes the stimulus (x), the response (y, as spike times), the + number (count) of sampling points it should cut out from the + stimulus and the sampling\_rate to convert from times to + indices. + \item \textbf{Beware:} sometimes the spike\_time may be too close + to the beginning or the end of the stimulus to cut out enough + data. + \item Calculate the STA for P-unt data. + \end{enumerate} +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Spike--Triggered--Average -- STA} + What does the \textbf{STA} tell us? + \begin{figure} + \centering + \includegraphics[width=0.25\columnwidth]{images/sta} + \end{figure} + \begin{enumerate} + \item Is there a relation between stimulus and response?\pause + \item Is there a lag between them and how large is it?\pause + \item How far in the past does a neuron encode?\pause + \item Can it see into the future? + \end{enumerate} +\end{frame} + + +\begin{frame}[plain] + \huge{3. Reverse reconstruction using the \textbf{STA}} +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Reverse reconstruction using the \textbf{STA}} + \textbf{Basic idea:}\newline + \begin{itemize} + \item The \textbf{STA} is the average stimlus that led to a spike. \pause + \item The other way round: Whenever we observe a spike, the stimulus + was likely to have looked like the \textbf{STA}.\pause + \item Thus, we should be able to reconstruct how the stimulus looked + like to evoke the observed response. + \end{itemize} +\end{frame} + + +\begin{frame}[fragile] + \frametitle{Relating stimulus and response} + \framesubtitle{Reverse reconstruction using the \textbf{STA}} + \textbf{Algorithm:}\newline + \begin{enumerate} + \item Estimate the \textbf{STA}. \pause + \item Create a binary representation of the spike train (vector of + \verb+zeros+ in which each \textbf{1} signals the occurence of a + spike).\pause + \item Use the STA as a Kernel (compare to the spike convolution + method to create the PSTH) and convolve (\verb+conv+) it with the + \textbf{STA}. + \end{enumerate} + \textbf{Exercise:} + \begin{enumerate} + \item Reconstruct the stimulus from the P-unit data. + \item Plot original and reconstructed stimulus into the same plot. + \item Are they similar? In which instances are they different? + \end{enumerate} +\end{frame} + +\end{document} + +\begin{frame}[plain] + \huge{4. Spike--Triggered--Covariance \textbf{STC}} +\end{frame} + + +\begin{frame} + \frametitle{Relating stimulus and response} + \framesubtitle{Spike--Triggered--Covariance \textbf{STC}} + \textbf{Problem:} + \begin{itemize} + \item The \textbf{STA} is a linear filter. + \item It can only reveal linear relationships between stimulus and response.\pause + \end{itemize} + \vspace{1em} + Further relations can be recovered using the Spike-triggered-covariance.\pause + \begin{equation} + STC = \frac{1}{N-1} \cdot \displaystyle\sum_{n=1}^{N}[\overrightarrow{s}(t_n) - STA] + [\overrightarrow{s}(t_n) - STA]^T + \end{equation} + where:\\ $N$ is the total number of spikes, + $\overrightarrow{s}(t_n)$ the stimulus snippet for the $n^{th}$ + spike, and $STA$ the Spike-Triggered-Average. +\end{frame} + + +\begin{frame} + \frametitle{Relating stimulus and response} + \framesubtitle{Spike--Triggered--Covariance} + \begin{itemize} + \item + \end{itemize} +\end{frame} + + +\end{document}