PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of springeropenLink to Publisher's site
Journal of Computational Neuroscience
 
J Comput Neurosci. 2010 August; 29(1-2): 127–148.
Published online 2009 June 5. doi:  10.1007/s10827-009-0163-5
PMCID: PMC2950077

An online spike detection and spike classification algorithm capable of instantaneous resolution of overlapping spikes

Abstract

For the analysis of neuronal cooperativity, simultaneously recorded extracellular signals from neighboring neurons need to be sorted reliably by a spike sorting method. Many algorithms have been developed to this end, however, to date, none of them manages to fulfill a set of demanding requirements. In particular, it is desirable to have an algorithm that operates online, detects and classifies overlapping spikes in real time, and that adapts to non-stationary data. Here, we present a combined spike detection and classification algorithm, which explicitly addresses these issues. Our approach makes use of linear filters to find a new representation of the data and to optimally enhance the signal-to-noise ratio. We introduce a method called “Deconfusion” which de-correlates the filter outputs and provides source separation. Finally, a set of well-defined thresholds is applied and leads to simultaneous spike detection and spike classification. By incorporating a direct feedback, the algorithm adapts to non-stationary data and is, therefore, well suited for acute recordings. We evaluate our method on simulated and experimental data, including simultaneous intra/extra-cellular recordings made in slices of a rat cortex and recordings from the prefrontal cortex of awake behaving macaques. We compare the results to existing spike detection as well as spike sorting methods. We conclude that our algorithm meets all of the mentioned requirements and outperforms other methods under realistic signal-to-noise ratios and in the presence of overlapping spikes.

Keywords: Realtime spike sorting, Extracellular multi electrode recordings, Tetrode recordings, FIR filters, Deconfusion

Introduction

In order to understand higher brain functions and the interactions between single neurons, an analysis of the simultaneous activity of a large number of individual neurons is essential. One common way to acquire the necessary amount of neuronal activity data is to use simultaneous extracellular recordings, either with single electrodes or, more recently, with multi electrodes like tetrodes (O’Keefe and Recce 1993). However, the recorded data does not directly provide the isolated activity of single neurons, but a mixture of neuronal activity from many neurons additionally corrupted by noise. The task of so called “spike sorting” algorithms is to reconstruct the single neuron signals (i.e. spike trains) from these recordings. Many approaches for analyzing the data after acquisition, i.e. offline spike sorting algorithms, have been developed in the last years; see for example Vargas-Irwin and Donoghue (2007), Delescluse and Pouzat (2006), Pouzat et al. (2004), Kim and Kim (2003), Takahashi et al. (2003), Shoham et al. (2003), Hulata et al. (2002), Lewicki (1998), Fee et al. (1996a). Although more methods are available in this category, there are several reasons to favor methods which provide results already during the recordings, termed realtime online sorting algorithms. For example, realtime online spike sorting techniques are indispensable for conducting “closed-loop” experiments and for brain-machine interfaces (Rutishauser et al. 2006; Obeid and Wolf 2004). The few existing approaches to realtime online sorting (Thakur et al. 2007; Rutishauser et al. 2006; Aksenova et al. 2003) are clustering based and have at least one of the following drawbacks: 1) They are not explicitly formulated for data acquired from multi electrodes, 2) they do not resolve overlapping spikes, 3) they do not perform well on data with a low signal-to-noise ratio 4) they are not able to adapt to non-stationarities of the data as caused by tissue drifts. We discuss the reasons and importance of these issues in the following:

  1. Multi electrodes (e.g. tetrodes) provide significantly more information about the local neuronal population than single electrodes (Harris et al. 2000; Rebrik et al. 1999). Having several recording electrodes closely spaced instead of one, the same action potential is present on more than one recording channel. The so called stereo-effect—a neuron specific amplitude distribution among the recording channels—allows for a better discrimination between action potentials from different neurons (Gray et al. 1995). This allows also for more a reliable resolution of overlapping spikes.
  2. With tetrodes recording an increased number of neurons compared to high impedance single electrodes, overlapping spikes are more likely to occur. Also, studies stress the relevance of ensemble coding, which translates into local synchronized firing and hence a raised occurrence frequency of overlapping spikes (Sakurai and Takahashi 2006). To identify such a code, the resolution of overlapping spikes is crucial and efforts have been made addressing this issue (Ding and Yuan 2008; Wang et al. 2006; Zhang et al. 2004; McGill 2002; Chandra and Optican 1997). However, the cited approaches are all computationally very expensive, making a realtime online implementation difficult. One of the reasons for this computational complexity is the implementation of separate sub-routines for the processing of overlapping spikes, which, additionally, are more complex than the processing steps for non-overlapping spikes.
  3. Most of the spike sorting approaches use a stand-alone standard spike detection technique (see for example Choi et al. 2006; Obeid and Wolf 2004; Rebrik et al. 1999 for commonly used spike detection techniques), and a separate classification procedure. Neither the shape of the waveforms nor their change over time or their amplitude distribution across the recording channels is taken into account by the spike detection method. This leads to a poor detection performance, in particular when the signal-to-noise ratio (SNR) is low. Further, the spikes are cut and aligned on some feature (e.g., peak position) as a preprocessing to the classification algorithm. However, overlapping spikes, which severely alter the spike waveform, are not identified as such. This leads to wrong alignments and false classifications by the sorting procedure.
  4. There are two general approaches to extracellular recording with electrodes, namely acute and chronic recording methods. In acute recordings, individual electrodes are advanced into tissue at the beginning of each recording session anew, causing a compression of the tissue (Cham et al. 2005). During the experiment the tissue relaxes and the distances between the electrodes and neurons change; an effect called tissue drift (Branchaud et al. 2006). As a consequence, the shape of the measured waveforms and the characteristic of the background noise changes. Sorting algorithms which do not take into account such variations will perform poorly on data from acute recordings.

An approach based on blind source separation (BSS) techniques and addressing primarily problems 1) and 4) was presented in Takahashi et al. (2002), in which independent component analysis (ICA) was applied to multichannel data recorded by tetrodes (4 channels). Later, the method was adopted to data recorded by dodecatrodes (12 channels) (Takahashi and Sakurai 2005). However, both approaches had to deal with several new problems: Amongst others, time delays between the channels were not considered, biologically meaningless independent components had to be discarded manually, and different neuronal signals with similar channel distributions could not be classified correctly. Furthermore, the methods can only be applied to data recorded with certain electrode types (i.e. tetrodes, dodecatrodes). The most severe problem, though, is the fact that the method cannot deal with data containing neuronal activity from a greater number of neurons than recording channels (over-completeness).

In this work, we present a realtime online spike sorting method based on the BSS idea, which explicitly addresses the four issues 1)–4), but also avoids the drawbacks of the method in Takahashi et al. (2002) and Takahashi and Sakurai (2005). In sum, a spike sorting algorithm for multi electrode data, which detects and resolves overlapping spikes with the same computational cost as non-overlapping spikes, is formulated. The method makes optimal use of an arbitrary number of simultaneously recorded channels and can even run on single channel data. Moreover, since spike detection, spike alignment, and spike classification are not separate parts, but are combined into a single algorithm, our method performs well on data with low SNR and containing many overlapping spikes. By incorporating a direct feedback, the algorithm adapts to varying spike shapes and to non-stationary noise characteristics. The algorithm is fully automatic and due to its linear and parallel computation steps it is ideally suited for realtime applications (see Fig. 4 for a summary of our method).

