MIEN Spike Sorter

The spikesort block provides a collection of tools for addressing the "spike sorting" problem: e.g. given a recording of electrical activity from a collection of neurons, how do we detect action potentials and classify them as the activity of particular individual Neurons. The toolkit includes a custom spike sorting GUI, built as a MIEN Dataviewer extension, and a set of DSP blocks and command line tools for batch processing. The spike sorting tools in this module focus on using phased array sorting methods, which are appropriate for use in systems where action potentials travel in a particular direction, and recordings can be made from an array of locations in space. In general these conditions occur in peripheral nerves.

Summary of the phased array method

Phased array spike sorting is a method for taking advantage of information inherent in the propagation rate of spikes. In tissues where spikes propagate over long distances (most notably peripheral nerves), The signals from different cells will, in general, propagate at different rates, due to differences in axon diameter, channel density, and ion concentrations. These differences in propagation rate can be used to assist in rapid identification of signals from a particular neuron. In order to take advantage of propagation rate information, the recording electrode array must be constructed so as to record signals at several stages of propagation, so that each signal is recorded earlier in time at more proximal electrodes, and again later in time at more distal electrodes.

Given such a recording array, signals with a particular propagation speed v will reach the first electrode e1 at time t1, and a second electrode e2 at distance x further away at t2=t1+x/v. Knowing x and v thus allows us to construct an optimal average signal by first offsetting the voltage time series from e2 by t2-t1 relative to the trace from e1 and then adding or averaging the two time series. This offset average will attenuate signals with propagation velocities other than v, but maintain the intensity of signals with propagation velocity v, and also improve the signal to noise ratio of these signals by signal averaging. Signals associated to units with propagation velocity v can then be extracted from the offset average (the discriminant function) by simple thresholding. We refer to this process as shift and sum discrimination. The process works for any array of electrodes, provided that a shift template specifying the time offsets between each adjacent pair of recording sites can be found for each signal of interest.

The shift and sum method can be implemented very quickly, and shift templates require substantially less human interaction to construct than more conventional waveform templates. The shift and sum method is also excellent at removing crosstalk and other electrical artifacts from recordings, since these artifacts often occur with 0 time delay between electrodes, and are thus tremendously attenuated by any shift template. In some systems, shift and sum alone is sufficient to separate spikes. In other systems, however, additional methods must be used to aid in separation.

Obviously, the simple shift-and-sum method will not work in systems where many units have very similar propagation velocities, or if individual units have highly variable propagation velocities. Also, even when the propagation velocity does provide sufficient information to separate spikes, more complex algorithms may be needed to account for additional complexities in the data, including:

In some cases, propagation velocities may be too similar (or too variable) to allow for discrimination at all. Even in these cases, the phased array method provides some advantages. Here, it is used as a pre-processor for additional sorting methods that use information in the spike waveforms. This multiple-pass sorting method has several advantages:

The spike sorter extension provides tools for all three types of phased array spike sorting: simple phased array separation, advanced phased array algorithms, and use of the phased array system as a preprocessor for waveform template, filter, and cluster-based spike separation algorithms.

Outline of Spike Sorting Procedure

The process of spike sorting using the phased array method follows the following outline:

The spike sorting package provides tools to assist in all of these steps except #1.

One goal of the package is, obviously, to implement the sorting procedure described above. The basic procedure is quite simple, however, and could be quickly implemented with custom software. The need to choose appropriate templates, discriminants, and thresholds interactively, and to evaluate the reliability of the results presents a more challenging problem. Two additional goals of the package are to provide intuitive user-interaction features for developing and testing sorting criteria for particular data sets, and to provide a large toolbox of discriminant functions, error functions, and spike waveform statistics to facilitate the development of effective spike sorting tool chains for a wide variety of different types of extracellular data.

Features of the Spike Sorting Package

The spike sorting package includes the following modules:

Discriminants

