Brain Research, 169 (1979) 455-479 © Elsevier/North-HollandBiomedicalPress

455

COMPUTER SEPARATION OF UNITARY SPIKES FROM WHOLE-NERVE RECORDINGS

CLAUDE CAMP and HAROLD PINSKER Marine Biomedical Institute, Departments of Physiology and Biophysics, and Psychiatry and Behavioral Sciences, University of Texas Medical Branch, Galveston, Texas 77550 (U.S.A.)

(Accepted October 19th, 1978)

SUMMARY A practical and efficient off-line computer technique is described for automatically separating unitary waveforms from multiunit whole-nerve spike train data with no a priori knowledge about the number of units or their waveforms. The procedure requires two recording electrodes which provide 3 measures on each putative unit: (1) peak-to-peak amplitude on the proximal channel, (2) amplitude on the distal channel, and (3) temporal offset (depends on conduction velocity) between proximal and distal spikes. On the basis of these 3 measurements, individual unitary spikes are automatically separated into clusters according to empirically-determined limits of variability. The results of the program are displayed in 3-D plots of the 3 measures on each unitary spike and in plots of superimposed waveforms from each cluster. These plots can be used to interactively correct clustering errors. The procedure is illustrated with a 1-min segment of spike train data recorded in vivo from the siphon nerve of a freely-behaving Aplysia. We routinely obtain about 10 relatively well-isolated units in such segments. By utilizing the average waveforms and conduction velocities for individual clusters, it may eventually be possible to separate unitary spikes from compound waveforms resulting from simultaneous firing of two or more units.

INTRODUCTION Studies of both vertebrates and invertebrates have classically utilized nerve fascicles to isolate and study single axonsl,L However, to analyze behavior it is necessary to examine concurrent activity in large populations of neurons s. Extracellular action potentials recorded from heterogeneous neuronal structures such as dendrites and somata can be complex, whereas axons generate simple waveforms that are easier to analyzeiv. Many functionally related axons can be monitored simultaneously with whole-nerve electrodes. Furthermore, such electrodes can be implanted to

456 record in vivo from peripheral nerves or central connectives in freely-behaving animals 3,5,z6. In all types of multicellular recording, the individual units must be separated to specify the activity in the neuronal network underlying the behavior. Although a number of approaches have been developed s, separation of action potentials from multiunit spike trains remains a basic problem in neurobiology. Most spike separation procedures are based on single channel recordings 7. The simplest technique is a window discriminator in which a single threshold can separate out a single unit on the basis of spike amplitude. Several windows can be used as discriminatorsl~, 25 if all units have dissimilar amplitudes. However, 2 units with similar amplitudes will not be separated accurately, even if their waveforms are different. This limits the usefulness of window discriminators except when a small (1-4) population of distinct units is present. A better single-channel technique is to take multiple measurements of individual waveforms, e.g., amplitudes, areas and durations of selected portions of the spikes. Two or more variables can be used as orthogonal axes to plot each occurrence of a spike as a single point. Two assumptions underlie this procedure: (1) repeated occurrences of the same unit will tend to cluster, and (2) the clusters corresponding to different units will not overlap when the right variables are selected. An automated procedure can then be used to detect spikes occurring within specified boundaries v,l t, 12,14.

Cross-correlation has also been used for spike separation from single channel recordings. Identification of a spike can be made by cross-correlating its waveform with an amplitude-normalized matched filter which is tuned to a known nerve impulse22. A large peak in the correlation affirms similarity. However, the filter contains only waveform information, and amplitude is not used as a discriminator. Hence, units of similar shape but different amplitudes can not be sorted by this method. Recent spike separation procedures have utilized multiple channel recordings. One correlation technique involves recording on two channels to take advantage of differencesin conduction velocityand has been used to separate up to 6 unitsI°. This is an effective and fast separation technique; however, units with similar conduction velocities cannot be discriminated even if they have differentamplitudes. Another technique involves recording on up to 6 channels2L If the waveformof each unit of interest is known in advance, then a digitallinear filtercan be tuned to any of the units, and when convolved with a segment of data will produce a sharp peak whenever the unit fires. If the number of recordingchannels approximatesthe number of units, then separation can be accomplished with great accuracy, regardless of superposition; however, as the number of recording channels decreases, the procedure becomes less effective. Although this procedure has been used successfullyin vitro where it is possible to record from a singlenerve at manydifferentsites, it is of limited practical value for in vivo recording. An efficientspike separation procedure would discriminatea maximumnumber of units for a minimum number of recording channels and measurements per spike. We have developed an automated, off-line, spike separation procedure than can

