“f’” LlBRARY
O l 0 Michigan State
’ University

 

 

 

This is to certify that the
thesis entitled

NEUROQUEST: A COMPREHENSIVE TOOL FOR LARGE-
SCALE NEURAL DATA ANALYSIS

presented by

Kl YONG KWON

has been accepted towards fulfillment
of the requirements for the

MS. degree in ELECTRICAL ENGINEERING

 

 

 

WSW Signature

”I-

Date

MSU is an Affirmative Action/Equal Opportunity Employer

_ _... -.._ — . _. . _.__ . _.—.—o—.-o—o—o---o-o-.-.—-—a—--- .—

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 ICIProj/Achres/CIRCIDateDueJndd

NEUROQUEST: A COMPREHENSIVE TOOL FOR LARGE-SCALE NEURAL
DATA ANALYSIS

By

Ki Yong Kwon

A THESIS
Submitted to
Michigan State University
In partial fulfillment of the requirements
for the degree of
MASTER OF SCIENCE
Electrical Engineering

2010

ABSTRACT

NEUROQUEST: A COMPREHENSIVE TOOL FOR LARGE SCALE NEURAL
DATA ANALYSIS

By
Ki Yong Kwon
Processing the massive amounts of neural data to extract biologically relevant
information from the activity of large ensembles of neurons in noisy recordings is
a major challenge. Efficient and practical software development is indispensable
to deal with these challenges, and many scientific findings will rest on the ability
to overcome these challenges.

We developed a comprehensive MATLAB-based software package -
entitled NeuroQuest - that bundles a number of advanced neural signal
processing algorithms together in a user-friendly environment.

In this thesis, we also proposed novel spike detection and feature
extraction algorithms that are fully integrated in NeuroQuest. These methods
capture a number of useful features in a sparse representation domain of the
neural signals. These sparse "footprints" permit efficient and reliable identification
of neural spike events, even in the presence of interfering signals. Results
demonstrate the efficiency and reliability of these methods compared to other
similar methods and software packages, and versatility to a wide range of

experimental conditions.

TABLE OF CONTENTS

LIST OF TABLES ............................................................................. vi
LIST OF FIGURES .......................................................................... vii
1. Introduction ................................................................................. 1
1.1 Study of Neurophysiological Signals ........................................ 1

1.2 Existing software packages and their limitations ......................... 3

1.3 Thesis contributions .............................................................. 4

1.4 Organization ....................................................................... 6

2. Background ................................................................................. 7
2.1 Neurophysiological signals and the characteristics ...................... 7

2.2 Extracting information from neural activity ..................................... 10

2.2.1 Pre-Conditioning ........................................................ 10

2.2.2 Spike Detection ......................................................... 11

2.2.3 Spike Sorting ............................................................ 13

2.2.4 Spike Train Analysis ................................................... 15

2.2.4.1 Single Neuron Properties .................................. 15

2.2.4.2 Multi Neuron Properties .................................... 16

3. Overview of NeuroQuest ............................................................. 18
3.1 Structure and Organization of NeuroQuest ............................... 18

3.2 The Contributions of NeuroQuest .......................................... 19

3.2.1 Comprehensive toolbox ................................................ 22

3.2.2 Simple GUI and Visualization tools ................................ 23

iii

3.3

Chapter 4.
4. 1
4.2

3.2.3 Capability of handling large-scale neural data .................. 24

3.2.4 Spike Detection and Sorting methods ............................ 26
3.2.5 Developer friendly environment .................................... 27
Summary of individual modules ............................................. 27
3.3.1 Data format ............................................................... 27
3.3.2 Pre-Conditioning ......................................................... 28
3.3.2.1 LFP Removal ................................................. 29
3.3.2.2 Artifact Removal ............................................. 30
3.3.2.3 Wavelet Denoising ........................................... 30
3.3.3 Spike Detection ........................................................... 32
3.3.4 Spike Sorting .............................................................. 33
3.3.4.1 Spike Extraction and alignment .......................... 35
3.3.4.2 Feature Extraction ........................................... 37
3.3.4.3 Clustering ...................................................... 38
3.3.4.4 Sub-Clustering ................................................ 39
3.3.4.5 Extra spike sorting tools .................................... 39
3.3.5 Spike train analysis ..................................................... 41
3.3.5.1 Single Unit Analysis tools .................................. 41
3.3.5.2 Multi Unit Analysis tools .................................... 43
3.3.5.3 Advanced Analysis tools ................................... 46
Spike detection and feature extraction ...................... 48
Introduction ....................................................................... 48
Theory ............................................................................. 51
4.2.1 Observation Model ...................................................... 51
4.2.2 Wavelet Transform and its Properties .............................. 52

iv

4.2.3 Optimal Threshold Selection for 3 x2 Distributed Signal ...... 57

4.2.4 Wavelet Footprint ........................................................ 60

4.3 Results ................................................................................. 62

4.4 Occlusion .............................................................................. 71
APPENDIX ..................................................................................... 73
BIBLIOGRAPHY......................................... .................................... 75

LIST OF TABLES

TABLE 1.1 EXISTING SOFTWARE PACKAGES ....................................... 3
TABLE 3.1 FEATURES OF THE SOFTWARE PACKAGES ........................ 23
TABLE 3.2 SPIKE DETECTION RESULTS OF NEUROQUEST AND

OSORT WITH FOUR SIMULATION DATASETS ............................... 25
TABLE 4.1 SPIKE DETECTION RESULT OF EBD AND WD ...................... 66

TABLE 4.2 CONFUSION MATRICES OF THE SPIKE SORTING
RESULTS .......................................................................................... 69

vi

LIST OF FIGURES

Figure 2—1 Types of neurophysiological signals ......................................... 8
Figure 2-2 Artifacts and action potentials in the extracellular recording ......... 11
Figure 2-3 Processes to obtain spike train from extracellular recording. . . . .....13
Figure 3-1 Organization of NeuroQuest ................................................. 20
Figure 3-2 Flowchart of NeuroQuest ..................................................... 21
Figure 3-3 Comparison of Spike Detection performance between

NeuroQuest and OSort ................................................................ 25
Figure 3-4 Pre-Conditioning GUI .......................................................... 29
Figure 3-5 Spike Detection GUI ............................................................ 33
Figure 3-6 Flowchart of Spike Sorting Algorithm in NeuroQuest .................. 34
Figure 3-7 Spike Sorting GUI ............................................................... 36
Figure 3-8 Clusters of detected spikes in temporal PCA feature space

with different alignment methods ..................................................... 36
Figure 3-9 Clusters of extracted spikes in different feature spaces .............. 38
Figure 3-10 Scheme of Sub-Clustering ................................................... 40
Figure 3-11 Example of Cluster Merging Tool ........................................... 40
Figure 3-12 Single Unit Analysis GUI ...................................................... 43

vii

Figure 3-13 Multi Unit Analysis GUI ........................................................ 45
Figure 3-14 Advanced Spike Train Analysis GUI ....................................... 47

Figure 4-1 Block Diagram of the proposed wavelet footprint detection
method .............................................................................................. 51

Figure 4-2 Sub-space reorientation of action potentials in wavelet domain....55

Figure 4-3 PDF and CDF of X2 distribution with different parameter v
and the equalized CDPs using the general histogram equalization (HE)

and the modified histogram equalization (MHE) ................................. 59
Figure 4-4 Histogram and modified histogram of sufficient statistic T............61
Figure 4»5 Wavelet footprint extraction method for the correlated

multi-channel signal ..................................................................... 63
Figure 46 R00 for different datasets with different SNR ........................... 65
Figure 4-7 Sufficient statistic T of a segment of dataset 3 and the different

threshold selections ..................................................................... 67
Figure 4-8 Cluster of temporal PCA and wavelet footprint .......................... 69
Figure 49 Spike sorting results of spontaneous recordings in

an anesthetized rat ....................................................................... 70
Figure 4-10 Spike sorting result of the muIti-channel dataset ........................ 71

viii

CHAPTER 1

Introduction

Recording and analyzing neurophysiological signals provide a window to
understand how the brain interprets the physical world. This thesis presents
NeuroQuest, a new software package for analyzing neurophysiological data, as
well as novel spike detection and feature extraction methods implemented in the
software. NeuroQuest is designed to assist researchers to translate recorded
neural activity into quantitative measurements by combining a number of
advanced neural signal processing algorithms in a unified Graphical User
Interface (GUI). NeuroQuest is equipped to play an import role in processing the
massive amounts of neural data to extract biologically relevant information about
brain function. The algorithms are designed to improve the quality of the analysis.
This chapter provides introduction of the study of neurophysiological signals,
motivation for this work, the contributions, and an overview of the following

chapters.

1.1 Study of Neurophysiological signals

Electrophysiology is the study of the electrical properties of biological cells and

tissues. It involves measurements of voltage or current on a wide variety of

scales from single ion channel proteins to whole organs. In neuroscience, the
study includes measurements of the electrical activity of neurons [1]. Recording
action potentials - or spikes - emitted by neurons to signal each other is key to
understanding the complex relationship between physical features of the world
and the brain’s interpretation of those features [2]. Neural data in the form of
spikes from a single neuron or from p0pulations of neurons have been
extensively used in various studies, for e.g., to measure connectivity between
different brain areas, to examine the dynamics of neuronal response and its
relationship to external stimuli, to design Brain-machine Interface (BMI) systems,
to clinically diagnose and treat brain disorders such as Parkinson’s disease (PD)
and certain types of epilepsy, and many other research areas. Despite the
significance of extracellular recordings, there are many technical challenges in
recording and processing the data [3-6]. For example, extracellular recordings
require a surgery to implant the recording electrodes which causes a risk of
complications and permanent tissue damages. Neurons can be easily injured,
and perturbation of their immediate environment can affect their firing patterns
[2]. In addition, recordings are typically contaminated by artifacts and other noise
components. Therefore, the data must go through multiple steps of processing
before any statistical analysis or information extraction can be done [7]. Despite
the risk of the surgical procedure and the commonly observed instability in the
recordings, studies of extracellular recordings seem indispensable to many basic

and translational neuroscience researches [2].

1.2 Existing software packages and their limitations

Many techniques and algorithms have been developed to analyze extracellular
recordings, and they can be categorized into two groups: spike sorting software
and spike train analysis software as shown in TABLE 1-1.

Spike sorting software is designed to extract spike trains from the
extracellular recordings, and it consists of signal processing tools, such as
filtering, spike detection, and spike sorting. Spike sorting is the most
computationally intense process in neurophysiological signal processing, and
despite significant improvements, spike sorting remains an imperfect process [7].
These packages only extract spike trains and provide very basic spike train
analysis tools. MClust and KlustaKwik are example tools that require detected
spikes as input data [8],[9].

The second category of software is designed to evaluate the statistical
characteristics of spike train data. The input data to these software packages is
spike train data that a user must obtain from extracellular recordings using other
spike sorting software.

TABLE 1-1. EXISTING SOFTWARE PACKAGES

 

Category Spike Sorting Spike Train Analysis

 

MClust, KlustaKwik, Chronux, MatOFF, Spike Train Analysis
Software NeuroMAX, OSort, Wave_clus, Toolkit, FIND, etc
etc

 

 

 

 

 

Very few academic software tools exist to date that integrate all the
needed processing steps in one comprehensive package. As a result, every lab
relies on custom built analysis tools in one form or another. This is coupled with a
significant lack of community-wide standardized set of tools that enables
replicating experimental data analysis across labs to facilitate objective
comparison between scientific findings.

Automated processing is necessary to process large scale datasets, but it
is inevitable to select some parameters manually in some situations such as low
Signal-to-Noise ratios (SNR). Manual selection of large number of parameter is
performed empirically, and this makes manual processing tedious and

inconsistent.

1.3 Thesis Contributions

NeuroQuest is a MATLAB graphic user interface (GUI) software package
developed in this thesis that contains a number of signal processing and analysis
tools exclusively designed for extracellular neural recordings. NeuroQuest is
distinguished from other software packages in several aspects:

First, the package includes algorithms for spike detection and sorting, as
well as spike train analysis. This creates a unified processing environment that

makes a sequence of processes and analyses more efficient and time saving.

Second, spike detection and sorting algorithms implemented in
NeuroQuest outperform those in other software packages such as KlustaKwik,
WAVE_CLUS, and OSort [8],[10],[11]. Robust spike detection and sorting are
crucial because virtually all-subsequent neural data analysis depends on the
outcome of these two steps. These two steps may vary in the way they are
implemented depending on the type of electrode array used and brain area of
recordings, and NeuroQuest offers multiple spike detection and sorting methods
designed for different situations.