Fig. 4
Schematic illustration of the way data is processed: The data is bandpass filtered and periods containing artifacts are excluded from further analysis (Section 2.7). During the initialization phase a conventional spike detection and clustering ...

This paper is organized as follows: In Section 2 we present our method step by step. First, we briefly introduce linear filters. These filters were used in radar applications (Turin 1960), geophysics (Robinson and Treitel 1980) as well as for spike detection (Thakur et al. 2007; Vollgraf et al. 2005), but to our knowledge have not been applied to spike sorting yet. Moreover, in contrast to those studies, we do not directly apply a threshold to the filter outputs, but consider them as a new representation of the data. In this representation the spike sorting task can be handled as a well defined BSS problem, which we solve with a un-mixing technique we will refer to as “Deconfusion”.

The evaluation of our method is done on two different datasets from real recordings and also on simulated data. The experimental setup, used equipment and the characteristic of recorded data are described in Section 3. The advantages and abilities of the method are demonstrated in Section 4. Evaluations of the spike detection performance are done using data from simultaneous intra- and extracellular recordings made in slices of rat visual cortex, and show that the proposed algorithm is superior to conventional spike detection methods. The noise robustness and the ability to successfully resolve overlapping spikes is evaluated systematically on synthetic data. Finally, the method is applied to data from extracellular recordings made in the prefrontal cortex of awake behaving macaques. This data is particularly challenging, because the tetrodes are not implanted chronically, but inserted before every experiment anew, leading to tissue drifts. We conclude that our method adopts to non-stationarities and also successfully resolves overlapping spikes in real data. A summary and a discussion of further improvements is given in Section 5.

Methods

Glossary of mathematical notation

We use a notation in which symbols for scalar quantities are represented by lower case letters, vectorial quantities are represented by bold lower case letters, and operators or matrices are represented by bold upper case letters. Matrices representing several vectorial quantities, but not linear transformations, are labeled with an additional bar. In Table 1 all important quantities are listed. The corresponding vectorial quantities are defined by concatenating all channel-wise defined vectors. As an example the vectorial template equation M1 of neuron i is given by

equation M2

where the superscript equation M3 means transpose. The vectors equation M4, equation M5, equation M6 are defined in the same way. Analogously, covariance matrices, e.g, the data covariance matrix equation M7, are defined as

equation M8

with equation M9. equation M10 is a symmetric N ·Tf by N ·Tf Toeplitz matrix. Alternatively, it can be expressed as

equation M11
1
Table 1
Definitions of important quantities and their meaning

Generative model

We assume an explicit model for the neuronal data recorded extracellularly. The underlying assumptions are:

  1. Each neuron generates a unique spike waveform equation M32 (called template), which is constant over a time period of length T.
  2. All time series equation M33 of spike times of neuron i (called spike trains) are statistically independent of the noise equation M34. Furthermore, these quantities sum up linearly.
  3. The noise statistic is entirely captured by a covariance matrix equation M35.

As discussed extensively in Pouzat et al. (2002), these assumptions are reasonable and are used explicitly or implicitly in most spike sorting techniques. Consequently the measured data equation M36 can be expressed as

equation M37
2

The measured data are a convolution of the mean waveforms with the corresponding intrinsic spike trains corrupted by colored Gaussian noise (see also Fig. 1(a)–(c)).