Successful detection of events in a phased array spike sorting tool chain depends on selecting an appropriate shift template for a particular unit. In most cases, however, detection depends equally strongly on the choice of an appropriate discriminant function. This function transforms the many channels of (appropriately shifted) input data into a single time series, and it must do so in such a way that action potentials in the unit being selected consistently generate larger (or smaller) values of of the discriminant than any other type of event in the data set. The spike sorting package uses aa Schmidt trigger to detect spikes. This detection method provides some protection against accidental detection of noise events and multiple detection of single events, but it is still a simple thresholding algorithm in the sense that if it is set to consistently detect some events with some magnitude A, it will also detect any other events with magnitude >=A. Consequently selection of a particular event type over another must be implemented by the discriminant function.

In a classical phased array, there is only one choice for discriminant function, and this is the superposition of all the (shifted) input channels . The spike sorting package uses the mean of the channels instead of the sum, but this has the same properties as the sum.

The mean discriminant has the advantage that it is extremely quick to calculate, and introduces few assumptions beyond the fundamental assumption that different signals have different propagation rates. This discriminant can generate good results if each signal has a unique propagation rate and a similar absolute amplitude. In communications this is often the case. For example, a phased array antenna used for directionalizing a signal is in fact comparing a signal to other recordings of the same signal, taken with different orientations relative to the antenna. In this case, the amplitude of all signal types is identical, and the effective propagation velocity is simply the speed of light (c) times the cosine of the angle between signal source and array (Theta). Here a shift template based on a velocity of c and the sum discriminant will reliably select the array orientation Theta=0.

For spike waveform data, however, the sum (or mean) discriminant is often insufficient. There are several common causes for this:

8Signals from several units can occur at the same time (superposition). Temporal alignment will partially resolve this situation, especially if the superposed units have different propagation speeds, but even so, a superposition of two nearly-aligned units may generate a sum greater than a single perfectly aligned unit. This is a sufficiently difficult and common problem that a separate section is devoted entirely to it.

Clearly, if difficulties 1 and 2 are too severe, the basic phased array method simply is not sufficient for spike separation, and the discriminant function must use additional information (beyond what is contained in the shift template) in order to separate spikes. In the systems we have tested so far, however, propagation rate from a particular unit appears to be highly consistent. The rates are so consistent that even though many units have quite similar propagation rates, there is still enough information in the the rate alone to separate the signals.

On the other hand, difficulty 3 has proven to be a major obstacle to rapid and accurate spike sorting. In the cricket abdominal nerve cord, extracellularly recorded action potential amplitudes vary by more than a factor of 10. Consequently, even when the shift operation attenuates non-optimal signals robustly, for example by a factor of two or three, a very strong, but attenuated, unit can still be stronger than a weak, but perfectly aligned unit. One possible solution to this difficulty is to detect the largest units first, and subtract out their signals before moving on to detect smaller units. This iterative process can be effective, but requires a a human operator to spend a great deal of time in the interactive phase of the sorting process. We have found that choosing a more appropriate discriminant function is an effective way to simplify this process.

Choosing a discriminant (or set of discriminants) appropriate to the features of a particular data set is often the best way to improve sorting efficiency and accuracy. To this end, the spike sorting package provides a range of discriminant functions, and is designed to facilitate rapid development of additional algorithms. The mean discriminant function has already been described. The following sections will cover the other functions currently provided by the package, followed by a section describing the processes for adding a function.

NOTE concerning sign conventions:

In general, extracellularly recorded spikes have a more robust negative phase than positive phases. Consequently, most of the functions in the spike sorting package operate on minima and negative-going threshold crossings, rather than maxima and positive threshold crossings. The included discriminants are all constructed so as to attain smaller (more negative) values when a probable event occurs. This allows them to be used in the same manner as the mean discriminant for typical extracellular data. For data sets that have positive-going events, it is probably easiest to invert the data during pre-conditioning, rather than construct positive threshold methods particular to these data.

