INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F APPENDIX G APPENDIX H REFERENCES

The functioning of the brain depends on the concurrent activation of large populations of neurons. Yet most contemporary neurophysiological theories still focus on the individual properties of single neurons without consideration of the potential role played by the interaction among large neuronal groups. Recently electrophysiological and computational studies suggest that the precise timing (1-2 ms) of spikes in pre- and postsynaptic neurons may be used in neural networks to decipher information encoded in spike timing (Abeles et al. 1993; Bi and Poo 1998; Gerstner et al. 1996; Riehle et al. 1997). Therefore a technique for detecting precise spike timings from closely neighboring neuronal ensembles is much needed to characterize the dynamic interactions in populations of neurons involved in processing and storing information.

Activities from two or more closely neighboring neurons can be recorded simultaneously by only one common electrode. In fact, one electrode for extracellular recording often contains signals from more than one neuron around the electrode. To separate the single-neuron components from such recorded multi-neuronal activities, the ordinary spike-sorting methods use spike waveform differences of the single neurons. However, there are many situations in which different neurons generate action potentials having very similar shapes in the recorded waveform. This happens when the neurons are similar in morphology and about equally distant from the recording microelectrode. To avoid this problem, methods for recording with a small tetrode (Wilson and McNaughton 1993), constructed from four microwires insulated to the tips, have been developed. The relative amplitudes and waveforms of neuronal signals on the four recording wires are used to obtain single-neuron isolation. The idea is that multiple neurons are less likely to be equidistant from four microwires, which provide the minimal number of recording points necessary to separate sources (single neurons) based on the relative spike amplitudes and waveforms on different microwires.

Techniques to separate the individual sources using the different spike amplitudes and waveforms are referred to as "spike sorting" (for a recent review, see Lewicki 1998). Most current spike-sorting techniques assume that the noise and action potentials are Gaussian (stationary) and that there are few neurons that generate their action potentials completely simultaneously around the electrode. Concerning the Gaussian assumption, however, it is known that action potentials may have varying shapes (Fee et al. 1996a). Moreover, concerning the nonsimultaneity assumption, we have no way of knowing whether or not some neighboring neurons generate simultaneous action potentials because spike waveforms overlap on a common electrode when two or more neurons fire simultaneously, and the overlapped waveforms cannot be separated by the ordinary spike-sorting methods. In addition to these problems, the tetrode may integrate both axonal action potentials and dendritic action potentials under some conditions (Buzsáki et al. 1996; Cohen and Miles 2000). These problems limit the ability to sort multi-neuronal activities recorded with tetrodes and to detect precise spike timings among closely neighboring neurons around the microwires of tetrodes.

To overcome some of these problems, our previous work (Takahashi et al. 2002) has been devoted to the development of a procedure to separate single-neuronal activities from multi-neuronal recordings using independent component analysis (ICA) (Hyvärinen 1999). ICA determines what spatially fixed and temporally independent component activations compose an observed time-varying response without attempting to specify where in the cortex these activations arise. However, ICA has a serious limitation that the number of single neurons must be less than the number of electrodes, whereas it is possible that a tetrode (4 microwires) can detect more than four neurons from one recording site (Gray et al. 1995).

The present study combines the ICA and the efficiency of the ordinary spike-sorting technique (k-means clustering) to solve the nonstationarity and the spike-overlapping problems with no limitation on the number of single neurons to be separated. Here we report an automatic procedure using an innovative method for spike sorting that is applicable to multi-neuronal data recorded with tetrodes in the working brain.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F APPENDIX G APPENDIX H REFERENCES

Overall procedure

Our procedure has three steps, i.e., clustering the spike waveforms by k-means clustering, decomposition of each clustered spike waveform by ICA, and aggregation of the decomposed clusters to reconstruct single neuronal activities. Before the first step, we divide the condition of the spike waveforms into stable and unstable forms (Fig. 1).

View larger version (27K): [in this window] [in a new window]   Fig. 1. An example of waveform variability on clustering of spike waveforms. The dots represent the distribution of spike waveforms of a real single neuron. The 2 coordinates are the maximum and minimum amplitude of spike waveforms. If the cluster is represented by a single mean, an example of cluster boundary produces the area outlined by the dotted ellipsoid. We define the spike waveforms in that area (red dots) as being in the stable condition. The residual spike waveforms (blue dots), which may be influenced by overlapping spikes, complex bursts and/or other nonstationary factors, are from the same neuron. We define the residual spike waveforms (blue dots) as being in the unstable condition.