Fig. 1
Sketch of the generative model (ac) and the processing stages of the algorithm (de). (a) Spike trains of two neurons. (b) Simulated waveforms of each neuron on a hypothetical multi electrode (two recording channels, without noise). ( ...

Calculation of linear filters

Spike sorting is achieved when the intrinsic spike trains equation M38 are reconstructed from the measured data equation M39. Since, according to the model assumptions, the data were generated by a convolution of intrinsic spike trains with fixed waveforms, the most straightforward procedure would be to apply a deconvolution on equation M40 in order to retrieve equation M41. For an exact deconvolution a filter with an infinite impulse response is necessary. In general, such a filter is not stable and would amplify noise (Robinson and Treitel 1980). Nevertheless, a noise robust approximation for an exact deconvolution can be achieved with finite impulse response filters, to which we will refer as linear filter.

Let us briefly summarize the idea of these filters: The goal is to construct a set of filters equation M42 such that each filter equation M43 has a well defined response of 1 to its matching template equation M44 at shift 0 (i.e. equation M45), but minimal response to the rest of the data. This means that the spikes of neuron i are the signal for filter equation M46 to detect but will be treated as noise by filter equation M47.

Incorporating these conditions leads to a constrained optimization problem

equation M48
3

to which the solution are the desired filters (see Appendix A for a more detailed derivation). A major advantage is the fact that the mentioned optimization problem can be solved analytically. In particular, the filters are given by the following expression:

equation M49
4

where equation M50 is the data covariance matrix defined in Section 2.1. Linear filters maximize the signal-to-noise ratio and minimize the sum of false negative and false positive detections, and are, therefore, optimal in this sense (Melvin 2004).

Filtering the data

Once the filters are calculated, they are cross-correlated with the measured signal, i.e. equation M51. Note that we do not have to pre-process the data with a whitening filter, but the filters can be applied directly to equation M52. This is because the noise statistics is already captured in the matrix equation M53.

From a different point of view, the filtering just changes the representation of the templates. While in the original space the template i was represented by equation M54, its representation in the filter output space is given by the vectors equation M55, j = 1,...,M, where equation M56, see also Fig. 2. This interpretation of filtering will be useful in the next section.

Fig. 2
This figure exemplary illustrates the representation of the templates in the filter output space and the calculation of the Deconfusion parameters. In this example, three templates (equation M57, top row of the figure) originating from tetrode recordings are used. ...

Deconfusion

The linear filters derived in Section 2.3 should suppress all signal components except their corresponding template with zero shift. Thus, the filter response to all templates (and their shifted variants) has to be minimal. This already leads to equation M69 minimization constraints; a number which is normally greater than the number of free variables of a filter which is Tf ·N. In addition, if the SNR is low, the noise covariance matrix equation M70 dominates Eq. (1).

The lower the SNR, the less spikes from other neurons a filter will suppress. Thresholding of every filter output equation M71 individually will, thus, lead to false positive detections. The idea is to de-correlated the filter output in order to achieve an improved spike detection and classification.

We have seen in the previous section that each template equation M72 can be represented in the filter output by M vectors equation M73, j = 1,...,M. Since the detection and classification of the spikes is based on the detection of high positive peak values in the filter output (by construction), all values below zero in the filter output are irrelevant, and thus, can be discarded. As a result, we ignore all values below zero by applying a half-wave rectification I(x) to the filter output equation M74, where

equation M75
5

The next step is to consider equation M76 as a linear mixture of different sources, where every source is the intrinsic spike train equation M77 of a neuron. Since there are as many filters as neurons, the dimension of the filter output space is equal to the number of neurons, and therefore, the detection and classification problem can be considered as a complete BSS problem. However, it is not guaranteed that the maximal response of filter equation M78 to spikes from neuron j will be at a shift of 0, i.e., when the filter and the template overlap entirely. This leads to the following model for the rectified filter output:

equation M79
6

with equation M80 being the mixture matrix, and τi,j being the shifts between the maximal response of filter equation M81 to template equation M82; i.e.,

equation M83
7

where equation M84 and τi,i = 0 ∀ i by construction. We want to reconstruct the sources equation M85 by solving the corresponding inverse problem:

equation M86
8

with equation M87. Here, the relation to ICA becomes clear, since this is a similar inverse problem ICA solves. In contrast to ICA, we do not have to estimate equation M88 and τi,j from the data, but can calculate them directly from the responses (i.e. cross-correlation functions) of all filters to all templates, as illustrated in Fig. 2.

All steps of these procedure are summarized under the term “Deconfusion” (see also Fig. 1(d)–(e) for a schematic illustration). After Deconfusion the false responses of the filters to non-matching templates are suppressed (see Fig. 3). In principle, it is possible that the inverse problem in Eq. (8) is not exactly solvable, if the shifts are not consistent. Consistent shifts have to satisfy the following equation:

equation M89
9

A derivation is given in Appendix B. For arbitrary templates and data covariance structures, Eq. (9) can in principle be violated. However, with templates from real experiments we did not observe this to be a problem.

Fig. 3
The figure shows the effect of Deconfusion on the filter outputs. The input for Deconfusion were the filter responses equation M90, i,j = 1,2,3 shown in Fig. 2. After Deconfusion the signal of neuron i is mainly present on the output channel ...

Spike detection and classification

In the final step, thresholding is applied to every row i of equation M92. Again, by construction we have only to consider positive peaks. All local maxima after a threshold crossing are identified as spiking times of neuron i. In this sense, spike detection and spike classification is performed simultaneously.

The threshold is set for each row of equation M93 individually such that the total error of false negative and false positive detections is minimal. Amongst others, the threshold depends on the variance of the noise, on the Deconfusion output, and on the firing frequencies of the neurons. A detailed derivation is given in Appendix C.

Artifact detection

Artifacts were removed from our data in two ways. First, all periods during which the animal had to perform a physical task (e.g., pressing a button) were not considered for further analysis. Secondly, for each period of length 10 ms the number of zero-crossings on each data channel was counted and summed up. All periods, in which this number was below 10% of the maximal number of possible zero crossings, were not considered for further analysis. This second type of heuristic removal aims at eliminating artifacts caused by oscillations of the electrode shaft inside the guiding tube (e.g., caused by movement of the animal).

Noise estimation

The noise covariance matrix equation M94 is determined by calculating the auto- and cross correlation functions of every channel. Only data points which were not part of any spike nor any artifact period, were used for the calculation. The noise covariance matrix is needed for the initialization phase, see Section 2.10, and for evaluation of the sorting result on real data, see Section 4.2.3.

Adaptation

Due to tissue relaxations the measured waveforms change over time as the relative distance between the multi electrode and the neurons change. In order to track these changes we re-estimate the templates as well as the data covariance matrix after every time period of length T. Each template equation M95 is re-estimated as the mean of the last 350 spikes (see Section 5 for a discussion of this value) detected from neuron i; whereas the spikes of neuron i are aligned on the maximal peak of the response of filter equation M96. For the re-estimation only spikes which were classified by our method as non-overlapping spikes are used. The data covariance matrix is re-estimated from the last 30 s of the recordings and the linear filters are re-calculated. Consequently, the Deconfusion and the thresholds are re-computed as well. In Section 4.2.3 we show that we can indeed track drifts with this approach.

Templates whose SNR decreases over time might be a concern. By constantly adapting the template, finally, there is a risk of getting a template which is very close to the noise signature, and the corresponding filter will detect pure noise. This can be prevented by removing filters at the appropriate moment. Consequently, we stop tracking templates whose SNR drops below 0.65. This value proved to be appropriate during simulations (see Section 4.2.2).

Initialization phase

Most of the analysis done in the precedent sections was based on the assumption of known initial templates. Hence, before applying our method, one needs an initialization phase during which the templates are found. In principle, any supervised or unsupervised learning method can be applied.

We want to emphasize that the initialization phase is only necessary at the beginning of a recording session (Fig. 4): Once the initial templates are estimated, the main algorithm runs online. Furthermore, because of the feedback described in Section 2.9, the initialization does not have to be very accurate, as the templates are re-estimated after every period of length T. Usually we used an initialization phase of about 30 s in our real recordings (Section 3.3). This time window is short enough so that the templates change only very slightly in time and can, therefore, be clustered reliably, but long enough to acquire enough spikes to estimate robustly the mean waveforms.

Initial spike detection and initial spike alignment

During the initialization phase spike detection can be done with any conventional technique. We used an energy based approach, since it usually delivers a better performance than other methods (Mtetwa and Smith 2006; Obeid and Wolf 2004).

In particular, we applied the MTEO detector (see Section 4.1 for definition) with k-values [1,3,5] to each recording channel separately and set the threshold to 3.5 times the median of its output. Spike periods were defined as intervals of length 1.5 ms, in which the output of the MTEO detector exceeded the threshold value at least once.

Correct spike alignment is crucial for a good clustering result. While in many studies an alignment based on the maximal and/or minimal peak value of a spike is used, again, methods based on the energy of a spike usually yield better results (Fee et al. 1996a). After cutting out all spikes around the peak of the detector, we used the following algorithm for alignment:

  1. Calculate the average template over all spikes
  2. Minimize the energy difference between every spike and the template by shifting the spikes
  3. Repeat until convergence or a maximal number of iterations is reached

In our experiments described in Section 3.3 the average number of spikes in the first 30 s of recordings is around 2500 and convergence is obtained after 15 to 20 iterations.

Initial clustering

Although a broad range of sophisticated clustering algorithms is available, we used a standard approach, since a very accurate initialization is not crucial for our method. The aligned spikes are whitened (e.g., see Pouzat et al. 2002) and projected into the space of the first 6 principle components. The clustering consists of a Gaussian mixture model in combination with the Expectation-Maximization algorithm (Xu and Wunsch 2005). For every number of cluster means between 1 and 15 the clustering procedure is executed 3 times with random initial means. The covariance matrices are fixed to 2.5 times the identity matrix. The run and the number of means with the highest score according to the Bayesian inference criterion (Xu and Wunsch 2005) are selected as initialization for the main algorithm.

Signal-to-noise ratio (SNR)

The SNR is a scalar value which is an indicator for the difficulty of detecting a signal in noisy data. In this sense, the SNR definition should be dependent on the method used for signal detection. Several definitions of the SNR are used in the spike sorting literature. A very common one is to define the SNR by some maximal value, e.g., the maximal amplitude, the maximal difference in amplitudes (peak to peak distance), or the maximum of the absolute value of the amplitude, divided by the variance of noise σ2, i.e.,

equation M97

(e.g. see Choi et al. 2006). Another current definition for the SNR is based on the energy of a signal, i.e.,

equation M98

(e.g. see Rutishauser et al. 2006). We introduce a definition of SNR which is based on the Mahalanobis distance of a template equation M99 to zero:

equation M100
10

In the special case of single electrode data and of 1-dimensional templates (Tf = 1), all SNR definitions are equivalent. To show that equation M101 is an appropriate SNR definition for the linear filters, while the other definitions are in contradiction with the meaning of signal-to-noise ratio, we simulated datasets containing a single neuron, which fired according to a Poisson statistic, and a noise covariance matrix equation M102, where equation M103 denotes the identity matrix, and equation M104 is a noise covariance matrix from one of the experiments described in Section 3.1, with equation M105 for all i. The used template was extracted from the same experiment. We simulated datasets for ten different α values between 0 and 1. The equation M106 decreased with increasing α, and consistently the detection performance of our method decreased, see Fig. 5. Note that equation M107 for all α values, which means that those definitions are inappropriate for the proposed method. Nevertheless, we always provide values for all three definitions of SNR in order to allow comparisons with other publications.

Fig. 5
(a) Template equation M108 (in arbitrary units) used for the simulations. (b) Noise autocorrelation function of the same experiment from which the template was extracted. This autocorrelation was used to calculate equation M109. (c) Plot of equation M110, of equation M111 and of equation M112 in dependence of α ...

Experiments and datasets

For the performance evaluation of our method, three different datasets were used. All experiments were performed in accordance with German law for the protection of experimental animals, approved by the local authorities (“Regierungspräsidium”), and are in full compliance with the guidelines of the European Community (EUVD 86/609/EEC) for the care and use of laboratory animals.

Simultaneous intra/extra-cellular recordings

The experiments were done in acute brain slices from Long Evans rats (P17–P25). In every experiment a pyramidal cell from visual cortex, Layer 3 or 5 depending on the experiment, was simultaneously recorded intracellularly and extracellularly. Extracellular spike waveforms were recorded using a 4-core-Multifiber Electrode (Tetrode) from Thomas RECORDING GmbH, Germany. The cell was intracellularly stimulated by a current injection (varying from experiment to experiment between 80 pA and 350 pA). Extracellular recordings were sampled at 28 kHz and filtered with a bandpass FIR filter (300 Hz to 5000 Hz).

The intracellularly recorded spikes were detected using a manually set threshold on the membrane potential. The threshold crossings in the membrane potential were used as triggers to cut out periods from the extracellular recordings (2 ms before and 5 ms after the trigger). In total, data was recorded from 6 different cells, which resulted in 9957 intracellularly detected spikes. For analysis only the recording channel with the highest SNR was considered. The SNR of the different experiments varied from equation M114 (equation M115, equation M116) to equation M117 (equation M118, equation M119). A short period of recordings with a moderate SNR (equation M120, equation M121, equation M122) is shown in Fig. 6, top row.

Fig. 6
A short piece (approx. 400 ms) of extracellularly recorded data from slices of rat visual cortex is processed with different spike detection techniques (rows 3–6). Data were recorded simultaneously intracellularly (first row) and extracellularly ...

Simulated data

The artificially generated data simulates a single channel recording of 15 s length at a sample frequency of 32 kHz containing activity from three neurons. Every dataset contained exactly 750 equidistantly distributed spikes of every neuron, which corresponds to a firing frequency of 50 Hz. The three used templates were extracted from the recordings described in Section 3.1 and had a length of 2.1 ms. The noise was generated by an ARMA model (Hayes 1996) approximating the noise characteristic shown in Fig. 5(b).

Dataset with overlapping spikes

The relative number of overlapping spikes was systematically varied from 1% up to 50%. 75% of all overlapping spikes consist of overlaps between two templates (25% for each combination), and 25% of all overlapping spikes consist of overlaps between all three templates. The amount of overlap, i.e., how much the templates overlap, is distributed according to a uniform distribution on the interval [1/3, 1]. The SNR was kept constant for all overlapping ratios, namely, all three templates were scaled to an equal SNR, which was equation M123. This corresponds to equation M124 and equation M125 (average values over the three templates).

Dataset with SNR variation

The equation M126 was systematically varied from 0.6 to 1.4 (which is equivalent to 2.71 to 6.32 average equation M127 and 1.06 to 2.48 average equation M128). The amount of overlapping spikes was constant and set to 7%, which is approximatively the overlap ratio resulting by chance under the assumption of independent spike trains.

The over-completeness, the equal SNR of all three templates, and the presence of overlapping spikes make these datasets particularly challenging.

Acute recordings

Tetrodes were placed in ventral prefrontal cortex for individual recording sessions, sampling data from the same region across experiments. Recordings were performed simultaneously from up to 16 adjacent sites with an array of individually movable fiber micro-tetrodes (Eckhorn and Thomas 1993). Recording positions of individual tetrodes were manually chosen to maximize the recorded activity and the signal quality. Data were sampled at 32 kHz and bandpass filtered between 0.5 kHz and 10 kHz.

Neuronal activity was recorded while 2 macaque monkeys performed a visual short-term memory task. The task required the monkeys to compare a test stimulus to a sample stimulus presented after a 3 s long delay and to decide by differential button press whether both stimuli were the same or not. Stimuli consisted of 20 different pictures of fruits and vegetables which were presented for 0.5 s (test stimulus) or for 2 s (sample stimulus). Correct responses were rewarded. Match and non-match trials were randomly presented with an equal probability. This experimental setup was presented in Wu et al. (2008).

Approximately, the monkeys perform 2000 trials per session, which is equivalent to almost 4 h of recording time. For the evaluation of our algorithm only the first 5 s of every trial were processed, as the remaining data might contain severe artifacts caused by the monkey’s movement.

Results and discussion

The performance of a spike sorting method depends on its capability to detect spikes and to assign every spike to a putative neuron. As described in Section 2.6, our method achieves both simultaneously. We evaluated the performance of our approach, first, as a pure detection method, and then, as a combined detection and classification technique. In both categories we compared it against techniques commonly used.

Spike detection performance

The evaluation was done on the in-vitro dataset described in Section 3.1. Although the extracellular signal was recorded with a tetrode, we used only one recording channel for further analysis, since most conventional spike detection methods are only defined for single channel data. The detectors used are:

  1. Mahalanobis distance: This method is described in Rebrik et al. (1999). In brief, periods having a greater Mahalanobis distance to zero than a certain threshold are identified as spikes. The noise covariance matrix was estimated from data pieces in which the neuron was not stimulated. The size of the matrix was chosen to match the observed length of spikes in the experiment and was then applied window-wise. Local maxima crossing the threshold are identified as spike times.
  2. Squaring: The raw data is squared and normalized. Local maxima crossing the threshold are identified as spike times. In case of an one-dimensional noise covariance matrix, this method is equivalent to the method “Mahalanobis distance”.
  3. Squaring smoothed: A Savitzky-Golay filter of span 5 and order 2 is additionally applied to the output of the method “squaring”. This method is very similar to the one used in Rutishauser et al. (2006).
  4. MTEO: This method is described in Choi et al. (2006). In brief, the data is smoothed with a Hamming window and a quantity (which depends on parameters k) related to the energy of the signal is computed. We used two parameter sets for this method, one with k-values of equation M129 and one with k-values of equation M130.
  5. Optimal filter: Since the occurrence of the spikes is known (due to the intra-cellular recording), the optimal filter is calculated using the average waveform of all spikes of the recorded neuron.
  6. Our method: In the case of a single neuron, our spike sorting method corresponds to a single “estimated filter” detector, i.e., the initial filter is calculated using the average waveform of all spikes found by the equation M131 with a threshold set to 3.5 times the median of its output.

A short piece of the recordings and some of the corresponding detector outputs are shown in Fig. 6.

We compared the performance of the different spike detection methods using receiver operating characteristic (ROC) curves. For every detector the threshold is systematically varied between 0, resulting in zero false negative detections (FN), and the minimal value which does not detect any spikes; i.e., zero true positive detections (TP). For every threshold the percentage of TP is plotted against the false positive (FP) rate. Such a curve is shown for one exemplary experiment in Fig. 7. The curve for the best possible detector (i.e. no FP, but 100% TP detections) would pass through the point (0,100). The area under such a curve (AUC) is, thus, a measure for the performance of a detector. The normalized AUC values for the area up to 30 Hz of FPs of all detectors averaged over all available datasets are shown in Table 2. Although only the average performance is presented, our method and the optimal linear filter also achieved higher scores on every individual dataset described in Section 3.1. In all experiments the optimal filter was superior to the other detectors, while our method scored second with a very similar performance. This shows that taking into account the full waveforms as well as the data statistic always greatly improves the detection performance. The optimal linear filter was included into the evaluation to provide an upper bound on the performance one can achieve with our method. Our method offers another advantage for the detection of spikes, namely a bigger robustness to threshold variations, see Fig. 8. This means that a deviation from the optimal threshold has a less drastic impact on the total error (FP + FN) than for the other methods.

Fig. 7
ROC curves for different spike detection methods. In this experiment equation M132, equation M133, equation M134
Table 2
Average normalized area under the curve for each evaluated spike detection method
Fig. 8
Number of errors for each spike detection method in respect to varying thresholds. The color coding is the same as in Fig. 7. For each detection method the maximal threshold was determined and normalized to 1. The maximal threshold is defined ...

Spike sorting performance

Resolution of overlapping spikes

We recall that the applied operations to the recorded data could be summarized in Eq. (8). The cross-correlation between the filters and the data is a linear operation. The following Deconfusion consists of a half-wave rectification, which is a non-linear operation, but affects only noise and not the action potentials (represented in the filter output), and the un-mixing, which is linear again. Hence, one can expect that if the superposition of spike waveforms is also linear, overlaps should be resolved successfully. We validated this assumption on the dataset described in Section 3.2.1. The algorithm was executed in the same way as described in Section 2. In order to allow the method to adapt (Section 2.9), the method was iterated 5 times on the same dataset. We also compared the performance of our method to those of two popular clustering based offline methods, one of them being the method described in Section 2.10.2, which will be abbreviated as “GMM”. Since this is also the method which is used for initialization of our algorithm, the comparison with GMM directly provides information about the improvements in sorting when our method is used.

The other algorithm, called “KlustaKwik”, was explicitly developed for clustering neuronal data and was first introduced in Harris et al. (2000). The clustering parameters were set to their default values. Spike detection and alignment was done in the same way as described in Section 2.10.1. To provide an upper bound on the performance our approach could achieve, we included the evaluation with the optimal filters calculated directly from the real templates. Note that other existing, purely clustering-based sorting methods, either in the PCA space or in the original data space, would perform similarly to GMM and KlustaKwik.

For the evaluation the relative number of TP was counted (Tables 3 and and4).4).