457 discriminate many different unitary spikes from whole nerve recordings. The use of two recording electrodes allows us to utilize conduction velocity, which is a highly reliable characteristic of a given axon. A high-speed algorithm is used to measure temporal offset between spikes on the two channels. This method requires no a priori knowledge of the number of units or their exact waveforms. The procedure relies heavily on differences of conduction velocity of neuronal units and, to a lesser extent, peak-to-peak amplitudes on the two recording sites. On the basis of these 3 measures, clusters of unitary spikes are automatically identified. A variety of interactive quality control checks and graphical displays can be used to detect and eliminate artifacts and errors. At present, this technique can only deal with unitary spikes of stable amplitude. However, the resulting clusters provide templates of the currently active units which might later be used to resolve more complex waveforms. A preliminary description of this procedure has appeared (Kanz, Camp and Pinsker27). METHODS To illustrate the spike separation procedure we shall utilize spike train data from the marine mollusc, Aplysia californica. The Aplysia nervous system is organized into ganglia which connect to the periphery and to other parts of the nervous system by means of long lengths of free nerves that are amenable to extracellular recording. Whole-nerve double cuff assemblies, consisting of two stainless steel monopolar electrodes~ separated by exactly 1 cm, were implanted on the siphon nerve of the abdominal ganglion in order to examine in vivo central neuronal activity underlying triggered and spontaneous movements of the siphonlS, 19 (Kanz, Cobbs and Pinsker is, 19,~s). This nerve consists of many hundreds of axons4 of different sizes and shapes varying in distance from the recording electrodes, and generating unitary triphasic action potentials of characteristic amplitudes and conduction velocities. When action potentials from two axons occur together, the currents summate algebraically. The proximal and distal channels (with respect to the ganglion) were filtered (60 Hz notch filter, Burr-Brown), amplified (WPI, DAM-5A), and recorded on an FM tape recorder (Ampex, SR-300). When playing the analog data into the computer, an adjustable high-pass filter (Krohn-Hite) was used to eliminate lowfrequency noise. Both channels were filtered at the same setting, to equate any resulting phase shifts. The filtering affects action potential amplitudes and waveforms but is necessary for DC stability. Voltage calibration pulses from the WPI amplifiers (d: 100 #V) were digitized and measured by an interactive procedure to obtain a conversion factor for each channel. Analysis of l-min segments of data was done by a DEC PDP 11/45 computer with 32K core memory, two R K l l movable-head disk drives, a Gould 5000 printer/plotter, and a DEC LPS11 with real time clock and 12-bit A/D and D/A converters. The analog records from each channel are first digitized. To avoid dealing with unreliable signals a noise threshold is established below which voltage fluctuations are ignored. Spike amplitudes are measured on each channel separately. Temporal offsets for spikes are measured between the two channels. Spikes with the same amplitudes

458 DIGITIZE 2 CHANNELS 5000 Hz

CALCULATE NOISE LEVELS

'-1""

AND TIMES OF OCCURRENCE (PARABOLICFIT)

1 GRAPHIC OUTPUTS

IDEN'[IFY CLUSTERS OF TRIPLETS

-,,(-

FORMULATE TRIPLET LIST

Fig. 1. Flow chart of spike separation procedure. (1) One-minute segments of data on two channels are digitized at 5000 Hz. (2) The noise level for each channel is then calculated and used as a threshold to locate individual spikes, whose (3) peak-to-peak amplitude and time of occurrence are measured with parabolic fit. (4) For each spike on the reference channel, the opposite channel is scanned for occurrences of spikes within a selectable time window. (5) Each pair of spikes (1 on each channel) whose offset tails within the time window is converted into an ordered triplet (x, y, z), where x -~ peak-to-peak amplitude on reference channel, y ~ peak-to-peak amplitude on opposite channel, and z = offset between the negative peaks. The list of ordered triplets is then scanned for well-isolated, compare clusters which correspond to different neuronal units. (6) Various 3-dimensional and superimposed-waveform plots are then used to inspect the results. The firing pattern and interspike interval histogram is plotted separately for each clustered unit. and conduction velocities cluster together and are identified. The results o f the automatic clustering p r o g r a m are displayed graphically to allow visual inspection a n d / o r interactive modification. A flow chart outlines the procedure (Fig. 1).

(1) Digitization Analog-to-digital conversion on two recording channels is done at 200 #sec intervals, allowing at least 25 sample points for even the fastest spikes, and providing an accurate representation o f the waveforms. Subsequent analysis is done digitally with F o r t r a n IV and assembly language programs.

(2) Noise calculation F o r each channel, the block o f 256 points (51.2 msec) with the smallest standard deviation (S.D.) is located a m o n g the first 30 blocks o f data. Based on the assumption that this block contains the fewest signals, the S.D. is used to estimate baseline noise. In the slow runouts (Fig. 3), the portion o f the record that is blanked out represents ~ 6 S.D., and only those spikes whose negative-going peak falls outside o f this range are included in the subsequent analysis. The negative peak is selected because there is only one per spike whereas there are usually two positive peaks 17.

(3) Spike measurements Independent lists o f spike measurements are created for each channel. Utilizing thresholds calculated in step 2, negative deflections are located subject to the constraint that if there are two or m o r e peaks within a 10 msec time interval, only the largest peak is selected. Two measures are obtained for each spike on each channel: peak-to-peak amplitude and temporal location. A 7-point least squares parabolic fit (Fig. 2A, inset) is used to estimate spike amplitude. The time of occurrence ot a spike is taken as the location o f the axis o f symmetry o f the p a r a b o l a fitting the negative peak.

459 A

t .........""-

1

""",..., ................................

1. PROXIMAL AMPLITUDE \ /

2 DISTAL " J "AMPLITUDE ~

PROX

1

1

Tt I

DIST 1

1

I i

2 ~CLUSTER

I 56 7

i!,

1

8 -- TRIPLET

~

CLUSTER

""-OF

Fig. 2. Schematic illustration of spike measurements and triplet formation. A: three measures are obtained for each putative unit. Peak-to-peak amplitude (in/~V) is measured from the first positive peak to negative peak for each spike on each channel. A least square parabolic fit (see inset) is used to interpolate between sampled points for more accurate readings. Offset (in msec) is measured between the axes of symmetry of the parabolas fitting the negative peaks. Inset: diagram of 7-point least squares fit of parabola superimposed on a negative peak. The axis of symmetry of the parabola is used to calculate the precise amplitude and time of occurrence of each positive and negative peak. B: The 3 measures obtained for each putative unit are used to form a list of triplets. Since there is not a perfect correlation between spike sizes on different channels, and since units can be either efferent or afferent, there is no way to tell at first which spike on the proximal is to be paired with which distal spike. Consequently, each proximal spike is paired (dotted lines) with every distal spike that occurs within a'reasonable' time window. In the diagram, 2 units (clusters 1 and 2) each fire 3 times, however 8 combinations (triplets) are selected. Two of the combinations (2 and 6) are incorrect, however repeated (3 or more) occurrences of the same unit serve to identify the clusters and eliminate the incorrect combinations.

(4) Triplet formation Using the measures just obtained, a list o f 'triplets' is created in which temporal offsets between spikes on opposite channels are measured. A triplet consists o f 3 pieces of information (Fig. 2A): (1) peak-to-peak amplitude on proximal (#V), (2) peak-topeak amplitude on distal (/iV) and (3) offset (msec) between the axis o f s y m m e t r y o f the parabolas fitting the negative peaks on proximal and distal channels. The channel with better signal/noise ratio is selected as the reference channel. F o r each spike in the reference channel (R spike), the opposite channel is searched over an empiricallydetermined time window b o t h before and after the R spike. F o r example, if the two electrodes are separated by 1 cm and the range o f conduction velocities is f r o m 0.5 to 1.0 m/see, then an afferent spike will occur on the distal electrode f r o m 10--20 msec before it appears on the proximal, whereas an efferent spike will appear on the distal electrode f r o m 10-20 msec after. A n y spike on the opposite channel that occurs f r o m 10-20 msec before or after the R spike is combined with the R spike as a 'triplet'. Therefore, a given spike f r o m either channel can appear in several triplets (Fig. 2B). It is also possible for a spike on the opposite channel to be excluded f r o m all triplets due to having no R spike within the time window. The present algorithm does n o t address itself to this situation, however if a R spike has no corresponding spike on the opposite channel within the specified time window, a d u m m y triplet is constructed whereby the peak-to-peak amplitude on the opposite channel and the offset are b o t h given a value o f zero. Thus, only one real measurement is taken, and the same procedure used to

460 distinguish spikes with 3 measurements (see below) can be used. However, since one measurement is usually insufficient to reliably identify a unit, only the largest spikes are considered. For ease of processing, two triplet lists are maintained and crossreferenced; one chronologically ordered, the other ordered by peak-to-peak amplitude on the reference channel. (5) Cluster identOqcation

The clustering of triplets associated with individual neuronal units has two stages: searching for a 'nucleus', and finding the rest of the cluster. We have established by intrasomatic stimulation of two giant neurons, RI and R2, that the coefficient of variation (s.d./mean) is less than 3 ~o for peak-to-peak extracellular amplitudes, and less than 1 ~ for temporal offsets across 50 spikes (Table I). Based on these findings, and using the amplitude-ordered triplet list to reduce search time, a nucleus of 3 triplets is sglected if (a) the coefficient of variation (c.v.) for each of the 3 variables is less than or equal to those found for R1 and R2, and (b) the 'distance' value, calculated as the sum of the c.v,'s for each variable (triple weight on offset) is minimal. The algorithm for scanning the ordered triplet list is a one-pass procedure which creates a list of nuclei with their associated distance values. For each triplet (except the last two) T(i), T(j) and T(k) and inspected such that i < j < k and k < i + 5 0 . If T(j) differs from T(i) by 5 ~ or more on P-P amplitude on the best channel, or if the offsets of T(i) and T(j) are of different sign, then no further permutations involving T(j) are considered, at which time j is incremented and k is reset to j + 1. For every i, the permutation T(i), T(j), T(k) with the smallest distance value becomes an element of the list of nuclei. From then on, the search for a nucleus only involves scanning the list for the minimum distance: If test (a) above fails, the selected list element is eliminated, and a new minimum distance is selected. To avoid reclustering problems, the list of nuclei is continually updated by removing elements containing clustered spikes. After a nucleus is selected, a sublist of triplets in the neighborhood of the nucleus is searched for the best candidate, i.e., the triplet whose sum of percent difference from the mean for each variable (triple weight on offset) is minimal. For example, i f a triplet TABLE I Coefficients of variation for two identified giant neurons in 3 preparations

In isolated abdominal ganglia, cells R1 and R2 were impaled and fired approximately 50 times. Extracellular spikes were recorded from R1 and R2 axons inright connective by two bipolar hook electrodes. Proximal and distal electrodes were separated by 1 cm (See Pinsker et al. 1976 for details). Cell R2

Preparation Proximal amplitude Offset Distal amplitude

No. 1 0.75* 0.98 1.10

Cell RI

No. 2 0.56 0.84 0.70

No. 3 )~(SD) No. 1 0.8 0.70 (0.13) 2.80 1.0 0.94 (0.09) 0.79 1.6 1.13 (0.45) 0.94

No. 2 2.2 0.6 1.9

No. 3 2.40 0.74 2.10

~(SD) 2.47 (0.30) 0.71 (0.10) 1.65 (0.62)

* Coef~cierit~'Ofvariation are calculated as standard deviation divided by mean, and are expressed as percentages.

461 has values (x,y,z) and the current mean is (R,S',z) then the sum l(x-~)/~l + I(y6')/?l + 31 (z[)/[1 is calculated. The candidate must then pass two closeness-of-fit tests: (1) individual variation, and (2) group variation. The allowable limits of each test were empirically determined, and can be tailored for use in different situations. These include a fixed parameter for offset and a variable parameter (based on signal-to-noise ratio) for the amplitudes. The reason is that offset will not be greatly affected by baseline noise whereas smaller amplitudes will. The first test is passed if the percent difference from the mean does not exceed 2~o for offset and [75/(9/n)] ~o for amplitudes, where 9 is the average amplitude on a single channel for the current cluster, and n is twice the noise threshold calculated during the 2nd step of the analysis procedure. For the second test, the triplet is added to the current cluster and the c.v. of each variable is recalculated. At this point the means of the variables may be subject to some drifting relative to the original nucleus. Thus, if the nucleus is off-centered, the algorithm makes continual adjustments to identify the entire cluster. If the c.v.'s do not exceed 1.2 ~o for offset and [50/(~/n)] ~o for amplitudes, the test is passed. The more stringent criteria of the second test are employed because there is always less group variation than individual variation. If both tests are passed, the candidate remains in the cluster and a new candidate is searched for. Otherwise, any changes made to the cluster are restored, and the candidate is temporarily flagged so that a search for a new candidate will not reinspect the old one. As soon as another candidate is accepted, however, all temporary flags are removed to allow reprocessing in the dynamic environment. Furthermore, since multiple triplets with the same spike are localized in the chronically-ordered triplet list, the latter is used to find and flag any triplets containing clustered spikes. Thus, as more clusters are identified less time is required to find new ones. This process is continued until no more candidates can be found which satisfy the above conditions, at which point the cluster, now complete, is recorded and a new nucleus is searched for. When no new nucleus can be found, all clusters are accounted for. The flagging procedure insures that the correct number of spikes are accounted for, regardless of the number of triplets to which a given spike belongs. This is done by maintaining a list of three elements for each triplet: (1) temporal location of proximal spike, (2) temporal location of distal spike, and (3) flag bit. As a triplet is given to a cluster, its flag bit is set, disabling any further processing involving that triplet. Then, all other triplets are flagged which contain a spike with the same proximal or distal temporal location as the clustered triplet. Thus, several triplets may be flagged (accounted for) while only one triplet is clustered. This technique eliminates the possibility that a spike might become associated with more than one cluster.

(6) Graphic output Once the automatic clustering is finished, several quality-control graphical outputs and tables can be obtained. These include (a) slow runouts of raw data, (b) 3dimensional plots, (c) fast runouts of raw data with labeled spikes, (d) superimposed individual and averaged waveforms of clustered units, (e) firing pattern for each unit and (f) summary tables. (a) Slow runouts are used to show gross multiunit activity and

462 size ranges of neurons, as well as revealing 'glitches' which can be eliminated interactively. Also, the slow runouts indicate the threshold beyond which further analysis will be carried out by blanking-out the appropriate central (i.e., noise) portion of the data from each channel. (b) The most flexible outputs are 3-D plots in which each of the 3 measurements for a triplet is expressed as one dimension. Clusters and/or triplets can be viewed from any angle on either side of the plane of the graph paper. The plots use vertical straight line segments to represent a point in 3 dimensions. The bases of the segments are affixed to a 2-dimensional scaled grid, and the lengths of the segments are proportional to the value of the third dimension. With angular perspective variation, well-isolated clusters can be distinguished from neighboring ones, and unclustered triplets within graphical boundaries can be plotted. Each clustered spike is plotted with an identification number. If unclustered triplets are included, they can be labeled with index numbers which can be used to locate the raw data. This allows inspection of unclustered triplets in the vicinity of clusters to determine the nature of the triplets, i.e., whether the spikes are compound or are different units. (c) Fast runouts with labeled spikes display raw data with labels on every spike involved in the analysis. All clustered spikes are indicated with the cluster number in parentheses, and the triplet number(s) are indicated for all spikes above threshold, thus identifying the distal spike(s) that were paired with each proximal spike. (d) In order to inspect homogeneity within a cluster, individual waveforms can be superimposed on top of one another. Both channels are DC adjusted using a parabolic fit on the negative peaks, and the proximal channel is time-triggered on the parabolic axis. The distal channel is plotted with time delays adjusted to reflect the actual offset as originally measured for each individual spike. Occasionally, with some of the smaller units, we can detect mistakenly clustered spikes which happened to have similar measurements as the rest of the cluster, but which differed in waveform. Differences in waveform and conduction velocity across units can be seen by plotting the averaged wavet0rms for the clusters. On the proximal channel, the waveforms displayed are at successive DC shifts and triggered on the parabolic peaks. (e) Two types of graphical output display firing pattern information for most of the clustered units. One type employs solid vertical bars placed along a horizontal axis at every point corresponding to the time-of-firing of the unit. Another indicates the interspike interval histogram for the individual units. (f) Some tables enumerate the clusters and lists size as well as the mean, standard deviation, and coefficient of variation for each variable. Others provide accounting information, including number ot peaks/channel, numbers of spikes with 3 measures, and number of clustered and unclustered spikes. RESULTS The 1-min segment of in vivo background spontaneous activity in the siphon nerve of a quiescent animal is shown in a slow runout of the digitized spike train data on two channels (Fig. 3). The procedures for displaying the slow runout are different from those used to separate out the individual spikes, and are intended solely to provide a visual record of the sample of spike train data to be analyzed. The height of

463 CC144 SPONT BCT,ITEH 2 R l , t - 2 7 - 7 7 POST-OP ORY # 1 . T E H P ( C ) = I S , B BLOCKS 1-12BB

,l,r,,

Jt ,J,

," i

,.I, ,i,, h , , , LIJ,I..I1,1, ,, tr,,,jr

15

26

25

'*"I 3'~----~-

,1

iI

1~'111111'r ,'o

58

--+

55

6~

TIHECSEC.)

Fig. 3. Slow runout of background spontaneous activity in siphon nerve. The blank portion of each channel represents 4- 6 S.D. of the noise level, calculated during the quietest 51.2 msec (256 sampled points) interval of the first 1.5 sec of data. The rest of the plot is drawn with line segments depicting the time, height and polarity of every inflection beyond the noise threshold. This method of plotting is unrelated to any analysis done on the data, and is used merely to obtain rapidly an overall picture of the l-rain segment. The information in the upper left corner identifies the experiment in which this segment was recorded (Kanz, Eberly and Pinsker, unpublished data).

each positive and negative inflection above the noise threshold is measured, and a straight line of that height and polarity is constructed at the appropriate point in the record. Note that there are many units with different spike amplitudes that contribute to the background spontaneous activity, and that a given unit can have a large spike on the proximal and a small spike on the distal electrode. The objective of the program is to account for all of the signals above the calculated noise level by separating them into individual clusters, each corresponding to activity of a different axon in the nerve. The clustering program makes quantitative comparisons among the lists of triplets (see Methods). The separation procedure is based on 2 amplitudes and 1 offset. Therefore, each spike is considered to be a point in a 3-dimensional space and displayed as a line segment extending from a plane. If the variability within units is less than that across units, then spikes belonging to different units will cluster in different locations in the 3-D plots. A total of 203 spikes were automatically separated into )2 clusters for this 1-min segment of background spontaneous activity (Fig. 4). In the 3-D plots, temporal offsets are converted into conduction velocities. All of the clustered

464 spikes in the present s a m p l e are efferent; however, the p r o g r a m can also cluster afferent spikes. Clusters are n u m b e r e d in o r d e r o f spike a m p l i t u d e on the channel with better signal-to-noise r a t i o ( p r o x i m a l in this p r e p a r a t i o n ) . A n u m b e r o f the largest CC144 SPONT

AET

2Al

ZI

CE .r-)

12

Fig. 4.3-D plots of clustered triplets. The 12 clusters (total spikes : 203) were automatically identified, numbered in order of size on the proximal channel, and plotted by computer. Each cluster represents repeated (3 or more) firings of the same neuronal unit, giving rise to approximately the same values for the three measurements: peak-to-peak amplitude on proximal channel, peak4o-peak amplitude on distal channel, and conduction velocity. Thus, the location of the base of each line segment indicates the proximal'amplitude of the spike and its speed of propagation, whereas the length of the line segment indicates the distal amplitude. The arrows with numbers indicate 9 well-isolated clusters in which we have the most confidence.

465 spontaneously active units in the nerve have been identified intracellularly as siphon m o t o n e u r o n s (D. Allen,.unpublished data). N o t e h o w well units that are large on at least one channel (e.g., clusters 1-4) separate in c o m p a r i s o n to signals that are small on both channels. Clearly then, the tightness o f the cluster and its degree o f isolation f r o m neighboring clusters depends primarily on signal-to-noise ratio. Although 12 putative units were clustered by the program, only 9 o f t h e m (Fig. 4, arrows) have sufficient spatial separation to convince us that they represent unique units.

1.

CLUSTER#1, 6 SPIKES gi MEANCOND. VEL.-e.607 MI5

~[pnox

o 9~.2 ±2.28uv

OFFSET- 16.5

± 0.e9 MS

~I~!OlST = 19.3 ± 1.35 UV

(2.5%) (O,S~) (7.0~)

u

i:i

|

2.

CLUSTER #2, 19 SPIKES gT MEAN coNo. V E L . - 9 . 6 0 6 MIS I

PROX : ~ SS. 4 FFSET= 16.5 , O I S T , - 9B.6

,a

s

~

i 1, 42 UV ± B.SS MS ± 1.4B UV

~s TZM~MS)

~s

(1.7%)

(8.2~) (I.S%)

31

~s

Fig. 5. Well-isolated clusters. A: the 3-D blowup of the two clusters with largest proximal amplitudes demonstrates the lack of variability for conduction velocity in larger units, as well as total 3-D separation. Note that only one variable, distal amplitude, is responsible for the separation; the other two measurements alone could not be used to distinguish the two units. B: real waveforms of the triplets in clusters 1 and 2 are displayed. Both channels are DC adjusted, and the proximal channel is timetriggered on the parabolic axis of each negative peak. The time lag on the distal channel corresponds precisely with the original raw data, and is the most important measurement used to cluster the spikes. In extraeellular axonal spikes, the early positivity is usually larger and faster than the late positivity, however the filters introduced some distortion of the waveforms. Numbers in parentheses refer to the coefficients of variation for the particular variable. For well-isolated clusters these values are usually less than our empirically determined values (Table I). The 'double spike' for cluster 1 on the distal record may reflect an axon that bifurcated after the proximal electrode.

4a

466 H u m a n judgement is still required to evaluate and correct the results of the automatic clustering p r o g r a m ; however, this can be extremely time consuming. For the present 1-min sample, the automatic clustering took approximately 15 rain. However, obtaining the various runouts for visual inspection added several hours, and interactive corrections could add m a n y more hours to the data analysis. The a m o u n t o f time that is w o r t h investing in interactive processing depends u p o n the degree o f accuracy required by the experimental hypothesis. We shall examine more closely the results o f the automatic p r o g r a m and illustrate several ways in which h u m a n interaction can improve these results. Two potential types o f clustering errors must be evaluated: inclusion o f incorrect spikes in a cluster (Figs. 5-7) and exclusion of correct spikes (Figs. 8, 9).

B1.

PROX OFFSET-

l.~t

15+.~

±

17.7

± 0.15

UV M5

(?,97.) (~.8%)

ca

2.

CLUSTER #10,

l e SPIKES

t.

• •O;,C,~p ,.'~f

++o+,.-,"++

OFFSETOIST •

1B.2

12.6

i 0,15 MS ± 1.06 UV

. . ,~,. .,'s Ttt.,";.S),E( z'+~

(O,BI) (B,41)

,;'

+'s '

Fig. 6. Partially overlapping clusters. A: clusters 9 and l0 illustrate two clusters similar in amplitude on both channels, but with only minimal overlap on conduction velocity. Because of the close proximity between the clusters, one or two misidentifications of triplets might have occurred in the area of overlap (at approx. 0.52 m/sec velocity). B: the superimposed wavvforms of clusters 9 and 10are triggered on the proximal channel. Note the significant increase in coefficients of variation for the amplitudes on the proximal channel as compared with clusters 1 a~d 2 (Fig. 5]3). The offsets also show more variability (0.8 ~), but are still within empirically determined limits.

;

467

&

B1.

CLUSTER a8. Z5 SPIKES I~ HEAN COND, VEL.-11.465 H/S

L =

QFF$ET-

.

21.5

.

.

= 11.13 HS

(11.6I)

A

N 0

2,

CLUSTER I | | ,

37 SPIKES R/S

~, MEAN CONO. VEL.-B.476

PROX • | 3 , 5 OFFSET= E l . 8

t 1.7E UV i | , 1 6 HS

(|Z.71) (i.6|)

0[5T

t

'1 ,

-

|

11.9

tl

1.97

UV

(15.71)

Is II~115)

['1

31

~s

Fig. 7. Overlapping clusters. A: the 3-D plot of clusters 8 and 11 shows two units that are difficult for the computer algorithm to separate accurately because of dissimilarity on only one variable (amplitude on distal channel). Since the signals are small, more variability in amplitudes is allowed, resulting in several misidentifications of triplets marked with '8'. With our present procedures, this problem is difficult to avoid with small amplitude units whose clusters overlap substantially. B: superimposed waveforms of clusters 8 and 11. The coefficients of variations are high (over 12 ~o) for amplitudes, but still less than 1 ~ofor offset. Note that cluster 8 contains 6 spikes on the distal channel that appear to be rnisidentified and should have been grouped into cluster 11. T o determine whether discrepant waveforms have been clustered inadvertently, all o f the spikes in each cluster are superimposed and inspected visually. Since the cluster identification p r o g r a m considers only one measurement per waveform (excluding offset for the moment), it is possible to have several discrepant waveforms with the same peak-to-peak amplitude. This possibility is minimized by the fact that the 3 measures (2 amplitudes and 1 offset) are relatively independent, since there is not a perfect correlation between spike sizes on the two recording channels. Furthermore, there is a relatively low correlation between spike amplitude and conduction velocity due, in part, to infolding o f axonal membranes in molluscan axons ~0 and to the fact that individual axons are different distances f r o m the recording electrodes. Thus, if two units are similar in amplitude but different in shape on one channel, it is unlikely

468 that they will have the same amplitudes on the other channel as well as identical conduction velocities. Equivalently, any well-isolated, compact cluster will most likely contain homogeneous waveforms. We have found this to hold true extremely well, with exceptions occurring only with clusters involved in partial overlap with other clusters. A pair of well isolated clusters are plotted in 3-D (Fig. 5A) along with superimposed spikes (Fig. 5B). Cluster 1 fired 6 times and cluster 2 fired 19 times during the l-rain sample. Cluster 1 has the largest spike on the proximal channel, but a relatively small spike on the distal channel, whereas cluster 2 has large spikes on both channels. Although there is some overlap between the two clusters on 2 measures (proximal amplitude and conduction velocity), there is no overlap on the third (distal amplitude). This permits the program to distinguish these two units reliably. Despite the fact that the clustering program does not utilize any information about the individual waveforms other than peak-to-peak amplitude, the spikes of well isolated clusters are remarkably superimposable (Fig. 5B). Only the proximal negative peaks are CC144

SPONT

FILI

. ~ - - - - - - ~

.......

Jl~.;'.~ 1 27 / 7 PuST--OP [n4Y ~ 1 IEi'4P ( t : ) ~ 1 !,. l~ EUFF b E P , I ~ Mt,1

IK

1

v

V~ t~

h

8

12

ll6

20

26

~8

OFFSET(H5)

Fig. 8. Averaged waveforms. The 9 units with best separation (Fig. 4) are plotted, triggering on the proximal channel. Each trace is the average waveform for all spikes in the cluster (cluster number indicated below negative peaks), Distal spikes are displayed at the average offset for the cluster. This plot is similar to an oscilloscope tracing with certain refinements: (1) the entire proximal waveforms can be seen, (2) the averaged proximal waveforms are time-synchronized on their respective parabolic axis, (3) only the average waveform from each unit is displayed, allowing easier observation of differences between units.

469

v

I-t-~

|

|

||

0

|

| b

|

a., Co~,f . ,~_

(

Fig. 9. Unclustered triplets. The 3-D plot shows all triplets left over after the clusters were automatically separated. The largest triplet is a compound spike, belonging to duster 2 and another smaller unit, which was not identified by our algorithm. The remaining triplets might be different units that fired less than 3 times or compound spikes, but without analysis of more data to reveal new clusters, or individual inspection of the triplets to locate complex waveforms, this question remains unanswered. There are 37 unclustered triplets, and therefore approximately 85 ~ of the total 3-valued triplets were accounted for in the 1-min segment of background spontaneous activity. synchronized in these plots; the location o f the distal peaks reflect the actual measured offsets. Therefore, the high degree of superimposability o f the distal spikes in these clusters also reflects the reliability o f the conduction velocity as a 'fingerprint' for a given axon in the nerve. The averages, s.d.'s and coefficients of variation (c.v.'s) for the three measures are displayed between the proximal and distal records. N o t e the extremely small c.v.'s for the offsets (_< 0.5 ~o). The large spikes (cluster 1 proximal, cluster 2 proximal and d i s t a l ) a l s o have very small c.v.'s for amplitudes ( < 2 . 5 ~ ) , whereas small spikes (cluster 1 distal) have proportionally larger c.v.'s. These statistics are s u m m a r i z e d for all clusters in Table II.

470 TABLE 1I Individual average cluster measurements

The 12 automatically separated clusters are ordered from the largest proximal peak-to-peak amplitude to the smallest. The number of spikes in each cluster is indicated in the column labeled 'N'. The man, standard deviation, and coefficient of variation for each of the three measurements is given. Note the small coefficients of variation for the offsets compared to those for the amplitudes. Proximal amplitude

Distal amplitude

Cluster No.

t~ V

SD

°o

tt V

SD

%

MS

SD

%

N

1 2 3 4 5 6 7 8 9 10 11 12

90.2 85.4 35.6 28.0 24.2 17.5 16.5 16.1 15.3 13.9 13.5 8.7

2.28 1.42 1.71 1.35 0.38 1.36 1.56 1.99 1.21 1.43 1.72 1.96

2.5 1.7 4.8 4.8 1.6 7.7 9.5 12.3 7.9 10.3 12.7 22.6

19.3 90.6 35.2 20.3 18.8 14.9 15.5 14.7 14.3 12.6 11.9 8.9

1.35 1.40 1.67 1.47 0.44 2.06 1.77 2.82 1.11 1.06 1.87 1.17

7.0 1.5 4.7 7.2 2.4 13.8 11.4 19.1 7.7 8.4 15.7 13.1

16.5 16.5 20.2 21.3 22.5 22.1 19.9 21.5 17.7 18.2 21.0 22.7

0.09 0.03 0.05 0.07 0.07 0.17 0.11 0.13 0.15 0.15 0.16 0.20

0.5 0.2 0.2 0.3 0.3 0.8 0.6 0.6 0.8 0.8 0.8 0.9

6 19 19 10 3 21 22 25 19 10 37 12

Offset

As the amplitudes o f the spikes decrease, the clusters typically show more overlap a m o n g measures and poorer visual isolation in the 3-D plots, and clustering errors o f b o t h types are more likely to occur. Nonetheless, some smaller units still show relatively g o o d isolation f r o m neighboring ones (e.g., Fig. 4, clusters 7 and 12). In some cases the overlap is not too extensive on at least one variable, and the probability o f including incorrect spikes can be evaluated from a close inspection o f the 3-D plots. F o r example, clusters 9 and 10 are completely overlapping on both proximal and distal amplitudes, but are only minimally overlapping on conduction velocity (Fig. 6). In this case, it seems reasonable that n o m o r e than 2 or 3 spikes (out o f a total o f 29) have been mistakenly clustered. The waveforms associated with the specific triplets that show overlap in the 3-D plots can be displayed separately, and c o m p a r e d with the remaining waveforms o f the 2 clusters. W h e n neighboring clusters overlap m o r e extensively on all 3 variables (Fig. 7), then confidence necessarily decreases. F o r example, clusters 8 and 11 show p o o r visual isolation in the 3-D plots (Fig. 7A) and substantial variability in waveforms (Fig. 7B). Closer inspection suggests that at least 6 spikes f r o m cluster 11 were mistakenly placed into cluster 8. Since cluster 6 is also involved in overlap with these two clusters (see Fig. 4), and since these signals all have relatively p o o r signal-to-noise ratios, our technique provides little basis for reliable discrimination a m o n g these 3 putative units. Thus, o f the 12 clusters that were separated automatically, we have confidence in 9, after considering compactness, isolation from other clusters, and homogeneity within a cluster. The average waveforms for the 9 clusters (Fig. 8) suggest that each has a unique waveform and conduction velocity, and belongs to a different axon. With interactive procedures, any waveforms that are suspected to be incorrect can be eliminated from a cluster, and other waveforms can be added. The effect o f

471 these procedures can then be evaluated on (1) the superimposed waveforms, (2) the visual isolation in the 3-D plots, and (3) the values for the 3 measures. The important point is that, with some investment of time, it is often possible to 'sanitize' clusters in which we previously had little confidence and thus obtain pure unitary waveforms of individual units. This information can be used to form 'templates' of these units that may eventually allow them to be recognized in compound action potentials when other units are firing simultaneously. Inspection of the waveforms of the clustered triplets can provide information concerning the first type of clustering error (i.e., inclusion of incorrect spikes). However, in order to evaluate the second type of clustering error (i.e., exclusion of correct spikes), a 3-D plot of all of the unclustered triplets (Fig. 9) is helpful. For the lrain sample in which 203 triplets were clustered (Fig. 4), 37 were not (Fig. 9), i.e., about 15 ~o of the total number of triplets remained unclustered. These triplets might represent unitary spikes that fired less than 3 times, or compound spikes. Note that there are relatively few unclustered spikes in the regions of well isolated clusters and more in the regions of poor isolation (compare with Fig. 4). In some cases, it may be necessary to guarantee that no triplets belonging to a particular cluster have been missed, i.e., to obtain an accurate picture of a given cell's firing pattern. In such cases, it is possible to display a 3-D blowup of one or more individual clusters, with all of the unclustered triplets in the general region (Fig. 10A). B

A

1.

52

2.

"~33

3.

118

33

°.4;,

~g

B.~4

B,eS

TIME(SEC)

, e.12

Fig. 10. Cluster with unclustered triplets. A: blowup of cluster 3 with all unclustered triplets in range of graph. This plot is used to examine how well the clusters were identified by the computer by inspecting any neighboring triplets that might have been excluded inadvertently from the cluster. Line segments with numbers other than 3 are unclustered triplets with an index number used for fast data retrival. Given.the typically small variability of conduction velocity for a unit, the only triplets that require inspection for similarity to cluster 3 are triplet numbers 33, 52 and 68. B: the waveforms of triplet numbers 33, 52 and 68 indicate superposition with another (smaller) unit. (1) On both channels, a rising phase from each unit adds to increase peak-to-peak amplitude. (2) On distal channel, a rising phase from each unit increases peak-to-peak amplitude. (3) On both channels, a falling phase of the smaller unit cuts into the initial rising phase of the large unit to decrease peak-to-peak amplitude. Thus, triplets 33, 52 and 68 are probably unit 3 spikes that the automated algorithm could not handle.

472 Specific unclustered triplets can be displayed in 100 msec time samples of raw data containing the spikes involved (Fig. 10B). To account for every such triplet would be time consuming. However, for all but the smallest units, the number of triplets that must be inspected is reduced by the fact that while superposition of units affects amplitude, the measured offset if relatively unaltered. F o r example, 3 unclustered triplets have the same conduction velocity as cluster 3 but are different in one or both amplitudes (Fig. 10A). Closer inspection (Fig. 10B) of the 3 spikes reveal superposition with another unit, creating complex waveforms in each case. Since our basic spike separation technique is automated, and designed to handle only unitary spike trains, we do not want to emphasize this interactive type of quality control, but make it clear that we can quantify and document our mistakes. To evaluate the computer results, we can also plot longer samples of raw data with triplet and cluster numbers indicated for each spike (Fig. 11). Plots of this kind allow us to see not only which spikes were clustered (number in parentheses), but also which were totally left out of the analysis (unnumbered spikes). In the sample shown, all 5 triplets having a spike on both channels were clustered, whereas the two proximal triplets (41 and 42) containing no corresponding distal spike were not. The accounting information for the l-min sample is shown in Table IlL Note that these values summarize the results of the automatic clustering program, and would be substantially improved if we had utilized our interactive procedures to correct clustering errors. Once the results of the clustering program have been evaluated, and the errors

u~ c~

TIME(MS)

Fig. 11. Fast runout with labeled spikes. Plotted are 1.5 sec of digitized data with labels attached to all spikes involved in the analysis. The numbers in parentheses are the cluster numbers, and those without parentheses are triplet numbers. Parentheses without a number indicate an unctustered spike. This plot shows (1) how many spikes from each channel were detected, (2) how many spikes were clustered, (3) which spike on proximal channel corresponds with which distal spike, and (4) which peaks in the record were left out of the analysis. In this particular case, 71% of the spikes on the proximal channel and 100 % on the distal channel were automatically clustered. The two unclustered proximal spikes have no distal spike with which to calculate an offset. Thus, only one measurement per spike was used to try to cluster them, ending in failure.

473 TABLE III Accounting information from spike train analysis

'Total spikes' refers to the absolute number of spikes detected beyond the noise level. 'Spikes with 3 measures' includes every spike belonging to a 3-valued triplet. Recall that ifa given spike on one channel has no corresponding record on the opposite channel, the spike is still placed in a triplet, however only one measure (amplitude) is made. 'Clustered spikes' normally include spikes which were part of 3-valued triplets, as in this case. Occasionally, however, some of the largest units have records on only one channel, resulting in clusters based on only one w~veform parameter. By visual inspection, these clusters are accepted or rejected, and the spikes are tallied in the above columns. 'Unclustered spikes' includes spikes in both 3-valued and 1-valued triplets, and therefore, this number (46) is greater than the number of unclustered triplets (37) displayed in Fig. 9. A good measure of per cent accounted for is made by dividing the lesser of 'clustered spikes' by the lesser of 'spikes with 3 measures'. In this case, 81.5 ~ of the total were accounted for. Proximal channel

Total Spikes 373 Spikes With 3 Measures 260 Clustered Spikes 203 Unclustered Spikes 46

.

Distal channel

278 249 203 46

corrected (if so desired) by means of interactive procedures, then the firing patterns of individual 'dusters can be displayed. For the 9 clusters in which we have good confidence, we have indicated the times of occurrence of the individual spikes for each putative unit across the 1-min sample (Fig. 12A), and also the corresponding interspike interval histograms (Fig. 12B). The 9 units shown have a variety of rates (322 spikes/minute) and patterns ~of spontaneous activity. The individual firing patterns (Fig. 12A) of the larger units can be compared with the slow runout of the multiunit spike train data (Fig. 3). The unit corresponding to cluster 1 only began to fire after about 30 secs. This unit had the largest proximal and a small distal spike (Fig. 5), and can easily be seen only on the proximal channel in the slow runout (Fig. 3). The firing pattern of cluster 2 can readily be identified as that o f the largest distal spike in the slow runout. Note that a window discriminator on the proximal channel would most likely have confused clusters 1 and 2. Some units(e.g., cluster 5) fire only rarely, whereas others (e.g., clusters 2 and 3) fire at much higher rates and with long stretches of highly periodic activity. Once the time of occurrence of individual spikes has been determined, it is possible to utilize higher order techniques (e.g., autocorrelation, cross-correlation) to characterize rhythmicity of a single unit and/or excitatory or inhibitory interactions between units2,6,13,16,zL The present procedure works well when most of the spikes in the record are unitary, as in the case of background spontaneous activity in the siphon nerve. We have analyzed successive 1-min samples of background spontaneous activity from the same nerve separated by 30-min intervals. Despite gradual changes in amplitude and conduction velocity, we were able to follow the same population of units. This indicates that the technique is sufficiently reliable to be used for longer-term neuronal analyses.

474

A

INOIVIDURL

B

FIRING PRTTERNS

INTERSPIKE INTERVRL HISTOGRRMS

CC144 SPONT RCTo 2R1

I

CC144 SPONT ACT, 2RI BIN WIDTH - ,5 SEC

TEMPORRL LOCRTION OF SPIKE

"

CLUSTER

CLUSTER

,

#t

:

.I

#2

I I1,111]111.111

#3

III

IIIIIII

I

#4

I

#7

#9

:1

#10

,I

#12

I.I

I, .

I J Ill

.I

II I.IIII

10

I

I ,I

I

I,I

11.1

I

IUI

Ill

26

11

I

I II

I

I.

I

I

I

I

TIME(SEC)

I

m,,~,

#2

i

~

#3

I

,

~,4

:

,I

,

#5

I

~

i ,I

~?

•

! 111

.I

I

I_~L

I I

3'6

#1,

_,11

I:.il

,

H

I il

: II

.I

I

~S

I .I

I a6

I

•

iL~.,

m,,

.m__

l,

m

-

i

ml

m

.m m

j

,I

,__ , 5'6

,,

~

1

,

#16

:

~

:1 60

,lz

:1 6

.J.

.

.

. ~

.

: B

TIME

.

I

;

I

12

16

210

INTERVflL(SEC)

Fig. 12. Individual firing patterns and interspike interval histograms. Only the 9 weU-isolateddusters in which we have good confidence (Fig. 4, arrows) are plotted. A: each axis corresponds to the cluster whose number is at the left, and the bars on the axis represent the time of firing of the particular unit. For some of the larger units (e.g., 1 and 2) the patterns can be directly compared with the slow runout of the raw data (Fig. 3). The regularity of unit 2 is in close agreement with the largest distal spike in the raw data, and the pattern for unit 1, whose proximal spike is slightly larger than that for unit 2, can also be visually compared with the slow runout. B: the interspike interval histograms of each of the units are plotted, based on the individual firing patterns. These are used to indicate regularly firing units (2, 3, 7, 9) and to display modal interspike intervals. This technique can also be used, in a limited way, under conditions of heightened activity when many compound spikes are present. A characteristic burst of efferent activity occasionally occurs in the siphon nerve that involves many of the same units that contribute to the background spontaneous activity, but can bring in other units as well. This burst is triggered by activity in the interneuron II network 2a. In vivo records during spontaneous interneuron I1 bursts show to phases: increased activity in smaller units during which time larger units are inhibited, followed by postinhibitory rebound firing of larger units (Kanz, Cobbs and Pinsker2S). A 1-min segment of data from the same preparation, but including a spontaneous interneuron II burst (slow runout, Fig. 13A), was automatically clustered (Fig. 13B). There were a substantial number of unclustered triplets (due to compound spikes) during the height of the spontaneous interneuron II burst. Nonetheless, 4 new well-isolated clusters (Fig. 13B, large arrows) appeared in the 3-D plots in addition to m a n y of the same clusters (large numbers) that contributed to the background spontaneous activity in the previous sample (compare with Fig. 4).

2~

475

B

A C¢lk& SP0NT INT II,ITEM ~1,1-E7-77 P0$1-0p 01~¥ I I , r E N p ( C ) . I S . B

BLOCI~S 1-12m|

Fig. 13. Spontaneous interneuron II burst in siphon nerve. A: slow runout of 1-min segment of increased activity from the same animal during a spontaneous firing (at approx. 30 sec) of interneuron II. The data were automatically analyzed by the program to discover new units which had not fired previously during background spontaneous activity. Note that the gain on the proximal channel is lower than that in Fig. 3, due to the occurrence of large proximal spikes that fired shortly after the high frequency burst of activity in the smaller units. B: 3-D plot of clusters reveal new units (large arrows) in addition to previously clustered units (large numbers). The small numbers on each line segment are unrelated to those of Fig. 4. However, by comparing spatial relations between clusters on the two graphs, many of the same units are easily identified. For ease of comparison, the large numbers are identical to the large numbers in Fig. 4. There were many unclustered triplets (due to superimposed spikes) during the 10-20 sec of heightened activity; nevertheless, enough unitary waveforms were present to allow us to generate templates for the new units. The four well-isolated new clusters indicate 4 units that did not fire during the previous segment of background spontaneous activity but were involved in the spontaneous interneuron II burst.

DISCUSSION

W e have described a n a u t o m a t e d off-line technique for s e p a r a t i o n o f u n i t a r y spikes f r o m m u l t i u n i t spike t r a i n d a t a . A s is true for o t h e r techniquesii,i2,i4, 23, the present a p p r o a c h involves clustering o f units on different w a v e f o r m measures. U n l i k e m o s t o t h e r techniques, which e m p l o y m u l t i p l e d e p e n d e n t m e a s u r e m e n t s o f a single w a v e f o r m r e c o r d e d on o n l y one channel, we utilize two r e c o r d i n g channels, allowing us to t a k e o n l y one m e a s u r e ( p e a k - t o - p e a k a m p l i t u d e ) for each o f two waveforms, a n d also to calculate the c o n d u c t i o n velocity. This provides us with t w o distinct a d v a n tages: (1) spike a m p l i t u d e s are u s u a l l y different on the t w o channels - - a unit which has a large p e a k on one channel can have a small one o n the other, a n d (2) c o n d u c t i o n velocity is m i n i m a l l y affected b y noise a n d is a stable p r o p e r t y (at c o n s t a n t t e m p e r a ture) o f a given unit. W i t h the a d d e d refinement o f the p a r a b o l i c fit, the coefficients o f

476 variation for offset are so small ( < 1%) that units of similar amplitude can sometimes be distinguished reliably just on the basis of offset. The present technique for spike separation is both practical and efficient, in that it requires only two recording channels and 3 measures to separate out reliably many different units, with no a priori knowledge of the number of units present or their waveforms. There is no theoretical limit to the number of units that can be discriminated, other than the ultimate limitations determined by signal detection theory7. Two types of mistakes that any spike separation routine must guard against are (1) to include spikes that do not belong, and (2) to fail to include spikes that do belong. We guard against the first possibility by examining the superimposed waveforms of the clustered units, and the second by examining unclustered triplets in the region of the clusters. If necessary, the appropriate spike(s) can then be added to or deleted from the cluster by an interactive procedure. Although we routinely inspect our data in these ways, 'mopping up' procedures are extremely time consuming and are not necessary for some experimental goals. Automated procedures would save time and improve the quality of the data first presented to the human observer. At present, we prefer to analyze more samples from the same experiment rather than invest more time in a single sample. A major limitation of the present technique is that it can only be utilized on unitary waveforms, since it depends upon stable spike amplitudes. To handle the problem of spike superposition we may be able to merge our technique with that of Roberts and Hartline 21, who have developed an automated procedure for separating out known unitary spikes from a background of simultaneous activity in other known units (see Introduction). Since our procedure automatically provides templates for many units (e.g., Fig. 8), these might be used to generate the required digital filters. In order to separate a large number of units from only a few channels, Roberts and Hartline 21 use an iterative procedure, which involves separating out the largest units first and then subtracting the templates from the record. An alternative procedure to analyze compound spikes will also be explored. As long as there are sufficient unitary waveforms in the records to allow clustering of different units, the offsets can be used to identify the units contributing to a compound spike. If the negative peak is large enough to be detected when it is firing in isolation, it should produce a detectable inflection in the record even when other units are firing at almost the same time. Inflections on the proximal and distal records can be used to measure temporal offset, and the library of clustered units can then be searched for a unit with that particular offset. If an appropriate unit is found, the average waveforms on proximal and distal records can be lined up and subtracted from the complex spike train to reveal the underlying spikes. A less serious limitation of the present approach is that it is conducted off-line, i.e., the spike train data are recorded on FM tape for subsequent processing. We have not explored the possibility of on-line analysis because most experimental questions do not require this refinement. Furthermore, when a permanent record of the original

477 data is stored, it is possible to analyze the same set of data repeatedly and thereby gain confidence in the results. To our knowledge, there is no other practical procedure than can separate as many individual units as the present approach. The success of our analysis results, in part, from our selection of appropriate electrophysiological parameters. (1) We record from a relatively homogeneous population of straight-through axons which transmit information in parallel down the nerve. (2) The unitary waveforms are homogeneous, and therefore, little information is lost by reducing them to a single number representing peak-to-peak amplitude. (3) We use two channels to obtain an additional measure (i.e., conduction velocity) that is relatively independent of noise and is an extremely reliable 'fingerprint' for an individual axon. Preliminary results indicate that we can follow the same units in our 1-min in vivo recordings over prolonged time intervals (i.e., several hours) despite the fact that both amplitudes and offsets might be changing gradually. This is done by comparing the geometrical relationships between clusters across 3-D plots. It is also possible to relate our electrophysiological measurements to the underlying anatomy of the nerve. We have now compared cross-sections of several different nerves with 3-D plots of the clustered units, and have seen clear correlations between electrophysiological (i.e., spike amplitudes and offsets) and anatomical (i.e., size and relative numbers of' axons, amount of infolding of surface membrane) measurements. The present technique for spike separation should be generally useful in increasing the yield of experiments which now use nerve fascicles or window discriminators to eliminate all but a single unit. Although it was developed to allow us to quantify in vivo whole nerve recording from Aplysia axons, this approach should be applicable (with minor program modifications) to any preparation in which it is possible to record from the same population of units at two different locations in the nervous system. When dealing with a large number of faster action potentials, as would typically be the case in a mammalian preparation, several problems may arise. First, depending on the computer system configuration, the need for a taster sampling rate may impose limitations on the amount of data which can be analyzed. Second, a very active nerve bundle might contain an abundance of superimposed waveforms, which the clustering technique could not handle effectively. Third, heigthened activity results in the formation of many more triplets which increases processing time, and also makes it difficult to separate the clusters from the remaining 'duplicate' triplets. Provided the record contains a majority of unitary spikes, however, this technique may help to answer questions regarding the correlation between activity in neuronal populations and behavior. ACKNOWLEDGEMENTS We thank Drs. W. Aspey and H. Bryant and Mr. W. Farr and R. Zears for helpful discussions; Dr. W. Willis, Jr. for support throughout this project; Dr. J. Kanz and Mr. L. Eberly for providing the spike train data; Drs. M. Correia, R. Glantz, E. Glaser, R. Grisell, and D. Hartline, and Mr. C. Taylor for critical comments on an

478 ear l i er version o f this m a n u s c r i p t ;

M r.

W. P o g u e and

Ms.

N.

Rothman

for

p h o t o g r a p h y , a n d Ms. E. Preslar f o r typing the manuscript. This research was s u p p o r t e d , in part, by R e s e a r c h G r a n t s N I H N S 12223 and B N S 76-17480 to H . M . P .

REFERENCES 1 Adrian, E. D. and Zotterman, Y., The impulses produced by sensory nerve-endings. Pt. 2. The response of a single end-organ, J. Physiol. (Lond.), 61 (1926) 151-171. 2 Bryant, H. L., Jr., Ruiz Marcos. A. and Segundo, J. P., Correlations of neuronal spike discharges produced by monosynaptic connections and by common inputs, J. Neurophysiol., 36 (1973) 205-225. 2a Byrne, J. H. and Koester, J., Respiratory pumping: neuronal control of a centrally commanded behavior in Aplysia, Brain Research, 143 (1978) 87-105. 3 Cobbs, J. and Pinsker, H., In vivo responses of paired giant mechanoreceptor neurons in Aplysia abdominal ganglion, J. Neurobiol., 9 (1978) 121-141. 4 Coggeshall, R. E., A light and electron microscope study of the abdominal ganglion of Aplysia californica, J. NeurophysioL, 30 (1967) 1263-1287. 5 Coulter, J. D., Sensory transmission through lemniscal pathway during voluntary movement of the cat, J. Neurophysiol., 37 (1974) 831-845. 6 Gerstein, G. L., Functional association of neurons: detection and interpretation. In G. C. Quarton, T. Melnechuk and G. A. Adelman (Eds.), The Neurosciences: Second Study Program, Rockefeller Univ. Press, 1970, pp. 648-660. 7 Glaser, E. M., Separation of neuronal activity by waveform analysis, Adv. Biol. Med. Eng., 1 (1970) 77-135. 8 Glaser, E. M. and Ruchkin, D. S., Principles of Neurobiological Signal Analysis, Academic Press, 1976, pp. 471. 9 Hartline, H. K. and Graham, C. H., Nerve impulses from single receptors in the eye, J. cell. comp. Physiol., 1 (1932) 277-295. 10 Heetderks, W. J. and Williams, W. J., Partition of gross peripheral nerve activity into single unit responses by correlation techniques, Science, 188 (1975) 373-374. 11 McCann, G. D., Interactive computer strategies for living nervous system research, IEEE Trans. Biomed. Engr., 20 (1973) 1-11. 12 Mishelevich, D. J., On-line real-time digital computer separation of extrace!lular neuroelectric signals, IEEE Trans. Biomed. Engr., 17 (1970) 147-150. 13 Moore, G. P., Segundo, J. P., Perkel, D. H. and Levitan, H., Statistical signs ofsynaptic interaction in neurons, Biophys. J., 10 (1970) 876-900. 14 O'Connell, R. J., Kocsis, W. A. and Schoenfeld, R. L., Minicomputer identification and timing of nerve impulses mixed in a single recording channel, Proc. 1EEE, 61 (1973) 1615-1621. 15 Olds, J., Disterhoft, J. F., Segal, M., Kornblith, C. I. and Hirsch, R., Learning centers of rat brain mapped by measuring latencies of conditioned unit responses, J. Neurophysiol., 35 (1972) 202-219. t6 Perkel, D. H., Spike trains as carriers of information. In G. C. Quarton, T. Melnechuk and G. A. Adelman (Eds.), The Neurosciences: Second Study Program, Rockefeller Univ. Pre~s, 1970, pp. 587596. 17 Phillips, M. I. (Ed.), Brain Unit Activity During Behavior, Charles C. Thomas, 1973, pp. 360. 18 Pinsker, H. M., Cobbs, J. and Kanz, J., Neural correlates of the siphon withdrawal reflex in intact, unrestrained Aplysia, F.A.S.E.B., 35 (1976) 769. 19 Pinsker, H. M., Cobbs, J. and Kanz, J., Neuronal correlates of short-term habituation and sensitization of the siphon withdrawal response in freely-behaving Aplysia, Soc. Neurosci. 6th Ann. Meet., 1976, 353. 20 Pinsker, H. M., Feinstein, R., Sawada, M. and Coggeshall, R., Anatomical basis for an apparent paradox concerning conduction velocities of two identified axons in Aplysia, J. Neurobiol., 7 (1976) 241-253. 21 Roberts, W. M. andHartline, D.K.,Separationofmulti-unitnerveimpulsetrainsbymulti-channel linear filter algorithm, Brain Research, 94 (1975) 141-149. 22 Roberts, W. M. and Hartline, H. K., Multi-channel nerve impulse separation techniques, In P. B. Brown (Ed.), Computer Technology in Neuroscience, John Wiley, New York, 1976, pp. 221-236.

479 23 Sim•n••.•Therea•-time••rting•fneur•-••ectricacti•np•t•ntia•sinmu•tip••unitstudies• Electroenceph, clin. Neurophysiol., 18 (1965) 192. 24 Soucek, B. and Carlson, A. D., Computers in Neurobiology and Behavior, John Wiley, New York, 1976, pp. 324. 25 Strong, D. R., Tatton, W. O. and Crapper, D. R., Solid-state amplitude discriminator of neural units, IEEE Trans. Biomed. Eng., 18 (1971) 237-240. 26 Wine, J. and Krasne, F. B., The organization of escape behavior in crayfish, J. exp. Biol., 56 (1972) 1-18. 27 Kanz, J. E., Camp, C. and Pinsker, H., Computer separation of unitary spikes from whole-nerve recordings, Fed. Proc., 37 (1978) 219. 28 Kanz, J. E., Eberly, L. B., Cobbs, J. S. and Pinsker, H. M., Neuronal correlates of siphon withdrawal in freely-behaving Aplysia, J. Neurophysiol., in press.