In the stable condition, spike waveforms show isotropic Gaussian variability, and their density in the feature clustering space is large enough that we expect to be able to recognize and sort them (Fig. 1, red dots). In the unstable condition, there are a small number of spikes whose waveforms are too irregular to cluster them by ordinary spike-sorting techniques (Fig. 1, blue dots). For this reason, our procedure requires two phases: detection and determination of the number of firing neurons and their relative position on different microwires under the assumption of the stable condition and detection of residual spikes under the assumption of the unstable condition. We give a detailed description of each phase and step in the algorithm. A flow chart summarizing these phases and steps is shown in Fig. 2.

View larger version (25K): [in this window] [in a new window]   Fig. 2. Flow chart of the automatic classification algorithm.

Spike detection

First of all, candidate events, where one neuron probably fired, are identified in the tetrode recordings. The root-mean-square (RMS) amplitude in each channel, an upper bound on the noise power, is calculated, and each threshold level is set to three times the RMS amplitude in each channel (50-80 µV under the conditions used in our laboratory). When a threshold crossing is detected on either microwire of a tetrode pair, 32 samples of the waveform from each of the four microwires of the microelectrode are saved. Both waveforms are centered with respect to the particular waveform that first crosses the threshold level. The voltage sample with the largest amplitude is set as the fifth sample in the stored waveform, and a time stamp saved with each waveform indicates the time of the waveform absolute peak with a 40-µs resolution. We deal with these segmented sets of spike waveform and its time stamp throughout the following procedure.

Detection and decomposition of the stable waveforms

In the first step, the segmented waveforms are sorted into a large number of clusters, typically two times the expected number of neurons, using a feature clustering technique, the so-called k-means algorithm (see APPENDIX A). For computational efficiency, the features used to perform the classification are maximum and minimum spike amplitude. All spike waveforms are then classified into these initial clusters. At this point, the problem of clustering the entire data set is reduced to one of aggregating an initial set of clusters, typically 20 in our examples, into a final set of clusters containing spike waveforms from the same single neurons.

Due to the possibility that the tetrode records not only axonal but also dendritic action potentials and the possible presence of dendritic backpropagating spikes, there may be many similar combinations of action potentials in a segment window. Hence, we must decompose the clustered spike waveforms into single neurons in this phase. ICA yields data decompositions since the spatially stable and sparsely active components sum to the observed multichannel response (see APPENDIX B). ICA determines what spatially fixed and temporally independent component activations compose an observed time-varying response without attempting to specify where in the cortex these activations arise from. Data from N microwires can be reconstructed as the weighted sum of the N independent components. Our previous work (Takahashi et al. 2002) indicated that FastICA (Hyvärinen 1999) requires at least as many microwires as firing neurons to separate single-neuronal activities. Although the number of neurons in an entire data set is unknown, we assume that each initial cluster sorted at this point has a putative single neuron. Thus in the second step, the spike waveforms of each initial cluster are successfully decomposed into single neurons by FastICA under the assumption described in the preceding text.

We reconstruct each of the four decomposed components of tetrode data to identify individual spike sequences (see APPENDIX C). This means that the tetrode data in one initial cluster are expanded to four tetrode data containing only one single-neuronal activity. The decomposed and reconstructed waveforms are upsampled in the segment over each threshold, and such segments are then extracted. The threshold used is twice the RMS amplitude of the entire segmented waveform. We assume that the extracted segments of each component are the same putative single neurons. Each of the decomposed data sets has an independent component basis vector (ICBV) (see APPENDIX D). If there are observed neurons at a close position, these ICBVs must be nearly equal. This means that if the decomposed data have similar ICBVs, they may be the action potentials generated from one neuron. Then in the third step, the decomposed data are automatically aggregated by means of the comparison of ICBVs between the decomposed clusters, M_ab (Eq. D1), as a measure of the similarities of the distribution of the spike waveform relative to the different microwires (see APPENDIX D). A threshold level is set at 0.99.

Detection and decomposition of the unstable waveforms

The unaggregated data of the first phase contain the unstable spike waveforms including overlapping spikes. In this phase, we deal with those data sets that are likely to give inaccurate results if spike waveforms were modeled under the assumption of Gaussian mixtures. In the first step, we take account of irregular spike waveforms caused by overlapping spikes. The larger the number of clusters that are sorted, the smaller the proportion of overlapping spikes in the clusters that must be considered. Therefore the segmented waveforms in the unaggregated data are sorted into a large number of clusters, typically four times the expected number of neurons, using k-means clustering. At this point, although we cannot sort the segmented waveforms into individual clusters containing only one single-neuronal activity, we can assume that each cluster contains a combination of similar single-neuronal activities. In the second step, the spike waveforms of each cluster are decomposed into several components containing a putative single neuron, using FastICA. As shown in Detection and decomposition of the stable waveforms, the decomposed and reconstructed waveforms are upsampled in the segment over each threshold, twice the RMS amplitude of the entire segmented waveform, and such segments are then extracted. We assume that the extracted segments of each component are the same putative single neurons. In the third step, the decomposed data are automatically aggregated by means of the comparison of ICBVs between the decomposed clusters, M_ab (Eq. D1), as a measure of the similarities of the distribution of the spike waveform relative to the different microwires (see APPENDIX D). A threshold level is set at 0.9. Then to reduce the combination of single neurons in the feature clustering at the first step, the decomposed and reconstructed spike waveforms identified as putative single neurons at the second step are removed from the original waveforms, and we return to step 1 until no more spikes are detected. Typically, 50-150 iterations are needed.