Table 3
Average performance of the proposed method for non overlapping and overlapping spikes
Table 4
Same evaluation as in Table 3, but for the method “GMM” described in Section 2.10

The simulations show that our method indeed resolves overlapping spikes and outperforms the clustering based methods; see Fig. 9. Our method works even for datasets with a large amount of overlapping spikes, and the performance is close to the theoretical bound of this approach. On the other hand, the performance of the purely clustering based methods rapidly decreases with an increasing amount of overlapping spikes. Overlapping spikes are mostly detected as single events by conventional spike detection techniques, which leads to a high FN rate. Furthermore, since the waveforms of overlapping spikes are distorted, their distances to the corresponding cluster means are large, making it difficult to assign them to a neuron. This results in a low TP score for clustering based methods.

Fig. 9
Average performance of the different spike sorting methods over 10 simulations. The x-axis indicates the overlap ratio, i.e. the relative number of overlapping spikes (see Section 3.2) while the y-axis represents the correct classifications in ...

Performance for various SNR

The evaluation on the dataset with a varying SNR (see Section 3.2.2) was done in the same way as in the previous section. The results are shown in Fig. 10. The performance of the clustering based methods is severely affected by a low SNR. The performance of the proposed method follows the one of the GMM algorithm, since it relies on its output for initialization. Nevertheless, our method is always superior to it. Because of the rapid decrease in performance from a SNR level of 0.7 to an SNR level of 0.6, we stopped the algorithm from detecting spikes for templates with a lower SNR than 0.65 in real recordings by deleting the corresponding templates and filters. In contrast, the optimal filter method is only slightly affected by a low SNR level, indicating that a more elaborate initialization would increase the performance of the proposed method on datasets with very low SNRs.