Third, the Graphical User Interface (GUI) of NeuroQuest reduces the
complexity of selecting the large number of parameters needed during the
analysis. This significantly reduces learning time for users. NeuroQuest also
provides a variety of graphical tools to help the user manually set the parameters
of choice by instantaneously illustrating the effect of parameter changes on the
analysis results.

Fourth, individual modules are designed as independent GUIs. This
design scheme allows a developer to integrate easily their own processing
modules into NeuroQuest.

Finally, NeuroQuest provides advanced spike train analysis tools that
enable identifying the functional and effective connectivity between
simultaneously recorded cells. This is fundamentally important when studying

cortical plasticity at the population level that may accompany learning and

memory formation - an important design consideration for neural decoding

algorithms in adaptive BMls.

1.4 Organization

The remainder of this thesis is organized as follows. In Chapter 2 the structure
and the organization of NeuroQuest are presented. It further discusses the
contribution of NeuroQuest to overcome the limitations of other software
packages and summary of each process module of NeuroQuest. Finally in
chapter 3, spike detection and spike sorting methods implemented in
NeuroQuest are presented. This unique spike detection method utilizes sparse
representation of the signal for better neural yields and more robust detection in

low SNR cases.

CHAPTER 2

Background

This chapter presents a brief summary of neurophysiological signal processing.
In section 2.1 different types of neurophysiological signal and their characteristics,
applications, and limitations are presented. Section 2.2 presents the complete
procedure of spike train analysis from extracellular recordings, as well as

previous work on spike detection and spike sorting methods.

2.1 Neurophysiological signals and the characteristics

Electrophysiology is the study of the electrical properties of biological cells and
tissues that involves measurements of electric potential and current changes. In
neuroscience, the study focuses on measurements of the electrical activity of
neurons, known as neurophysiological signals [1]. The commonly known types of
neurophysiological signal are Electroencephalogram (EEG), Electrocorticogram
(ECoG), and extracellular recordings, depending on the source of the signal and
the recording sites [2]. EEG, recorded from the scalp, is the summation of the
synchronous activity of large populations of cortical neurons. EEG has been

extensively used in various fields such as neuroscience, cognitive science,

clinical applications, and Brain Machine Interface (BMI) mainly because of the
advantage of being non-invasively recorded. However EEG suffers from poor
spatial and temporal resolutions caused by the filtering effect of the skull and the

artifacts from muscle movements [2].

 

I: ' Scal
T Extracellular *‘A- p
E ! EEG/I Recording \ ECoG ‘
| , ,-
8 | , .ECQG Non-mvaswe Skull
8 I LFP Semi-invasive __.____‘ .. . - _ . ._\___
(D l “c Invasive ‘
E Am?“ 3 Extracellular «- _ -.
=‘ Porem'als Recording “’ "
, ' Cortex
msec Time sec 1

Figure 2-1. Types of neurophysiological signals adapted from [12]

The second type of signal is ECoGs and is directly recorded from the
cortical surface to circumvent the rapid signal attenuation effect of the skull [2].
ECoG has higher spatial resolution, broader bandwidth, and higher amplitude
than EEG, because the electrodes record the neural activity from a smaller brain
area [13]. They are also free of muscle, eye-movement, and other artifacts
consistently present in the scalp EEG of waking, moving patients. ECoG has
been used to localize epileptogenic foci during pre-surgical planning, to map out
cortical functions, and to predict the success of epileptic surgical re-sectioning

[14],[15]. Although ECoG delivers these superior features, they are semi

invasively recorded, requiring a craniotomy, and therefore cannot be collected in
healthy humans [2].

Extracellular recording - the observation of action potentials generated by
a single or ensembles of neuron — can be obtained by directly implanting voltage-
sensing microelectrodes. Extracellular recording typically contains two types of
signals: Local Field Potential (LFP), and Action Potential (AP). LFP is believed to
be the sum of action potentials generated by cells within approximately 50-
350um from the tip of the electrode [16]. LFP is obtained by low-pass filtering the
recordingswith cutoff frequency at ~300Hz.

Action potentials - or spikes — constitutes the main communication method
among neurons, and are key to understanding the complex relationship between
physical features of the world and the brain’s representation and interpretation of
those features [2]. APs from a single or population of neuron have been used in
various studies such as measuring connectivity between different areas of the
brain, examining dynamics of neuronal response and their relationship to
behavior, summarizing experimental data, BMI applications, clinical diagnosis of
the brain disorders such as Parkinson’s disease (PD) and certain types of
epilepsy, the application of Deep Brain Stimulation (DB8), and many other
research areas [3-6]. Despite the significance of studying extracellular
recordings, there are technical challenges in recording and processing the data.
Extracellular recordings require a surgery that causes a risk of complications and

permanent tissue damage to implant the recording electrodes [2]. Since the

recordings are contaminated by artifacts and other noise components, the data
must go through multiple steps of signal processing before the statistical analysis
[7]. Despite the risk of surgical procedure and invasive and unstable recordings,
studies of extracellular recordings seem indispensable to biomedical applications

and basic research.

2.2 Extracting information from neural activity

Extracellular recordings must go through multiple steps of signal processing, Pre-
Conditioning, Spike Detection and Spike Sorting, before examining the statistical

property of the neural responses obtained.

2.2.1 Pre-Conditioning

Pre-Conditioning is a combination of processes to enhance the quality of
the signal to make it suitable for further information extraction. The raw
extracellular recordings are contaminated by many noise components such as 1)
background neural activity attributed to neural sources far from the electrode
array, 2) thermal and electrical noise produced by electrical devices in the data
acquisition system, 3) other artifacts from muscle movement and mechanical

noise [17]. The raw extracellular recording is filtered to remove low frequency

10

components such as LFP to be able to detect spikes. Artifacts that have similar
spectral properties with APs as shown in figure 22 must be removed, because
these artifacts could be incorrectly classified as spikes causing a decrease in
accuracy of further analysis results. If the recordings suffer from low signal to

noise ratio (SNR), various noise reduction techniques can be used to improve the

SNR.
A. a ( i-‘ ., -‘ . u I‘ o. - l ' «H I. " 3°. A 49L 4‘ f“ ' I ‘I . ‘l ‘-, I -.i IV‘ .I '\ p In?" -- A- . ' t) '4 " D r ‘ '1 l ; ‘ ’
'l‘ \' L v' L." "—<.~.'~a‘ V. rd}: / .7 ,1..- r. ‘A ' I I} 7‘ I, Y" ‘J ‘l _ l "f. «V in , , 4‘ .“’ rt 3 I. .‘ v :43 L 'I "1 "".... l w. ‘ y. I. U‘ .( "I _. 4 . .', '
V Spik
O Artifact
, A
. _. A
u, ‘I d "A. I ‘ Ll 'l ~fl
' ' r i I: ’ ' I
5
.. ..>~- '1 l‘. .- 2 ~ I . ‘ I
,A f 1 3|. ‘ I .9“ ', t a. , .I “I L,» I‘“\ u. ‘V n .
I

Figure 2-2. Artifacts and action potentials In the extracellular recording

2.2.2 Spike Detection

Detecting the presence of action potentials from the extracellular recording is
referred to as Spike Detection. It is believed that spike arrival times carry all the

information about information processing and not the actual waveform shape.

11

Therefore extracting the temporal information of spikes is a fundamental step in
the analysis of neural recordings [2]. The main challenges in spike detection are
the noisy environment of neural recordings and the similarity between the spike
waveforms and the background noise. Many spike detection algorithms are
available to overcome these challenges. The simplest technique for spike
detection is amplitude thresholding that searches for an event that crosses user-
specified amplitude thresholds [7]. Although its computational simplicity that
makes easy to implement in hardware, the performance of the method declines
rapidly in low SNR. The instantaneous energy of an extracellular signal has been
used to emphasize the spike peak [18],[19]. This method is more robust to the
noise than the amplitude threshold method, but the energy computation without
any noise reduction technique causes a poor performance in low SNR. In
template matching method, templates, obtained during a training session, are
used to locate possible events in the signal by measuring the resemblance
between the template and the segments of the signal [18]. This method relies on
a priori knowledge of the spike shape to form the template. The performance
again decreases in low SNR, since the automatic selection of a template in a
noisy signal is very difficult. Transformation and decomposition of signals are
also used in spike detection methods. Spike features, observed in the wavelet
domain, have an advantage that separates signals from noise by thresholding
the wavelet coefficients [10],[20],[21]. The wavelet based spike detection

methods suggest that sparsity plays an important role in spike detection by

12

increasing the SNR. The main drawback of this technique is the assumption of a

similarity between a family of wavelet base and a spike shape.

Extracellular Recording Action Potential Extraction Noise Reduction

_, _ ,_ - _, _ __-,_,,__ _, I“
. } [lwv'v.u E

  
  

 

 

 

 

    

 

 

F'Identified Spikes Feature Extraction Spike Extraction Spike Detection
and Spike Sorting and Alignment

T I III I 1
"fl-“ ____- -l_- l l l I ll

Raster Plot of Spike Train

 

 

 

 

Figure 2-3. Processes to obtain spike train from extracellular recording

2.2.3 Spike Sorting

As the tip of the electrode is surrounded by many neurons, the
extracellular recordings can simultaneously pick up spikes of an unknown
number of neurons. The detected spikes can be assigned into a particular
neuron. This process is called Spike Sorting, and this is valid under the

assumption that each neuron generates a distinguished spike shape from others

13

[7]. Spike sorting methods are categorized into two approaches: 1) The pattern
recognition approach: operates on the individual spikes extracted after the Spike
detection step; 2) The Blind Source Separation (BSS) approach: operates on the
raw data without the need to perform spike detection a priori [12]. The Pattern
recognition approaches typically consist of feature extraction from the detected
spike waveforms and the assignment of the spikes to the individual neurons by
clustering the features. Many feature extraction methods were proposed such as
feature of the shape, such as spike height and width or peak-to-peak amplitude,
principal component analysis (PCA), a set of orthogonal basis vectors that
capture the largest variation in the data, and Wavelet Transform [7],[10],[22].
These extracted features are clustered using unsupervised classification
methods such as hierarchical clustering, k-means and fuzzy c-means, Gaussian
mixture model and the t-distribution mixture model, and self organizing map
[7],[23-26]. K-means clustering is relatively sensitive to outliers, and mixture
model based clustering methods usually need assumption of certain distribution
to yield accurate clustering results. The clustering is data-driven process which
means that there is no absolutely superior clustering method than others. Each of
clustering method is designed for the certain type of data, and a proper clustering
method must be selected regarding the nature of the data to yield an optimal
clustering result [27].

In 888 approach, the observation signal is considered to be a mixture of

multiple signals, and the original signals can be estimated from the mixture by

14

finding the unmixing weights [7]. Independent Component Analysis (ICA) and
Multiresolution Analysis of Signal Subspace Invariance Technique (MASSIT) are
examples of the BBS approach [28],[29]. Spike sorting has become crucial in
multiple spike train analysis, because the accuracy of spike sorting influences the

validity of subsequent analysis [6].

2.2.4 Spike Train Analysis

Statistical properties of population of neuron are estimated using the
temporal arrival time information of identified neurons [6]. Action potentials can
be represented a stream of binary events, or called a spike train, where ‘1’
represents the arrival of spike and ‘0’ is not. Spike trains are widely used in the
study of neurophysiology especially neural coding, characterizing the relationship
between the stimulus and the individual or ensemble neuronal responses and the
relationship among electrical activity of the neurons in the ensemble [2],[5],[6].

Spike train analysis is the attempt to find patterns in spike trains that
reflect some aspect of neural functioning. This could include relating neural
activity to stimuli, finding functional interactions among neurons, and estimating

codes distributed across a neuron population [30].

2.2.4.1 Single Neuron Properties

15

One of the simplest ways to study the patterning of spike activity of a
neuron is to construct an lnterspike interval histogram (ISlH). This is simply a plot
of the distribution of the observed times between spikes collected in fixed width
[5]. Peristimulus Time Histograms (PSTH) are histograms of the times at which
neurons fire. These histograms are used to visualize the rate and timing of

neuronal spike discharges in relation to an external stimulus or event [31].

2.2.4.2 Multi Neuron Properties

The cross-correlogram is a function that indicates the firing rate of one
neuron versus another. Cross-correlograms give a measure of the firing rate of
one neuron around the time that another neuron fires which indicate the
dependency of the pair neurons [32].

Joint Peristimulus Time Histogram (JPSTH) gives not only the ability of
cross-correlograms, but also displays the stimulus-related dynamics of the
relationship. The temporal information that is given by the two-dimensional nature
of the JPSTH is important in studying neural connectivity [33].

Advanced analysis tools consist of identifying relationships between the
observed neurons from spike train ensembles. NeuroQuest offers two algorithms
to achieve that goal: multiscale clustering and Dynamic Bayesian Network (DBN).
The first algorithm identifies any potential statistical dependency between spike