Aggregating the stable and unstable waveforms into single-neuron clusters

The remainder of the procedure deals with the aggregation of these clusters in the stable and unstable conditions into the more global structure of single-neuron clusters. Each of the clusters detected in the second phase is assigned to a cluster whose ICBV is the nearest one among all of the clusters detected in the first phase, using the ICBV of clusters in both stable and unstable conditions, M_ab (Eq. D1), as a measure of the similarities of the distribution of the spike waveform relative to different microwires (see APPENDIX D). A threshold level is set at 0.9. When all possible pairs of clusters have been rejected, the algorithm is terminated. Due to the anisotropic non-Gaussian variability of spike waveforms and overlapping spikes, the aggregation will often not be adequate after the first calculation of the process. The unaggregated data are usually re-estimated by two to five iterations of the present procedure. Each of the aggregated clusters are simply added to the previous clusters in the first phase using ICBVs. Using the M_ab derived by ICA, which is a common measure through the first and second phases, we can detect precise single neurons in the presence of overlapping spikes and non-Gaussian variability of spike waveforms.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F APPENDIX G APPENDIX H REFERENCES

Examples of sorted single neurons

We focus on a data set which contains roughly 10⁵ spike waveforms collected from a tetrode in 1 h (see APPENDIX E). The average signal-to-noise ratio of the data set is around 4.1. The entire procedure is under control of the host computer (AMD Athlon 1.2 GHz processor and 1.5 GB Memory) and takes about 12 h. The program was implemented in MATLAB (Mathworks, Natick, MA) and the C programming language.

Figure 3 shows results from the procedure in Detection and decomposition of the stable waveforms. The superimposed spike waveforms, the auto-correlation functions, and two three-dimensional (3-D) scatter plots of maximum-to-minimum spike amplitude for 4 of the 10 aggregated clusters in the first phase for this data set are shown to be single neurons. As our goal in this phase is to sort the stable waveforms into single neurons, most of the superimposed waveforms must be similar. We identify all clusters as single neuronal activities, whose waveforms are similar and whose auto-correlation functions indicate a clear refractory period (1-2 ms). The scatter plots reveal clearly separated clusters reflecting the different amplitude ratios of the different neurons recorded on each of the four microwires. These results demonstrate that the present procedure in this phase can successfully sort the stable spike waveforms into single neurons.

View larger version (40K): [in this window] [in a new window]   Fig. 3. Performance in the stable condition. A: representative spike waveforms for 4 of 10 clusters identified in the stable condition. The waveforms recorded on each microwire of the tetrode are shown. B: corresponding auto-correlation functions. The bin width is 0.1 ms. Note that the clear refractory period (1-2 ms) of all of the auto-correlation functions indicates that they are all single neurons. C and D: 2 aspects of the spatial projection of the same sets of clusters for the same neurons as in A and B. In C, V1 is the minimum amplitude of channel 1 of the tetrode, P2 is the maximum amplitude of channel 2 of the tetrode, and V3 is the minimum amplitude of channel 3 of the tetrode. In D, P1 is the maximum amplitude of channel 1 of the tetrode, P2 is the maximum amplitude of channel 2 of the tetrode, and P3 is the maximum amplitude of channel 3 of the tetrode.

Figure 4 shows results from the procedure in Detection and decomposition of the unstable waveforms. Two typical segments having overlapping spikes are shown [Fig. 4, A (O1) and B (O2)]. The present procedure estimates that both O1 and O2 have a combination of two putative single-neurons, cluster 2 (S1) and cluster 4 (S2; Fig. 4C), and identifies the decomposed and reconstructed waveforms, E1_1 and E2_1, and, E1_2 and E2_2 (Fig. 4, A and B), as S1 and S2, respectively.

