Automatic online spike sorting with singular value decomposition and fuzzy C-mean clustering
© Oliynyk et al.; licensee BioMed Central Ltd. 2012
Received: 1 February 2012
Accepted: 16 July 2012
Published: 8 August 2012
Skip to main content
© Oliynyk et al.; licensee BioMed Central Ltd. 2012
Received: 1 February 2012
Accepted: 16 July 2012
Published: 8 August 2012
Understanding how neurons contribute to perception, motor functions and cognition requires the reliable detection of spiking activity of individual neurons during a number of different experimental conditions. An important problem in computational neuroscience is thus to develop algorithms to automatically detect and sort the spiking activity of individual neurons from extracellular recordings. While many algorithms for spike sorting exist, the problem of accurate and fast online sorting still remains a challenging issue.
Here we present a novel software tool, called FSPS (Fuzzy SPike Sorting), which is designed to optimize: (i) fast and accurate detection, (ii) offline sorting and (iii) online classification of neuronal spikes with very limited or null human intervention. The method is based on a combination of Singular Value Decomposition for fast and highly accurate pre-processing of spike shapes, unsupervised Fuzzy C-mean, high-resolution alignment of extracted spike waveforms, optimal selection of the number of features to retain, automatic identification the number of clusters, and quantitative quality assessment of resulting clusters independent on their size. After being trained on a short testing data stream, the method can reliably perform supervised online classification and monitoring of single neuron activity. The generalized procedure has been implemented in our FSPS spike sorting software (available free for non-commercial academic applications at the address: http://www.spikesorting.com) using LabVIEW (National Instruments, USA). We evaluated the performance of our algorithm both on benchmark simulated datasets with different levels of background noise and on real extracellular recordings from premotor cortex of Macaque monkeys. The results of these tests showed an excellent accuracy in discriminating low-amplitude and overlapping spikes under strong background noise. The performance of our method is competitive with respect to other robust spike sorting algorithms.
This new software provides neuroscience laboratories with a new tool for fast and robust online classification of single neuron activity. This feature could become crucial in situations when online spike detection from multiple electrodes is paramount, such as in human clinical recordings or in brain-computer interfaces.
Electrophysiological recording of single neuron activity represents a fundamental tool for investigating brain functions. Since a recording electrode often picks-up spikes from more than one neuron, a spike sorting technique is needed to identify and separate spikes of different neurons [1, 2]. Most currently available computational procedures provide accurate sorting and classification, but are often highly interactive and time-consuming and require specific experience and subjective judgments. Fast automatic methods are available [3–5], but they are usually not as accurate as the offline ones and often suffer from problems such as false match or double match, spike overlap and errors in classification [1, 6]. Solving the tradeoff between automation, speed and accuracy of spike sorting is thus a crucial challenge in systems neuroscience.
Here, we aim at contributing to the progress of the field by achieving accurate, fast and fully automated spike sorting. To this purpose, we present a new method (and a software package) based on the Fuzzy C-mean (FCM) classification of spike waveforms in the low-dimensional feature space of Principal Components (PCs). Many currently used offline spike sorting algorithms already employ Principal Component Analysis (PCA) as a preprocessing step to compress the dimensionality of the patterns to be clustered [1, 7]. However, its practical application for fast and automatic separation of neurons is still limited for many reasons. In fact, it is commonly reported that the performance of the PCA heavily depends on the accuracy of the waveform alignment , thus requiring strong computations and human supervision of the clustering results . In addition, eigenvectors accounting for the largest variance of the data are not necessarily providing the best separation of the spike classes . Finally, it is often pointed out that PCA, in its basic configuration, is a static technique, not suitable for monitoring of non-stationary behaviour , while in vivo single unit activity represents mostly non-stationary system with nontrivial dynamics .
The article is organized as follows. First, we describe the spike sorting method used in our approach. We then evaluate its accuracy using simulated spike datasets at different background noise levels proposed by Quiroga et. al. (2004). Then we demonstrate the performance of our method by applying it to the analysis of real extracellular recordings from the macaque premotor cortex, and we investigate the robustness of the algorithm to sample size and inhomogeneities between cluster size to deal with the non-stationarity of the data during the recording session. Finally, we illustrate how our method of spike sorting can be implemented to monitor online the activity of single neurons during electrophysiological recordings.
Simulated extracellular recordings were used to test the spike sorting procedures. These simulated data were the ones used in  and are available online at http://vis.caltech.edu/~rodri/Wave_clus/Wave_clus_home.htm. Briefly, simulated signals consisted of spike shapes of three neurons compiled from recordings in the neocortex and basal ganglia. For generating background noise, spikes randomly selected from the database were superimposed at random times and amplitudes. Next, a train of three distinct spike shapes was superimposed on the noise signal at random times. The amplitude of the three spike classes was normalized to have a peak value of 1. The noise level was determined from its standard deviation, which was equal to 0.05, 0.1, 0.15, and 0.2 relative to the amplitude of the spike classes. There were four different example simulations, number from one to four in order of increasing sorting difficulty (see Ref ).
Electrophysiological recordings were made from freely behaving, partially restrained, macaque monkey (Macaca fascicularis). All experimental protocols were approved by the Veterinarian Animal Care and Use Committee of the University of Ferrara, by the Italian Ministry of Health and complied with the European laws on the use of laboratory animals. The surgical procedure were the same as previously described . Multi-unit recordings were performed by using varnished (Sivamid 595, ELANTAS Deatech S.r.l., Italy) tungsten microelectrodes with impedance 0.15–1.5MΩ (measured at 1 kHz), slowly inserted in the cortex by a hydraulic microdrive (Kopf Instruments, CA, USA; step resolution, 10 μm). Recorded signal was amplified 10,000 times (BAK Electronics, Germantown MD, USA), filtered by a dual variable filter VBF-8 (KEMO Ltd., Backenham, UK) (bandwidth 300–5000 Hz) and digitized (USB-6229, National Instruments, USA) at 10 kHz. During online classification, isolated spike shapes of selected unit are flashed on computer display and reproduced by a sound device to provide experimenters with audio feedback on neuron response.
where A is a real-valued data matrix A (n×24) containing spike shapes, n is a vector holding n th row average in the n th element, 1 is a n-vector of ones and X is a matrix holding the centering data.
where U is a n × 24 left orthogonal matrix of singular vectors giving the principal component scores which represent the spike waveforms in term of the PCs; V is a 24 × 24 orthogonal matrix detailing the spike profile and mapping the vector space; and S is a 24 × 24 diagonal matrix where the diagonal elements are the singular values of X.
A crucial problem in multivariate data analysis is the number of components to retain when applying PCA. This should be determined considering the tradeoff between dimensionality and the loss of data information . Both underextraction and overextraction may have consequences that adversely impact the efficiency and meaning of PCA , resulting in classification errors. Our software solves this problem by automatically using the Scree Test Optimal Coordinates (n oc ) method, as numerical approach to the Cattell’s scree test introduced in Ref . For online classification we maintain the parameter obtained during test acquisition phase until our measures of goodness of clustering detect that a new training phase is needed because the data have changed their properties too much (see Measures of cluster quality section).
where σ is the standard deviation of the distribution and N is the number of available samples, which corresponds to the total number of spikes or singular vectors in matrix and their -norm values in our specific case. The histogram was than thresholded to eliminate the noise content (low-amplitude peaks in the histogram which do not correspond to any distinct cluster) by using the value of first lower valley after the mode of -norm distribution. Thus, the number of peaks on histogram above the selected threshold was used as indication of the number of objects to be used for FCM clustering technique. However, in the interactive mode, FSPS leaves the possibility to choose the number of clusters manually, according to visual examination of the clustering results and expert judgment.
To perform FCM clustering (whose details are described in Additional file 2), one has to specify also an exponent m (m > 1.0), which determines the degree of fuzziness of the resulting clustering process. As m→1 the fuzziness of the clustering result tends to the results derived with the well known ISODATA method . As m→∞ the membership values of all the objects to each cluster tend to the reciprocal of the number of classes 1/c. The analysis of our data obtained during recordings, the distribution of PCs and measurements of cluster quality showed that, with m = 1.1 the FCM algorithm is performing clustering correctly on both real and simulated data (see Results). Besides, this value was consistently found to lead to good results with the well- and poor- separated classes and we were able to classify the cells with acceptable accuracy. Thus, we set this default value of m in our program, although it can be changed by the user if needed.
Since clustering algorithms define clusters that are not known a priori, it is fundamental to define a performance criterion to quantify the goodness of the resulting partition. In our FSPS software we have implemented the standard and popular figures of merit associated with FCM introduced by J.C.Bezdek (1981), such as the Partition Coefficient (pc), the Partition Entropy (pe) and Proportion Exponent (pex) that make use only of membership values and have the advantage of being easy to compute (see Additional file 2). However, these measures are often subject to numerical instability during the quantification of overlapping clusters of unequal size . Therefore, we also included in our software one recently introduced objective validity index L ratio . This parameter allows to obtain stable cluster evaluation from a particular recording site, and to take into account the clusters with a larger number of spikes. The evaluation of L ratio changes during single unit recordings was also used as measure of classification performance to monitor the stability of data acquisition. The threshold L ratio value is set to 5 as default, and can be modified by the user (we recommend a choice in the range 3–6). Once L ratio becomes bigger than its threshold value, the program alerts the user that it is advised to recompute full SVD and to update FCM prototype because classification is deteriorating.
Finally, for each isolated unit the FSPS software allows to compute a number of standard quantities and statistics of interest to the neurophysiologist, such as Peri-Stimulus Time Histograms (PSTHs), raster plots and inter-spike interval (ISI) and some clinically important indices that measured tonic and phasic activity, including burst index (BI) and pause index (PI) .
In order to validate our spike sorting approach and to compare it with other known algorithms, we tested it on simulated datasets described by Quiroga et. al. (2004) and compared to their already published results obtained by superparamagnetic clustering (SPC) and K-mean clustering techniques applied to different spike features (wavelets, PCA, using the first three PCs, and the whole spike shape) . The dataset contains two types of spike shapes: noisy “non-overlapping” spike shapes, which were generated by taking the target waveform and adding noise, and “overlapping spikes” which were generated overlapping spikes with a latency shorter than 0.7 ms and then adding noise. Performance was quantified in terms of number of classification errors.
Number of classification errors and noise levels, obtained using FSPS, SPC and K-means, in all simulated examples
Number of noisy spikes
Number of classification errors for all simulated examples and overlapping spike shapes
Number of overlapping spikes
We then tested our spike sorting algorithm and its quality of separation using real polyspikes recorded from premotor area F5 of a Macaque monkey. To demonstrate the functionality of our method in a realistic situation we present spike sorting results of two different datasets, obtained from the same recording site, but in different experimental conditions, thus introducing the complexity of non-stationarity in the data. First dataset, referred to as A, was acquired when monkey performed goal-directed grasping movements in full vision, while other recordings (referred to as dataset B) were done 30 min after A during the same grasping in the dark. Each dataset contains raw signals of twelve 3 s trials of each of the two conditions with a 10-12 s intertrial interval. The total duration of acquisition of one dataset was therefore 2–3 min. The structure of the two datasets and the result of spike identification are shown in Additional file 2: Table S1 (dataset A) and Additional file 2: Table S2 (dataset B). In the waveform extraction phase, 475 and 295 false positive peaks exceeding the threshold value were automatically removed from datasets A and B, respectively. Thus, 6191 spikes out of 6666 identified in the dataset A and 6379 out of 6674 in the dataset B were processed. Dataset A was used at first as test acquisition for the unsupervised spike sorting, creation of SVD model (U k , S k , and V k ) and to obtain the parameters that were necessary to further supervised spikes classification of the dataset B.
Performance of FSPS™ clustering/classification of datasets A and B with their respective quality measures
Number of spikes processed
Average processing time (at our system), ms
684 ± 21
44 ± 1
Number of Steps
Partition Coefficient (PC)
Partition Entropy (PE)
Class 1 (blue)
Class 2 (green)
Class 3 (pink)
Class 4 (yellow)
It is known that FCM algorithms have a tendency to find clusters of comparable size because they use a sum of squared errors objective function and approximately equal cluster populations result in smaller values of this objective function [17, 39]. This might become a serious problem in online applications, where the algorithm may be applied to relatively small stretches of data and so random fluctuations of spike rates may make the relative size of clusters at a given time very disproportionate.
Figure 14B shows the results obtained with real dataset A. The classifier is still performing well and true positive rate is higher 80% if it remains at least 20,1%, 4,9%, 1,1% and 5,1% spikes, respectively to Class 1, Class 2 Class 3 and Class 4. Figure 14C shows the outcome of the classification real dataset B. The true positive rate for each class remains still at higher than 80%, if classes 1, 2, 3 or 4 contain at least 18,9%, 16,1%, 13,3% and 12,5% spikes, respectively.
An additional exploratory test in which all clusters were modified in a balanced and “Uniform” manner showed excellent performance in simulated, real A and real B datasets (black dotted line in Figure 14A-C). In all these cases the classification accuracy was completely independent of the density of the clusters and the true positive rate was at its maximal initial value until the each cluster contained at least one spike.
The present work is devoted at addressing the challenges involved in balancing the different needs for accuracy, speed and automation in spike sorting.
The first point of discussion regards the selection of the optimal spike-shape features to be used for sorting [1, 40–42]. Here we chose the multivariate techniques for selection of spike features. The multivariate approach has proven successful in many industrial online applications [15, 43]. In this work, it is shown that neurophysiological research does not constitute an exception. We improved the performance of PC based classification by using a careful alignment of spike shapes and PSVD to reduce noise and computational time and to select the optimal number of components to be used, thereby choosing a reduced variable sets as inputs for the clustering algorithms. In agreement with other studies [8, 44–46], we have indeed found that when careful alignment and PCA were implemented, the results of classification of simulated datasets showed high tolerance to the noise, with better performance in comparison with similar methods.
Another goal of this paper was to demonstrate the applicability and efficacy of our technique for online isolation of single neurons. Indeed, a straight PCA approach is not sufficient to describe the changing pattern of neural activity adequately [11, 13]. There are mainly two reasons for this. First, the neuronal activity displays a non-stationary behaviour . Second, PCA is not an optimal method for feature extraction when the features are used in a supervised classifier . Finally, PCA is a computationally intensive preprocessing technique, making hard its use in real-time processing. Despite these conceptual difficulties our way of implementation of the mentioned algorithms into LabVIEW environment allowed us to run online classification with small and acceptable delay as far as our experimental conditions are concerned. There are several examples of applications of multivariate statistical online monitoring (and modeling) associated with PCA and FCM overcoming problems associated with non-stationarity [48–50].
A main focus of this article was the accuracy assessment of the FCM, which have been incorporated into FSPS software to produce both crisp and fuzzy classifications. The main advantage of FCM implementation is that fuzzy classification gives faster detection and smoother control than crisp classification. The importance of this becomes more evident during online classification. The tests reported in this article suggest that the FCM clustering copes well with the problems generated by the non-stationarity of the real data. It is interesting to note that the FCM procedure is based on an iterative clustering algorithm and can thus be regarded as an essentially unsupervised classifier. However, we also implemented a partially supervised mode, benefitting from partition and membership function previously obtained from training dataset. Consistently with this, the results of our tests performed on datasets prove that the implementation of FCM overcomes the problem of sensitivity for unequal cluster sizes, which is crucial for correct online classification. Thus, together with PCs extracted by PSVD from accurately preprocessed spike waveforms, FCM becomes a versatile noise tolerant technique for the sorting of neuronal action potentials having even small variation in their discharge.
Any kind of clustering or classification needs an objective measure of its quality. Although we implemented in the FSPS software most conventional indices associated with FCM, including partition coefficient, partition entropy and proportion exponent, our tests showed that the index L ratio was superior to classic FCM indexes and so was implemented as the default quality measure in our software. In terms of real recordings we found L ratio useful not only to determine whether the quality of a cluster is within acceptable limits, but also to control the stability of recordings, to predict the future behavior of the neuron and check whether the SVD model is going out of control.
An important feature of the FSPS software is that it is implemented entirely within LabVIEW. The latter constitutes one of the most frequently used programming languages for the data acquisition, analysis, control and visualization. LabVIEW is often faster than many other high-level programming languages used in neuroscience, such as MATLAB , and it is far better equipped for the development of experimental and clinically-oriented spike sorting applications . The upshot is that the entire spike time acquisition process can be run within a single environment, which has the all-important added benefit of simplifying experimental procedures.
In recent years there have been successful attempts at creating a cross-platform GUI for data visualization, navigation and spike sorting features within another software environment using the Python framework [52, 53]. Python applications, like “SpikeSort”, “Spikepy”, “spyke” and “OpenElectrophy”, provide adequate tools for the exploration of data and offline spike sorting, while “NeurOnline” provides the means for online spike sorting. However, the LabVIEW code we used is far more convenient because it rarely calls the Operating System (OS) directly, so it can be used with different OSs without the need for major modifications. Moreover, LabVIEW supports thousands of hardware devices and, in addition to the popular desktop OSs (Windows, Mac, and Linux), it can target several embedded real-time controllers, ARM microprocessors, and field-programmable gate arrays (FPGAs), allowing the deployment of our FSPS code with the most appropriate hardware platform without the need to learn new toolchains.
All mentioned properties contributed to the creation of a fast, powerful, user-friendly and stand-alone multi-platform software, designed for clustering/classification of neural data. Moreover, our online approach will help physiologist to overcome new challenges in experimental electrophysiological research.
We believe that the software developed here complements existing spike sorting toolboxes and will be a useful tool for fast on- and offline sorting of spike trains with limited supervision or fully automated. Because of these properties, our tool will be particularly useful for the analysis of large parallel recordings (where human supervision is practically impossible or inconvenient) and will therefore be important for improving our understanding of population codes [54–56] and for online applications such as the decoding of neural ensembles to control Brain Computer Interfaces or for clinical applications.
Project name: Neurolab
Project home page: http://www.spikesorting.com
Operating system(s): It was tested on Windows XP, Windows Vista, Windows 7
Programming language: NI LabVIEW 2009
Other requirements: for running in the online mode, the requirements are as follows. Hardware: Digital acquisition board from National Instruments (PCI or USB). Additional software: 1. LabVIEW Run-Time Engine 2009 for Windows 2000/7/7 x64/Vista/Vista x64/XP - (32-bit Standard RTE) - available free at: http://www.ni.com. 2. NI-DAQmx Run-Time Engine 9.3 or higher - (Core) for Windows 7 64 bit/7 x86/Server 2003 R2 (32-bit)/XP x86/Vista x64/Vista x86/Server 2008 R2 (64-bit) - available free at: http://www.ni.com;
License: FSPS software is distributed under Creative Commons Public License (CCPL BY-NC-ND) and can be used for non-commercial academic applications providing they properly reference this work in any publication that uses results generated by FSPS software.
Any restrictions to use by non-academics: Commercial License needed.
Field-programmable Gate Array(s)
Fuzzy SPike Sorting
Graphical User Interface
Principal Component Analysis
Peri-Stimulus Time Histogram(s)
Partial Single Value Decomposition
Single Value Decomposition
This work has been supported by RobotCub (ROBotic Open-Architecture Technology for Cognition, Understanding and Behaviour), IST-004370, and by the BMI Project of RBCS Department at the Italian Institute of Technology. The authors thank Stefano Panzeri for his feedback on software development and manuscript writing.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.