trains, often referred to as functional connectivity [34]. The second algorithm

16

infers the type of connection (excitatory/inhibitory) and directions between the
functionally-interdependent neurons, often referred to as effective connectivity
[35]. This is achieved through DBN. This two-stage framework can efficiently
identify neural circuits in large neuronal populations. It can also be utilized to

track plastic changes associated with learning and memory.

17

CHAPTER 3

Overview of NeuroQuest

This chapter provides context for the contributions of the thesis. First, the
structure of NeuroQuest and the organization of individual processing modules
are presented in Section 3.1. Second, Section 3.2 discusses the contributions of
NeuroQuest that overcome the limitations of current existing software packages.

Finally, summary of individual process modules are presented in Section 3.3.

3. 1 Structure and Organization of NeuroQuest

NeuroQuest is designed to process the large-scale raw extracellular recordings
to obtain dynamics of neural population with statistical tools. NeuroQuest bundles
multiple processing modules designed as Graphical User Interface (GUI), and
they are connected to each other through the main workspace as shown in figure
3-1. Total 8 processing modules provided in the software are classified into two
groups: spike sorting group and spike train analysis group. Spike sorting modules
require the raw or pre-processed extracellular recordings, while spike analysis
tools handle single or multiple spike trains. Once the input data is loaded, the
corresponding group of modules becomes available on the main menu. Each

module contains sub-modules that assist to yield more accurate analysis results.

18

Figure 3-1 illustrates the organization of the individual processing modules and
their sub-modules of NeuroQuest. Details of each process are discussed in
section 3.3.

The flowchart in Figure 3-2 demonstrates the complete neural data
processing and analysis steps in NeuroQuest. After the extracellular recording
data is loaded, the group of spike sorting tools is activated for further process.
The first stage is to denoise the data to enhance the neural yield. After denoising,
spikes are detected from the recording. Detected spikes are extracted from the
data and sent to the spike sorting algorithm to obtain spike trains. The obtained
spike trains are then further analyzed using the primary spike train analysis tools
such as lnterspike Interval Histogram(lSlH), Peristimulus Time Histogram(PSTH),
Joint Peristimulus Time Histogram(JPSTH), Cross-Correlogram (CC), and the
advanced spike train analysis tools such as functional and effective connectivity
estimation among observed neural population. Since individual modules are
designed as an independent GUI, the analysis results of each module can be

saved and loaded separately for the later data analysis.

3.2 The contributions of NeuroQuest

The main contribution of NeuroQuest is to provide a tool that assists a researcher

to translate neural dynamics into statistical measurements. Many

19

neurophysiological signal processing software packages were developed for the

purpose , but they did not fulfill the needs with their limited capabilities to deal

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

       
   
 

 

 

 

 
  

Spike Sorting fl
_ — O) .—
m +- tu *5 c to 8 *" 8 713 75
o. > 3 > '55 2‘2) " '= 9 =- 9 C E C
5% “at: as $28. at :9: a:
a: < t: 3 8 ,9 8 3 8 m o 0
Pre- Spike
Conditioning Detection Sorting

 

t

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

L

Figure 3-1 Organization of NeuroQuest

20

 

 

 

 

\\ i
' LI \v’
L . Neural
[ LMaIn GUI W Data
, r 7\\"‘
/’7’7 ‘1‘ §
.2/ J \.
l‘y/f‘ V . ‘3' '4\
\‘v 7 .
Single Unit Multi Unit Advanced
Analysis Analysis Analysis
E — a 2'
E I m E ‘1’ E
g E g 8’ l— .5 ‘3 2% ‘6
(D (D h B (D 8 C a) QC)
- 0- 0 t % c c 3: c
0 LE 0 UJ o
O 0 0
Spiketrain Analysis

 

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Multi Unit
Data
save the
processed
data
save the
processed .
data Detection
Spike
Sorting Single
{k Unit
Data
’ Spike Train
Analysis JL
Single Spike
Train Analysis
(ISIH, PSTH)
V
Clustering Circuit Inference
(Functional M (Effectivel
connectivity) connectivity)
l
l
Spike Train ¢

 

 

Analysis Result

 

 

 

Figure 3-2 Flow Chart of NeuroQuest

21

with the experimental data. NeuroQuest is designed to overcome the limitations
of the existing software packages, and following aspects make NeuroQuest
distinguished from the others: 1) comprehensive toolbox, 2) simple GUI and
visualization tools, 3) capability of handling large-scale neural data, 4) varieties of

spike detection and sorting methods, and 5) developer friendly environment.

3.2.1 Comprehensive toolbox

As mentioned in chapter 1, most of the current existing software packages fall
into two groups, spike sorting software and spike train analysis tools. Since both
groups are required to perform a complete neural data analysis, a combination of
multiple software packages is used to analyze experimental data. However this
scheme has several drawbacks. First of all, input data must be converted into
suitable format for the different software packages, since each software package
requires different input data format. It is also time consuming and inefficient
process due to the stream of non-unified processes for the single analysis.
TABLE 3.1 illustrates the provided features of several software packages
available. Among spike sorting software, only NeuroMAX and OSort offer limited
spike train analysis tools [11]. Among spike train analysis tools FIND is the only
software equipped spike detection and simple spike sorting tools [36].
NeuroQuest offers comprehensive tools from pre-processing of

extracellular recordings to analyzing multiple spike train to achieve seamless

22

TABLE 3-1. FEATURES OF THE SOFTWARE PACKAGES

      
 

k
0,08

0x
w;

r°°llrit

Ve