View larger version (33K): [in this window] [in a new window]   Fig. 4. Performance in the unstable condition. A and B: O1 and O2 are overlapping spikes recorded on each microwire of the tetrode. E1_1 and E1_2 are waveforms decomposed and reconstructed from O1. E2_1 and E2_2 are waveforms decomposed and reconstructed from O2. E1_1 and E2_1 are aggregated to cluster 2 in Table 1. E1_2 and E2_2 are aggregated to cluster 4 in Table 1. C: the average waveforms of clusters S1 and S2, which were recorded on each microwire of the tetrode. D: one aspect of the projections of clusters S1 (red dots) and S2 (green dots) from the decomposed and reconstructed waveforms of E1_1, E1_2, E2_1, and E2_2. The two coordinates are the maximum spike amplitude of channel 2 (P2) and the minimum spike amplitude of channel 1 (V1). Note that the estimated and reconstructed waveforms E1_1 and E2_1, and E1_2 and E2_2 are located approximately on the projections of cluster S1 (red dots) and cluster S2 (green dots), respectively.

Study of simultaneous paired intracellular and extracellular recording (Harris et al. 2000; Henze et al. 2000; Wehr et al. 1999) suggests that the electrical linearity of the extracellular medium is valid. From this point of view, the decomposed and reconstructed waveforms were expected to accurately predict the shape of source spikes. The decomposed and reconstructed waveforms, E1_1 and E2_1, and, E1_2 and E2_2 approximately match the shape of the average waveforms, S1 and S2, respectively. Their relative spike amplitudes on different microwires are similar. Moreover, Fig. 4D shows that the spatial projection of clusters 2 and 4 contain approximately the decomposed and reconstructed waveforms, E1_1 and E2_1 and E1_2 and E2_2, respectively. These results demonstrate that the procedure in the second phase can decompose overlapping spikes accurately. Note that the absolute peak of channel 2 of the decomposed waveform E2_1 is appropriately time shifted (Fig. 4B, blue arrow). Because such delay must be added to the time stamp of the spike segments, we reset the time stamp of the overlapping spike segments using the time of the absolute peaks of the decomposed and reconstructed waveforms.

Proportions of overlapping spikes

Table 1 shows the final result classified using the present procedure. All of the 10 estimated clusters contain overlapping spikes, and there are various proportions of overlapping and total spikes in each estimated cluster. Neurons with very low firing rates such as clusters 9 and 10 may be silent cells (Thompson and Best 1989) and often tend to be deleted from the database in studies using ordinary spike-sorting techniques (Henze et al. 2000). However, although the proportions of the overlapping and total spikes of clusters 9 and 10 are higher than the other estimated clusters, we did not delete these clusters from the database. We discuss the meaning of the high overlapping rate of these clusters in the DISCUSSION.

The procedure presented here has been used to classify spike waveforms from 20 sets of data obtained from tetrodes implanted in the prefrontal cortex of two monkeys during delayed matching-to-sample tasks (Sakurai 2001; Sakurai et al. 2001). An average of 4.9 ± 0.7 (mean ± SE) putative single-neuronal activities were isolated per tetrode. There were various proportions of overlapping spikes in all the data sets as shown in Table 1.

Examples of variable sorted clusters

The present procedure can deal with data of varying quality and reveals variable types of sorted clusters. Figure 5 shows examples of variation in the quality of the sorted clusters, i.e., the scatter plots of the maximum and minimum amplitude of spike waveforms in the feature clustering space. Figure 5, A-C, is from different neuronal data obtained at different recording sessions. Figure 5, left, shows all triggered wave functions including both spikes and noise. Figure 5, right, shows clusters sorted from the data on the left. Figure 5A is from the data shown in Table 1 and shows four clusters (Clusters 1, 3, 5 and 8) on the right. Figure 5, B and C, right, shows two clusters separated from the left. These clusters in Fig. 5, right, are only clusters that can be clearly plotted in the present two-dimensional clustering space. The present procedure detects more clusters from the data than those shown on the right (e.g., see Table 1 for Fig. 5A) but all separated clusters are in the higher dimensional spaces using ICA and cannot be clearly plotted in a two-dimensional space. This is a difference from the ordinary cluster-cutting procedure. The separated clusters in Fig. 5A, right, have some outliers that may be assigned to the different clusters or may be discarded as noise if ordinary manual or automatic spike-sorting methods were used. We will discuss the meaning of such spikes with irregular amplitudes in the DISCUSSION. The other examples of sorted clusters obtained from different neuronal data (Fig. 5, B and C, right) are well-separated in the feature clustering space. There are various types and qualities of neuronal data showing various distributions of these spherical and elliptic clusters on the feature clustering space.