Mean/Std

This method returns the mean discriminant divided by he across-channel standard deviation of the input data. Consequently, minima only occur if the mean is strongly negative, but they also only occur if the negative values across the various recording channels are consistent. This can work well for normalized data, if the real spike events have a consistent waveform across channels. In this case, the standard deviation "weighting" will reduce the signal due to anomalies, spike superpositions, and poorly aligned high amplitude spikes.

Mean/OptMin

This is a relatively new discriminant method that has proven to be very successful for dealing with data that contains a range of spike amplitudes. The function returns the mean discriminant divided by the "optimal minimum" value of the data. The optimal minimum is a parametric function that takes a single parameter: the template width. For each time point in the data, this function calculates the best possible alignment within the template width. This is the collection of shift values that would result in the most negative possible mean, under the constraint that no shift value is greater than the template width. The function then returns the resulting most negative possible mean. When used as a discriminant function, the mean is divided by the absolute value of the OptMin, so the function attains a minimum of -1 when a perfect alignment is detected.

The optimal minimum weighting is particularly useful for avoiding detection of poorly aligned but powerful spikes. A strong spike with bad alignment will still return a very strong (highly negative) mean, but it will return an even stronger optimal minimum (since that will correspond to the mean that would be attained if the strong spike was also perfectly aligned). This results in a weak value for the Mean/OptMin ratio. In our artificial data, and also in cricket abdominal nerve cord data, the OptMin algorithm has been shown to detect spikes with less than half the amplitude of the largest spikes in the data, without first removing the signals from the large spikes. This even works when the smaller spikes propagation velocities that differ form the large spikes by only ~10%.

In addition to rejecting poorly aligned large spikes, OptMin algorithms will tend to reject most spike superpositions. This occurs since, in a complex superposition waveform, there is often some way to line up the waves so as to obtain a stronger mean than would be present for any of the individual component spikes taken alone, even if they were to be ideally aligned. This property is both a strength and a weakness for the OptMin algorithms. Superpositions are often hard to classify accurately, so these algorithms will often trade potential misclassification errors for detection failures, and this is often desirable. Also, if the goal is to calculate clean waveform templates, superpositions should not be included in the detection in any case, so the OptMin algorithms are ideal. However, if the purpose of an experiment is to generate finalized event IDs for use in, e.g. a study of spike time correlations, a systematic failure to detect superposed spikes would be a major weakness, and a second round of detection after the OptMin algorithm will be needed to resolve the superpositions (see this section).

One flaw in the Mean/Optmin ratio is that it can generate strong peaks when there are no spikes, or when there are artifacts, because, although the mean is low, the optimal minimum is also extremely low. The OptMin algorithm included in the spike sorter package is protected against division by zero, so it won't fail (or even generate a strong peak) if, for example, the input channels go dead at some time, resulting in an OptMin of zero. There are other more subtle glitches that can confuse the algorithm however. Some examples include a spike that fails after the first recording site, or an electrical artifact that occurs on only one channel, during an otherwise quiet stretch of recording. In these cases, the Single event is, trivially, perfectly aligned (with itself), so the mean/optmin can attain a value of -1.

Several modifications of the OptMin algorithm are provided to avoid this sort of problem. The first of these is Mean squared/OptMin. This is simply the Mean/OptMin ratio times the mean, and it consequently has a minimum value of the mean minimum, rather than -1. This avoids detection of events with a weak mean, but it is less powerful than the raw Mean/OptMin for avoiding detection of poorly aligned strong events. The second variation is Mean/OptMin with Threshold. This variant takes two parameters, the template width and a threshold. The value returned is the same as Mean/OptMin if the mean is stronger than the threshold, and zero otherwise. This avoids detection of anomalies with very weak means, while maintaining the strength of the basic algorithm for avoiding detection of poorly aligned strong spikes. For our cricket nerve cord data, Mean/Optmin with Threshold has proven to be the most effective discriminant function.