ex
[”81
'04,“

x K K X K >< \ FIND

x K K x x x K ””057:

x K K x x x K 87-4

x x x K K x K We

x x K K K K K 0.30”

K K K K K K K Mamba”
98:

”ea

GUI

 

Pre-Conditioning

 

Spike Detection

 

Spike Sorting

 

Single Unit
Analysis

 

Multi Unit Analysis

 

Connectivity
Estimation

X X X ‘ X X ‘ Chro"
x x x K x x K lat,”
x x x K x x K ”C

><><\\\><K

 

 

 

 

 

 

 

 

 

 

 

 

 

neural data analysis. This seamless and unified process eliminates the burden of

converting data and speed up the analysis.

3.2.2 Simple GUI and Visualization Tools

A large number of parameters selections are required in each step of the analysis
and the complexity of the parameter selection increases the learning time of the
software for the users. The simple GUI of NeuroQuest reduces the complexity of

the parameter setting and shortens the learning time. Various graphical tools help

23

the manual parameter selection by instantaneously illustrating the effect of the
parameter changes on the analysis results. The extensive visualization of the
manual process helps to improve the accuracy of the analysis results. We
demonstrate the effectiveness of the GUI by comparing the spike detection result
against OSort, a spike sorting software package, using the same energy-based
spike detection method with three simulation data sets obtained from the OSort

package [1 1].

Detection Rate is defined as a total number of positive detection over a total
number of the actual spikes. As Figure 3-1 shows, NeuroQuest yields more
reliable positive spike detection results with various noise levels while
maintaining the same false detection rate. Details of the spike detection
comparison between NeuroQuest and OSort results are found in TABLE 3-2. The
improved detection performance is due to the threshold selection tool in the spike

detection GUI that assists a user to select a proper threshold.

3.2.3 Capability of handling large-scale neural data

Complex behavioral related experiments require hours of data recordings, and
processing these data become another challenge in neuroscience and BMI
development [7]. Most of existing spike sorting tools are missing the capability of

handling large-scale neural data mainly due to the MATLAB environment

24

Sim3
100 f2 _;:._ _'_', '7 '— —"
‘1 NeuroQuest ,
I; 080"

     

80

 

Detection Rate

 

 

 

 

Figure 3-3. Comparison of Spike Detection performance between
NeuroQuest and OSort

TABLE 3-2 SPIKE DETECTION RESULTS OF NEUROQUEST AND OSORT
WITH FOUR SIMULATION DATASETS

 

 

 

 

 

 

SNR 4 3 2 1

TP:1535 TP:1523 TP:1463 TP:1316
A (97.40%) (96.57 %) (92.83%) (83.50%)

E E FP: O FP: O FP: 13 FP: 75

t5 8 TP:1513 TP:1503 TP:1407 1094
(96%) (95.37%) (89.28%) (69.42%)
FP: O FP: 2 FP: 46 FP: 163

TP:1544 TP:1420 TP:967 TP:857
(98.47%) (90.56%) 61 67%) (57.27%)
g {-35 FP: 0 FP: 39 FP: 85 FP: 225
0-, e TP:1460 TP:1251 TP:762 TP:312
V (93.11%) (79.78%) (48.60%) (19.90%)
FP: 0 FP: 14 FP: 202 FP: 134
T132870 T192278 TP:14OO TP:1378
A (96.12%) (76.29%) (46.89%) (46.15%)
(E) g FP: 6 FP: 60 FP: 321 FP: 348
“(75 g; TP:2531 TP:2145 TP:1193 TP:820
(84.76%) (71 84%) (39.95%) (27.46%)
FP:7 FP:317 FP:27O FP:308

 

The shaded area indicates the results of NeuroQuest and the white background

is the result of OSort.

25

 

,known for poor handling of large dataset, and this is a significant drawback in
the process of the experimental data. Although NeuroQuest also runs on
MATLAB environment, we circumvent the memory issue by allocating large data
set into small segments. This segmentation helps to speed up the process and
resolves the ‘out of memory’ problem. The segmented process also helps to
improve the accuracy of the analysis results by estimating the important

parameters for the process locally.
3.2.4 Spike Detection and Sorting methods

Some characteristics of extracellular recordings, caused by the different
geometry of the recording electrode array, require specially designed processing
technique. Most of the software packages only offer single spike sorting method
that is not suitable for processing the different types of data. NeuroQuest is the
only software that offers array detection and sorting algorithms that are specially
designed to analyze the data recorded using the closely spaced electrode array

such as stereotrodes, tetrodes and polytrodes.
3.2.5 Developer Friendly Environment

Very few software tools exist to date that integrate all the needed processing

steps in one comprehensive package. As a result, every lab relies on custom

26

built analysis tools in one form or another. This is coupled with a significant lack
of community—wide standardized set of tools that enables replicating
experimental data analysis across labs to facilitate objective comparison between
scientific findings. In NeuroQuest individual modules are designed as an
independent GUI, and this design scheme allows a developer to integrate their

own processing module into NeuroQuest easily for the performance comparison.

3.3 Summary of individual modules

3.3.1 Data format

Since most of the commercial extra-cellular recording systems store the
recordings with their own data structure, proper data conversion is the inevitable
step to use the academic software packages. Without a standard data format,
data conversion is one of the difficulties in the development of neurophysiological
signal processing tools. There have been a consensus on unifying data structure
of neurophysiological recordings, and NeuroShare (http://neuroshare.orgl) is
one of the many efforts. This data structure is designed to efficiently store
neurophysiological signals, and several commercial and academic software
packages are supporting this format. However the NeuroShare data converter

requires mediate levels of understanding in computer language C++.

27

To make usage of the software more universal, a generic MATLAB data
file is used for the input data of NeuroQuest. NeuroQuest works with a MAT file
that contains a specific data structure (see Appendix A) that can be created
easily with a beginner level of understandings in MATLAB script language. Once
the following components are saved in a MAT file, NeuroQuest automatically

enables the group of modules that can be applied to the data.

3.3.2 Pre-Conditioning

Pre-Conditioning tools are designed to improve the quality of data and preserve
the information at the different frequency range other than APs (frequency range
of 300Hz - 5000Hz) such as LFP. Properly pre-conditioned data improve the SNR
and fewer artifacts that increase the accuracy of further processes such as spike
detection and spike sorting. Three pre-conditioning tools in NeuroQuest, LFP

Removal, Artifact Removal, and Wavelet Denoising are presented in this section.

28

DECO n 51.4.. .1, A L - .‘L _.. .L- :23 .-.r - 0.0015109_GUI A -- 4; E

 

 

 

 

    

 

 

 

I‘m PMWM‘ -‘ — I) 1 r T 1 r r T r T
[ I TFP'Rémoval J I -
, _ -1 3, .. .
I r: Artifact Removal 1 ' ‘5 fI ~ , ‘ , ' ‘ I |
~— V . . . , . I I I I | . ' I ‘
]_ LSpike Data Denoisl ng J I 1.1 ' '
T I I l l
l l l I I | l
I Advanced Settings '1 I I I I I I I! I I
I_.g. "1:_--__‘.-__I_Ag--'_. J It I I ‘ I
F8' .~. ~ 1' 4 . 1 .
,.’ § ’1' .2 . 1 I |
G 3
I, ’5 e L I I I l I
, b . I— ~' I I I I I
a: 2 f
8 , ‘ ' a zvl
{ . 3.; I j . '
, .. 1 ' ‘ 'i , Scale 0511751117 . l 1 l
'1. g . r . . '1 .
.‘9. j " .. I 1" l I ’ II
. - t , "i I 1
,. ”I: a I * I I 1 I I I i
:IPN. I .. .5 - _- A 4— .}: ‘
Tot-Imam 1938 s'oc . 1 T 1 1 . L l I ~ I , ,
w ,. - 0.1 ‘ 0.2 0.3 0.4 0.s-~ . 0.0 0.7. ,.. 0.0-. 0.0
I v . O f 1 . w. 5' )“I MI." I . 05—;
:va .FI—A‘ Save A I 3 I Ii ‘ I 6% . ' Jul ‘1’ -' N ' ‘ .
l .1 T” ‘ '7 _ .L I "u ’. . . . . I II '3‘” "-‘-’ I 3."
J I 'J' I“ I .3 3 .i , n“ .4?“ I I t

Figure 3-4 Pre-Condltloning GUI

3.3.2.1 LFP Removal

LFP is low frequency oscillation, range of OHz - 300Hz, that measures the activity
in a population of neurons. LFP is robust over time comparing with AP, while it.
provides highly specific information [2]. Conventionally LFP was considered as a
noise component and filtered out from the recordings [7]. However current
studies show that LFP has excellent decoding properties for Brain Machine
Interface and other studies [2]. NeuroQuest extracts LFP data from the unfiltered

extra-cellular recording and displays it using a spectrogram [37]. 5th order

29

Butterworth filter is implemented in NeuroQuest for LFP extraction, and a user

can specify upper and lower band limits of the filter [38].

3.3.2.2 Artifact Removal

One of the common noise components in extra-cellular recording is the artifacts
from different sources such as electromyogram (EMG) from muscles in the scalp,
jaws, neck, and body, and many types of electrical artifacts as the animal moves
[39]. These artifacts are highly correlated through neighboring electrodes, and
the artifacts have similar spectral content to the desired spikes. This similarity
causes cluster overlap in feature space that make difficult to define clean
boundaries for spike sorting. These artifacts can be identified using principle
components analysis (PCA). First two or three PCs, the best representation of
the data, reflect the correlated noise components and the artifacts can be
removed from the recordings by eliminating these dominant PCs. Details of the

method are found in [39].

3.3.2.3 Wavelet Denoising

Multi-electrode arrays are intended to simultaneously monitor large number of
neurons and therefore detection of as many neurons as possible from the

recorded signals is of fundamental importance. Effective denoising improves the

30

results of further signal processing steps, such as spike detection and spike
sorting [9].

Spatially-correlated noise across channels is one of the fundamental
problems faced during spike detection. It arises due to the synchronized firing of
large cell assemblies that are not close enough to the recording array to be
individually detected. Suppressing this noise component is necessary for optimal
spike waveform detection. Conventional linear filtering methods, such as Fourier
transform-based filters, are not suitable to discriminate this noise component
from the signal of interest [9]. The Discrete Wavelet Transform (DWT) is a
nonlinear transform that possesses many desirable properties not present in
linear filtering methods. The DWT is a powerful signal processing tool that has
been demonstrated to have excellent properties in many applications of
denoising and compression with relatively simple mathematical computations
[40],[41],[17],[20],[10].

For example, in the DWT domain, the signal is concentrated in very few
coefficients with large amplitude while small amplitude noise coefficients are
widely spread. This characteristic of DWT permits to suppress the noise by
thresholding small coefficients. The performance of DWT denoising depends on
the choice of the threshold. Six different thresholding methods are provided in
NeuroQuest: Heursure, Rigrsure, VisuShrink, SureShrink, BayesShrink, and

Minimaxi. Details of these methods are found in [40],[42],[43].

31

3. 3.3 Spike Detection

The compactness property of the wavelet transform facilitates the detection of
neural spikes in the wavelet domain. NeuroQuest performs multi-level stationary
wavelet packet decomposition (SWT) of the recorded data. It uses a likelihood-
ratio test (LRT) for detection, in which a sufficient statistic is computed depending
on the array geometry [44],[45].

If the array is closely spaced, the statistic is computed from a snapshot of
the entire array to minimize the effect of the spatially correlated noise component.
If the array is not closely spaced, then the sufficient statistic is computed from
snapshots of individual channels and therefore improves spike detection in both
single and multi channel scenarios [41]. Not only wavelet detection, NeuroQuest
also provides three time domain detection methods: single amplitude, absolute
amplitude, and energy-based. Single amplitude detection method identifies any
signal that crosses the threshold as a spike, while absolute amplitude detection
applies both positive and negative thresholds simultaneously [7]. Energy-based
Spike detectors compare the local power of the signal with a threshold estimated
from the power of noise [19]. This method is more robust to the noise than the
amplitude threshold methods.

Several detection tools help to maximize detection performance. Manual
detection allows a user to select a threshold value from the selected data

segments. Different data segments, the lowest SNR, the smallest noise variance,

32

and a large data segment (30 sec of data segment is selected empirically) can be
selected for the threshold value estimation. In wavelet detection, choices of
different wavelet bases help to maximize the compactness of wavelet

coefficients.

3. 3.4 Spike Sorting

Spike sorting is a step where spikes are assigned to the individual neurons that
have emitted them [7]. Overall spike sorting algorithm in NeuroQuest is described

in figure 3-6.

 

 
 
 
 
 
 
 
   
   

 

.‘ * .4. 5 . 3
. Wavelet Detection
it.
Time Detection
,III‘BRGAJISHIIM—T ' '5
”4.2.... I.’—-.. I
MW - '.
(h .I‘fi] " . _I II
If- . .‘J. - ‘
.fim‘ -"-
. “'1
”It Symon :I '
l 2 : I

 

 

 

 

 

.~ I I I ' I I I l . I I I I I I I
.. l I I I I I II I . I I I I I I I I II I I I I I
III‘, I I . I ' ‘I CI l I I I ' I. I

OI
.. II .
'I '. :Ir

I

I'II

l I ,
I.
I. . I

..-.-'IIIi? I’I
.3 '

 

 

 

- . --°-' , a , ..
. . v L.__“"'Ic""""..-.._ fimwnmm

3.‘..... .. ......- -..._ --.;_ :4;

Figure 3-5. Spike Detection GUI

33

 

 

 

 

I

Multi Channel Single Channel
Spike Sorting Spike Sorting

 

 

 

 

 

 

 

 

 

 

 

 

 

 

'— MASSIT

 

 

 

 

Wavelet __
Domain PCA

 

 

 

 

 

C #3:: )

Figure 3-6. Flowchart of Spike Sorting Algorithm in NeuroQuest

First, the specific type of the recording electrode must be specified to choose the
proper spike sorting method. Geometrical differences of the recording electrode
array require different spike sorting approaches. lf spacing between neighboring
electrodes is close such as stereotrodes, and tetrodes, there is high chance to
record the same action potential with a group of electrodes. Recordings from
these electrodes allow additional information to be used for more accurate spike

sorting [15]. Two types of spike sorting techniques, single and multi channel

34

sorting, are embedded in NeuroQuest. Single channel sorting method does not
consider any correlation between the recordings from the adjacent electrodes
and processes individual signal independently, while multi channel sorting
captures the correlation among multiple electrodes and use the information to
identify the source of the spikes.

Second, if two adjacent events are close enough, the event is considered
to be an overlap. In the case, Multiresolution Analysis of Signal Subspace
Invariance Technique (MASSIT) resolves the spike overlap based on an
augmented representation of the observation space to simultaneously
incorporate the spectral, temporal and spatial information of the spike waveform
[44].

The algorithm verifies whether the detected events consist of a simple spike or
overlap of two or more spikes. In the non-overlapping case, wavelet features of
the spike waveform are compared to previously extracted features during
training. If matching occurs, the spike is classified to belong to the matching
class. If not, a new class is created. More details are provided in [44]. Spike
sorting in NeuroQuest consists with three steps: spike extraction, feature

extraction and clustering. Details of each step are discussed below.

3.3.4.1 Spike Extraction and Alignment

35

Detected spikes are extracted and aligned in this step. Since duration of spikes
may vary, proper spike extraction is crucial to preserve the information for spike

sorting. Spike alignment is also important factor in accurate classification [7].

 

 

       
   

 

 

 

“A LA- A J A .3... A 5 MA A __ 4
mm 22. _ - u . . ,9.“ M: .. .r ;.2
' 09* O” "1 W __3___,__ ‘ I J
...t-_--. -., " m ,
um - 1 - V -. 2' c . I
t.-?°-té'.'::‘°:, of! -' ' . ‘ ~=""~:....=':'~”’
‘ Featur- Extraction J' '2 I 05 I mm
Gull-.Cuttng I? l
“.4. -0 4" .---l A * o I . I
1h stoowomq I ‘ I If _ ._,
- l wags .
. ' ' a Q I
. 4 ‘ g?‘ , ‘_
3 .
I I"? I
3i
t- .-_J
,v
~ ,-..~
I 'rOunuErm I

I

11". :ML—IS—L—J “

L.;‘__ _. ' l

 

   

(a) Negative Peak Alignment (b) Positive Peak Alignment (c) Energy-Based Alignment

Figure 3-8. Clusters of detected spikes in temporal PCA feature space with
different alignment methods.

36

3.3.4.2 Feature Extraction

Next step is feature extraction from the properly extracted and aligned spikes.
Extracting a feature set that represents a group of spikes emitted from the same
neurons is the key to maximizing spike sorting accuracy. Typical feature sets for
spike sorting are features of the shape, such as spike height or width, and
dominant PCs of a spike. Recent studies show that DWT is a powerful feature
extraction method, and with properly selected wavelet base DWT outperforms
the PCA—based feature extraction [46]. NeuroQuest provides a flexible feature
extraction method that allows a user to visualize the extracted feature set in the
feature space instantaneously; different combinations of features or wavelet
bases can be adjusted to achieve the best cluster separation in the selected
feature space. Figure 3-9 shows the clusters of the same set of extracted spikes
with the different feature sets. Wavelet footprint is a wavelet feature set that
represents the transient of the spike by extract the largest wavelet coefficients
through the nodes within a certain range or called the cone of influence [47].
Details of the method are discussed in chapter 4. As mentioned in chapter 2,
properly selected wavelet base DWI“, in this example Symlet wavelet base with
scale 4, outperforms the PCA-based feature extraction by clearly forming 4

clusters while only three clusters appear in the temporal PCA domain.

37

   

a) Wavelet footprint with a) Wavelet footprint with a) Temporal PCA

sym4 wavelet base db4 wavelet base

Figure 3-9. Clusters of extracted spikes in different feature spaces.

3.3.4.3 Clustering

Clustering is a method for finding clusters in multidimensional data sets and
classifying data based on those clusters. There are many methods for clustering
and the best clustering method is driven by the nature of the data [27]. Estimating
the number of cluster is one of the difficulties in clustering, and the automated
class selection can yield a reasonable estimation only in certain condition such
as sufficient cluster separation and accuracy of assumption of the distribution
model [27]. Among many class estimation methods NeuroQuest employs the
modified Expectation-Maximization (EM) algorithm for the class estimation
method. This specific algorithm can be applied to any type of parametric mixture
model for which it is possible to write an EM algorithm, while most of the literature
on finite mixtures focuses on Gaussian mixtures [48]. Clustering of the extracted
features is achieved through four clustering methods: Fuzzy c-means, EM, k-

mean and manual cluster cutting [27],[25].

38

3.3.4.4 Sub-Clustering

Once the spike sorting is performed, initially classified classes can be further
sorted into smaller classes using sub-clustering algorithm. Sub-clustering
algorithm allows selecting one of the classified classes, and projects them into
the different angle of the original feature space or a different feature space. If
there is any improvement in cluster separation in the new feature space, new
clusters are classified using the different clustering methods. Details of sub-

clustering are illustrated in figure 3-10.

3.3.4.5 Extra spike sorting tools

NeuroQuest also provides two graphical tools to evaluate and enhance the spike
sorting result. lnterspike Interval Histogram (ISIH) is used to validate the sorting
result, since a significant number of interspike interval within the refractory period
indicates an error of the sorting results. Another assistant tool is class merging

GUI that allows merging multiple classified classes into single class.

39

 

Cluster Extraction Feature Extraction Clustering

Figure 3-10. Scheme of Sub-Clustering.

Cluster Merging

   

Initial Clustering Final Clustering
Result : 10 clusters Result : 5 clusters

Figure3-11. Example of Cluster Merging Tool

40

Merging tool helps to correct the clustering error without re-clustering entire data.
Merging tool can be also used as a manual clustering tool by initially clustering
data into many classes and combining them. In figure 3-11 initially clustered 10

classes using Fuzzy-c mean are merged into 5 clusters.

3. 3.5 Spike train analysis

Since-it is believed that the time of arrival of the spikes carries all the information,
the Spike waveform is abstracted into a stream of binary events where an
isolated '1' represents an action potential [6]. The binary waveform is referred to
as a spike train. The patterns of spike activity are influenced by three things: a)
the intrinsic properties of the neuron, especially the properties of its membrane,
b) network interactions, because spike activity in one neuron might have
feedback effects on that neuron because of the changes that it produces in
reciprocally connected neurons and c) the nature of the inputs to that neuron [5].
Spike train analysis is the attempt to find patterns in spike trains. NeuroQuest
provides a number of spike train analysis tools categorized into three groups:
single unit analysis (SUA), multi unit analysis (MUA), and advanced analysis

tools.

3.3.5.1 Single Unit Analysis tools.

41

One of the simplest ways to study the patterning of spike activity in a neuron is to
construct an lSIH. This is a plot of the distribution of the observed times between
spikes collected in ‘bins’ of fixed width.

PSTH, a histogram of the times at which neurons fire, is used to visualize

the rate and timing of neuronal spike discharges in relation to an external

(k)
stimulus or event. Let "t (u) be the spike count of neuron iin the uth bin of the

kth trial. Then, the average spike count,

K
1
<n§k><u>> = E Znikkm
19:1 (3.1)

represents the probability of the occurrence of the spike event of neuron iat bin

(’9)
u, where K being a total number of trials. The histogram of (”i (71)) over index

u represents the PSTH [49].

SUA GUI displays the raster plot of multiple spike trains and the stimulus,
if presented. Number of parameter selections for ISIH and PSTH such as Bin
Size and pre and post stimulus duration for the PSTH display, are available in
SUA GUI. NeuroQuest also provides a curve fitting tool for the ISIH that validates
the quality of spike trains. ISTH of a typical neuron demonstrates a single
negative exponential distribution with the empty first few bins of the histogram
due to the refractory period, and the curve fitting tool verifies the quality of the

spike train by fitting the distribution of the ISIH with exponential distribution [6],[5].

42

 

 

 

 

 

m- ~ 5—- .. s
'5 m": ' ” «U
l . U’ - ' ..r n
5 "l I:

 

 

 

    

 

 

.43- .- a
9“ " 'r ‘ " 3
' am I 3
'0 mug {If - ‘ q
mm I, ‘ .
fl" II2 in“, ‘ '. 3‘ n:
J» 9'1"“? A " li~
9"“ g " ‘E .. ' :3 . _ . .‘ ' _ , ' up, 0.4 -
‘ .a.3.7 ".p‘ F n ..T l h ‘ l 'V k' . I. ..." t ”0., I 4; \ . .. .
‘3 ... l “I ‘I IYI‘I rtuj ""1 Y 1 I41" 1 I IIIIII IIII Ill Jinn
3’ I - 3", lli Ill n'fiIII'IHIII'I'ii'II”)l "Iii ' I "I I"
S p P" I III“ I r III II Ilnll| l."l ' ' ”6' I' It. 'I. In.“ 'I II n. 'l iLI .
' I 3 II i I III I I I ' I II
- ; 5.; Y: n . a»? I":h' H | III) II '4 II I il [I I'll hlul'il |III1”l|'IIII:’I."I f’l lqll I" .4
..I ' ’ V . I ll ii I II I II I II I I I I
1'” fl .. l *2: "' 1” ”ltl'ti i” ”ill ill;l l’llm'l III! I II I ”I | u Iil'b' I'l'il I II" If;
a It: I: .‘ - 1’ . ... ”v ‘.. "J .... I“ 'Lll'r‘l'.":' "".1.l'."l"". :Ih. 'I... l'l'l'” minim.
"I. L .‘ '1 " II llllll II ‘I H“ II ill" I I I""' II I II I
i - u ‘ “' '_ ”I“ a; 5 III I III (”I’ll ‘H I it‘l I Ifi Iii I".I 4" 4"
fi: ‘ ., ‘u ' A.’ ”W: I '“ PIII J I I II III I It“ I H I IIIIIi I. I” III" I :"l'yl‘ I‘InlulilllI 1"“ I. “ II | I. . '4
"1" _ .~ 3': . " a " 3 I ' . II it! ' I ”if“ ' "PIMP. I II I If: I '0 IIIII I IL 'II. II II 'l'l II

 

 

 

. l
.-- .-f --'_ -."..' I ‘ 'I'I "ll' .. ""‘l
1'},- .{ - ,- - 3L 1 is lmll ill! ”H H I I III lIIHl I "S'HI' ' '1' I ... II I I ‘4
i, - “ .,' g I h
4" o ‘ , 0 ,~ .‘IT— 1 .
. , . I _ . ”w TT'T

 

 

 

 

' mam~ 8 " — In. .- F f . _
gr. , . IIIII ’ ‘ _ u ,.
”47923...- » :4 * III” ~’~ _
. ‘ MN" . i °
I .g]" 7’ . -_ (- q

. ‘~.. . 3‘ .. 1‘ 9" t:
(‘23. ..3...‘ ... " * -5. .5 Q J" " ‘
' L “1' J( 3'" Ji“" 1

with-mu? i aw i -

 

 

  
 

 

-.- ‘ -1 1 /.-
-Plfl*"' “ ’ 2‘
- “I 3 :4. -
1 h
mm * . -
I '~ -<
.1 ' i 5.x
A“, ' -‘ V ' l
‘J- M I ‘ g '7‘ .5 ; mgr -r.sQ‘I.-» In In In «II III 'III an
...... . .. mfim

 

Figure 3-12 Single Unit Analysis GUI

3.3.5.2 Multi Unit Analysis tools.

The cross-correlogram (CC) is a function which indicates the firing rate of the

target neuron versus the reference neuron. CC give some measure of the firing

43

 

rate or firing probability of the target neuron around the time that the reference
neuron fires. Therefore, the CC provides some indication of the dependence of

the two neurons [50]. The cross-correlation function can be defined as
Qi,j(7') = Et[(32‘(n)3j(n + 7)] (3.2)

where EtH denotes expected value over time n, and 330 is a sum of Dirac
delta functions of neuron j at the time of firing events [51].

The JPSTH provides not only the ability of the CC, but also displays the
stimulus-related dynamics of the relationship. The temporal information which is
given by the two-dimensional nature of the JPSTH can be quite important in
studying neural connectivity [50]. In the computation of the JPSTH, a square
matrix with a size of the trial duration T is prepared. The two time axes are locked
to the stimulus onset at the left-bottom corner. The square matrix is divided into
the bins, each of which is specified by a pair of two integer indices (u, v) [49]. At
each compartment, we assign the matrix element that quantifies the probability of
the occurrence of a joint event, where neuron 1 has a spike event at bin u and

neuron 2 has one at bin v,

K
(n12( u, ’0» M:
K1k= (3.3)
where mm) (u v)- — ni’“ (207220002) (3.4)

Similar to SUA GUI, MUA GUI also provides the number of parameter selections
for both JPSTH and CC such as Bin Size, Window Size for CC, and pre and post

stimulus duration for the JPSTH display.

44

 

 

CCJPSTHI GUI

 

fl

I Y

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. - ~ . ‘4; N.. ”“001 'l f‘
I 3_ '" WW” II I lYlnl I”: mII ” IHIIIHH II IIIII II I llll ITT'HH I I ' u 1' ' "I ”I'm Illl’fI‘q
,IJQ II i I' |I IImIL IIIII ‘ llIII II|"IIIIII ”III '"III u IIHII'I lAII Illl II I III I I II II I
— ”Ii II II I I III l IllI'I III ' ' "I ”h II'I' II I I I III
30+.“ iI II l IIIIII I I, I II ll ll l iii—I
' IIIIII III ' I I I II |lI ll" II"I II III " ' I'I‘\III'I|'I 'l'l '
I #IIII I I II I III! II II I l I '1 ”II I ll I II |
. 25 ’- I 'llllIH.I ”I III I” l I II III J 1|; II:I lll“ 'FI " lI II lll‘l| III“ I I I I'I I Ill" 1' I'll H H L‘
.. "Illl.”“ll| I" II III I III Ill I I II ’IIIII‘IIIIIIQ I II I” I‘I’ I'll H“ " I Illl I
gm‘h- II [III H‘lll I I I IIII I ‘IHIH :Iull II III I I II' I II I) 'I II I I I II II—
4, ... I lIl|lll 'I ' II“ III ' ”It. IHII II llII'I'I i i I I H
1ST I". Ill ‘ lI’ll l I l I l I II lullI "I I! 'Ilfll II “an”. :I'IIH‘ I'll‘III: lIll ll I'IVI' I I i I ' C]
I I I II I II lIlll I llI I I I I I III I H | ' Ill .|' IIH' IIIII I lIIIIII I. ll $I.' "”ll 1"
‘10: III III I III Illl II” 'I Il IIIII IIlIll II IIIIIII I|l| IIIII IIII'IIIIIII IIIIIII IIIII II IIIIIII lL
Il I I”H'll‘lll II I I II III I I | l I III [II II I I lllfl" I ' II .I"' III “III'I IIIII I II III‘III
II I l I II II IIIIIIII I I I II I
”- g_.t I I Illllll I I llll I II I
42‘ I I L l l L J I
Iii-sl- on 012 _013 4 I15 m In as l
~o 0.1 £10.: 3» 0.3 can WI. , :03 0.0 III I
efimdwslt * - ~ — ‘- 0.8 H
I I. I
Caro... @JPSTH . I I i f
P“ - -. T p 'h‘ 7 '107_
. Neuront Nouw i ' . - ll -_
.2 e] 2 g I : -. .. _ . ":—
. I 0A, ' .:. : -' '_
“' s ' - -I q " I "
R50 . 1.. y - 5-! n
0.5 > - '11.
II : ,, ‘- I- 110
I" . 1 0.4, I I. - I
,r - 1 - - -.l ' .
Panel ' . 18
aInquI : 0.3~
‘0 muc I
,sIuIIng'PohI - - I 0‘2
106 I
m mm ‘1 0.1
“onto. 20 '
786 "'_'°° ' {I Save Plot 520"

 

 

. 0.2

0.3 0.4

Matron 2

0.5

 

Figure 3-13 Multi Unit Analysis GUI

45

 

3.3.5.3 Advanced Analysis tools

Advanced analysis tools consist of identifying relationships between the
observed neurons from spike train ensembles. NeuroQuest offers two algorithms
to achieve that goal: multiscale clustering and Dynamic Bayesian Network (DBN).
The first algorithm identifies any potential statistical dependency between spike
trains, often referred to as functional connectivity [34]. The second algorithm
infers the type of connection (excitatory/inhibitory) and directions between the
functionally-interdependent neurons, often referred to as effective connectivity
[35]. This is achieved through DBN. This two-stage framework can efficiently
identify neural circuits in large neuronal populations. It can also be utilized to

track plastic changes associated with learning and memory.

46

 

 

.ClustertngResult— —_- -77 ”we ,____ __

 

Best Temporal Depth: 3.072 sec

Dunn Index forTrIal 1 Clustering Result for Trial 1

&

 

L

3 4

 

M (.0
o
I

K 1
0‘ J a"
-o’ I l Save Plot

513151715 (90
~ ‘9

 

 

Dunn index
\Q
Q

Q CD a

 

Temporal Depth (5) ‘9 I
Augmented Correlation Matirx Trial 1 <1 1 ®

 

 

 

 

19 18
c:
g I
a) I
Z i

5 10 15 20 I
Neuron

 

(a) Functional connectivity Estimation GUI

 

 

 

Network Inference Result—-w-.- -————~w——~-~-—r ___-_--._____. -—~-n --____,, _._,_-._ _, -——~- ~

 

 

Connectivity Matrix for Trial 1

Graph Style Ineeto ”v

2 4 8 8 10
From

Save Plot

 

 

 

 

 

 

(b) Effectiveness Connectivity Estimation GUI

 

Figure 3-14 Advanced Spike Train Analysis GUI

47

 

 

 

CHAPTER 4

Spike Detection and Feature
Extraction

4.1 Introduction

Spike detection refers to the identification of the arrival time of action potential
waveforms produced by a neuron during active communication with other
neurons in the nervous system. It is believed that spike arrival times carry all the
information about information processing and not the actual waveform shape.
Therefore extracting this temporal information is the first step to analyze neural
recordings [7].

Extracellular recordings are corrupted by many noise components such as
thermal and electrical noise, caused by signal amplifiers and other components
of the data acquisition system, background neural activity, and the occasional
similarity between the spike waveforms and the background noise [44]. Spikes
are non-stationary over a long period of time due to many reasons such as
electrode shifting, cell migrate, etc, and non-stationary background noise makes
the presence of spikes unclear. The goal of spike detection is to obtain temporal
information believed to be the most important parameter of neural activity from

the noisy extracellular recordings [7].

48

Many spike detection algorithms were proposed to overcome these
difficulties and they can be categorized as supervised or unsupervised, manual
or automated methods [7],[19],[18],[44].

In this study, we focus only on automated and unsupervised detection
methods. The simplest detection method is the amplitude threshold detection
method which identifies any signal that crosses the threshold as a spike. Even
though this method is simple, the detection performance is sensitive to noise.
Energy-based spike detectors compare the local power of the signal with a
threshold estimated from the power of noise [19]. This method is more robust to
noise than the amplitude threshold method.

The wavelet method has been motivated as an alternative to threshold or
energy based detection methods. Wavelet-based spike detection methods
transform the extracellular recordings into multiple sparse representation spaces
and compare them to a threshold. Wavelet-based spike detection methods were
suggested due to the sparsity they introduce and therefore plays an important
role in spike detection by increasing the signal-to-noise ratio (SNR) and localizing
transients in extracellular recordings [18],[44].

However the sparsity in wavelet detection methods requires a rule for
threshold selection because the distribution of the test statistic is not trivial to
derive. This causes a complexity in automating the detection methods, and the
threshold selection inevitably relies on an empirical approach.

Another problem of current wavelet detection methods is the relatively high

49

computational complexity due to the redundant detection routines in entire or
selected wavelet subspaces. In wavelet domain, spikes are represented in
relevant subspaces, and the general wavelet methods detect the spikes from
each subspace separately and combine the detection results either by majority
voting or logically OR them [44],[21]. Under this detection scheme, a single spike
is detected multiple times creating a redundant detection.

in conventional spike sorting algorithm, detected spikes are extracted and
properly aligned to obtain the feature set for spike sorting. These time consuming
steps require to store entire recording or segments of detected spikes which is
undesirable in large-scale neural data process.

In this chapter we propose a novel spike detection method using sparse
representation in wavelet domain. Combining information scattered across
multiple wavelet subspaces can reduce redundant detection. The proposed
threshold selection method modifies the distribution of the observation to
maximize separability between the noise and the signal distributions. We also
propose the feature extraction method that captures the compact representation
of the transient of the signal, called wavelet footprint. Under the proposed feature
extraction scheme, a compact feature set is obtained at the same time detecting
spikes that eliminates spike extraction and alignment steps as shown in figure 4-
1.

The proposed wavelet detection method incorporated with the threshold

selection was tested and the results were compared to several commonly used

50

spike detection methods. The proposed methods outperformed other methods in
low SNR cases, and the wavelet footprint extraction method reduced the several
steps in conventional spike sorting method: spike extraction, alignment and
feature extraction, while keeping the desired separability of clusters for the spike
sorting. Under the proposed wavelet footprint detection, significant amount of

data reduction can be achieved for offline processing.

 

Conventional Spike Sorting

 

 

Spike p Spike I Spike __ Feature

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Extra Detection Extraction Ali nment Extraction
-cellulaf _ Clustering i" Spike
Record: Tram
ng Wavelet
4‘ Footprint

 

 

 

 

 

 

 

Detection

 

 

LProposed Method

4

Figure 44. Block Diagram of the proposed wavelet footprint detection
method.

4.2 Theory

4.2.1 Observation Model

We formulate spike detection as a binary hypothesis test problem, where

51

under the null hypothesis Ho the signal is not present, while under hypothesis 7711

the signal is present

H0:y[n]=z[n] n=0,1,...,N—1
H1 :y[n]=a:[n]+z[n] n=0,1,...,N——1 (4.1)

1><N lxN . .
where y 6 ER is the observation, 117 E §R is action potentials

§RIXN

from single source, and Z 6 is a white Gaussian noise with

Z N N“), 02) where N is the number of samples in time domain [44].

Sparse representation of a transient signal has advantages in spike
detection mainly because it maximizes the SNR by characterizing non-transient
samples of the signal with large number of small coefficients, while discontinuities
of the signal with small number of large coefficients. Second, the sparse
representation tends to be highly localized in time domain, and this characteristic
makes the detection robust to variations in spike durations. Transforms that yield
the most compact representation of the signal are also desired to maximize the
spike detection performance, and wavelet-based techniques demonstrate

compact representation with precise frequency and time localizations [44].

4.2.2 Wavelet Transform and its Properties

52

Wavelet LIJ, also called a mother wavelet, is a function of oscillation with
finite energy and zero mean. A family of wavelet can be obtained by scaling and

translating the mother wavelet [40]

wa,b(t) = 713% (7) a,b e R

(4.2)

where a > O is the scale and b is the translation.
The projection of a random signal y(t) onto the subspace of scale a has the
form

y<a><t> = [Rutiwakodt

(4.3)

where 9(a)“) is the sparse representation of the signal y(t), a wavelet coefficient,
at the ath wavelet subspace [40].

The wavelet transform has several attractive properties for signal and image
processing [47]. Locality allows the wavelet atom localized simultaneously in time
and frequency, and Compression transforms real-world signals to sparse
representation. Persistence is a propagating tendency of wavelet coefficients
across scales, and it is the key property supporting the proposed method [47].
Since a sparse representation of a spike propagates across multiple subspaces,

there will be a redundancy across the relevant subspaces. Measuring the

53

correlation among the subspaces can capture the redundant representation of a
single spike, and this can reduce the number of detection routines.

To maximize the correlation, wavelet coefficients must be invariant among
different subspaces, and Stationary Wavelet Transform (SWT), designed to
achieve the translation-invariance, satisfies our need for this study [40]. Figure 4-
2 shows the sparse representation of spikes in wavelet domain and the
properties of wavelet are clearly illustrated.

Now we derive the model in wavelet domain. By the linearity of SWT at

scale j, observation y can be expressed as

ye) ___ We)

. 1 N .
where 31(3) 6 R x is wavelet coefficients at jth subspace and 20(3) is a
corresponding wavelet base [44]. With L level wavelet domain transform, the two

hypotheses in (4.1) have the following form:

H0:g[n]=_z_[n]n=0,1,~- ,N—1
’Hl:_y_[n]

alnl+élnl n=0,1,~-,N—1 (4.5)

54

100 uV
m
5
(I:

 

 

 

 

 

; DI
52222
:‘l 2 :l il 3; 103
:l :llt 21;:
till: lll
_l-§_‘flff

Figure 4-2. Sub-space representation of action potentials in wavelet
domain The signal is transformed into 4 level wavelet domains using SWT and
sym4 wavelet base. D and A denote detail and approximate nodes of SWT,
respectively, and the number indicates level.

 

 

 

 

 

gin. = .y(1)in]:- . - :y(L)[nll
gin = [sum nj , . . . ,33(L) [72]]
érni = :3“) 772, - ~ - i 23¢) [71]] (4.6)

where y[n], x[n], and z[n] are a snapshot of wavelet coefficients of the
observation, the spikes, and the noise, respectively, across multiple subspaces
mflmen

Let RWol’g) and 73011 L11) be the conditional risks associated with accepting

and rejecting the hypothesis 7'10 given the evidence y respectively. These risks

can be expressed as

55

73(H1ly) = A00P(H0|y)+/\01P(H1|y_)
73(Holy) = AtoPlffolyHAtthtly) (47)

’\(ii‘) = Amimj) is the loss incurred for deciding Hi when the true state of
nature is Hi. Since correct decisions are not penalized, /\00 and /\11 are zero.
This leads new expression of conditional risks described as

73(Hlly) = AOIPWIIQ) and Rmolai = AIOPWOIQ). The fundamental rule is to

decide 711, imelifl) < RWOIE) and vice versa. After invoking the Bayes rule

sz'ly) = PQI'HJ/PQ), the decision rule becomes

P(QI’H1) >71, /\01 Pail) é
)

Home) <“0 2:: Pitta ’7 (4.8)

 

Note that n represents the acceptance threshold for 711. Under the unsupervised
detection scheme, the covariance matrices of the signal and the noise are
unknown and they can be estimated from the observation. Under the assumption
that the noise is Gaussian we have PQIHO) "’ N“): 2) and PQIHI) ~ N(i#t2)
where pi is the L-component mean vector of the signal x and Z is the L-by-L
covariance matrix of the observation. The generalized likelihood ratio test

(GLRT) can be expressed as [45]

56

T —1 H
E 2 2222321097 (4.9)

This is essentially a blind energy detector [44]. The sufficient statistic T for spike

detection is

(4.10)

4.2.3 Optimal Threshold Selection for a x2 Distributed Signal

The threshold n in (4.9) is estimated considering the costs of false detection and
a priori probability of Ho and 'Hi [27]. Since information about the true action
potentials and the noise is unknown in unsupervised detection, n is estimated as
mean and variance of the observation [52]. Median Absolute Deviation (MAD) is
used to estimate the noise variance, but it is only valid with a Gaussian noise
assumption [52]. Since the square of a random variable with a normal distribution
is x2 distribution, the sufficient statistic T in (4.10) is x2 distributed. x2 distribution

is a special case of a gamma distribution expressed as

(4.11)

57

with a scale parameter [3 = 2, where Na) = (a — W, a shape parameter c =
v/2, and v being the degree of freedom of the random variable x.

An optimal threshold must be selected to maximize the positive detection
rate and to minimize the false detection rate [27]. Figure 4-3.a and 4-3.b illustrate
the probability density function (PDF) and the cumulative distribution function
(CDF) of x2 distribution for different v. Since T is a combination of the x2
distributed signal and noise with the parameter v close to 1, discriminating one
from another is not easy. Figure 4-4.a illustrates the x2 distribution of the
sufficient statistic 7'. Red bars represent histogram of noise and black bars
correspond to the signal. There is no clear separation between the two
distributions.

A modification of the histogram of T helps to estimate the threshold.
Histogram equalization (HE) is a method of contrast adjustment in image
processing using the histogram of the image [53]. This method increases the
global contrast of the image presented by close contrast value by linearizing the
cumulative distribution function (CDF) [53]. Linearization of the CDF to lead an
increase in the parameter v. If the HE shifts v from 1 to 3 or larger values, the
distributions of the noise and the signal are separable as shown in figure 4.3.a.

The general HE at bin index k is expressed as

_ cdf(k) — cdfmm
has) — Cdfmam '— Cdfmin X (L - 1) (4.12)

 

58

 

0.2

 

 

 

%fl_m_ 5'mvflmw _w

 

 

 

 

 

 

 

 

 

 

 

 

a) PDF of x2 distribution a) CDF of x2 distribution
1l—‘_ " T T f".
(7 l
0.8 ...:"i’ 5
5
0.6 ,........J 5
0‘ :¥cos_~w ii
CDF of MHE, .
----CDF of HE 5' l
0.2 5
00 of 0:4 " comma?“ __ "i

c) CDP of the three histograms
Figure 4-3. PDF and CDF of X2 distribution with different parameter v and
>the equalized CDPs using the general histogram equalization (HE) and the
modified histogram equalization (MHE)

where cdfmin being the minimum of CDF, and cdfmax being the maximum of

CDF, L being a total number of bins. However the regular HE method creates a
significant quantization distortion in the case of x2 distribution with v=1 as shown

in figure 4-3.c. The modified HE (MHE) for x2 distribution is expressed as

59

, (cdf(k)—cdf(k—1))xk k>1
h (k) :{ cdf(k) x k k =1 (“3)

The MHE yields much smoother linearization of the CDF of x2 distributed signal
with v=1. Once the histogram is equalized using the MHE, the histogram is
transformed into the modified histogram (MH) with v greater than 1 where the
separation between the noise and signal is more favorable.

Figure. 4-4.b shows the modified histogram of T that has a clear separation
between the noise and the signal distributions. The estimated threshold is located
at a local minimum between the two peaks. This is similar to the gray-level

histogram threshold selection in image processing [54].

4.2.4 Wavelet Footprint

The discontinuous structures of a signal often carry critical information, thus
efficient characterization of the discontinuity of the signal is a central task in
signal processing [47]. In general, larger wavelet coefficients tend to be around
the edge of a signal and these wavelet coefficients collect most of the energy of
the original signal. A wavelet footprint is defined as scale space vectors obtained
by gathering all the wavelet coefficients around the discontinuities of the signal

that model discontinuities in piecewise polynomial signal exactly [47].

60

15X10"4 Histogram of T

 

 

 

 

 

1O
:2 i
8 I
0 t.
5 ltiI -
8 2 o o 2 o 4 o o o a 1 1 2
Amplitude of T
(a) Histogram of sufficient statistic T
o Modified Histogram of T
m ' l - 1
m P
N
o
a, N
E
g :2
2
“’ l
O ' ' . ~ . . 1i I i IIIIIL
-0.2 0 0 5 1.0

Amplitude of T
(b) The Modified Histogram of sufficient statistic T

Figure 4-4. Histogram and modified histogram of sufficient statistic T (a)
The histogram illustrates the sparsity of T. Bright colored bars represent noise
while dark ones represent signals. (b) The modified histogram of T is shown.

61

Given a piecewise constant signal y with only one discontinuity at position k,

wavelet footprint 1km) is the scale-space vector obtained by gathering together

the largest wavelet coefficient from each node in the cone of influence of k.
This footprint 1km) can be written as

D} D1

D
13:0) _—_ [max(y(1)[k — ? : k + 7 L '

D
(L) __ __I:
], ,maa:(y [k 2 .k+ 2 I] (4.14)

. .th
where Di is duration of the cone of influence of l subspace at k. The cone of

influence (COI) is the region of the wavelet spectrum in which edge effects
become important and details of COl are found in [40].

We can expand the feature extraction method from the single channel
case to multi-channel data set when high correlations between neighboring
channels is present by concatenating a footprint from single channel together
with ones from other channels into one event vector, a long feature vector can
consolidate all the information in three domains: time, scale and space. Figure 4-

5 demonstrates how the method can be applied to multi- channel data.

4.3 Results

To evaluate the performance of the presented spike detection and the threshold

selection methods, simulated datasets with different SNRs were used. For

62

comparison, the datasets were obtained from OSort, a spike detection and

sorting package [11]. Stimulated datasets were generated by using a database of

'0': '00- (nor... r..- r...-
I I O O I
U 0 I u

| O . ‘ .
v o ,
' .6. -. ' ’ 2 . Q -. --' I -Q ‘- I.
. I. .' e ,2 I t I
I . O U

. O
I -'_ , ,. . J
'.r, . ' .

channel

4 i 5
2 5 ' at;
1 .

a
I

...........

channel
.' """""""" I' """" .' ' '7 """""" I- """""" I
o o . I l 9
Events ___________ i[”_“,n_g£ ____________ 5"-u_"_45
D'I — D3 A4
-- 02 — D4
Wavelet :' r r r '

footprint§___L_]__I IIIIIL WIWii

(b)

Figure 4—5 Wavelet footprint extraction method for the correlated multi-
channel signal. (a) Multi-channel data. A spike event, a snapshot of neural
action potential across the multi channels, is extracted and forms a vector of
event. (b) Wavelet footprint extraction from an event vector. Different colors
indicate specific nodes of SWT.

150 mean waveforms taken from well-separated neurons recorded in previous
experiments. The random background noise is generated by selecting randomly
scaled spike waveforms from the database and added to the noise traces.

Identifiable spikes are added by simulating a number of neurons (3 in dataset 1

63

and2, and 5 in dataset 3), with a renewal Poisson process with a refractory
period 3ms and a fixed firing rate between 1 to 10 HZ [11]. The proposed
detection algorithm was fully implemented in NeuroQuest, a software package for
neural data analysis, and all the results were obtained using NeuroQuest [55].

First, receiver operating characteristic (ROC) curves were plotted to
examine the detection accuracy of each method. The SNR was calculated as the
root-mean square (RMS) of the spike divided by the standard deviation of the
observation [11]. Three commonly used spike detection methods, single
amplitude detection (SAD) with a single threshold, absolute amplitude detection
(AAD) with both positive and negative thresholds, energy-based detection (EBD),
and the proposed wavelet detection (WD) were used for the comparison. For the
WD, the observation signal was transformed into wavelet domain using 5-level
SWT with a symlet 4 wavelet base. The sufficient statistic T was estimated from
D2 to D4 where D denotes a detail node of wavelet transform, because typical
action potentials have frequency range between 1 KHz to 5 KHz. From Figure 4-
6, SAD performed poorly, while WD performed better than the other methods in
the lowest SNR cases.

We compared two different threshold selection methods, the proposed
modified histogram equalization (MHE), and the MAD. For the comparison, we
used 3 and 5 times of the MAD, because other wavelet and energy based

detection methods typically estimate the detection threshold with 3 or 5 times of

64

the MAD [11]. The results shown in figure 4-7 illustrate the robustness of MHE
method in low SNR cases.

We compared the spike detection performance with EBD. EBD was selected for
the comparison because of its robust detection performance from the ROC test

and the x2 distributed EBD samples that the MHE method can be applied to.

 

 

    
   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Datasetf Dataset 2 Dataset 3 I-O-AAD l
1 1 5 . . 5 - 1 5 +SAD
5 ,g—r-o- —~b ----- c: ----- -$.\~°r-*"*"“""'""+ A EBD

O. ( 5 WE)

o. ‘5
Q.
l- l

O. I

0.

ROC Of SNR 2.2
0.5
FD
55..w.+.__ ...: 1 15—" '———"“"— 5
o. ; o 85 I 0.81 l/
: ”W‘”'*“~fi 5 ’1’
5 )

Q 0. Q o 655‘ Q 0.555 ”We
l- i— i—

0.4 0 0.4

O. 0. 0

HOC Of SNR 1.2 ROC Of SNR 1.3 HOC Of SNR 1.3
0 0.5 6 o5 ofl—w ”—m as
Fp Fp Fp

Figure 4-6. ROC for different datasets with different SNR AAD: absolute
amplitude detection, SAD: single amplitude detection, EBD: energy-based
detection, WD: the proposed wavelet detection. Tp and Fp denote True Positive
and False Positive respectively. Each column represents ROC of different
dataset and each raw represents different SNR.

65

TABLE 4-1. SPIKE DETECTION RESULTS OF EBD AND WD

 

 

 

Dataset3 (2986) SNR:5.4 SNR:2.7 SNR:1.8 SNR:1 .3
EBD + TP:100% TP:98.6% TP:91 .3% TP:76.1%

5 x MAD FP: 0% FP: 6% FP: 12% FP: 33%
EBD + TP:98.6% TP:93.5% TP:64.1% TP:54.4%

MH FP: 0% FP: 2% FP: 1% FP: 2%

WD + TP:96.4% TP:95.3% TP:93.5% TP:80%

5 x MAD FP: 1% FP: 5% FP: 48% FP: 54%
W0 + TP:96.4% TP:96.7% TP:75.4% TP:58.3%

MH FP: 1% FP: 1% FP: 1% FP: 3%

 

 

Dataset 3 was used for the comparison. There are 2986 spikes in the dataset 3.
TP and FP represent true positive and false positive of detection, respectively.
FP rate is calculated by the number of false detection divided by the maximum
false detection with a detection threshold value at 0.

Four different combinations of two detection methods and two threshold
selection methods were tested. From the result shown in Table 4-1, we conclude
that W0 is robust to the noise and MH yields low false detection rate for high
positive detection rate. In the highest SNR case, however, the performance of
WD was lower than EBD. The reason was that a small number of overlapped
spikes were represented as a single spike in the wavelet domain.

We examined the spike sorting result using the wavelet footprint feature
set. Once spikes were detected, the wavelet footprint was obtained by extracting
a local maximum of each subspace within the cone of influence of the detected
spike. Since the typical duration of single spike is 1-2 ms, the range of the cone
of influence is K :i: D/2, where K is being a position of the detected local

maximum and D being total samples equivalent to 2ms. D in this study was 50

with 25kHz of the sampling rate.

66

 

(a) Raw dataset3 with SNR 5.4 in time domain

88
+Histogram
83 x Median

- 5 x Median

 

 

 

 

 

_ab.._

 

 

Li}—

I")

 

 

 

 

 

 

   

Ml ' . . "I '
(b) Corresponding Sufficient Statistic of subspace redundancy and different
threshold selection methods

(c) Flaw dataset3 with SNR 1.8 in time domain

SS
->l<-Histogram
B3 x Median
213:: 5 x Median

 

 

 

 

 

 

    
 

 

    
    

 

 

  
     
   

2.7.... a- '
l' III-III II IIIIIl
l lllll'll‘ I

» :ll'l'llll (inhibit!

(d) Corresponding Sufficient Statistic in subspace redundancy and different
threshold selection methods

Figure 4-7. Sufficient statistic T of the segment of dataset 3 and the

different threshold selections Markers indicated the detection points with the
corresponding thresholds.

   

\\
\,u
_-!E 1"]

 

  

 

 

 

 

 

67

The spike sorting result of two feature sets - the wavelet footprint and the
entire spike waveform - were compared. First, we tested the method with
simulation data set that consists of two independent recordings. For display
purposes, we selected only the largest two principle components of each feature
set. All the spike sorting results were obtained by fuzzy c-mean clustering
method. Both Channel 1 and channel 2 were recorded from three neurons. As
the confusion matrices in figure 4-8 demonstrated, both the wavelet footprint and
the temporal PCA achieved the similar accuracy of spike sorting results.

We applied the method to extracellular recordings from the barrel cortex
of an anesthetized rat. The sampling frequency of the data was 25 KHz. Total
1849 spikes were detected, and the spike sorting results with the temporal PCA
and the wavelet footprints were compared. The results of spike sorting are shown
in figure 4—9. Three clear clusters were identified in temporal PCA, while four
clusters were identified in the wavelet footprint. The identified clusters in the
wavelet footprint were more scattered than those found in the temporal PCA
which is not desirable. However it helped to separate two clusters which were not
separable in the temporal PCA.

Finally, we applied the method to the multi-channel data recorded in the
dorsal cochlear nucleus of an anesthetized Guinea pig. The wavelet footprint
feature set was obtained as described in section 4.2.4, and the temporal PCA
feature set was extracted from the concatenated spike waveform extracted from

the each channel of the array.

68

TABLE 4-2. CONFUSION MATRICES OF THE SPIKE SORTING RESULTS

 

  
 

 

 

 

 

 

 
  
 
  

 

 

Channel Neuron Neuron Neuron Channel 2 Neuron Neuron Neuron
1 1 (181) 2 (159) 3 (25) 1 (173) 2 (190) 3 (32)
Classified 97.5% 1 (0) O Classified 100% 2 (7) 0
Neuron1 (95.1%) Neuron1 (100%)

Classified 2 (4) 98.3% 0 Classified 0 98% 0
Neuron2 (100%) Neuron2 (93%)

Classified 0 0 1 00% Classified O 0 100%
Neuron3 Q 00%) Neuron3 (100%)

 

 

( ) indicates the sorting result using the wavelet footprint feature set. Confusion
matrices show the similar spike sorting results of both feature sets.

1 00V

   

   
      

0.1Sec

r

Temporal PCA

Temporal PCA

Channel 1

5 5 Channel2 5

  
 

  

Wavelet PCA

5:

Wavelet PC

Figure 4-8 Clusters of temporal PCA and the wavelet footprint Raw signals,
feature spaces in temporal PCA and wavelet footprint, and templates of sorted
spikes. Channel 1 and 2 are recorded from three different neurons. Color
indicates links between spike templates and clusters.

69

 

   

Temporal PCA Wavelet Footprint

Spike1 (108) Spike2 (262) Spike1 (107) SpikeZ (271)
3.2451 to w w
Spike?’ (1479) Spike3 (894) Spike4 (577)

Figure 4-9 Spike sorting results of spontaneous recordings from an
anesthetized rat. 3 units are identified in Temporal PCA, whereas 4 units are
found in wavelet footprint domain. Spike 3 in temporal PCA can be further
separated into two units, spike 3 and spike 4 in the wavelet footprint domain.

For comparison, we fixed the number of clusters 6 and used the fuzzy-c
mean clustering method. As figure 4—10 illustrated, the two different feature
spaces displayed similar distribution. This implies that wavelet footprint feature
set shares the similar variation with the actual spike waveform, because PCA is a
projection of the data onto the eigenvector that has the direction of the largest
variation in the data [27]. This is an indication of the effective feature extraction

performance of wavelet footprint that can achieve a significant amount of data

reduction, 90% of compression rate per spike sampled at 25KHz.

70

 

 

   
  

       

."

Wavelet Footprint

 

 

 

 

 

 

 

Temporal PCA
Class 1 Class 2 Class 2
Class 4

\-

Y

Class 6

 

Figure 4-10 Spike sorting result of the multi-channel dataset The color-
coded clusters match with the corresponding spike waveforms. The width of the
shaded area in the spike waveform indicates the variance. Clusters and the
sorted spike waveforms in two different feature spaces, temporal PCA and
wavelet footprint, have similar distribution. This result implies that the wavelet
footprint is the effective sparse representation of the data.

4.4 Conclusion

Sparse representation of the extracellular recordings is key to achieve reliable
spike detection in low SNR, and the wavelet transform exhibits advantages in
characterizing these signals. We presented two contributions: 1) a novel spike
detection method, capturing the correlations across multiple subspaces, and 2)

the threshold selection, modifying the distribution of the observation to maximize

71

the separability between the noise and the signal. The proposed detection
method demonstrated improved detection performance in low SNFis.

We also have shown the simultaneous feature extraction method using
wavelet footprints, a sparse feature set that captures the information across
scales. Many different data sets including actual recorded data from anesthetized
animals were tested to verify the effectiveness of the method. The entire spike
waveform and the wavelet footprint feature set demonstrated very similar
distributions in PCA domain among many different data sets. Since the size of
wavelet footprint feature set is one tenth of the entire spike waveform, a massive
reduction in data transmission and storage can be achieved. The proposed
method also eliminates multiple processing steps in the conventional spike
sorting algorithm which is crucial step to implement the algorithm into a realtme

online sorting.

72

APPENDIX

73

Appendix A

NeuroQuest Input Data Structure

NeuroQuest works with a MAT file that contains a specific data structure.

The data file must have

data: A cell array that each cell contains 1 sec raw data.

BinWidth: Inverse of the sampling rate. 1/(sampling rate)

chanlnfo: channel labels. it is a raw vector. Ex) Four channel data with
channel index 1,3,5,7. chanlnfo = [1 3 5 7];

plotOption: a structure array that contains seven components (raw,
denoise, detection, stimulus, LFP, Spiketrain, and Trials). plotOption is a
set of flags that indicates types of data the file contains. If a data file has
only raw data, plotOption.raw is 1 and other values are 0.

fileDescription: a cell array that contains a short description of the data.

All the labels are case-sensitive and missing components will cause an error.

74

BIBLIOGRAPHY

75

[1]

[2]
[3]

[4]

[5]

l6]

[7]

[8]
[9]

BIBLIOGRAPHY

M. Scanziani and M. Hausser, “Electrophysiology in the age of light,” Nature, vol.
461, Oct. 2009, pp. 930-939.

G. Buzsaki, Rhythms of the Brain, Oxford University Press, USA, 2006.

“MEASUREMENTS OF VISUAL EVOKED POTENTIALS IN PARKINSON'S
DISEASE,” Dec. 1978.

A.M. Kuncel and W.M. Grill, “Selection of stimulus parameters for deep brain
stimulation,” Clinical Neurophysiology, vol. 115, Nov. 2004, pp. 2431-2441.

R.E. Kass, V. Ventura, and EN. Brown, “Statistical Issues in the Analysis of
Neuronal Data,” J Neurophysiol, vol. 94, Jul. 2005, pp. 8-25.

EN. Brown, R.E. Kass, and RP. Mitt-a, “Multiple neural spike train data analysis:
state-of-the-ait and future challenges,” Nature Neuroscience, vol. 7, May. 2004, pp.
456-461.

M.S. Lewicki, “A review of methods for spike sorting: the detection and
classification of neural action potentials,” Network: Computation in Neural Systems,
vol. 9, 1998, p. 53.

K.D. Harris, KlustaKwik. Automatic Cluster Analysis, version 1.0, 2002.

C. Fraley and AB. Raftery, “MCLUST: Software for model-based cluster analysis,”
Journal of Classification, vol. 16, 1999, pp. 297-306.

[10] KO. Quiroga, Z. Nadasdy, and Y. Ben-Shaul, “Unsupervised spike detection and

sorting with wavelets and superparamagnetic clustering,” Neural computation, vol.
16, 2004, pp. 1661-1687.

[11] U. Rutishauser, E.M. Schuman, and AN. Mamelak, “Online detection and sorting

of extracellularly recorded action potentials in human medial temporal lobe
recordings, in vivo,” Journal of neuroscience methods, vol. 154, 2006, pp. 204-224.

[12] KG. Oweiss, Statistical Signal Processing for Neuroscience and Neurotechnology,

Elsevier, in Press.

[13] EC. Leuthardt, G. Schalk, J .R. Wolpaw, J .G. Ojemann, and D.W. Moran, “A brain—

computer interface using electrocorticographic signals in humans,” Journal of
Neural Engineering, vol. 1, 2004, pp. 63-71.

76

[14] W.H. Pilcher and W.G. Rusyniak, “Complications of epilepsy surgery,”
Neurosurgery Clinics of North America, vol. 4, Apr. 1993, pp. 311-325.

[15] K.J. Miller, M. denNijs, P. Shenoy, J .W. Miller, R.P. Rao, and J .G. Ojemann,
“Real-time functional brain mapping using electrocorticography,” Neurolmage,
vol. 37, Aug. 2007, pp. 504-507.

[16] AD. Legatt, J. Arezzo, and H.G. Vaughan, “Averaged multiple unit activity as an
estimate of phasic changes in local neuronal activity: effects of volume-conducted
potentials,” Journal of Neuroscience Methods, vol. 2, Apr. 1980, pp. 203-217.

[17] K.G. Oweiss, “A systems approach for data compression and latency reduction in
cortically controlled brain machine interfaces,” IEEE Transactions on Bio-Medical
Engineering, vol. 53, Jul. 2006, pp. 1364-1377.

[18] X. Yang and SA. Sharnma, “A totally automated system for the detection and
classification ofneural spikes,” IEEE Transactions on Biomedical Engineering, vol.

35, 1988, pp. 806-816.

[19] Kyung Hwan Kim and Sung June Kim, “Neural spike sorting under nearly O-dB
signal-to-noise ratio using nonlinear energy operator and artificial neural-network
classifier,” IEEE Transactions on Biomedical Engineering, vol. 47, Oct. 2000, pp.
1406-1411.

[20] KG. Oweiss and D]. Anderson, “Tracking signal subspace invariance for blind
separation and classification of nonorthogonal sources in correlated noise,”
E URASIP Journal on Advances in Signal Processing, vol. 2007, 2007.

[21] Z. Nenadic and J .W. Burdick, “Spike detection using the continuous wavelet
transform,” IEEE Transactions on Biomedical Engineering, vol. 52, 2005, pp. 74-
87.

[22] MS. Fee, P.P. Mitra, and D. Kleinfeld, “Variability of extracellular spike
waveforms of cortical neurons,” J Neurophysiol, vol. 76, Dec. 1996, pp. 3823-
3833.

[23] MS. Fee, P.P. Mitra, and D. Kleinfeld, “Automatic sorting of multiple unit neuronal
signals in the presence of anisotropic and non-Gaussian variability,” Journal of
Neuroscience Methods, vol. 69, 1996, pp. 175-188.

[24] T. Hermle, C. Schwarz, and M. Bogdan, “Employing ICA and SOM for spike
sorting of multielectrode recordings from CNS,” Journal of Physiology-Paris, vol.
98, 2004, pp. 349-356.

[25] J.C. Bezdek and R. Ehrlich, “FCM: The fuzzy c-means clustering algorithm,”
Computers & Geosciences, vol. 10, 1984, pp. 191-203.

77

[26] S. Shoham, M.R. Fellows, and RA. Normann, “Robust, automatic spike sorting
using mixtures of multivariate t-distributions,” Journal of neuroscience methods,
vol. 127, 2003, pp. 111-122.

[27] RD. Duda, P.E. Hart, and DO. Stork, Pattern Classification, Wiley-Interscience,
2000.

[28] S. Takahashi, Y. Anzai, and Y. Sakurai, “A new approach to spike sorting for multi-
neuronal activities recorded with a tetrode--how ICA can be practical,”
Neuroscience research, vol. 46, 2003, pp. 265-272.

[29] K.G. Oweiss and D]. Anderson, “MASSIT-Multiresolution Analysis of Signal
Subspace InvarianceTechnique: a novel algorithm for blind source separation,”

Signals, Systems and Computers, 2001 . Conference Record of the Thirty-Fifth
Asilomar Conference on, 2001.

[30] JD. Victor and KP. Purpura, “Metric-space analysis of spike trains: theory,
algorithms and application,” Network: Computation in Neural Systems, vol. 8,
1997,p.127.

[31] G. Palm, A. Aertsen, and G.L. Gerstein, “On the significance of correlations among
neuronal spike trains,” Biological Cybernetics, vol. 59, 1988, pp. 1-1 1.

[32] G.E.P. Box, G.M. Jenkins, and GO Reinsel, Time series analysis: forecasting and
control.

[33] A.M. Aertsen, G.L. Gerstein, M.K. Habib, and G. Palm, “Dynamics of neuronal
firing correlation: modulation of "effective connectivity",” J Neurophysiol, vol. 61,
May. 1989, pp. 900-917.

[34] S. Eldawlatly, R. Jin, and K.G. Oweiss, “Identifying functional connectivity in

large-scale neural ensemble recordings: A multiscale data mining approach,” Neural
computation, vol. 21, 2009, pp. 450-477.

[3 5] S. Eldawlatly, Y. Zhou, R. Jin, and K.G. Oweiss, “On the use of dynamic bayesian
networks in reconstructing functional neuronal networks from spike train
ensembles,” Neural computation, vol. 22, 2010, pp. 158-189.

[36] R Meier, U. Egert, A. Aertsen, and M.P. Nawrot, “F IND--A unified framework for
neural data analysis,” Neural Networks, vol. 21, 2008, pp. 1085-1093.

[37] H. Bokil, P. Andrews, H. Maniar, B. Pesaran, J. Kulkami, C. Loader, and P. Mitra,

“Chronux: a platform for analyzing neural signals,” BMC Neuroscience, vol. 10,
2009, p. SB.

78

[3 8] AS. Sedra and PO. Brackett, Filter theory and design: active and passive, Matrix
Pub, 1978.

[39] P.G. Musial, S.N. Baker, G.L. Gerstein, E.A. King, and J .G. Keating, “Signal-to-
noise ratio improvement in multiple electrode recording,” Journal of Neuroscience
Methods, vol. 115, Mar. 2002, pp. 29-43.

[40] S.G. Mallat, A wavelet tour of signal processing, Academic Pr, 1999.

[41] K.G. Oweiss, “Multiresolution analysis of multichannel neural recordings in the
context of signal detection, estimation, classification and noise suppression,” Rice
University, 2002.

[42] S.G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image
denoising and compression,” IEEE Transactions on Image Processing, vol. 9,
2000, pp. 1532-1546.

[43] BL. Donoho, “De-noising by soft-thresholding,” IEEE transactions on information
theory, vol. 41, 1995, pp. 613-627.

[44] K.G. Oweiss, “Multiresolution analysis of multichannel neural recordings in the
context of signal detection, estimation, classification and noise suppression,” Rice
University, 2002.

[45] Y. Suhail, “Methods for neural signal processing and analysis,” MICHIGAN
STATE UNIVERSITY, 2005.

[46] A. Pavlov, V.A. Makarov, I. Makarova, and F. Panetsos, “Separation of
Extracellular Spikes: When Wavelet Based Methods Outperform the Principle

Component Analysis,” Mechanisms, Symbols, and Models Underlying Cognition,
2005, pp. 123-132.

[47] PL. Dragotti and M. Vetterli, “Wavelet footprints: Theory, algorithms, and
applications,” IEEE Transactions on Signal Processing, vol. 51, 2003, pp. 1306-
1323.

[48] M.A.T. Figueiredo and AK. Jain, “Unsupervised Learning of Finite Mixture
Models,” IEEE TRANSACTIONS ON PA TTERN ANALYSIS AND MACHINE
INTELLIGENCE, vol. 24, 2000, pp. 381--396.

[49] H. Ito and S. Tsuji, “Model dependence in quantification of spike interdependence
by joint pen-stimulus time histogram,” Neural computation, vol. 12, 2000, pp. 195-
217.

[50] P. Dayan, L.F. Abbott, and L. Abbott, Theoretical neuroscience: Computational
and mathematical modeling of neural systems, MIT Press, 2001.

79

[51] C.K. Knox, “Cross-correlation functions for a neuronal model,” Biophysical
Journal, vol. 14, 1974, pp. 567-582.

[52] D.L. Donoho and J. JOHNSTONE, “Ideal spatial adaptation by wavelet shrinkage,”
Biometrika, vol. 81, 1994, p. 425.

[53] SM. Pizer, E.P. Ambum, J.D. Austin, R. Cromartie, A. Geselowitz, T. Greer, B. ter
Haar Romeny, J .B. Zimmerman, and K. Zuiderveld, “Adaptive histogram
equalization and its variations,” Computer vision, graphics, and image processing,
vol. 39, 1987, pp. 355-368.

[54] N. Otsu, “A threshold selection method from gray-level histograms,” Automatica,
vol. 1 1, 1975, pp. 285-296.

[55] K.Y. Kwon, S. Eldawlatly, and K.G. Oweiss, “NeuroQuest: A Comprehensive Tool

for Large Scale Neural Data Processing and Analysis,” 4th Int IEEE EMBS,
Turkey: 2009, pp. 622-625.

80

MICHIGAN STATE UNIVERSITY LIB

III I III II lllllllllllms
3

4

l l

3 1293 03063 77