View larger version (44K): [in this window] [in a new window]   Fig. 5. Examples of sorted clusters. These data (A-C) are obtained from different neuronal data at different recording sessions. Left: the dots represent the distributions of all triggered wave functions including spikes and noise. Right: the dots represent the distributions of spike-waveform functions of separated clusters sorted by the present procedure. The dots of different colors represent separated different clusters. The two coordinates are the maximum and minimum amplitude of spike waveforms. Comparing the left with the right indicates that how clearly the clusters are away from noise. A, right: 4 clusters sorted from the data shown in Table 1 by the present procedure. Clusters 1 (red dots), 3 (green dots), 5 (blue dots), and 8 (black dots) are plotted. The outliers (solid and red dotted arrows) can be seen. B and C, right: 2 spherical or elliptic clusters in each clustering space. The clusters in A-C, right, are only clusters that can be clearly plotted in the present 2-dimensional clustering space. The present procedure detects more clusters from the data but all separated clusters are in the higher dimensional spaces using ICA and cannot be clearly plotted in a 2-dimensional space.

Examples of functional connectivity between single neurons

To test the functional connectivity among neighboring neurons, we calculated the shuffled-corrected cross-correlogram (see APPENDIX F) between two putative pyramidal neurons (Fig. 6). Clusters 3 and 4 in Fig. 6 are determined using the combination of the parameters (firing rate, spike shape, and the mean of the auto-correlograms) (Csicsvari et al. 1998). The large peak of the cross-correlogram (Fig. 6B), which cannot be verified by the ordinary spike-sorting techniques without overlap decomposition, indicates excitatory shared inputs to the neurons (Fig. 6A), demonstrating that the present procedure can detect and decompose the overlapping spikes for real data sets.

View larger version (12K): [in this window] [in a new window]   Fig. 6. Examples of shuffling-corrected cross-correlograms (see APPENDIX F for details of the calculation) from a pair of pyramidal neurons. A shuffling-corrected cross-correlogram represents the number of spikes per bin that occurred in one neuron before and after spikes in another neuron. The ordinate values are averaged numbers of spikes per bin. The bin width is 0.5 ms. To test that peaks in the shuffling-corrected cross-correlograms are statistically significant, a band of 99.5% confidence limits (Abeles 1982) for the equivalent independent Poisson processes is shown (· · ·) for each cross-correlogram. A: the estimated synaptic connections functioning between the pyramidal neurons. B: a shuffling-corrected cross-correlogram between the pyramidal neurons (clusters 3 and 4) with overlapping spikes. C: a shuffling-corrected cross-correlogram between the pyramidal neurons (clusters 3 and 4) without overlapping spikes.

To understand the effect of overlap decomposition using the present procedure for detecting the precise spike timings among neighboring neurons, we artificially made the cross-correlogram shown in Fig. 6C that subtracted the overlapping spikes from the cross-correlogram shown in Fig. 6B. Consequently, the cross-correlogram shown in Fig. 6C is constructed by the spike-sorting technique without overlap decomposition, and its classification for the nonoverlapping spikes is as accurate as the present procedure. Figure 6, B and C, indicates excitatory shared inputs such as that shown in Fig. 6A. However, the cross-correlogram shown in Fig. 6C misses 104 spikes as compared with the cross-correlogram shown in Fig. 6B, and the missing spikes are overlapping ones. Consequently, the results of Fig. 6 demonstrate that, unlike the spike-sorting techniques without overlap decomposition, the present procedure can be used to examine the precise functional relationship among neighboring neurons for real data sets.

Detection of monosynaptic connections between pyramidal neurons and interneurons

Recent in vitro and in vivo experiments in several brain regions have suggested that short-latency peaks may reflect a monosynaptic connection between the pyramidal neuron and the target interneuron (Ali et al. 1998; Cohen and Miles 2000; Csicsvari et al. 1998; Galarreta and Hestrin 2001; Krimer and Goldman-Rakic 2001). Note that Buzsáki and colleagues have confirmed such monosynaptic connections by making simultaneous intracellular and extracellular measurements (Buzsáki et al. 1996; Henze et al. 2000; Kamondi et al. 1998). To examine such monosynaptic connection, we divide the estimated clusters into pyramidal neurons and interneurons by using the combination of the parameters described in the preceding text.

The shuffling-corrected cross-correlogram between the putative pyramidal neuron (cluster 1, Fig. 7A, left) and the putative interneuron (cluster 6, Fig. 7B, right) shows a significant peak at monosynaptic latency (2 ms; Fig. 7B). This latency corresponds to that found in in vivo and in vitro studies (Cohen and Miles 2000; Csicsvari et al. 1998; Galarreta and Hestrin 2001; Krimer and Goldman-Rakic 2001), demonstrating that the present procedure can detect precise spike timings among neighboring neurons not only for the overlapping spikes but also for the nonoverlapping spikes.