The spike sorting package includes an algorithm for calculating the optimal minimum that is efficient, and implemented as a C binary extension, so the OptMin based discriminants can be calculated relatively quickly (though they are still noticeably slower than the simple mean).

Template Match Algorithms

The previous discriminants require only the information in the shift template and the input data to calculate. Consequently, they can be calculated and used in single-pass spike detection. However, in a situation where the propagation velocities (and therefore the shift templates) don't contain sufficient information to separate spikes, no discriminant function of that sort can possibly succeed. In these cases, a multi-pass sorting strategy is required. In this case, the initial pass generates a waveform template for a particular unit, and the second pass uses the information in this waveform template in addition to the input data and shift template to calculate a discriminant function. In the most general case, the waveform template is an ensemble of raw input data segments associated to events from a particular unit. To reduce the problem dimensionality, the spike sorting package actually extracts data segments from the shifted input data rather than the raw data. These segments can be shorter, since they don't need to contain samples representing the activity at a particular electrode before the spike of interest has reached there, or after it has passed. In addition, all the current waveform-based discriminants use only the first or first and second statistical moments (calculated across events, not across channels) of the waveform ensemble.

The simplest matching algorithm is template match. This function uses only the mean of the waveform template, and, for a sliding widow in the input data, calculates the euclidean distance of the data segment from the template mean. The function thus returns a minimum possible value of zero, if a given data segment exactly matches the template, and returns more and more positive values as the data segment differs more from the template.

An alternative to template match is spike match. This algorithm functions the same way, but uses a particular spike waveform in place of the mean of a template. In general, this is less reliable, but it can be useful for short data sets that contain many superpositions and few clean isolated spikes. It also requires less human interaction time to use than template match, so it can be handy for rapid interactive testing of new tool chains.

A refinement of the template match is the Mahalanobis match. This version uses the first two moments of the waveform template, and calculates the Mahalanobis distance of data segments from the template. Predictably, this generates better results for noisy data.

All of the template matching algorithms are implemented as C extensions for improved runtime speed.

Adding discriminant functions

Like the rest of MIEN, the spike sorter package has a modular design that facilitates rapid development of additional algorithms with little to no "wrapper" or interface code writing. Adding a discriminant function to the package requires modifying only one module, and can require as few as 3 lines of new code.

The module to edit is the discriminants.py module, and there are two steps:

The discriminant function API is even simpler than the general DSP function API. Every discriminant function takes exactly three arguments: dat, pars, and template, and returns a 1D float array.

Dat is a 2D numpy array of floats, containing the shifted input data. This array has shape (N, M), where M is the number of recording channels and N is the number of samples. Your function must return a 1D array of length N specifying the value of your discriminant at each time sample.

Pars is a tuple of parameters, used by parametric discriminant functions. It has a default value of an empty tuple, if no parameters are passed in. Discriminant functions can safely ignore this argument if they don't require parameters. If your function does require parameters, you should supply default values, and check to make sure that pars has an appropriate length before assigning parameters from it. If it is too short, use default values. In batch and command line scripts, pars should always be fully specified, but the GUIs may call your function without specifying any parameters, or with a user supplied tuple that might be of the wrong length. Regardless of the supplied parameters (including when there are none), the pars argument is always a tuple, so it is safe to call, e.g. "len(pars)" or "if pars:"

Template has a default value of None, but if it is a true value, it is a MIEN Data instance containing the spike sorter template that is currently being evaluated. Single pass discriminant functions can ignore this argument. Multiple pass functions will need to extract, for example, a waveform template from it. The following:

tem=template.getElements("Data", "template", depth=1)[0] td=tem.getData()

will result in td referencing a 2D float array containing the first the mean and variance of the waveform template, or will raise an IndexError if no template is specified. Td contains 2M rows, such that row 2i specifies the mean waveform on channel i, and row 2i+1 specifies the waveform variance of channel i. More detailed information can be accessed form the template object, but this process is beyond the scope of this overview, and requires understanding the MIEN Data API.

