Abstract
Experimental and clinical applications of extracellular recordings of spiking cell activity frequently are used to relate the activity of a cell to externally measurable signals such as surface potentials, sensory stimuli, or movement measurements. When the external signal is timevarying, correlation methods have traditionally been used to quantify the degree of relation with the neural firing. However, in some circumstances correlation methods can give misleading results. A new algorithm is described that estimates the extent to which a spike train is related to a continuous timevarying signal. The technique calculates the probability of generating a spike train with Poisson statistics if the timevarying signal determines the Poisson rate. This is accomplished by successive division of the signal and the spike train into halves and recursive calculation of the probability of each halfsignal. The performance of the new algorithm is compared with the performance of correlation methods on simulated data.
INTRODUCTION
When making extracellular recordings of spike trains for electrophysiology experiments, it is usually not possible to state preciselywhat is coded by the neuron or to what degree the neuron encodes any particular experimental variable. As a result, it is often more fruitful to compare several variables to different neural signals to determine which is most closely related to each neuron's activity. For example, extracellular recording in deep structures of the brain may be compared with the surface potentials recorded at many locations on the scalp to determine the topography of cortical projections. In such a case, it becomes important to determine which signals from scalp electrodes are most closely related to the signal from the deep electrode. Other applications commonly occur when using surface electromyographic (EMG) or quantitative movement kinematics measurements to determine the motor effects of a neuron, or when using complex visual or somatosensory stimuli to discover the pattern or modality sensitivity of a neuron.
Most current techniques to address such questions are based on computing the correlation between a smooth timevarying signal and a numerical representation of the neural spike train (Fu et al. 1995; Glaser and Ruckin 1976; Perkel et al. 1967a,b; Werner et al. 1997). The spike train may be represented by the number of spikes in a sequence of short time bins (the binned rate) or by the reciprocal of the interval between successive spikes (the instantaneous rate). In some cases, the numerical spike sequences are smoothed with a lowpass filter (Miller and Houk 1995; Miller et al. 1993). For short bin widths, the use of a lowpass filter applied to the binned rate is equivalent to assigning a smooth “bump” to the location of each spike (the kernel method). Correlation methods are the basis for attempts at realtime “decoding” of neural spike trains (Bialek et al. 1991), and when performed on subbands of the frequency decomposition of the signal lead to coherence analyses (Rosenberg et al. 1989). However, since the original spike train is a sequence of discrete events that occur at specific times, the correlation operation may not be appropriate (Moore et al. 1966; Shao and Chen 1987; Shao and Tsau 1996). By assigning the data to arbitrary bins or by lowpass filtering, the numerical representation may lose information in the signal.
A more important issue depends on the mathematical foundation of the correlation operation. Correlation normalizes the power in the two signals by dividing by the square root of the sum of the squared values at each time point (often called rootmeansquare or rms normalization). Normalization compensates for arbitrary changes in the gain or amplitude of one signal compared with another. Normalization serves a very important purpose in the correlation operation, since it favors signals that are spread out in time. In other words, given two signals with equal integral, the more sharply peaked signal will have a higher meansquared power and thus will be penalized in the correlation when divided by its power. Therefore one interpretation of the correlation between two signals is that it is maximized when1) the two signals have high values at the same time, and2) each signal is spread out as much as possible.
For a binned neural spike train, rms normalization may not be a meaningful operation, since the amplitude is set arbitrarily by the experimenter (by assigning a value of 1 to each spike). In this case, rms normalization divides the number of spikes in each bin by the square root of the total number of spikes. For small bin widths, this operation does not favor signals that are spread out in time. One possible effect of this is shown in Fig.1. In this figure, two smooth signals are used to generate two random Poisson spike trains with timevarying rate (see methods below). Because all the spikes in Fig.1 A occur near the highest peak of the signal in Fig.1 B, the correlation between the smooth signal inB and the spike train in A (using 1ms time bins) is higher than the correlation between the smooth signal inB and the spikes it generated itself. This effect occurs because the normalization of the spike train does not favor the more spread out spike train in Fig. 1 B, since there is no penalty for placing all the spikes at the maximum of the signal. This problem also occurs when using the instantaneous spike rate (Carvell et al. 1996). It can be alleviated by using a larger bin width so that there is a large number of spikes in each bin, but in this case the operation loses temporal resolution and may have difficulty with rapidly changing signals. Performing the correlation in shorter timewindows (Schwartz and Adams 1995) worsens the problem. The alternative technique of using a peristimulus time histogram (PSTH) (Gerstein and Kiang 1960;Gerstein and Perkel 1969; Murphy et al. 1974) is only applicable for external signals related to discrete events which can be easily repeated for averaging.
The new technique proposed here is based on probabilistic inference, rather than correlation. The approach is to determine which of several external signals is related to each possible spike train by calculating the conditional probability that each spike train could have been generated by a Poisson process with timevarying rate given by the external signal. This method is applicable to choosing between a set of candidate signals, but the absolute probability values are so small that they are difficult to interpret outside this context. To avoid the necessity of choosing arbitrary bin sizes, the probability is computed using a recursive multiscale approach, with the timewindow being successively divided in half and probabilities computed over smaller and smaller segments. Once the conditional probability has been calculated, choosing the external signal with the largest probability yields the maximumlikelihood (ML) estimate of the best match between continuous signals and spike trains.
METHODS
A neural spike train n(t) is modeled as a timevarying Poisson process, with the probability of firingn times in a small interval Δt given by
To calculate the probabilityP[n(t)‖r(t)], first adjust the baseline of r(t) by adding a constant so that r(t) is always >0 and normalizer(t) by dividing by its integral so that the total integral is 1. This must be done so thatr(t) can be interpreted as a probability density. Then split n(t) and r(t) into left and right halves, so thatn_{L}
(t) is the left half ofn(t) from time 0 to T/2,n_{R}
(t) is the right half from time T/2 to T, andr_{L}
(t) andr_{R}
(t) are defined similarly. Count the number of spikes inn_{L}
(t) andn_{R}
(t) and call thesen_{L}
and n_{R}
. Integrate the continuous signalsr_{L}
(t) andr_{R}
(t) and call the integralsr_{L}
and r_{R}
. Because of the independent increment property of Poisson processes, the probability of obtaining exactly n_{L}
spikes over interval (0, T/2) is determined by the Poisson rater_{L}
. To calculate the probability of “splitting” the total number of spikes into the two halfsignals, note that since r_{L}
+r_{R}
= 1, the probability of the split is the probability of dividing n spikes into two bins (the left and the right halves) of exactly n_{L}
andn_{R}
spikes. This is given by the standard equation for the binomial distribution
We now continue this process by normalizingr_{L} (t) andr_{R} (t) to each have integral 1, splitting each halfinterval again in half (yielding quarterintervals) and calculating the probability of the two new splits. A bin that contains only zero or one spikes is not split. The process continues for each interval until there are only zero or one spikes in each interval. The probability of the spike train n(t) given the rate r(t) is then the combined probability of making all the specific splits, and since the splits are independent this probability is given by the product of the probabilities of all the individual splits. The resulting number will be very small (the probability of any particular spike train will always be small), but it allows comparison of the relative probability that different rate functions r(t) might have led to the particular spike trainn(t). Given a series of candidate rate functionsr_{i} (t), we perform this procedure on each and select the rate function with the largest probabilityP[n(t)‖r_{i} (t)]. This gives the ML estimate.
The computation of the probabilities is very rapid, but it can be made even more rapid by noting that, since the final result is a comparison, we do not need to calculate an accurate value of the probability so long as we are certain that the rankings of the candidate rate functions r_{i}
(t) will remain in order. Note that the term
Function prob [n(t),r(t)]:
1. If n(t) = 0 return
2. Normalize: r(t) →r(t)/∫ r(t)
3. Split n(t) and r(t) atT/2 into n_{L} (t),n_{R} (t),r_{L} (t)_{,}and r_{R} (t)
4. Count spikes and integrate the signalsn_{L} = ∑n_{L} (t),n_{R} = ∑n_{R} (t),r_{L} = ∫r_{L} (t)_{,}and r_{R} = ∫r_{R} (t)
5. Calculate the log split probability factor P =n_{L} log r_{L} +n_{R} logr_{R}
6. Call recursively p_{L} = prob [n_{L} (t),r_{L} (t)] andp_{R} = prob [n_{R} (t),r_{R} (t)]
7. Return p + p_{L} +p_{R}
(Matlab code for the algorithm and a demonstration using simulated data can be downloaded from http://www.stanford.edu/∼sanger/demo1.zip.)
For each simulation, 50 continuous rate functionsr(t) were each generated by adding a series of five cosine functions with random amplitude and phase, and periods varying from one cycle per total time T to five cycles per total time T. A constant was added so that the minimum value of each r(t) was greater than zero. For each rate function r(t), a single Poisson spike trainn(t) was generated by placing spikes pseudorandomly within each time interval, with probability of firing proportional to the average value of the rate function in that interval. Thus for each comparison, 50 different 1s sets of spike counts n(t) were generated from the 50 generating rate functions r(t), with bins sampled at 1000 Hz. Bin sizes are therefore 1 ms, for a total of 1000 bins. The probability algorithm was run on all 2500 possible comparisons [n(t), r(t)], and for each spike train n(t), the rate functionr(t) with the largest resulting log probability was chosen as the best estimate. The best estimate is correct ifn(t) was generated from ther(t) with highest log probability, and otherwise incorrect. The number of correct associations ofn(t) and r(t) for each simulation is therefore between zero and 50. The entire procedure was performed 50 times, and the mean and SD of the number of correct associations were calculated.
For comparison, two correlationbased methods were tested on the same rate function and spike data. In the “bin correlation” method, each spike train n(t) is represented as a series of spike counts in each of the one thousand 1ms time bins. Each rate function is then correlated against the spike train by calculating
RESULTS
Figure 2 A shows a comparison between the probability method and the two different correlationbased methods for estimating the correct spike train associated with a given smooth signal. In this case, the average spike rate is 20 Hz. The mean number correctly categorized by the probability method is 26 (SD 3.6); the mean for the bin correlation method is 22 (SD 4.0), and the mean for the interval correlation method is 21 (SD 3.3). The probability method result is significantly different from both correlation methods by pairedsamples ttest (P < 0.001 in both cases). Although none of the methods was able to categorize more than about half the examples correctly, it should be noted that with an average spike rate of 20 Hz there are only about 20 spikes in the timewindow, so that this is a very difficult categorization problem with very little data to distinguish between possibilities. Also note that all algorithms do much better than chance performance, for which the expected number correctly categorized is 1.
Figure 2 B shows a comparison for an average spike rate of 100 Hz. In this case, the mean number correctly categorized by the probability method is 49 (SD 0.96), the mean for the bin correlation method is 48 (SD 1.2), and the mean for the interval correlation method is 48 (SD 1.1). The probability method is significantly different from both correlation methods by pairedsamples ttest (P < 0.001 in both cases), although in this case the magnitude of the difference in performance is negligible. The two correlation methods are not significantly different from each other.
DISCUSSION
The probability technique is a new method for determining which of several smooth signals is most likely to be related to a train of neural spikes. It is based on the assumption that the spikes can be described by a Poisson process with rate given by the correct smooth signal. Under this assumption, an algorithm has been developed that recursively calculates the probability of producing a particular spike train, and the signal with maximum probability determines the maximum likelihood estimate. This is the best possible estimate under these assumptions and with no additional available information. There is no need to choose bin widths a priori, since the recursive splitting algorithm dynamically selects bin widths by splitting until there is no more than one spike per bin. If there is an unknown lag between the spike train n(t) and the continuous signalr(t), the probability algorithm can be used to compare the signal at varying lags to determine the maximum likelihood solution. This is illustrated in Fig.3 B, which shows that the peak of the probability as a function of time lag occurs near the simulated correct value of 500 ms. If multiple spike trains are tested at multiple lags, it is possible that different spike trains would match a particular continuous signal at different lags, and such a result would need to be interpreted with caution.
As shown above, the probability algorithm outperforms the bin correlation method and the interval correlation method on simulated data that conforms to its underlying assumption. The difference is most noticeable for lower spike rates. It is not known whether the probability algorithm would perform better on real data, but this is likely to require verification on a casebycase basis, since not all neural spike trains will be well described by timevarying Poisson processes. In practice, it may be necessary to verify that a spike train has nearly Poisson statistics (with variance in spike counts between bins equal to the mean spike rate) to be confident that this algorithm will apply.
The main disadvantage of the new technique is the relative complexity of the equation and the moderately increased computational requirements compared with the correlation method. Another disadvantage is that the algorithm cannot be run as a realtime filter, and it is therefore not useful for neural control applications. The algorithm analyzes only the information in the spike train that is carried by the spike rate, and it may therefore underestimate relations due to other coding strategies (Sherry and Klemm 1984). The rate function is adjusted to be positive and its integral is normalized, so information in the DC baseline or the signal magnitude is not used.
It is important to realize that the absolute magnitude of the calculated probabilities is a small number that may not have an easily interpretable meaning. The probabilityP[n(t)‖r(t)] is therefore used to compare the relative probabilities between different signals r(t). For example, a spike train n(t) recorded from a single cell using an extracellular electrode can be compared with a set of surface potentials r(t) on the scalp to determine the spatial location on the scalp that is most closely related to the cell's firing rate. For each scalp electrode signalr(t), the value ofP[n(t)‖r(t)] is calculated, and the signal r(t) that maximizes this value can be interpreted as the most closely related scalp surface signal to the spike train n(t). A similar experiment could be performed to test which of several rectified and filtered surface EMG signals is related to an extracellular recording. Note that although the Poisson model treats r(t) as the “generator” of n(t), this does not imply causality. In other words, it is not possible to conclude thatr(t) drives the spike trainn(t), but only that the two are related.
There are many possible variations on the basic structure of the algorithm. Rather than splitting at onehalf the time window, the splits could be placed at experimentally relevant time points. Alternatively, the split could be chosen to place an equal number of spikes on both sides, or each split could form a division of the smooth signal into two orthogonal waveletlike components with a smooth boundary.
The utility and validity of this technique need to be tested in particular experimental situations. Nevertheless, under this model of spike generation it is possible to determine the probabilistically correct solution for associating a spike train with a continuous signal.
Acknowledgments
I thank S. Giszter and P. Cisek for helpful discussions and comments.
This work was supported in part by the Stanford University Department of Child Neurology.
Footnotes

Address for reprint requests: Dept. of Child Neurology, Stanford University Medical Center, 300 Pasteur Dr., MS:5235, Stanford, CA 94305.
 Copyright © 2002 The American Physiological Society