IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 5, MAY 1997 403
Detection, Classification, and Superposition
Resolution of Action Potentials in Multiunit
Single-Channel Recordings by an On-Line
Real-Time Neural Network
Rishi Chandra and Lance M. Optican*
Abstract—Determination of single-unit spike trains from mul-
tiunit recordings obtained during extracellular recording has
been the focus of many studies over the last two decades. In
multiunit recordings, superpositions can occur with high fre-
quency if the firing rates of the neurons are high or correlated,
making superposition resolution imperative for accurate spike
train determination. In this work, a connectionist neural net-
work (NN) was applied to the spike sorting challenge. A novel
training scheme was developed which enabled the NN to resolve
some superpositions using single-channel recordings. Simulated
multiunit spike trains were constructed from templates and noise
segments that were extracted from real extracellular recordings.
The simulations were used to determine the performances of the
NN and a simple matched template filter (MTF), which was used
as a basis for comparison. The network performed as well as the
MTF in identifying nonoverlapping spikes, and was significantly
better in resolving superpositions and rejecting noise. An on-line,
real-time implementation of the NN discriminator, using a high-
speed digital signal processor mounted inside an IBM-PC, is now
in use in six laboratories.
Index Terms— Multiunit spike sorting, neural networks, on-
line, real-time discrimination.
I. INTRODUCTION
T
HE analysis of simultaneous activity from several
neurons can lead to a better understanding of their
functional connectivity. Extracellular recordings with low-
impedance electrodes are capable of recording such activity
from several neurons near the tip of the electrode. The
shape manifested in the recording of an action potential,
or spike, emitted from a particular neuron is a function of
the position of that neuron relative to the electrode. Thus,
by sorting the different shapes in an extracellular recording,
the interleaved spike trains of the individual neurons can be
separated. There are two significant difficulties in analyzing
extracellular recordings: high levels of correlated neural noise
and superpositions of waveforms. Background noise consists
mainly of the activity of distant neurons. These noise levels
can be relatively high and since most neurons have similar
time constants for potential changes, the spectral contents
Manuscript received July 24, 1995; revised January 7, 1996. Asterisk
indicates corresponding author.
R. Chandra is with the Laboratory of Sensorimotor Research, National Eye
Institute, NIH, Bethesda, MD 20892 USA.
*L. M. Optican is with the Laboratory of Sensorimotor Research, National
Eye Institute, NIH, Bethesda, MD 20892 USA (e-mail: lmo@lsr.nei.nih.gov).
Publisher Item Identifier S 0018-9294(97)02959-5.
of the noise and signal overlap considerably. Superpositions
occur when two neurons fire within 1 ms of each other,
resulting in the addition of their individual waveforms. In the
case where firing patterns are not correlated, the frequency of
superpositions is determined by the number of neurons and
their firing rates. If there are
neurons, each firing with an
average frequency of
Hz, and having a spike duration of
s, the average superposition rate is given by the expression
[1]. For example, if there are three neurons each
firing at 50 Hz, with waveforms of 1-ms duration, then the
average number of superpositions would be 15 s
. Doubling
the common firing rate to 100 Hz causes the overlap rate to
increase fourfold to 60 s
. Thus, superposition resolution is
an important consideration in spike sorting. When the firing
patterns are correlated to each other or a particular stimulus,
the rate of superpositions can be much higher. This paper
presents a fully automated system for sorting spikes and their
superpositions on-line in real-time.
Three criteria governed the development of our spike sorter,
1) good isolation of single units in multiunit single-channel
recordings; 2) an on-line system, because details of experi-
ments are often determined by the selectivity of the identified
neurons; and 3) resolution of superpositions, because many
neurons burst together when a stimulus is presented, which
results in many superpositions.
Many signal-processing techniques have been applied to
the challenge of sorting multiunit recordings (see the review
paper by Schmidt [2], and a comparison of various sorters
by Wheeler and Heetderks [3]). One of the most popular
methods is matched template filtering, commonly referred to
as template matching [4]–[13]. This technique uses templates
that represent the typical waveform shape of each neuron
to filter the data. Variants of template matching have been
implemented in the form of reduced feature sets, such as prin-
cipal component analysis or simple waveform characteristics
(e.g., spike height and width), to provide faster processing
[14]–[16]. A learning phase is required for the template
generation and determination of classification thresholds for
all these sorters. Some investigators have considered the
problem of superpositions [8], [13]. The general scheme
for superposition resolution in template matching relies on
iterative subtraction of all possible template combinations from
unidentified waveforms. This approach is computationally
U.S. Government work not protected by U.S. copyright.
404 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 5, MAY 1997
intensive, making on-line implementation difficult, especially
when the rate of superpositions is high.
Matched template filtering is a time domain technique;
frequency domain methods have also been applied to spike
sorting, the most common technique being optimal linear
filtering [17]–[21]. This approach is based on deriving optimal
filters that respond to one unique template and reject the
others and the background noise. The learning phase generates
templates which are then used to derive their respective
optimal linear filters. Although linear optimal filtering can
resolve linear superpositions, it is dependent on multichannel
recordings. Roberts and Hartline [17] stated that the number
of channels should be greater than or equal to the number
of units to be identified. The optimal filtering technique is
ideal for multichannel recordings where the units appear on
at least two of the channels. All the investigators using this
technique [17]–[21] concluded that the discrimination perfor-
mance would deteriorate rapidly if this condition were not met,
although this dependency may be reduced by adding nonlinear
properties to the filters, e.g., an iterative subtraction scheme
that removes the largest remaining units on each iteration.
An off-line Bayesian modeling scheme has been proposed
by Atiya [22] and Lewicki [23]. This method uses both
the waveform and firing rate information to minimize the
probability of error during classification. To resolve super-
positions, Atiya implements a comprehensive search of the
space of overlapping spike shapes and event times to find the
sequence of maximum probability. This is a very computa-
tionally demanding approach. Lewicki uses the same principle
but restricts the search by using a Poisson model for the
firing of the neurons, then uses data structures and dynamic
programming to efficiently sort the spikes. Lewicki’s technique
can be implemented in a pseudo-on-line system by processing
one set of data while collecting another. If results of the spike
classification are needed before the end of an experimental
trial, a true on-line system would still be needed.
Artificial neural networks (ANN’s) have been used for
the classification of units in a multiunit recording with
fully connected, feed-forward networks [24]–[26]. An
implementation of the adaptive resonance theory (ART)-2
neural network (NN) algorithm was also applied to this
problem [27]. Although these applications have been
successful in isolating single units and are better than template
matching [25], none were trained to recognize, or were tested
on, superpositions of waveforms. A Hopfield-type network
has been implemented in real-time to separate units and their
superpositions [28], although the investigators concluded that
template matching performed better than the Hopfield-type
network. A unique pseudo-unsupervised training method
based on network relaxation has been implemented off-line
[29], [30] that learns the template shapes. This application
was developed to resolve motor unit action potentials
(MUAP’s), which are significantly longer than neuronal
action potentials. Superposition resolution was possible when
the MUAP peaks were sufficiently separated in time.
An NN was used in the work presented here because
it combines linear filtering (by the weights from the input
layer to the hidden layer) with nonlinear classification (by
the nonlinear hidden units), which should enable it to resolve
superpositions of units in a multiunit single-channel recording.
Another advantage of using an NN for sorting is that in
multidimensional feature space, decision regions are hyper-
ellipses, because the background neural noise is colored-
Gaussian-type noise (if the noise were white, then the decision
regions would be hyper-spheres, i.e., equal variance in all
dimensions), and NN’s are capable of generating nonspherical
decision boundaries. The final advantage of using an NN lies
in the ability to implement it on high-speed digital processors,
allowing on-line real-time multiunit sorting with superposition
resolution. This paper presents a method for training an NN to
sort spikes and resolve superpositions. This training method
has three novel components. First, the training set consists
of both individual templates and superimposed templates.
Second, the training set includes simulated noise with the
same spectral characteristics as the neuronal noise. Third, the
NN is used as a nonrecursive [finite impulse response (FIR)]
digital filter, simultaneously accomplishing both detection
and classification of spikes. The performance of this NN is
compared to that of an MTF on simulated data sets with
different templates and signal-to-noise ratios (SNR’s).
The method developed here has four characteristic features:
automatic determination of a detection threshold based on
noise analysis, automatic clustering of waveforms to obtain
templates, supervised training of an NN, and operation of the
network as a simultaneous detector/classifier in real-time. This
method has been successfully implemented on a PC-host with
a digital signal processor (DSP) co-processor system.
II. M
ETHODS
A. Data Collection and Spike Detection
Data were acquired during extracellular recording with
tungsten microelectrodes from the medial superior temporal
area in the cortex of a monkey during a visual experiment
with moving stimuli (see Duffy and Wurtz [31] for methods).
The system has also performed well on recordings from
primate superior colliculus. Data sets were collected from three
different recording sessions and stored on analog tape. From
each session, a section of the tape with low background noise
levels was selected and digitized at a rate of 24 KHz using
a 16-b analog–digital (A/D) converter (ADC-DBCS5339-50,
Communication Automation and Control, Allentown, PA).
The ADC utilizes a front end antialiasing filter, delta-sigma
modulation, 64-times oversampling, and a three-stage digital
FIR filter to achieve a 94-dB SNR.
Spikes were detected when the amplitude of the recorded
signal exceeded a positive or negative threshold. The threshold
was determined by breaking the record into pure neuronal
noise segments and setting the threshold at three standard de-
viations of that noise [32]. The neuronal noise was segmented
by an algorithm which separates signal from noise based
on the Gaussian characteristics of the noise [33]. Detected
spikes were clustered together to form noise-free templates
using a recently developed simultaneous clustering algorithm
[34]–[36]. This algorithm exploits the advantages of simul-
CHANDRA AND OPTICAN: DETECTION, CLASSIFICATION, AND SUPERPOSITION RESOLUTION OF ACTION POTENTIALS 405
Fig. 1. Artificial data sets based on three separate neurophysiological record-
ing sessions. Three templates that represented the most dense spike clusters
with high SNR were chosen from each recording session. The templates are
shown embedded in a short segment of their corresponding neural noise. The
SNR’s for each template in the three sets were: 1) 7.72, 3.46, and 3.26; 2)
7.51, 3.46, and 4.72; and 3) 2.98, 2.40, and 4.93.
taneous clustering and accounts for detection jitter when
segmenting waveforms from recordings. The algorithm is
described briefly in the Appendix. Three templates, which
represented dense spike clusters with high SNR’s, and the
segmented neural noise, were stored from each of the three
data sets. Each template consisted of 24 samples (1.0 ms in
duration at 24 KHz). Fig. 1 shows the templates embedded in
a short segment of their respective noise traces.
B. Discriminator Details
This study compares the spike discrimination performance
of a connectionist NN to that of an MTF. Fig. 2 illustrates the
structure of the two discriminators. The MTF [Fig. 2(a)] sets
circular decision regions [13] in multidimensional space. The
incoming data stream takes the form of vectors in this space,
i.e., values for each dimension (feature) are recorded over the
length of an analysis window. In our implementation, a 24-
sample window was used, where each sample in the window
constitutes one feature. The templates are represented in the
same manner and the sum squared difference between the
incoming data and the templates is calculated. If the error
is a local minima and less than threshold, then the incoming
pattern is said to match the template. Since the noise is colored,
the threshold was determined empirically from a pure noise
segment. This was done by building a histogram of the power
(sum of squared samples in a 24-sample window) of the noise
record. The threshold was set at a value that was greater
than 99.9% of the population. This method approximates the
squared Euclidean distance threshold for classification and
avoids the errors that can be caused by outliers.
(a)
(b)
Fig. 2. The structures of the discriminators. The input to each discriminator
is the digitized signal, which is shifted one sample at a time. In (a) the
MTF, the template is subtracted from the signal, the differences squared and
summed. The output of each template filter is compared to a threshold deter-
mined from the statistics of the noise. (b) The fully connected feed-forward
NN is used as a nonlinear FIR filter. The outputs are compared to a threshold
which is empirically determined. Input layer units have a linear function, the
hidden layer units have a Gaussian function, and the output layer units have
a sigmoid function.
The other discriminator [Fig. 2(b)] was an NN used as an
FIR filter [37], [38]. A fully connected, feed-forward, three-
layer NN was used. The network had twenty-four input units,
eight hidden units, and three output units. Linear activation
functions were used for the input nodes, Gaussian activation
functions for the hidden layer, and logistic activation functions
(range from zero to one) for the output layer. The last 24
samples of the incoming data formed the input to the network.
An event constituted a strong response (
0.95) in any of the
output units. In simulations we tested different numbers of
hidden units. Below eight hidden units, the classification per-
formance dropped drastically. Although eight hidden units is
the maximum we could implement in the real-time system, we
tested networks with ten and twelve hidden units and no signif-
icant improvement in classification performance was observed.
The network was trained with a modification of the back-
propagation algorithm [39]. Modifications that were incorpo-
406 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 5, MAY 1997
rated in the training were an adaptive learning rate, weight
momentum, weight decay, and weight annealing. The training
for each test consisted of each template, with added noise, at all
possible shifts along the input buffer of the NN. The weighting
for the error at the output layer was an exponential function of
the shift. The training input–output set was defined as follows.
If
is the total number of templates, is the number
of possible shifts (
, in this training set),
template (
) stands for template with shift , and
represents the weighting of the output error, then, the inputs
are
template noise, for
and
the outputs are
for and
otherwise
and the weights for the error are
for
The training set with superpositions adds the following inputs:
template template noise,
for
and
the corresponding outputs and weights are
for
,
for all and
for ,
for and
otherwise.
for all
This input–output set trains the network to respond only when
a spike is placed exactly in the middle of the input buffer. Any
shifts cause the network to reject the waveform. In the case
of superpositions that are
samples apart, ideally the network
would respond to the first spike and then respond again n
samples later to the second spike. Fig. 3 shows the output
of the two discriminators for a short sample of data. Data is
the input,
is the corresponding output of the network, and
is the squared difference between the waveform and MTF
(
).
C. Simulations and Testing
Software simulations were used to determine the perfor-
mance of the MTF and NN discriminators. As a simple model
of the recording process, assume that when a neuron fires it al-
ways produces the same waveform shape in the recording, but
is corrupted by colored Gaussian noise. Then each template,
which is an average of many action potentials determined
to be from one neuron, represents the “noise-free” shape
Fig. 3. A 10-ms segment of a simulated spike train with five spikes, and
the corresponding outputs of the NN and MTF. The top trace shows the data,
which has an SNR for the three templates of 5.74, 3.61, and 2.65. The numbers
show the position of each template. The outputs of the two discriminators are
shown in the traces below.
Ni
is the
i
th output of the NN,
Mi
is the output of
the
i
th MTF. “
C
” represents a correct classification, “
M
” a missed detection,
and “
F
” a false detection. The NN correctly classifies five out of the six
waveforms (it misses one of type 3) and has one false positive. The MTF
detects every waveform but also has two false positives. Note that the output
traces of the two detectors have been shifted forward by 12 samples to align
the outputs with the middle of the spikes for illustration purposes.
manifested in a recording by the respective neuron [1]. Thus
simulated multispike trains were constructed by embedding
the templates at various locations in extended noise traces.
Neural noise for each data set was modeled by generating
white Gaussian noise [40] and digitally filtering it with a
second-order recursive (IIR) filter with three nonrecursive
coefficients and two recursive coefficients to match the spectral
characteristics of the real neural noise. The coefficients of the
recursive filter were determined by fitting a single exponential
to the decay of the autocorrelation function [41].
Simulations allowed complete control of test data and exact
evaluation of each discriminator’s performance. For each
simulated spike train, the variables in the record that could
be manipulated were: 1) the firing rate of each neuron; 2)
superpositions of waveforms; and 3) the background noise
level. Three different types of simulated spike trains were
constructed for each data set; the first type had an approximate
firing rate of 120 Hz for each neuron with superpositions
occurring at a rate of 80 Hz, i.e., about 2/3 of the spikes present
were the superposition of two or three individual spikes. Fig. 4
shows two examples of superimposed waveforms from the
first type of simulations. The second type had a firing rate of
90 Hz for each neuron and no superpositions. For these two
CHANDRA AND OPTICAN: DETECTION, CLASSIFICATION, AND SUPERPOSITION RESOLUTION OF ACTION POTENTIALS 407
(a)
(b)
Fig. 4. Examples of the superpositions of (a) three and (b) two templates.
The bottom traces in each example show the superpositions of the templates
shown above with additive noise. The average SNR of the spike train these
segments were taken from was 5.0. In (a), the overlap is 0.71 ms between A
and B, 0.67 ms between B and C, and 0.38 ms between A and C. The NN
was able to resolve the earliest and largest (C) waveform and the matched
filter rejected all three. In (b), the overlap between A and B is 0.88 ms. The
NN was able to identify both constituent waveforms. The MTF rejected both.
simulation types, the SNR was varied from 1.5 to 10.0 by
simply scaling the background noise. It is important to note
that the noise had the same autocorrelation characteristics as
the real noise from the recording sessions, and only its variance
differed. At each SNR, 4 s of data were tested for each data
set. The third and last simulation type was designed to test the
performance of the two discriminators when the amplitude of
the spike changed, because spike amplitude commonly drops
during a high-frequency burst. The noise level was fixed to
give an average SNR of 5.0 (relative to the original template
sizes) and the magnitude of the template was varied between
0.5 and 2.0 times its original size. The original template sizes
were used for the matched filter and training of the NN to
test the robustness of the two discriminators in the presence
of amplitude variation.
D. Real-Time Implementation
The matched filter and NN discriminators were realized
using an 80-MHz AT&T DSP32C digital signal processor
(Communications Automation & Control) mounted inside an
IBM Pentium P5-133 MHz PC (Gateway). The DSP was fast
enough to run either discriminator in real-time. The maximum
size of the network that could be implemented with a 24-
sample input buffer and a 24-KHz sampling rate was eight
hidden units and three output units. The output of the chosen
TABLE I
AS
AMPLE OF THE CLASSIFICATION SCHEME.ANASTERISK INDICATES A
SUPERPOSITION, OR OVERLAP OF TWO TEMPLATES.THE LOCATIONS IN THE
TABLE ARE THE SAMPLE NUMBERS OF THE WAVEFORMS.WHEN SOMETHING IS
DETECTED, IT IS REGISTERED AS A CORRECT CLASSIFICATION IF THE LOCATION
IS EXACTLY RIGHT AND THE IDENTIFIED CLASS MATCHES THE TRUE CLASS.
A
NY
OTHER DETECTION IS CLASSIFIED AS A FALSE POSITIVE.ERRONEOUS
DETECTIONS OF NOISE AS UNITS ARE ALSO CLASSIFIED AS
FALSE POSITIVES
discriminator was accumulated and relayed to the host PC
every millisecond for storage and on-line display. The learning
stage consisted of two phases, template generation, and NN
training. Fifteen seconds of data buffers were uploaded to
the PC for noise modeling, spike detection, and template
generation. The required processing time on an IBM P5-133 is
approximately 3–4 s. Network training was performed on the
DSP, with a typical training time of 60 s. Thus the total time
to acquire new data, generate templates, and train the network
is just under 2 min for this system. The code for the host PC
was written in the C programming language under MS-DOS;
the code for the DSP was written in assembly language.
III. R
ESULTS
The simulations described above were designed to compare
the discrimination performance of the NN to that of the MTF.
Two measures of performance were considered, the percent
correctly identified and the rate of false positives. A correct
identification means that the waveform was detected at the
correct location and was classified correctly. An incorrect
identification means either that the waveform was not detected
at all, or was detected at the correct location but was misclassi-
fied. A false positive is any detection where no waveform had
been placed. Table I gives a sample of the evaluation process.
In this scheme, correct classification of both constituent
waveforms in a superposition counts as two correct classifi-
cations and resolution of two superimposed waveforms. Thus,
a superposition resolution of 50% could mean that only one
waveform was resolved from every superposition, or that both
waveforms from half the superpositions were resolved, or a
combination of the two extreme cases. In the examples given
in Table I, the superposition resolution is 75%, as three out of
four waveforms were correctly resolved. The results presented
are the average performances of the two discriminators over
the three simulated data sets. The SNR’s quoted are the
average SNR’s of the spike trains. The SNR is defined as
the root mean square value of the waveform divided by the
standard deviation of the noise. The SNR’s for a data set are
simply the average SNR’s of their constituent templates.
