Spike sorting for large, dense electrode arrays

Abstract

Developments in microfabrication technology have enabled the production of neural electrode arrays with hundreds of closely spaced recording sites, and electrodes with thousands of sites are under development. These probes in principle allow the simultaneous recording of very large numbers of neurons. However, use of this technology requires the development of techniques for decoding the spike times of the recorded neurons from the raw data captured from the probes. Here we present a set of tools to solve this problem, implemented in a suite of practical, user-friendly, open-source software. We validate these methods on data from the cortex, hippocampus and thalamus of rat, mouse, macaque and marmoset, demonstrating error rates as low as 5%.

Figure 1. High-count silicon probe recording

(a), Layout of the 32-site electrode array used to collect test data. (b), Short segment of data recorded in rat neocortex with this array. Color of traces indicates recording from the corresponding colored site in (a). Black rectangles highlight action potential waveforms; note the frequent occurrence of temporally overlapping spikes on separate recording channels.

Figure 2. Local spike detection algorithm

(a), Adjacency graph for the 32-channel probe. (b), Segment of raw data showing two simultaneous action potentials on spatially separated channels (scale bars indicate 0.5mV / 10 samples). (c), High-pass filtered data shown in pseudocolor format (units of standard deviation). Vertical lines on the colorbar indicate strong and weak thresholds, θ_s and θ_w (respectively 4 and 2 times standard deviation). (d), Gray-scale representation showing samples which cross the weak threshold (gray), and the strong threshold (white). (e), Results of two-threshold flood fill algorithm, showing connected components corresponding to the two spikes in orange and brown. Note that isolated weak threshold crossings resulting from noise are removed. White lines indicate alignment times of the two spikes. (f), Pseudocolor representation of feature vectors for the two detected spikes (top and bottom). Each set of three dots represents three principal components computed for the corresponding channel (arbitrary units). Note the similarity of the feature vectors for these two simultaneous spikes (top and bottom). (g), Mask vectors obtained for the two detected spikes (top and bottom; 0 represents completely masked, 1 completely unmasked). Unlike the feature vectors, the mask vectors for the two spikes differ. Each set of three dots represents the three identical components of the mask vector for the corresponding channel.

Figure 3. Evaluation of spike detection performance

(a), Waveforms of the 10 donor cells used to test spike detection performance, in order of increasing peak amplitude (left to right). (b), Fraction of correctly detected spikes as a function of strong threshold θ_s (left), weak threshold θ_w (center), and power parameter p (right). Colored lines indicate performance for the correspondingly colored donor cell waveform shown in A; black line indicates mean over all donor cells. (c-e), Dependence of the total number of detected events, timing jitter, and mask accuracy on the same three parameters.

Figure 4. Evaluation of automatic clustering performance

(a), Receiver-Operating Characteristic (ROC) Curve showing the performance of the Masked EM algorithm (blue) and Classical EM algorithm (red) on one of the 10 hybrid datasets; each dot represents performance for a different value of the penalty parameter. The cyan curve shows a theoretical upper bound for performance, the best ellipsoid error rate (BEER) measure obtained by cross-validated supervised learning. (b), Mean and standard error of the total error (false discovery plus false positive) over all 10 hybrid datasets for theoretical optimum (BEER measure), Masked EM and Classical EM algorithms. For each dataset and measure, the parameter setting leading to best performance was used. (c), Effect of varying the penalty parameter (as a multiple of the AIC penalty) on the total error for both algorithms. The dotted line indicates the parameter value corresponding to BIC. Note that the Masked EM algorithm performed well for all penalty values. (d), The number of clusters returned by the Masked EM algorithm as a function of the penalty parameter.

Figure 5. The “Wizard” for computer-guided manual correction

(a), Illustration of the measure used to quantify cluster similarity. p_ij represents the posterior probability with which the EM algorithm would assign of the mean of cluster i to cluster j. (b), To test this measure, the clusters corresponding to hybrid spikes were artificially cut into halves of high and low amplitude. In each case, the similarity measure identified the second half as the closest merge candidate. (c), The Wizard identifies the best unsorted cluster as the one with highest quality (top), and finds the closest match to it using the similarity matrix. (d), The Wizard algorithm. The best unsorted cluster and closest match are identified. The operator can choose merge the closest match into the best unsorted, ignore the closest match, or delete it by marking it as multiunit activity or noise; the wizard then presents the next closest match to the operator (blue arrows). After a sufficient number of matches have been presented, the operator can decide that no further potential matches could have come from the same neuron, and either accept the best unsorted cluster as a well-isolated neuron, or delete it as multiunit activity or noise. The wizard then finds the next best unsorted cluster to present to the operator (orange arrows).

Figure 6. Screenshot of the KlustaViewa graphical user interface

In order to make the decisions presented by the Wizard, the operator has access to information including waveforms (center panel; gray waveforms correspond to masked channels), principal component features (top right), auto- and cross-correlograms (bottom right), and an automatically computed similarity metric for each pair of clusters (inset). To enable rapid navigation, all views are integrated; for example, clicking on a particular channel in the Waveform View will update other views to show the selected channels or clusters.

Figure 7. Consistency of manual curation across operators

(a), Performance of 8 human operators (5 experts, 3 novices) on a “drifty” hybrid cell requiring manual curation (see supplementary figure 13b). A tick indicates correct merging of the split hybrid cell, a cross indicates this merge was not performed. (b-d), consistency of assignments of multiple operators over all cells in this dataset. Each submatrix shows the conditional probability of the first operator’s cluster assignments given the assignments of the second operator (color scale at bottom of (d)). (b), consistency of cluster assignments for spikes marked as well-isolated by all operators; (c), consistency of cluster assignments for spikes marked as well-isolated by at least one operator; (d), consistency of whether spikes were marked as well-isolated by different operators. (e-g): Operator consistency for the analyses of (b-d) was quantified using the Fowlkes-Mallows index, for which 1 represents complete agreement and 0 complete disagreement. Note that while cluster assignments were highly consistent between all expert operators, the operators were often inconsistent in their judgements of which units were well-isolated.