Quality metrics of spike sorting using neighborhood components analysis.

While an electrode has allowed for simultaneously recording the activity of many neurons in microelectrode extracellular recording techniques, quantitative metrics of cluster quality after sorting to identify clusters suited for single unit analysis are lacking. In this paper, an objective measure based on the idea of neighborhood component analysis was described for evaluating cluster quality of spikes. The proposed method was tested with experimental and simulated extracellular recordings as well as compared to isolation distance and L ratio. The results of simulation and real data from the rodent primary visual cortex have shown that values of the proposed method were related to the accuracy of spike sorting, which could discriminate well- and poorly-separated clusters. It can apply on any study based on the activity of single neurons.

INTRODUCTION

Studies of neuronal activity throughout the brain predominantly rely on extracellular recordings of spiking activity [1]. Since an electrode often simultaneously records spiking activity more than one neuron, spikes are always processed by sorting algorithms. Nevertheless, spike sorting algorithms are not perfect and classification errors may occur because of a number of reasons. Such as the recording hardware, electromagnetic interferences, the superimposed activity of multiple neurons, the spatially averaged activity of distant neurons and so on. Therefore, it is imperative to quantify unit quality after sorting for evaluating objectively results of spike sorting as well as improving the stability and reliability of neuron coding or decoding.

In most cases, the quality assessment was done subjectively by a human observer [2, 3], but the subjective methods don’t permit comparison of data quality across different studies, and unfortunately they are predisposed to personal bias. Therefore, several objective methods have recently been proposed for measuring unit quality. These methods can be classified into four categories: (a) Methods based on signal-to-noise ratio (SNR). These approaches have often served as a quantitative measure of unit quality in traditional microelectrode techniques because their significance is clear and can be easily calculated [4]. But they are not available when multiple units are recorded by the same electrode. (b) Methods based on refractory period. If spikes are from the same neuron, then there should be no interspike intervals less than the refractory period of neuron, which for most neurons is no less than 1 ms [5]. But these methods are only suited for checking the quality of isolation over a long collection period. (c) Methods based on statistic distribution. For example, Pouzat et al. [6] assumed that the noise obeys a Gaussian distribution, which used to evaluate the variability of spike waveforms. The Lratio proposed by Schmitzer-Torbert et al. [4] used a Chi-square distribution as a distance measure of the noise events in a feature space. But, there is no general agreement about the statistic distribution of noise. (d) Methods based on distance measure. For instance, Harris et al. [7] and Schmitzer-Torbert et al. [4] introduced isolation distance to quantify the quality of clustering by the minimal distance where the number of spike events and noise events is equal. The drawback of these methods is that they do not have a global scale to differentiate between well and poorly isolated units. Joshua et al. [8] described isolation score algorithm based on Euclidian distance to evaluate the quality of unit. The heavy computational burden makes it be sensitive to the number of spikes.

In order to circumvent these problems, a method based on neighborhood components analysis (NCA) was proposed for evaluating cluster quality in this paper. The NCA is a learning method for clustering multivariate data into distinct classes according to a given distance metric over the data, which is always used to dimensionality reduction [9] and fast classification [10]. Here, the idea to use the NCA is not to discover new variables which could be of interest, but to use the idea of NCA to estimate the probability that every event is correctly classified, and further estimate spike’s cluster quality. Functionally, the NCA serves the same purposed as the K-nearest neighbors algorithm, which makes no assumptions about the shape of the distribution of signals. So it can improve the estimate performance of cluster quality of spikes to some extent. In this method, each spike was firstly preprocessed through up-sampled, aligned and normalized; Then, the features of spike waveform were extracted by wavelet transform and Kolmogorov-Smirnov test; Finally, the quality of cluster was estimated through the NCA. Simulations and real data from primary visual cortex (V1) of rats were used to demonstrate the utility of this method.

MATERIALS AND METHODS

Real Data