Fig. 10
Average performance of the different spike sorting methods over 10 simulations in respect to various SNR levels. Note that the performance of the proposed method degrades with the performance of the GMM algorithm. This is because the output of the GMM ...

Performance on experimental data

We applied our method to data recorded in the prefrontal cortex of monkeys performing a short-term memory task as described in Section 3.3. For illustrative purposes, we show the results obtained by processing data from one tetrode, since the qualitative outcomes from processing other tetrodes and different recording sessions are similar.

For the initialization phase we used the first 7 trials of the recording. The initial spike detection and clustering was done as described in Section 2.10, resulting in a total of 3219 detected spikes, which were assigned to 8 clusters. This basic clustering was used as an initialization for the main algorithm, which was executed in the same way and with the same parameters as described in Section 2 (see also Fig. 4 for a summarization). The 7 trials used for initialization were also processed with the main method in order to improve the sorting quality.

The templates after the first 90 trials are shown in Fig. 11, and seem to be reasonable by visual inspection of an expert. In total, our method found almost 200000 spikes (57111, 18060, 50724, 51709, 3974, 7057, 444, 10915 for each template). Two well-established tests to quantitatively asses the sorting quality of a method performing on real data are the inter spike interval distribution and the projection test (Rutishauser et al. 2006; Pouzat et al. 2002); the evaluation of our sorting with both tests is shown in Fig. 11. The relative number of spikes during the first 3 ms is smaller than 1.5% for all neurons, implying that the refractory period is respected. On the other hand, the projection test verifies that the spikes of a single neuron have not been artificially split by the sorting algorithm into multiple clusters or that spikes from multiple neurons are assigned to the same cluster. The sorting of our method also passes the projection test since the cluster distributions do not overlap and are close to the theoretical prediction of a normal distribution with a variance 1. In sum, the good results of these two tests imply that the found clusters are well separated and indeed correspond to single neurons, as well as that the assumptions made in Section 2.2 are justified.

Fig. 11
(a) Plot of the concatenated templates and their standard deviation. For the averaging all detected spikes from trial 50 to trial 90 were used. The vertical lines indicate the concatenation points, while the colored dots on the right serve as a label. ...

Since we inserted the tetrodes before every experiment anew, our algorithm has to deal with the variability in the data caused by tissue drifts. The adaption procedure described in Section 2.9 was executed after every trial and adapted the algorithm correspondingly. The time period over which the templates were assumed to be constant was set to T = 5 s.1 As a result, 2 neurons could be tracked from the beginning to the very end of the experiment, see Fig. 12. The other templates were deleted earlier, since their equation M142 dropped below 0.65. The importance of taking temporal variations for sorting into account is demonstrated in Fig. 13. If the drift is not accounted for, the clusters are elongated and their spread is larger, making any classification more difficult.