You can use the existing functions in the discriminants.py module as examples when building your new functions. As the simplest possible example, here is an implementation of the simple sum discriminant (just like the mean, except that the amplitude is scaled by a factor of M):

def calc_sum(dat, pars=(), template=None): return dat.sum(1)

The discriminants module imports all symbols from the errors.py module, and error functions require the same API as discriminant functions, so your functions can get the value of any defined error function (for example for generating a new mean-to-error-ratio discriminant) simply by calling that function and passing allong the same arguments (some complex functions may need to first remove some entries from the "pars" tuple). You will need to call the error function using its Python function name, not the GUI string name. You can find out what this is by looking at the ERROR_FUNCTIONS dictionary defined at the end of the errors.py module (this dictionary acts exactly like the DISC_FUNCTIONS dictionary mentioned above).

Error Functions

Error functions are functions that increase in value for input data indexes that are, by some measure, poor matches for the desired unit template. They can be used in several ways. The most common useser are to include them in the denominator of discriminant functions (as described in the last section) or to display them in parallel with the discriminant function in the interactive GUIs in order to quickly detect false positive event detections, or poor choices of discriminant function.

The spike sorting package provides a set of predefined error functions in the errors.py module. Error function operate similarly to discriminant functions, and you can add custom functions of your own using the same API described in the previous section. Since discriminant functions should attain minimum values when there is a good match (and low error), any discriminant function can also be used as an error function, and discriminants will appear in GUI menus where an error function is called for. The reverse is not true, since some error functions can take minimal values when there is no event match (e.g. Std), and/or are not well defined over the whole set of input data (e.g. Residual).

The existing error functions are:

Std

The across-channels standard deviation of the input data.

OptMin

This quantity is described in the previous section as the denominator of Mean/OptMin. It can also be used by itself as an error measure (although some of its properties are not intuitive for an "error").

Residual

This function requires a previously calculated waveform template, and is only well defined in the vicinity of spike events. It is the value of the across-channels mean of the imput data, if the spike template was first subtracted from the data at the time of the putative spike event. This value is useful for determining if a particular event is similar to the template or not. In particular, for a putative superposition of two events, one of which is from the selected template, the residual should look like the mean discriminant for the other event.

Post-Hoc Analysis

One difficulty with constructing an effective spike sorting algorithm for a particular type of data is that any choice of discriminant, parameters, and shift templates will always "succeed" in the sense that a set of (supposedly) classified events will always be generated. Normal completion of the algorithm therefore can't be used as an indication that the algorithm is a good one (in the sense that it generates accurate classifications).

When a user is developing tool chains interactively, error functions (particularly residual) and the standard deviation view of the waveform template help to provide fast feedback about the consistency of events detected by a particular template. Consistency is not a certain proxy for accuracy, however, and even if a template is accurate over the interactive test data, it may not remain accurate over an entire data set. Various forms of non-stationarity can lead to errors during batch processing. The worst case is a change in spike waveform, recording quality or location, or spike propagation rate (the latter may be caused by drying or ion accumulation in the prep). Even if these conditions remain constant, changes in spike response rates or synchrony can lead to complex superposition waveforms that were rare or absent in the training data and will be miss-classified.

In order to trust the classifications from any experimental spike-sorting algorithm, some post-hoc analysis is needed to test the reliability of classification.

In general, the real event IDs are not known (otherwise why would one need to spike sort?), so direct validation of a sorting algorithm is impossible. Confidence in a sorting method must therefore be gained through less direct testing.