The real data were collected using a 16-channel microelectrode array (Microprobe Inc., 2×8 platinum-iridium microwire array, 250 µm inter-electrode spacing, electrode impedance at 1 kHz equals 0.5～l.0 MΩ) and recorded from experiments performed on V1 anesthetized Long Evans (LE) rats. All experimental procedures were approved by Life Science Ethical Review Committee of Zhengzhou University. Neural signals were amplified and filtered with a computer-controlled multichannel amplifier (CerebusTM, Blackrock Microsystem Inc.). Spike signals were band-pass filtered at 0.25～5 kHz by a digital two pole Butterworth filter and sampled at 30 kHz. The implant procedure and recording conditions were described elsewhere [11].

Threshold amplitude crossing was used to detect the spikes and the threshold was defined as [12]., where x is the bandpass filtered signals and is an estimate of the standard deviation of background noise. And each spike was represented by using 2.1 ms (64 sampling points). The spikes were sorted by WaveClus approach. The WaveClus, which was proposed by Quiroga et al. [12], is an unsupervised spike sorting algorithm that combines wavelet transform and superparamagnetic clustering.

Synthetic Data

To test effectively the performance of proposed method, the simulated data from previously published data were used (Quiroga et al., publicly available online, http://www.vis.caltech.edu/～rodri/publications) [13]. The simulated data include four sets of simulated extracellular recording signals with a sampling rate of 24 kHz and the noise level determined by background noise’s standard deviation set to 0.05, 0.1, 0.15 and 0.2. Each simulation contains three distinct dominant single-unit spikes with normalized amplitude. In all simulations, three distinct spike trains have a Poisson distribution of inter-spike intervals with the mean firing rate of 20 Hz and the refractory period between spikes within the same category was set to 2 ms [12].

Neighborhood Components Analysis

Neighborhood components analysis proposed by Goldberger and colleges [10], which is a supervised learning method for clustering multivariate data into distinct classes according to a given distance metric over the data. The NCA aims at finding a distance metric that maximizes the leave-one-out error for a stochastic variant of K-nearest neighbors by finding a linear transformation of the input data. The key insight to the algorithm is to find a matrix A that corresponds to the transformation by differentiable objective function.

Setting consists of n real-valued input vectors in Rand corresponding class labels . And each vector xi selects another vector xj as its neighbor with some probability pij which was defined using a softmax over Euclidean distances in the transformed space,

where denotes Euclidean distances. Under this stochastic selection rule, we can compute the probability pi that vector xj will be correctly classified,

where denote the set of points in the same class as i.

The objective that we maximize is the expected number of points correctly classified under this scheme,

which is the differentiable cost function. Then A can be obtained using an iterative solver, such as conjugate gradient descent. The NCA used a probabilistic way to obtain the differentiable cost function, in other words, the cost function can be calculated by the average probability that every point in sample set belongs to the same set.

Quality Metrics Based on Neighborhood Components Analysis

In extracellular recordings, spikes recorded on a microelectrode represent a mixture of spikes obtained from one or more neurons and some non-neural spikes, such as noise events from motion artifacts, electromyographic activity, electric field pickup and so on [14]. In this study, a quantitative metrics based on NCA was proposed for evaluating cluster quality of spikes. Note, the aim of NCA is not only to discover new features which could be of interest, but to also use the idea of NCA to estimate the probability that every event is correctly classified, and further to quantify cluster quality. The purpose of using NCA is always the data visualization and fast classification [14] but here we are discussing it from the cluster quality evaluation perspective.

Setting k spikes from m neurons are recorded by a microelectrode and corresponding class labels in . In other words, k spikes are sorted to m clusters . To reduce the variability of spike waveforms from the same neuron due to discrete sampling of analog traces, we firstly up-sampled the data using cubic spline interpolation, the factor was fixed at 2. Then, all events were aligned by the point of maximum amplitude and normalized by standardized z-scores. All spikes after being prepressed were simply sorted two clusters, spike cluster and noise cluster. The spike cluster is a cluster classified as a single unit and the noise cluster is a cluster of events not classified with that unit in multi-units.

For each spike, a wavelet transform was calculated. A four-level multiresolution decomposition was implemented using Haar wavelets. The Haar wavelets were chosen due to their compact support and orthogonality, which allow the discriminative features of spikes to be expressed with a few wavelet coefficients and without a priori assumptions on the spike shapes. Thus, the Lilliefors modification of a Kolmogorov-Smirnov (KS) test [15] was used to select the optimal coefficients that best separate the different spike classes. The goal of this step is to reduce the differences among spikes from the same clusters and increase the differences among the different clusters, and improve estimation precision. In our method, the first 10 coefficients with the largest deviation from normality were used. For specific details please see reference [12].

For a given class, Ck contain Nk spikes, firstly computing the similarity between each event sx in spike cluster, , to all other events sy in noise cluster, ,

where is the Euclidean distance between vectors sx and sy. The parameter λ is a gain contrast () that sets the gain of this stretch. For a given event sx, when , all events are similar and similarity . On the other hand, when , all events are dissimilar and similarity . Here, , is the average Euclidian distance in spike cluster.

In order to turn the above similarity index into a probability-like quantity (positive values that sum to 1), the normalization was performed by the following formula,

where . Note that here the similarity index was normalized by the sum of the similarity index between sx and all other events from multiunits were recorded simultaneously.

For each event in spike cluster, sx, all the normalized similarity values were summed for all the Sy in the Ck cluster,

where is the probability that spike Sx belongs to this cluster. The calculation of was illustrated in Fig. ().

Finally, quality of the cluster is computed as,

and can be intuitively considered as the average probability that an event classified as a spike belongs to the cluster. The range of the NCA is from 0 to 1. If , which means a better cluster quality, in other words, with a small distance among spike events, and a large distance between them and noise events; If , it means a worse cluster quality, i.e. events from spike cluster are surrounded by the noise events.

Fig. () depicts an illustration of the calculation of . Correctly categorized events are shown in Fig. (), false positive and false negative errors are shown in Fig. ( and ), respectively. If noise events are classified as spikes by sorting algorithms (false positive errors), that are nonetheless close to noise cluster. The value would be reduced, due to its proximity to the other noise events, thus will reduce the overall NCA value (see Fig. ()). Likewise, when spikes are missed (false negative errors), they are close to spike cluster. As a result, events in the cluster that are close to such misses will be reduced , which in turn again reduce the overall value of NCA (see Fig. (1C)).

Quality Metrics Based on Isolation Distance and L

This subsection described two measures used to compare with the NCA, namely isolation distance and Lratio, and a complete description of the two methods can be found in [4]. The isolation distance is a measure of how well-separated a cluster is from the rest of the data set and Lratio indicates the distribution of non-cluster spikes around a cluster.

For cluster C, containing nc spikes, the isolation distance is defined as the squared Mahalanobis distance of the closest non-c spike to the center of C. The squared Mahalanobis distance is calculated as,

where xi is the vector containing features for spike i , and µc is the mean feature vector for cluster C. is the covariance matrix of spikes in cluster C. The isolation distance is not defined when nc is greater than the number of non-cluster spikes. A higher value indicates that non-cluster spikes are located farther away. Note that the value of isolation distance is not normalized against cluster size, so that clusters with a large number of spikes will tend to have a higher isolation distance.

L ratio is calculated as follows for cluster C,

where is the squared Mahalanobis distance between non-c spike i and the center of C;

is the Chi-squared cumulative distribution function describing the distribution of spikes in cluster C. The number of degrees of freedom is equal to the number of features used in the cluster space. A low Lratio indicates that there is a relatively empty space between the cluster and other spikes in the data set. Therefore, a lower Lratio and a higher isolation distance together indicate better cluster quality.

In addition to isolation distance and Lratio, the SNR was calculated for cluster . The calculating formula of the SNR is as follows [3]:

where Si is a waveform vector of spike i ; std(εi) is the standard deviation of the noise εi, is the average spike waveform in cluster C. The SNR is a traditional method for evaluating quality of a cluster, a low SNR indicates worse cluster quality.

RESULTS

Validation of Algorithm by Simulation Data

To evaluate the performance of the proposed algorithm, a signal with the noise level of 0.1 was selected, which contains 3 well-sorted spike clusters and 521 spikes. The events which were included in a circle of radius around center of a given cluster composed the spike cluster and the rest of the events composed the noise cluster. The false negative and false positive errors could be simulated as the radius increase or decrease. The results are shown in Fig. (). The principal component (PC) projections, spike waveforms and two-dimensional histogram of waveforms are shown in Fig. (, and ), respectively. In Fig. (), the red and black points represent events in the spike cluster and the noise cluster respectively. Note that the PC features only for representing conveniently representing the signals, other features, such as peak amplitude, energy, wavelet coefficients and so on, can also be selected.

Fig. () shows the performance of the NCA and ACC methods. The ACC was defined as [16] ACC = TP/(TP + FN + FP), where TP (True Positive) indicates the number of spikes which were classified as the spike cluster; FN (False Negative) indicates the number of spikes which were classified as the noise cluster, FP (False Positive) indicates the number of noises which were classified as the spike cluster. The ACC provides a quantitative metrics of the false negative and false positive errors. The value closer to 1, the smaller the errors and the higher the cluster quality. As shown in Fig. (), as the cluster was moved out of the distribution of spikes (see the black circle in (A)), both measures improved. When the radius , the value of NCA is consistent with the ACC. It is clear that the NCA can effectively measure the quality of spike cluster for cases in which the number of cluster spikes is not far greater than the number of noise events.

To further evaluate the performance of the NCA method, it was compared with Lratio (LRO), isolation distance (ISD), and SNR. Fig. () shows a comparison of the cluster quality among the NCA, the SNR, the ACC, the ISD and the LRO. For an easy comparison among them, the ISD, the LRO and the SNR were normalized by their maximum respectively. As shown in Fig. (), with the increase of radius r, the values of NCA, ISD and ACC increased at first and then decreased, the LRO, by contrast, descend at first and then increased, and the SNR increased all along. All methods attain the stability value when the ACC reaches its maximum value except the SNR. This demonstrates that these methods can be successfully used for measuring the cluster quality. Compared to other methods, however, the NCA is the nearest to the ACC. Although the performance of the LRO and the ISD is similar to the NCA, they don’t have a global scale to differentiate between good and poor separation. The SNR on spike cluster also increased with the radius r, but did not provide an unambiguous measure of cluster quality.

Validation of Algorithm by Real Data

Although the NCA performed well on simulated data, it was also important to demonstrate its performance using real data. Figs. (-) show three signals taken from micro-electrode array collected in rat V1 area. Fig. () shows data from the same electrode recording in which a well-separated cluster of spikes was observed. Histograms of the projection of each cluster onto Fisher’s linear discriminant between the spike cluster and the noise cluster are also shown (the linear discriminant, H, between the features of spike event S and noise event n were given by where µs and µn are the mean vectors and are the covariance matrices of the features, respectively. The projection of an event V onto the linear discriminant was given by ). The values of LRO, ISD, SNR and NCA for spike cluster were indicated through the black text on the right top corner of figure. Dots are the measured distribution and the solid line is the fitted curve using Gaussian function. For both clusters, a very little overlap of the projection of the spike and noise cluster was observed, and values of the NCA indicated that the clusters were well separated.

As can be seen in the PC projections and the histograms of the Fisher’s projection, the spike cluster in Figs. ( and ) are not as well separated as the examples shown in Fig. (). The worst separation was observed in Fig. (), while Fig. () had an intermediate separation. The values of NCA agreed with this subjective categorization of cluster quality. The clusters shown in Fig. () had the best quantitative separation for both measures, Fig. () had an intermediate value, and Fig. () had the worst values of NCA. These examples support the use of these cluster quality measures as an objective method for evaluating cluster quality.

DISCUSSIONS

In this study, we presented a quantitative measure of cluster quality for evaluating the separation of a cluster from other spikes recorded on the same electrode. It is based on the idea of NCA to estimate the probability that every event is correctly classified and further estimate the cluster quality. In this method, the spikes were firstly preprocessed, as being up-sampled, aligned and normalized; Then, they were sorted as the spike cluster and the noise cluster; Thirdly, the wavelet transform and the KS test were used to extract the features of spike waveforms; Lastly, the quality of spike cluster was estimated by the NCA. The best advantage of our method is its range from 0 to 1, which enables easy comparison of units recorded at different times and even by different research groups.

Closer to our proposed method, the isolation score algorithm proposed by Joshua and co-workers [8] applied the idea of NCA to measure the overlap between the noise (non-spike) and spike clusters, and used a nearest-neighbors algorithm to estimate the proportion of false positive and false negative classifications. However, there are two drawbacks of this approach: first, this method was applied directly to the spike waveform which reduces the difference of spike waveform from different neurons, and the performance was much more variable and less reliable; second, it applies the time stamps of spike trains and the entire analog signals as the initial input, largely increases the computation burden. In particular, we showed that the method proposed here overcomes these problems, giving a more robust estimation and a better performance.

In our method, the wavelet-based features were used as the input, which did not require the entire analog signals. The advantage of using the wavelet transform as a feature extractor is that very localized shape differences of the different units can be discerned, and the interference of noise can be greatly reduced. The information about the shape of spikes is distributed in several wavelet coefficients, whereas with PCA or other features [4, 17], most of the information about the spike shapes is captured only by very few variables, such as the first three principal components, which are not necessarily optimal to represent signals [12]. Since wavelets are a linear transform, therefore all the wavelet coefficients yield nearly the same results as the entire spike shape, the dimensionality reduction achieved by combining wavelets with the KS test may have a broad range of interest [12].

The performance of the proposed method was examined with simulated signals and it was then applied to real data from V1 of rats. With the synthetic signals, the performance evaluation was better than the one obtained using isolation distance and Lratio. The NCA is consistent with the ACC which reflects the real proportion of false positive and false negative errors (see Fig. () and Fig. ()), so this shows that it correlates with both the false positive and false negative errors. For the real data, the NCA can effectively discriminate well- and poorly-separated clusters, and its value is related to the correct identification of spikes (see Figs. (-)). Although isolation distance and Lratio are similar to the method in this way, they do not have a global scale to differentiate between well and poorly isolated units. Following Joshua et al. [8], we suggested to exclude units with the NCA value below 0.8 in studies whose conclusions may be influenced by the quality of the recorded units.

Although, the spike detection and sorting technology have been rapid development at present, the quality of single-unit after clustering is not uniform due to the interference of noise with larger amplitude or other neuron spikes. In this paper, we used the idea of NCA to provide a quantitative method for evaluating cluster quality. Unfortunately, the proposed method is invalid for cases in which the number of spike events is far greater than the number of noise events (see Fig. () and Fig. ()). In future, we will conduct research on this question. In most of the cases however, we found that the number of cluster spikes is lower than the number of the rest.

CONCLUSION

In this paper, a quantitative measure to evaluate spike cluster quality was presented and compared with isolation distance, Lratio and so on. The proposed method is based on the idea of neighborhood components analysis to quantify cluster quality. In this method, we used the wavelet-based features to represent the events for reducing the effects of noise. Its performance was tested with simulated data and real data, and the results proved that it provides better metrics for spike as there are no assumptions about the shape of the distribution of spikes and insensitivity to the number of spikes. Quantitative metrics of cluster quality is a necessary step in interpreting studies based on extracellular recording, so it is suggested that wider use of quantitative measures of cluster quality would likely improve the reproducibility of results across laboratories.