Fig. 12
Effects of tissue drift on the amplitudes of all templates over the whole experiment. For every recording channel of the tetrode and for every trial a peak amplitude histogram was calculated. In every trial the number of local extrema of 50 samples width ...
Fig. 13
Effects of tissue drift on templates and cluster distributions in the PCA space. (a) For three filters which detected spikes nearly throughout the whole experiment, the corresponding templates of the initialization and at trials 50, 1000 and 2001 are ...

The disappearance of neurons from the recording volume is a common phenomenon in our recordings. However, the opposite, i.e., the appearance of new neurons during recordings, is rarely observed. This might be explained by the fact, that at the beginning of the experiments, the tetrodes are explicitly placed at a position where a lot of neuronal activity is measured. Therefore, it is more probable that during the tissue drifts the high activity population of neurons disappears than that new, highly active neurons appear. We discuss this problem also in Section 4.4.

In Section 4.2.1 we have already demonstrated on simulated data the ability of our method to resolve overlapping spikes instantaneously. This is also the case for real data, see Fig. 14. The same figure also shows, that it would be very difficult to classify correctly these overlapping spikes with a purely clustering based algorithm.

Fig. 14
Ability of our algorithm to resolve overlapping spikes in real data. (a): Projection of all detected and whitened spikes from the trials 50–90 into the space of the first two principle components (the solid bars correspond to 3 standard deviations ...

The evaluation in Fig. 11 and Fig. 13 shows that the clustered spikes, although whitened, are not perfectly Gaussian distributed. This deviation is caused by overlapping spikes, but it is also due to an intrinsic waveform variability, as it is observed for example during bursts (Fee et al. 1996b). In this sense, the generative model assumed in Section 2.2 is not strictly valid anymore. Nevertheless, our method achieves a good performance, even for datasets containing bursting neurons identified by visual inspection. This can be explained by the fact that the scaling of the waveform during burst is close to linear (Rutishauser et al. 2006). Because of the linear character of our method (e.g. see Section 4.2.1), the response to a linearly scaled waveform will also only be scaled by the same factor. Hence, the algorithm classifies spikes from bursting neurons correctly as long as the amplitude degradation of the spikes is not too strong.

Limitations of our method

We have shown that our method is of great potential for spike detection and classification applications. However, there is a principle limitation: Since the filtering and the Deconfusion are linear operations, it is impossible to discriminate waveforms which are strictly linear dependent, i.e., when the spike waveform of one neuron is a multiple of the waveform of another neuron. A possible way to solve this problem is to sort the templates according to their SNR. Spikes with the highest SNR are detected first. Whenever a spike is found, the corresponding template is subtracted from the data and all other filter outputs are re-calculated for the affected period. This procedure is repeated for templates with a lower SNR. Further, if the sum of the waveforms of two different neurons with a certain shift is nearly identical to another neurons spike waveform, it is impossible to judge whether a spike is an overlap or not. Only probabilistic methods or soft clustering could give a hint at where the waveform came from.

Newly appearing neurons

We have not addressed the problem of neurons which are not detected during the initialization phase. As we observe spikes from neurons whose SNR decreases due to tissue drifts, and finally disappear completely from the recorded data, the opposite might also happen; i.e., neurons, previously undetected, slowly appear in the recording volume. A possible solution would be to run a conventional spike detection method in parallel to our method. All spikes detected by the conventional spike detection technique, but not by our method, could be collected, aligned and clustered. Respecting the newly found clusters, corresponding filters could be initialized and the Deconfusion procedure adapted accordingly.

Implementation and computational complexity

Especially for a real-time implementation the runtime of an algorithm is crucial. After the initialization phase, the proposed method consists mainly of linear operations. The adaptation of the covariance matrix, of the templates and of the Deconfusion parameters need only to be computed every few seconds. Therefore, the computational burden lies in the application of the linear filters and the Deconfusion to a new sample of recorded (multichannel) data. The current implementation was done in Matlab, however the source code is not ready for publication yet. We will make the method available e.g. on ModelDB as soon as the implementation is finished.

Runtime analysis

If a new multichannel sample of data is recorded, first the cross-correlation between the filters and the data has to be calculated and afterwards Deconfusion is applied. The number of operations needed for the cross-correlation of a filter (the number of filters equals the number of neurons M) and the data is directly proportional to the product of the length of the filter Tf and the number of recording channels N. The Deconfusion procedure consists of a half-wave rectification, which is just a sample wise trivial non-linearity, and a matrix-vector multiplication between the square matrix W of dimension M×M and the shifted and half-wave rectified filter outputs. To sum up, the computational complexity for a newly arriving data sample is equation M143. Since we can assume the number of filters to be higher than the number of recording channels, the resulting complexity is equation M144. This means the runtime complexity mainly depends on the number of filters and the filter length.2

Parallel computing

It is important to note that the cross-correlation for every filter—even for every channel of every filter—are independent of each other and can, thus, be computed in parallel as simple vector-matrix multiplications. For a so called vector processor such a multiplication would be one single operation only. E.g this could be implemented on a modern consumer computer-graphics hardware or on programmable digital signal processors.

Conclusion and outlook

An automatic method for simultaneous spike detection and spike classification was presented, having several advantages which were demonstrated on various datasets. Explicitly, the method makes use of the additional information provided by multi electrodes and has no constraints concerning the number of recording channels or the number of neurons present in the data. It resolves overlapping spikes instantaneously, performs well on datasets with a low SNR, and it adapts to non-stationarities present in the data. Moreover, the method operates online and is well suited for a realtime implementation.

In the first step of our algorithm, optimal linear filters were used to enhance the SNR. Linear filters, being an approximation to an exact deconvolution, account for the noise statistics as well as for the full, multi-channel template, and are, therefore, superior to other methods in detecting spikes of a specific neuron. An evaluation on simultaneous intra/extra-cellular recordings in slices of rat visual cortex and on realistic synthetic data shows that the difference in performance is considerable.

Further, we used the output of the linear filters as a new representation of the data. The advantage of the filter output space is that its dimension is equal to the number of neurons, whereas this was not the case in the original data space. This allowed us to treat the spike sorting problem as a well defined source separation problem and solve it by Deconfusion.

In the final step, a channel specific threshold was applied providing simultaneous spike detection and classification. Unlike in many other methods, the thresholds need not to be set manually by a human supervisor but are determined automatically in an optimal way. The advantage of a combined spike detection and classification, in contrast to existing spike sorting methods, was demonstrated on simulated datasets. Especially in the presence of overlapping spike and low SNR, our method achieved better performances. We showed that, in the case of linear filters, a proper definition of the signal-to-noise ratio is based on the Mahalanobis distance, whereas other commonly used definitions do not reflect the difficulty in detecting the signal.

By iteratively updating all quantities, namely the linear filters, the Deconfusion parameters, and the thresholds, the algorithm adopts to non-stationarities present in the data. As such, the method is also suitable for recordings made in acute experiments in which the multi electrodes are inserted each time anew. The number of spikes detected by a filter which were used for the calculation of the template, was set manually to a fixed value, equal for all filters. Instead, one could develop a model for the tissue drift and derive an optimal value which depends on the estimated drifting velocity, the firing rate of the neurons, on the SNR, and on the error tolerance. This is the aim of a future study.

Two drawbacks of the proposed method were discussed, namely the incapability to detect newly appearing neurons and the problem of strictly linear dependent templates. However, for both problems a possible solution was sketched. The detailed study and realization of these solutions will be the scope a future study.

By qualitative arguments, systematic runs on realistically simulated data and on real data from awake behaving macaques, we have shown that the algorithm is capable of resolving overlapping spikes; without additional computing time. However, for the acute recordings in awake behaving monkeys we cannot proof that the found solution is correct, since the ground truth is unknown. Only massive simultaneous intra- and extracellular recordings in vivo could be used to asses the quality of the sorting in real experiments. Due to technical limitations, such a dataset is currently not available.

The algorithm mainly consist of linear, independent operations, which can be executed in parallel and implemented in hardware. Therefore, the algorithm can be used for realtime implementations, making it an potential spike sorting method for brain-machine interfaces and for the execution of closed-loop experiments.

Acknowledgements

This research was supported by the Federal Ministry of Education and Research (BMBF) with the grants 01GQ0743 and 01GQ0410. We thank Sven Dähne for technical support.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Appendix A: Derivation of optimal linear filters

Filter equation M145 should respond with a peak to its matching template equation M146, but should have minimal response to the rest of the data. In particular, one demands that the response to the matching template is 1, i.e. equation M147. The response of the filter to the data is equation M148, where equation M149. Using the third assumption of Section 2.2 the response of a filter to equation M150 will be small (and therefore well distinguishable from the peak response of 1 to the matching template) if the variance of the filter output is small, i.e., one has to minimize equation M151. In summary, the constrained minimization problem is stated as

equation M152
11

A short calculation shows that

equation M153
12

Thus, the Lagrangian L of this minimization problem is given by

equation M154
13

where λ is the Lagrange multiplier. Since the objective function is convex in equation M155, there exists a single minimum, which can be found by solving equation M156. In fact, the minimum is attained at

equation M157
14

Often, linear filters are derived in the frequency domain instead, but linear filter defined in the time domain have several advantages, see Vollgraf and Obermayer (2006).

Appendix B: Derivation of Deconfusion

equation M158 can be expressed as a linear combination of the sources υj at shifts τi,j:

equation M159
15

We show that

equation M160
16

with equation M161 is the corresponding inverse problem. By inserting the expression in Eq. (15) into Eq. (16) one obtains

equation M162
17

Hence,

equation M163
18

This is true, if

equation M164
19

Note that this condition is always satisfied for k = i.

Appendix C: Derivation of the optimal threshold

If we assume that the noise in the Deconfusion output is still a mixture of Gaussians (as an approximation for a mixture of truncated Gaussians), it follows for its variance

equation M165
20

where equation M166 are shifted covariance matrices, i.e. taking temporal correlations into account of order Tf + |τk,j  τk,i|.

The optimal threshold for the detection and classification of spikes from neuron k is chosen such that the overlap between the distribution of the spikes from neuron k and the distribution of the other spikes (from neurons j, j = 1,...,M, j  k) is minimal. We assume the distributions to be Gaussian, with means μk,j and variance σj2. The μk,j are given by the maximal response values of filter j to template k after Deconfusion, i.e.

equation M167
21

whereas the variance is given by Eq. (20). One has only to consider the maximal false response and not the whole response, because the refractory period is in general longer than the length of the template. Thus the optimal threshold εk is given by

equation M168

where erfc denotes the complementary error function, and βj is a normalized weight proportional to the firing frequency of neuron j in order to minimize the total error. Note that the threshold must lie in the interval [0,1], hence this minimization problem can be solved numerically with a line search algorithm, for example using the “fminbnd” command of MATLAB.

Footnotes

1The value of T was set to 5 s just for convenience of implementation, since the first 5 s of each trial were processed.

2In principle, the cross-correlation can be calculated with the help of the fast Fourier transform more efficiently. However, this pays off for long data pieces only, and thus would require to buffer the data first, spoiling the real-time idea.

Contributor Information

Felix Franke, ed.nilreb-ut.sc@ff.

Michal Natora, ed.nilreb-ut.sc@arotan.

Clemens Boucsein, ed.grubierf-inu.eigoloib@niescuob.

Matthias H. J. Munk, ed.gpm.negnibeut@knum.saihttam.

Klaus Obermayer, ed.nilreb-ut.sc@ybo.

References

  • Aksenova TI, Chibirova OK, Dryga OA, Tetko IV, Benabid AL, Villa AEP. An unsupervised automatic method for sorting neuronal spike waveforms in awake and freely moving animals. Methods. 2003;30(2):178–187. doi: 10.1016/S1046-2023(03)00079-3. [PubMed] [Cross Ref]
  • Branchaud, E., Burdick, J., & Andersen, R. (2006). An algorithm for autonomous isolation of neurons in extracellular recordings. In Proc. first IEEE/RAS-EMBS international conference on biomedical robotics and biomechatronics BioRob 2006 (pp. 939–945). doi:10.1109/BIOROB.2006.1639212.
  • Cham JG, Branchaud EA, Nenadic Z, Greger B, Andersen RA, Burdick JW. Semi-chronic motorized microdrive and control algorithm for autonomously isolating and maintaining optimal extracellular action potentials. Journal of Neurophysiology. 2005;93(1):570–579. doi: 10.1152/jn.00369.2004. [PubMed] [Cross Ref]
  • Chandra R, Optican LM. Detection, classification, and superposition resolution of action potentials in multiunit single-channel recordings by an on-line real-time neural network. IEEE Transactions on Biomedical Engineering. 1997;44(5):403–412. doi: 10.1109/10.568916. [PubMed] [Cross Ref]
  • Choi JH, Jung HK, Kim T. A new action potential detector using the mteo and its effects on spike sorting systems at low signal-to-noise ratios. IEEE Transactions on Biomedical Engineering. 2006;53(4):738–746. doi: 10.1109/TBME.2006.870239. [PubMed] [Cross Ref]
  • Delescluse M, Pouzat C. Efficient spike-sorting of multi-state neurons using inter-spike intervals information. Journal of Neuroscience Methods. 2006;150(1):16–29. doi: 10.1016/j.jneumeth.2005.05.023. [PubMed] [Cross Ref]
  • Ding W, Yuan J. Spike sorting based on multi-class support vector machine with superposition resolution. Medical and Biological Engineering and Computing. 2008;46(2):139–145. doi: 10.1007/s11517-007-0248-0. [PubMed] [Cross Ref]
  • Eckhorn R, Thomas U. A new method for the insertion of multiple microprobes into neural and muscular tissue, including fiber electrodes, fine wires, needles and microsensors. Journal of Neuroscience Methods. 1993;49(3):175–179. doi: 10.1016/0165-0270(93)90121-7. [PubMed] [Cross Ref]
  • Fee MS, Mitra. PP, Kleinfeld D. Automatic sorting of multiple unit neuronal signals in the presence of anisotropic and non-gaussian variability. Journal of Neuroscience Methods. 1996;69(2):175–188. doi: 10.1016/S0165-0270(96)00050-7. [PubMed] [Cross Ref]
  • Fee MS, Mitra PP, Kleinfeld D. Variability of extracellular spike waveforms of cortical neurons. Journal of Neurophysiology. 1996;76(6):3823–3833. [PubMed]
  • Gray CM, Maldonado PE, Wilson M, McNaughton B. Tetrodes markedly improve the reliability and yield of multiple single-unit isolation from multi-unit recordings in cat striate cortex. Journal of Neuroscience Methods. 1995;63(1–2):43–54. doi: 10.1016/0165-0270(95)00085-2. [PubMed] [Cross Ref]
  • Harris KD, Henze DA, Csicsvari J, Hirase H, Buzsáki G. Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements. Journal of Neurophysiology. 2000;84(1):401–414. [PubMed]
  • Hayes MH. Statistical digital signal processing and modeling. New York: Wiley; 1996.
  • Hulata E, Segev R, Ben-Jacob E. A method for spike sorting and detection based on wavelet packets and shannon’s mutual information. Journal of Neuroscience Methods. 2002;117(1):1–12. doi: 10.1016/S0165-0270(02)00032-8. [PubMed] [Cross Ref]
  • Kim KH, Kim SJ. Method for unsupervised classification of multiunit neural signal recording under low signal-to-noise ratio. IEEE Transactions on Biomedical Engineering. 2003;50(4):421–431. doi: 10.1109/TBME.2003.809503. [PubMed] [Cross Ref]
  • Lewicki M. A review of methods for spike sorting: The detection and classification of neural action potentials. Network: Computation in Neural Systems. 1998;9(4):53–78. doi: 10.1088/0954-898X/9/4/001. [PubMed] [Cross Ref]
  • McGill KC. Optimal resolution of superimposed action potentials. IEEE Transactions on Biomedical Engineering. 2002;49(7):640–650. doi: 10.1109/TBME.2002.1010847. [PubMed] [Cross Ref]
  • Melvin W. A stap overview. IEEE Aerospace and Electronic Systems Magazine. 2004;19(1):19–35. doi: 10.1109/MAES.2004.1263229. [Cross Ref]
  • Mtetwa N, Smith L. Smoothing and thresholding in neuronal spike detection. Neurocomputing. 2006;69(10–12):1366–1370. doi: 10.1016/j.neucom.2005.12.108. [Cross Ref]
  • Obeid I, Wolf PD. Evaluation of spike-detection algorithms for a brain-machine interface application. IEEE Transactions on Biomedical Engineering. 2004;51(6):905–911. doi: 10.1109/TBME.2004.826683. [PubMed] [Cross Ref]
  • O’Keefe J, Recce ML. Phase relationship between hippocampal place units and the eeg theta rhythm. Hippocampus. 1993;3(3):317–330. doi: 10.1002/hipo.450030307. [PubMed] [Cross Ref]
  • Pouzat C, Mazor O, Laurent G. Using noise signature to optimize spike-sorting and to assess neuronal classification quality. Journal of Neuroscience. 2002;122(1):43–57. [PubMed]
  • Pouzat C, Delescluse M, Viot P, Diebolt J. Improved spike-sorting by modeling firing statistics and burst-dependent spike amplitude attenuation: A markov chain monte carlo approach. Journal of Neurophysiology. 2004;91(6):2910–2928. doi: 10.1152/jn.00227.2003. [PubMed] [Cross Ref]
  • Rebrik S, Wright B, Emondi A, Miller KD. Cross-channel correlations in tetrode recordings: Implications for spike-sorting. Neurocomputing. 1999;2627:1033–1038. doi: 10.1016/S0925-2312(99)00101-0. [Cross Ref]
  • Robinson EA, Treitel S. Geophysical signal analysis. Englewood Cliffs: Prentice Hall; 1980.
  • Rutishauser U, Schuman EM, Mamelak AN. Online detection and sorting of extracellularly recorded action potentials in human medial temporal lobe recordings, in vivo. Journal of Neuroscience Methods. 2006;154(1–2):204–224. doi: 10.1016/j.jneumeth.2005.12.033. [PubMed] [Cross Ref]
  • Sakurai Y, Takahashi S. Dynamic synchrony of firing in the monkey prefrontal cortex during working-memory tasks. Journal of Neuroscience. 2006;26(40):10141–10153. doi: 10.1523/JNEUROSCI.2423-06.2006. [PubMed] [Cross Ref]
  • Shoham S, Fellows MR, Normann RA. Robust, automatic spike sorting using mixtures of multivariate t-distributions. Journal of Neuroscience Methods. 2003;127(2):111–122. doi: 10.1016/S0165-0270(03)00120-1. [PubMed] [Cross Ref]
  • Takahashi S, Anzai Y, Sakurai Y. Automatic sorting for multi-neuronal activity recorded with tetrodes in the presence of overlapping spikes. Journal of Neurophysiology. 2003;89(4):2245–2258. doi: 10.1152/jn.00827.2002. [PubMed] [Cross Ref]
  • Takahashi S, Sakurai Y. Real-time and automatic sorting of multi-neuronal activity for sub-millisecond interactions in vivo. Neuroscience. 2005;134(1):301–315. doi: 10.1016/j.neuroscience.2005.03.031. [PubMed] [Cross Ref]
  • Takahashi S, Sakurai Y, Tsukuda M, Anzai Y. Classification of neural activities using independent component analysis. Neurocomputing. 2002;49:289–298. doi: 10.1016/S0925-2312(02)00528-3. [Cross Ref]
  • Thakur PH, Lu H, Hsiao SS, Johnson KO. Automated optimal detection and classification of neural action potentials in extra-cellular recordings. Journal of Neuroscience Methods. 2007;162(1–2):364–376. doi: 10.1016/j.jneumeth.2007.01.023. [PubMed] [Cross Ref]
  • Turin G. An introduction to matched filters. IRE Transactions on Information Theory. 1960;6(3):311–329. doi: 10.1109/TIT.1960.1057571. [Cross Ref]
  • Vargas-Irwin C, Donoghue JP. Automated spike sorting using density grid contour clustering and subtractive waveform decomposition. Journal of Neuroscience Methods. 2007;164(1):1–18. doi: 10.1016/j.jneumeth.2007.03.025. [PMC free article] [PubMed] [Cross Ref]
  • Vollgraf R, Munk M, Obermayer K. Optimal filtering for spike sorting of multi-site electrode recordings. Network. 2005;16(1):85–113. doi: 10.1080/09548980500275378. [PubMed] [Cross Ref]
  • Vollgraf R, Obermayer K. Improved optimal linear filters for the discrimination of multichannel waveform templates for spike-sorting applications. IEEE Signal Processing Letters. 2006;13(3):121–124. doi: 10.1109/LSP.2005.862621. [Cross Ref]
  • Wang GL, Zhou Y, Chen AH, Zhang PM, Liang PJ. A robust method for spike sorting with automatic overlap decomposition. IEEE Transactions on Biomedical Engineering. 2006;53(6):1195–1198. doi: 10.1109/TBME.2006.873397. [PubMed] [Cross Ref]
  • Wu W, Wheeler DW, Staedtler ES, Munk MHJ, Pipa G. Behavioral performance modulates spike field coherence in monkey prefrontal cortex. Neuroreport. 2008;19(2):235–238. doi: 10.1097/WNR.0b013e3282f49b29. [PubMed] [Cross Ref]
  • Xu R, Wunsch ID. Survey of clustering algorithms. IEEE Transactions on Neural Networks. 2005;16(3):645–678. doi: 10.1109/TNN.2005.845141. [PubMed] [Cross Ref]
  • Zhang PM, Wu JY, Zhou Y, Liang PJ, Yuan JQ. Spike sorting based on automatic template reconstruction with a partial solution to the overlapping problem. Journal of Neuroscience Methods. 2004;135(1–2):55–65. doi: 10.1016/j.jneumeth.2003.12.001. [PubMed] [Cross Ref]

Articles from Springer Open Choice are provided here courtesy of Springer