The best possible way to validate a sorting method is to jointly record extra- and intracellular data in the same prep. In this case, the exact event times for one unit are known (from the intracellular record). The sorted extracellular data is thus expected to include a unit ID with the same event times as the known unit. With several such data sets from similar preps, including intracellular records from several cells, it is possible to directly validate sorting algorithms for this type of prep. A reliable algorithm is one that produces (at least to good approximation) a unit ID corresponding to the known unit in each data set. As the number of matched data sets grows, confidence in the algorithm, and in the stability of data recorded from this type of prep, also grows. Eventually, experimenters can discontinue use of the intracellular recording, and continue to use the extracellular electrodes and proven sorting algorithm in the same type of prep.

The method in the previous paragraph is an ideal situation. Unfortunately, dual extra- intracellular experiments can be difficult and time consuming, even in accessible preps, and may be entirely impossible in less accessible ones (including the very relevant case of chronic implants in mammals). In the absence of this type of data, even less direct methods of validation are needed. None of these indirect methods is entirely certain, but they can increase confidence in a sorting procedure, particularly if the same procedure performs well by several different metrics.

Indirect methods of validation supported by the spike sorting package include several broad categories:

Synthetic Data

These methods center around the development of synthetic spike data with similar statistics to the natural data that is being analysed. Because these data are syntheticlly generated the actual event IDs are known, allowing for direct evaluation of sorting performance on these data sets.

The weekness of these methods, of course, lies in the possibility that the synthetic data sets won't capture the actual range (or structure) of variability present in the real data. In particular, there is a bootstrapping problem, since estimating the statistics of the data accurately may require that the spike waveforms are already identified.

The strength of the synthetic method is that properties of the data, such as noise, propagation jitter, and rate of spike superposition, can be systematically varried, and the results of sorting the data for each choice of these praameters can be exactly quantified. This approach can be very useful for establishing bounds within which a particular algorithm is effective. It my then be possible to estimate equivilent parameters for a real data set and determine whether they fall within the acceptable bounds.

The spike sorting package supports synthetic data methods by providing tools to estimate statistics of real waveforms for use in the development of synthetic data, and by providing the test_resuts script for checking the accuracy of a sorting attempt on synthetic data. Test_results also provides preliminary analysis of the mode of failure when a sorting attempt doesn't succeed. We are working to add tools for automatically generating synthetic test data.

Shift Jitter

Usually, the shift template for a particular unit is originally generated by calculating the optimal alignment for a particular event waveform. Once this shift template is used in the detection of a collection of events, it is possible to calculate how well the template matches the shifts of the detected units. Statistics of the shift distribution, including the variance of the optimal shifts for each detected unit, provide a measure of how accurately aligned the detected units are. For some discriminant functions, this can be a good indicator of detection errors. For example, for the mean discriminant, detection of strong but poorly aligned units will result in a dramatic increase in template shift jitter.

In addition to calculating shift jitter across events detected by a template, the set of events can be used to calculate a new shift template (for example the mean value across detected events, or the mean amoung events falling in one mode of a multimodal shift distribution). This "dejittered" template can be used to redetect spikes. Convergence of this detect-dejitter-redetect cycle indicates the ilimination of one class of detection error.

The spike sorting package supports calculation, display, and removal (dejittering) of shift jitter in interactive mode, can use dejittered shift templates in batch mode, and can display shift jitter statistics in the post-hoc analysis GUI

Waveform Analysis

Working ...

Custom Stimuli

In some preparations, the problem of spike separation can be simplified by using stimuli that selectively invoke responses from only one (or only a few) of the units, with known identity. This is possible if the prep contains, for example, cells that tile a stimulus space with their receptive fields, and the receptive fields of these cells are known (e.g. from previous intracellular experiments). Even if this much control is not possible, sending simplified or low amplitude (or low contrast) stimuli can reduce the number of active cells and the spike rates of these cells, leading to response data with few superpositions, which is comparatively easy to sort.

When this sort of recording is possible, the recorded responses to the rarified stimuli can be used in two possible ways:

Gui Features

Working ...

Dealing with superposition

Working

 

Last edit: 05/29/09

Index