Robust, automatic spike sorting using mixtures of multivariate t-distributions

Abstract

A number of recent methods developed for automatic classification of multiunit neural activity rely on a Gaussian model of the variability of individual waveforms and the statistical methods of Gaussian mixture decomposition. Recent evidence has shown that the Gaussian model does not accurately capture the multivariate statistics of the waveform samples’ distribution. We present further data demonstrating non-Gaussian statistics, and show that the multivariate t-distribution, a wide-tailed family of distributions, provides a significantly better fit to the true statistics. We introduce an adaptation of a new expectation-maximization based competitive mixture decomposition algorithm and show that it efficiently and reliably performs mixture decomposition of t-distributions. Our algorithm determines the number of units in multiunit neural recordings, even in the presence of significant noise contamination resulting from random threshold crossings and overlapping spikes.

Introduction

Extracellular recordings of neural activity provide a noisy measurement of action potentials produced by a number of neurons adjacent to the recording electrode. Automatic and semiautomatic approaches to the reconstruction of the underlying neural activity, or ‘spike-sorting’ have been the subject of extensive development over the past 4 decades and reviews of early and recent efforts can be found in the literature (Schmidt, 1984, Lewicki, 1998). It is generally assumed that each neuron produces a distinct, reproducible shape, which is then contaminated by noise that is primarily additive. Identified sources for noise include: Johnson noise in the electrode and electronics, background activity of distant neurons (Fee et al., 1996b), waveform misalignment (Lewicki, 1994), electrode micromovement (Snider and Bonds, 1998) and the variation of the action potential shape as a function of recent firing history (Fee et al., 1996b, Quirk and Wilson, 1999). Given this signal+noise structure, the problem of automatically classifying the different shapes is a clustering problem and can be addressed either in the context of the full time-sampled spike-shape or of a reduced feature set, such as the principal components or a wavelet basis (Hulata et al., 2002).

While the application of general clustering methods such as k-means (Salganicoff et al., 1988), fuzzy c-means (Zouridakis and Tam, 2000), a variety of neural-network based unsupervised classification schemes (Ohberg et al., 1996, Garcia et al., 1998, Kim and Kim, 2000) and ad-hoc procedures (Fee et al., 1996a, Snider and Bonds, 1998) have been pursued by some authors, a number of other studies (Lewicki, 1994, Lewicki, 1998, Sahani et al., 1997, Sahani, 1999), attempting to provide statistically plausible, complete and efficient solutions to the waveform clustering problem, have focused their attention on clustering based on a Gaussian mixture model. The assumption underlying the latter approach is that after accounting for non-additive noise sources (e.g. misalignment, changes during neural bursts), the additive noise component is Gaussian-distributed. As a result, the waveforms resulting from each neuron are samples from a multi-dimensional Gaussian distribution with a certain mean and covariance matrix. Given this statistical structure, it is possible to construct an appropriate statistical model of the data and apply the powerful method of Gaussian mixture decomposition to solve the clustering problem (Jain et al., 2000, McLachlan and Peel, 2000). This allows estimation of model parameters such as the shape of the individual waveforms and the noise characteristics. The estimated model parameters are used to classify each ‘spike’ to one of several mixture components that correspond to different neural units (or possibly noise).

Although the statistical framework resulting from the multivariate Gaussian model is powerful and well studied, recent evidence suggests that it may provide an inaccurate description of the spike statistics (Harris et al., 2000). Examination of the distribution of Mahalanobis squared distances of spikes produced by a single unit reveals a discrepancy between the expected χ²-distribution and the empirical distribution, which exhibits wider tails. Algorithms based on the Gaussian assumption may therefore be ill suited for the task of automatic spike sorting, in particular as it is well known that they are not robust against a significant proportion of outliers. In this study, we provide additional evidence for the non-Gaussian nature of spike-shape statistics and demonstrate that an alternative model, one using multivariate t-distributions instead of multivariate Gaussians is better suited to model the observed statistics. Multivariate t-distributions have attracted some recent attention in the applied statistics literature (Lange et al., 1989), and a mixture decomposition algorithm for multivariate t-distributions was developed (Peel and McLachlan, 2000), based on the expectation-maximization (EM) algorithm. This algorithm requires computation of twice as many hidden variables as in Gaussian mixture decomposition algorithms, and involves an additional computational step for adapting the ‘degrees of freedom’ (DOF) parameter.

In addition to the choice of a statistical model for the mixture components, practical EM-based mixture decomposition algorithms need to address a number of issues including the determination of the number of components, the choice of an initialization procedure and avoiding convergence to local likelihood maxima or parameter singularities. Determination of the number of components in a mixture model has been the subject of extensive research (reviewed by Sahani, 1999, McLachlan and Peel, 2000, Figueiredo and Jain, 2002). The methods most widely used for this task are based on selecting the best mixture models from a set of candidates with different numbers of components. After fitting the parameters of the different models (using the EM algorithm) the different models are compared using a penalized-likelihood function, which penalizes the likelihood for ‘complexity’ (i.e. a larger number of components) and an ‘optimal’ model is found. This class of methods has the disadvantage of requiring estimation of the parameters of multiple mixture models. Other approaches include the use of stochastic model estimation using model-switching Markov-Chain–Monte-Carlo methods (Richardson and Green, 1997), and deterministic annealing based approaches (Sahani, 1999), which we have recently adapted to the case of the multivariate t-mixture model (Shoham, 2002). These approaches suffer from significant computational complexity, and, in addition, annealing approaches are quite sensitive to the specific choice of an annealing schedule. A recently introduced algorithm (Figueiredo and Jain, 2002), provides a new strategy where a process involving competitive elimination of mixture components drives a modified EM algorithm towards the optimal model size, simultaneously with the model parameter estimation. This approach appears currently to offer the best overall profile in terms of computational simplicity, efficiency and selection accuracy, and tends to avoid the usual difficulties of initialization sensitivity and convergence to singularities associated with the EM algorithm. We provide an adaptation of this algorithm for the case of multivariate t-distributed components. Our final algorithm is statistically plausible, simple and well-behaved and can effectively deal with many real data sets.