View larger version (17K): [in this window] [in a new window]   Fig. 7. An example of auto-correlograms and shuffling-corrected cross-correlograms from a pyramidal neuron and an interneuron. All parameters and symbols are as in Fig. 6. A: auto-correlograms of the pyramidal neuron (cluster 1) and the interneuron (cluster 6). The lower portion shows the estimated synaptic connections functioning between the neurons. B: a shuffling-corrected cross-correlogram between the pyramidal neuron (cluster 1) and the interneuron (cluster 6). Note that a large and sharp peak is shown at 2 ms.

Performance on synthesized data set I

To explore the strengths and limitations of the procedure, we ran two mixtures of synthesized data sets composed of the six template spike waveforms extracted from the real data set, including burst firing, plus Gaussian noise (see APPENDIX G). One of the two synthesized data sets contains no overlapping spikes, and the other contains overlapping spikes, each proportion of which in true clusters 1-6 is about 5, 10, 10, 10, 50, and 90%, respectively. The average signal-to-noise ratio is around 4.3. The synthesized waveforms were automatically classified by the Expectation Maximization (EM) algorithm (KlustaKwik) (Harris et al. 2000) and the present procedure, respectively. KlustaKwik generates a list of occurrence times, which were compared with the known positions of neurons that generate action potentials, for each cluster. For KlustaKwik, the features used to perform the classification were maximum spike amplitude, minimum spike amplitude and first and second principal components derived by PCA (Abeles and Goldstein 1977). Taking account of the influence of the classification result, we used six as the maximum number of clusters which must be set by a human operator. The present procedure was run under the same conditions as in the preceding text.

The results of classifying the synthesized data set are shown for the nonoverlapping spikes in Table 2 and for the overlapping spikes in Table 3. Perfect performance would have all zeros in the off-diagonal entries and no undetected events. A spike can be missed if it is not detected in an overlap sequence or if all its sample values fall below the threshold for spike detection.

For the nonoverlapping spikes, the performance of the present procedure is nearly perfect (Table 2, right). For the performance of KlustaKwik, true clusters 4 and 6 are collapsed into an estimated cluster 4 (Table 2, left). Due to the synthesized factor of complex bursts, these clusters include some spikes, which have similar relative spike amplitudes. For practical use, because KlustaKwik consists of a semi-automatic process followed by examination and reassignment by a human operator using the auto-correlogram and the cross-correlogram (Harris et al. 2000), this estimated cluster 4 may be detected and reclustered to true clusters. The result of Table 2 indicates that the present procedure can classify multi-neuronal recordings into single-neuronal activities accurately and is valid for reducing the inspection process by a human operator.

Table 3 indicates that for the overlapping spikes, true clusters 1 and 6 are accurately identified and classified by KlustaKwik. However, true clusters 2 and 4 are collapsed into an estimated cluster 3 and half of the true cluster 5 is identified as cluster 5 (Table 3, left). As these incorrect clusters are not well separated in the feature space due to the overlapping spikes, it is difficult to choose the correct clusters without overlap decomposition. In contrast, the accuracy of the present procedure is more than 90% (Table 3, right). For the two worst estimated clusters 2 and 5, there are some samples that would be expected to exceed the noise level. As the threshold for spike detection is lowered, there is a trade off between the number of real spikes missed and the number of false positives resulting from common instances when the noise contains a spike-like shape. Therefore for the real data sets, the threshold for spike detection must be chosen based on a careful examination of many data sets. The result of KlustaKwik suggests that there are significantly more missed spikes, and only one spike in a segment window without the technique of overlap decomposition. Thus for real data sets including some overlapping spikes, the present procedure is needed to detect precise spike timings among neighboring neurons.

Performance on synthesized data set II

To quantify the performance of detecting overlapping spikes from more realistic and complex data, we made another synthesized data set composed of two linearly mixed signals each of which is respectively derived from one microwire of two different tetrodes obtained at different recording sessions, plus noise (i.e., signals having no spikes, recorded from another tetrode; see APPENDIX H). Because, in the synthesized data set, we deal with a signal derived from one microwire of a tetrode that must have multiple spike waveforms as a signal from a single neuron, two single neurons identified from the synthesized data set have several spike shapes, that is, irregular spike waveforms. Moreover, provided that overlapping spikes are correctly detected two single neurons detected from this synthesized data set containing two independent neurons would produce, due to chance, a flat shuffled-corrected cross-correlogram.

The synthesized waveforms were automatically classified by KlustaKwik and the present procedure run under the same conditions, as shown in the previous subsection, except that for KlustaKwik we used two as the maximum number of clusters that must be set by a human operator, taking account of the influence of the classification result.

