\documentclass[10pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\usepackage[round]{natbib}
\let\cite\citep
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\begin{document}
\title{Spike sorting using Hidden Markov Models}
\author[1]{Roger Herikstad}%
\author[2]{shihcheng}%
\affil[1]{National University of Singapore}%
\affil[2]{Affiliation not available}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\selectlanguage{english}
\begin{abstract}
Recent years have seen a steady increase in the potential number of
neurons that can be recorded simultaneously. Yet, the techniques for
extracting single unit spike trains from multi-channel extra-cellular
recordings are still being established. Recently~\cite{Herbst:2008kn},
introduced a spike sorting algorithm based on fitting a Hidden Markov
Model to the joint firing of neurons producing a recorded signal. Here,
we implement their model the Julia programming language, and investigate
the advantages of this technique over previous algorithms.%
\end{abstract}%
\sloppy
\par\null
\section*{Introduction}
{\label{631996}}
\section*{}
{\label{631996}}
\section*{Methods}
{\label{631996}}
The methods are based on \citet{Herbst2008}~and we refer to their paper
for in-depth description of the ideas behind this algorithm. Here, we
present details of our implementation of their algorithm in the Julia
programming language.
\par\null
\subsection*{Template identification}
{\label{376154}}
The first step of the algorithm was to identify the characteristic
waveform shape of all putative single units comprising a recorded
signal. To that effect, we employed the Baum-Welch algorithm (cite),
which uses an expectation-maximization (EM) approach to estimate
emission probabilities for observations emanating from a set of hidden
states.~
The state transition matrix for each template is given by,
\par\null
\[
P(s_k|s_{k-1}) = \begin{bmatrix}
1-p & p & 0 \cdots & 0\\
0 & 0 & 1 \cdots & 0 \\
\vdots & \vdots & \cdots & \vdots \\
1 & 0 & 0 \cdots & 0
\end{bmatrix}
\]
where~\emph{p} is the probability of transitioning from the silent to
the active state. Once the template enters the active state, it will go
through each of its subsequent~\emph{K}-1 states with probability 1
before returning to the silent state. In each
state~\emph{s\textsubscript{k}}, the probability of emitting an observed
value \(Y_t\)~\emph{~}at \emph{} time\emph{~t} from \emph{}
state \emph{k~}is given by a Gaussian distribution with
mean~\(\mu_k\)~and standard deviation~\(\sigma_k\).~
To estimate these probabilities, i.e. the state
transitions~\(P(s_k|s_{k-1})\) and the emission
probabilities~\(P_G(Y_t;\mu_k,\sigma)\), the algorithm uses the well-known
forward-backward algorithm to
estimate~\(\alpha_i(t) = P(Y_1,\cdots,Y_t,X_t=i | p,\mu,\sigma)\)~and~\(\beta_i(t) = P(Y_{t+1},\cdots,Y_T,X_t=i|p,\mu,\sigma)\), given recursively
by,
\par\null
\[
\begin{aligned}
\alpha_i(1) & = & \pi P_G(Y_1;\mu_i,\sigma)\\
\alpha_i(t+1) & = & P_G(Y_{t+1};\mu_i, \sigma)\sum_{j=1}^K\alpha_j(t)P(s_i|s_j)\\
\beta_i(T) & = & 1\\
\beta_i(t) & = & \sum_{j=1}^K\beta_j(t+1)P(s_j|s_i)P_G(Y_{t+1};\mu_j,\sigma)
\end{aligned}
\]
\par\null
\section*{Results~}
{\label{631996}}
\subsection*{High SNR condition}
{\label{186701}}
We first tested our algorithm on a simple problem. We created a
synthetic signal by randomly initiating two separate templates against
of background of gaussian noise. The amplitude of the templates was
chosen so that the signal-to-noise ratio (SNR) was 10. Figure ?? A shows
the resulting signal, while the two templates used to create it is shown
in Figure ?? B. We modelled a template as a damped sine function, with
amplitude a,~
To check the performance of our algorithm on these synthetic data, we
ran it 100 times, each time with a set of 7 randomly generated
templates. These templates were generated using the same function as we
used for generating the ground truth templates, but with random
parameters. Figure ?? shows an example of a random set generated this
way.~
\par\null
\subsection*{Low SNR condition}
{\label{884873}}
To test the performance in a low SNR condition, we generated a signal
consisting of the same two templates as above, but now with double the
noise. That is, we used a sigma of 0.8, with the same firing
probabilities of 0.003 and 0.001 for the two templates. The resulting
signal is shown in Figure {\ref{353854}}. We then ran
the template finding algorithm on this signal 100 times, each time
starting from 7 random ~generated initial templates. We then calculated
the fraction of resulting templates that matched the original set, and
found that only 46\% did so. In other words, in 54\% of the runs, the
templates found by the algorithm did not match the ground truth. We
further found that in 12\% of the runs, no templates were returned by
algorithm, indicating that all the initial seven templates for these
runs coalesced to noise.
The algorithm did return a single template in 37\% of the runs, and this
single template appeared to be a mix of the two ground truth templates.~
\par\null
~Overall, we found that the template finding algorithm was able to
converge to the correct templates as long as the SNR was larger than
about 4.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/figure2/figure2}
\caption{{Performance of the algorithm for when SNR is low. A) 500 ms of simulated
data. B) Ground-truth templates (solid lines) and the variance (shaded
areas) in the estimated templates from 100 runs, using only runs in
which the algorithm returned 2 templates (47\%). C) Variance around the
mean for those runs where the algorithm returned a single template
(37\%). This single template was a mean of the two ground-truth
templates.~
{\label{353854}}%
}}
\end{center}
\end{figure}
\subsection*{}
{\label{345134}}
\subsection*{High firing rate
condition}
{\label{345134}}
As the probability of state transition from silence to activity
increases, the number of overlapping spikes also increases. We tested
our algorithm for such a condition by simulating two neurons with firing
rates of 100 Hz. This meant that, for a simulated 1 second recording
with sampling rate 30kHz, we set the transition probabilities to ~
\par\null
\subsection*{Low firing rate condition}
{\label{407078}}
\subsection*{More than two cells}
{\label{481596}}
We also tested the algorithm for a signal consisting of more than two
cells. We created created three ground truth templates
(Figure~{\ref{981184}} B) with firing probabilities
0.003 (blue trace) , 0.001 (orange trace) and 0.0075 (green trace). We
then ran the template finding algorithm 100 times, each time with a
different initial templates. For 77\% of the runs, the algorithm
returned 3 templates which matched the 3 ground truth templates
(Figure~{\ref{981184}}~B and D), while for 23\% of the
runs, the algorithm returned either 2 or 4 templates
(Figure~{\ref{981184}}~C). As we can see from Figure
{\ref{981184}}~B, most of the differences between the
estimated an ground truth templates comes from relative shifts.
\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/figure-3-templates/figure-3-templates}
\caption{{Signal composed from 3 templates. A) 500 ms of simulated high pass
filtered data. B) The three templates (bold colours) with templates sets
from 77 runs (weak colours) superimposed. Each individual template was
colour-coded according to the ground truth template to which it was
closest. C) The same three ground truth ~templates as in B), but with
templates from 23 sets where the number of templates found was either 4
(7 runs) or 2 (16 runs) superimposed. D) The distribution of mean
relative square difference between the estimated and the grout truth
template from B).~
{\label{981184}}%
}}
\end{center}
\end{figure}
\subsection*{Overlap resolution}
{\label{345134}}
One of the advantages of the current spike sorting algorithm over others
is that it does overlap resolution in a seamless manner.~
\par\null
\subsection*{Hyper parameters}
{\label{345134}}
A number of hyper parameters needed to specified in our algorithm. Among
them was the probability \emph{p}\textsubscript{condense} associated
with merging templates. Setting this parameter too high resulting in
over clustering, typically because very similar templates were kept
distinct. Conversely, if set too low, we ended up with fewer templates
than the ground truth. We ran a number of simulations where we varied
this parameter. Figure ?? shows the result, where we plotted the number
of false positives (over-clustering) as well as the number of false
rejections (under-clustering) as a function of
\emph{p}\textsubscript{condense}~
\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/figure3/figure3}
\caption{{False positive rate (rate of over clustering) and false rejection rate
(under clustering) as a function of~\emph{p\textsubscript{condense}}.
While the false rejection rate did not change much for the two values
shown here, there was 30\% reduction in the number of false positives
when using~\(p_{\mathrm condense} = 0.001\)~compared to~\(p_{condense} = 0.05\)~. We saw
minor improvements when further decreasing p to 0.0001.~
{\label{829979}}%
}}
\end{center}
\end{figure}
Similarly, the algorithm has a parameter ~\(p_{small}\) which
decides when a template is considered too close to noise to represent a
single unit. Thus, a larger value of ~ resulted in final templates that
were close to the noise threshold, while a smaller value resulted in
fewer such templates. It is worth noting that we
consider~\(p_{small}\)~a somewhat less important parameter than
~\(p_{condense}\). This is because while it is easy to manually reject
noisy templates after the algorithm has completed, manually selecting
templates that should in fact be merged required re-running parts of
algorithm. ~
Nevertheless, we also ran simulations in which we
changed~\(p_{small}\)~and tracked, as in
Figure~{\ref{829979}}, the false positive and false
rejection rates. The results of these simulation are shown in
Figure~{\ref{355530}}. We found that, as expected,
decreasing~\(p_{small}\)~decreased the false positive rate, but at
the expense of increasing the false rejection rate. Indeed, we found
that decreasing the value of \(p_{small}\)~from 0.05 to 0.001
caused the false positive rate to go from 0.16 down to 0.02, while the
false rejection rate went from 0.11 up to 0.16.~
\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/figure-small-p-analysis/figure-small-p-analysis}
\caption{{False positive rate (orange) and false rejection rate (blue) as a
function of~\(p_{small}\), the threshold for rejecting a template
as noise. As ~\(p_{small}\) decreased, the false positive rate
decreased, while the false rejection rate increased.
{\label{355530}}%
}}
\end{center}
\end{figure}
\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plainnat}
\bibliography{bibliography/converted_to_latex.bib%
}
\end{document}