The result of classifying the synthesized data set is shown in Table 4. For the performance of KlustaKwik, true clusters 1 and 2 are collapsed into estimated clusters 1 and 2, and only half of the overlapping spikes are identified as estimated cluster 2 (Table 4, left). In contrast, the accuracy of the present procedure is more than 99% (Table 4, right). The shuffling-corrected cross-correlogram between two single neurons detected by the present procedure shows no peaks near the 0-ms bin (Fig. 8A); this corresponds to our expectation. Two error rates, proportions of undetected overlapping spikes to the total overlapping spikes, indicate 0.90 and 0.75, respectively. In contrast, the result of KlustaKwik, a spike-sorting technique without spike decomposition, shows a negative peak around 0 ms bin (Fig. 8B, arrow). These results demonstrate that the present procedure can detect overlapping spikes and provide more precise functional connectivity than spike-sorting techniques without overlap decomposition in realistic, complex cases.

View larger version (16K): [in this window] [in a new window]   Fig. 8. Shuffling-corrected cross-correlograms from a pair of the estimated clusters from the synthesized data set II. All parameters and symbols are as in Fig. 6. A: a shuffling-corrected cross-correlogram between clusters 1 and 2 estimated by the present procedure. Note that no peak is found near 0 ms. B: a shuffling-corrected cross-correlogram between clusters 1 and 2 estimated by KlustaKwik. Note that a sharp and negative peak is obtained at 0 ms.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F APPENDIX G APPENDIX H REFERENCES

Significance of the present procedure

The procedure reported in this paper significantly reduces the two largest problems of spike sorting, i.e., overlapping spikes and the non-Gaussian variability of spike waveforms. These problems limit the usefulness of multi-neuronal recording for detection of precise spike timing interactions among neighboring neurons. Although some techniques have been developed to permit sorting in the presence of overlapping spikes (Lewicki 1994) and non-Gaussian variability of spike waveforms (Fee et al. 1996b), no previous technique can completely overcome both problems. The present procedure completely solves both problems and provides automatic and reliable spike sorting for multi-neuronal data recorded with tetrodes.

There are two distinguishing features of our procedure and its implementation. First, we make no a priori assumption about the form of the variability in waveforms, allowing us to effectively sort spike waveforms into putative single-neuron clusters in the presence of overlapping spikes and realistic distributions of neural noise. Second, we incorporate ICA directly into the algorithm. This is particularly valuable for the identification of waveforms from a putative single neuron that produces overlaps and bursts of spikes with significantly different shapes. These features of the present procedure should overcome the limitations of multi-neuronal recordings previously encountered. For example, the present procedure can precisely detect brief bursts of high-frequency firing, which may reliably transmit signals to postsynaptic targets (Lisman 1997) and may increase the possibilities of generating overlapping spikes and non-Gaussian variability of spike waveforms. Possibilities concerning the significance of overlapping spikes are also indicated by the theory of synfire chains (Abeles 1991), spike timing-dependent plasticity (Bi and Poo 1998), and the occurrence of dendritic spikes, which may have special importance in neuronal networks.

The present procedure finds immediate application in multi-neuronal recordings from awake animals performing behavioral tasks. This approach is crucial to further understanding of actual working modes of the functional relationships among neighboring neurons and their dendrites, such as those suggested in the theory of "cell assemblies" (Hebb 1949; Sakurai 1996, 1999), which have previously been very difficult to observe. Some previous recording studies (Sakurai 1993, 1996, 2003) suggested that the shuffled-corrected cross-correlogram between two pyramidal neurons recorded from two different electrodes about 200 µm apart sometimes showed functional connections between the neurons; moreover, those connections often change when rats perform different behavioral tasks. It is, however, difficult to detect such task-dependent dynamic connectivity among closely neighboring neurons located within about 100 µm because of the presence of overlapping spikes in the same electrode. Similarly, even if we could record neurons from two different electrodes less than 100 µm apart, then the same action potential of one neuron might be detected from the two different electrodes. Even if we could take account of simultaneous paired intracellular and extracellular recording from two different electrodes (Cohen and Miles 2000; Henze et al. 2000), then axonal and dendritic action potentials of the same neuron might be confused in the recording. Therefore a sorting procedure with overlap decomposition, such as that presented here, is absolutely necessary to detect precise functional connectivity among closely neighboring neurons.

Outliers in the clustering space

The ordinary spike-sorting methods discard outliers, i.e., data not assigned to spherical clusters in the feature clustering space, but the present technique identifies some outliers as overlapping spikes. For instance, the sparse and spherical cluster located at the top right of Fig. 5A, right (indicated by an arrow), is identified as overlapping spikes of two clusters (clusters 1 and 3). Spike sorting with Gaussian mixture models cannot partition the overlapped clusters in the feature clustering space into single-neuron clusters, but the present procedure identified the overlapped cluster located at the top left of Fig. 5A, right (indicated by a red dotted arrow), as overlapping spikes of two clusters (clusters 5 and 8). The present study indicates that overlapped clusters and outliers have a high probability of being due to overlapping spikes.

The ICA, a core thread of the present procedure, finds independent components that potentially generate the neuronal data, but this does not necessarily mean that these components form distinct clusters such as shown in Fig. 5A, right. The other sorted clusters obtained at different recording sessions (Fig. 5, B and C, right) indicate that there are well-separated clusters in the feature clustering space because these clusters include extremely fewer proportions of overlapping spikes than the data shown in Fig. 5A, right. These results demonstrate that when the data have nonoverlapped and spherical clusters in the feature clustering space, the present procedure can sort single-neuronal components as clearly as the ordinary spike-sorting methods. The present procedure has the advantage that it can provide insights into the data that the ordinary spike-sorting methods discard as outliers or noise and can detect overlapping spikes from the data. It should be noted, however, that outliers can be due to the nonstationarity of spike waveforms in some neuronal data and, as shown in Fig. 5A, right, the present procedure identifies some of the data as overlapping spikes.

Tetrode may integrate both axonal spikes and dendritic spikes

In the experiments using the combination of intradendritic, intracellular, and extracellular recordings in vivo, Buzsáki and colleagues suggested that the tetrode might be used to record action potentials not only from axons but also from dendrites (Buzsáki et al. 1996; Henze et al. 2000; Kamondi et al. 1998). In addition, some vitro experiments have suggested that an axonal action potential was followed by a backpropagating action potential in neocortical pyramidal neurons (Stuart et al. 1997). Due to the propagation delay in the dendritic tree, the waveforms of the axonal action potential and its backpropagation, overlapped in the extracellular tetrode in a segment window, may result in overlapping spikes that the present procedure can detect and decompose. Moreover, some in vitro experiments (Golding and Spruston 1998) have also shown that dendritic action potentials can trigger axonal action potentials. This relationship of dendritic and axonal action potentials may also result in overlapping spikes in a segment window. From these experimental facts and interpretations, it is possible that cluster 9, of which about 92% of the spikes were overlapped with only the spikes of cluster 4 (Table 1), had action potentials that occurred from a part of a dendrite of the neuron of cluster 4.

Importance of assessing the roles of precise spike timings

Recent studies on the precise timing of spikes have suggested that there might be no special role of spikes with small variations in their firing rate (Baker and Lemon 2000; Brody 1999a,b; Oram et al. 1999). However, the data of spike trains used for those measurements were produced by the ordinary spike-sorting technique without overlap decomposition. The results from the synthesized data sets including overlapping spikes and complex spikes in the present study (Table 3) imply that the ordinary spike-sorting technique without overlap decomposition often missed significantly more spikes than the present procedure when real data include many overlapping spikes. Moreover, the functional connectivity shown in Fig. 6A can be driven from a single bursting neuron where some portions of spikes of the burst are artificially placed in different clusters separated by the ordinary spike-sorting techniques without overlap decomposition. The present procedure can detect overlapping spikes that do not occur from a single neuron, suggesting that the functional connectivity shown in Fig. 6B indicates excitatory shared inputs to the two different neurons.

The single-neuronal activities detected by our procedure from real data sets suggest that the second largest proportion of overlapping spikes to total spikes in the monkey prefrontal cortex during task performance is around 45% (Table 1). Such a high rate of overlapping spikes, produced in a few of the single neurons, cannot be due to dendritic spikes, as discussed in the preceding text, and have not been detected in previous studies that have examined spike timings. Recent electrophysiological and computational studies suggest that the precise spike timings may be used to decipher neuronal information (Abeles et al. 1993; Bi and Poo 1998; Gerstner et al. 1996; Riehle et al. 1997). It is, therefore possible that there are special roles of precise spike timings among closely neighboring neurons, as detected by the present procedure.

Further improvement of the present procedure

Due to the limitations of ICA, we assumed in the present study that no more than four neurons generated action potentials simultaneously in a segment window. Because there is the theoretical possibility that each segment of tetrode recording contains more than four neurons that fire simultaneously, in future studies, we should consider algorithms that find overcomplete representations (Amari 1999; Lewicki and Sejnowski 2000), although such techniques have not yet been successfully applied to real neurobiological data. On the other hand, our procedure can be extended to the analysis of records with larger numbers of simultaneously acquired channels than that of the tetrode (4 channels). In such cases, a given neuron can only contribute to the output of a small fraction of the total number of microwires. This will be useful to increase the number of decomposable action potentials in a segment window in the present procedure and improve the accuracy of spike sorting.