Abstract
The crosscorrelation coefficient between neural spike trains is a commonly used tool in the study of neural interactions. Two wellknown complications that arise in its interpretation are1) modulations in the correlation coefficient may result solely from changes in the mean firing rate of the cells and2) the mean firing rates of the neurons impose upper and lower bounds on the correlation coefficient whose absolute values differ by an order of magnitude or more. Here, we propose a modelbased approach to the interpretation of spike train correlations that circumvents these problems. The basic idea of our proposal is to estimate the crosscorrelation coefficient between the membrane voltages of two cells from their extracellular spike trains and use the resulting value as the degree of correlation (or association) of neural activity. This is done in the context of a model that assumes the membrane voltages of the cells have a joint normal distribution and spikes are generated by a simple thresholding operation. We show that, under these assumptions, the estimation of the correlation coefficient between the membrane voltages reduces to the calculation of a tetrachoric correlation coefficient (a measure of association in nominal data introduced by Karl Pearson) on a contingency table calculated from the spike data. Simulations of conductancebased leaky integrateandfire neurons indicate that, despite its simplicity, the technique yields very good estimates of the intracellular membrane voltage correlation from the extracellular spike trains in biologically realistic models.
INTRODUCTION
Sensory information and motor plans in many organisms are represented by the joint activity of neurons in a population (see e.g., Abbott and Sejnowski 1999). To study the association between the activity of neuronal pairs, investigators have often relied on the joint peristimulus time histogram (JPSTH) (Aertsen et al. 1989) or the crosscorrelogram (Perkel et al. 1967). Both measures become a correlation coefficient with appropriate normalization. The result of such normalizations are usually called the normalized JPSTH (or nJPSTH) and the normalized shuffledcorrected crosscorrelogram, also called the normalized crosscovariogram (Brody 1999). The correlation coefficient in the nJPSTH is a function of two time variables,t _{1} andt _{2}, representing the times relative to the onset of the stimulus at which the activities of the cells are considered. The correlation coefficient obtained in the normalized crosscovariogram is a function of a single time variable, τ, representing a relative time lag between the activity of the neurons.
The interpretation of the crosscorrelation coefficients obtained from either the nJPSTH or the normalized crosscovariogram is complicated by two wellknown facts (Palm et al. 1988). First, due to the binary nature of the spike train representation, the crosscorrelation coefficients are not bounded by –1 and +1 (representing perfect anticorrelation and perfect correlation) as is normally the case, but tighter bounds apply, which depend on the mean firing rates of the cells. Second, the upper and lower bounds are not equal; the absolute value of the minimum (negative) possible correlation is much smaller than the maximum (positive) possible correlation. The theoretical bounds of the correlation coefficient are dictated by the equations (see
)
Some of the consequences of the dependence of the correlation bounds on the firing rates of the cells are exemplified by Table1, which shows the theoretical minimum and maximum correlations for four pairs of cells with different mean firing rates. A direct comparison of J_{N} values across pairs can be difficult to interpret. For example, if we calculate J_{N} = +0.4962 for one pair andJ_{N} = +0.7053 for another pair, we might conclude that pair 2 is “more correlated” than pair 1. But if the mean firing rates are 5 and 20 spikes/s for pair 1 and 10 and 20 spikes/s for pair 2, both correlations actually represent their maximum theoretically possible values (see Table 1). Under these circumstances, it is more appropriate to consider the cell pairs as having the same degree of correlated activity.
A second consequence of Eqs. 1 and 2 can be seen in the disparity between theoretically possible maximum and minimum correlation for the same pair calculated for different time bins in the nJPSTH or different time lags in the normalized crosscovariogram, even when the mean firing rates are constant (see Table 1). Comparing the absolute magnitude of a negative correlation and a positive correlation in the same pair is problematic; given the range of firing rates in the cortex, the absolute value of the negative correlation will be at least an order of magnitude smaller than the positive correlation (Aertsen and Gerstein 1985; Palm et al. 1988).
The general dependence of the upper and lower bounds on the crosscorrelation coefficient as a function of the firing rates of the cells is shown in Fig. 1. In Fig.1 A, the upper bounds, calculated from Eq. 1, are represented as a pseudocolor image for pairs of neurons with firing rates ranging from 1 to 100 spikes/s. The parallel diagonal structure in this graph indicates that the upper bound of the correlation coefficient is determined approximately by the ratio between the firing rates (see derivation in the ). Figure1 B shows the lower bounds of the correlation coefficient that result from Eq. 2; the asymmetry between the absolute value of the upper and lower bounds is clear when comparing the scales of A and B. The parallel diagonal structure in this graph indicates that the lower bounds are determined approximately by the product of the mean firing rates (see derivation in the ).
In this report, we propose a new measure of association between neural pairs that circumvents some of the problems associated with the interpretation of the correlation coefficients in the nJPTH and the normalized crosscovariogram. Our method is based on a simple model of the joint firing of neurons in which the correlation of the intracellular membrane potentials can be estimated from their extracellular (spike) activity. The key idea of our proposal is to use the estimated intracellular correlation between the membrane potentials as a measure of association between the activity of the cells instead of the correlation coefficient that results from using the spike trains.
In the remainder of this manuscript, we first describe the model that underlies the proposed analysis. Then we show that the technique provides very good of estimates of the membrane voltage correlation in a more biologically realistic situation when many of the model's assumptions are violated. This was done by simulating a pair of conductancebased, leaky integrateandfire neurons. Details of the mathematical derivation of the method are provided in the .
RESULTS
Tetrachoric correlation—underlying model assumptions
We propose a new measure of association between neural spike trains within the context of a model of their joint activity. First, we assume that the membrane potentials of the cells follow a bivariate normal distribution with correlation coefficient ρ. Second, we assume that each cell fires only if its voltage exceeds a firing threshold (θ′
Figure 2 shows an example of a joint normal distribution of the membrane voltages in a pair of neurons. In this diagram, p
_{11} denotes the probability of both cells firing together in the time bins under consideration, which is obtained by integrating the joint distribution of the membrane potentials over the quadrant in which both thresholds are exceeded. p
_{1·} denotes the probability of cell 1 firing and is obtained by integrating the distribution over the area in which θ′
The data obtained in an extracellular recording experiment can be summarized in a contingency table representing the number of times each of the possible firing patterns occurred out of a total of Ninstances (assumed independent) in the time bins under consideration. These numbers are denoted by N _{00},N _{01},N _{10}, andN _{11}, and they sum to N. With these data, it is appropriate to estimate all three parameters of the model simultaneously using maximum likelihood estimation (Tallis 1962). The resulting estimate of ρ is called the “Gaussian fourfold correlation” or “tetrachoric correlation” coefficient (Pearson 1900; Pearson and Heron 1913). Our proposal is to use ρ as the measure of association between the activity of two cells instead of the crosscorrelation coefficient obtained from the spikes trains themselves. As discussed in the following text, this measure avoids many of the problems associated with the correlation coefficient calculated from the spike trains. A Matlab function that performs maximum likelihood estimation of the parameters of the model given the counts in the contingency table N _{00},N _{01},N _{10}, andN _{11} can be downloaded fromhttp://manuelita.psych.ucla.edu/∼dario/.
Testing the method on leaky integrateandfire neurons
The tetrachoric correlation coefficient is only expected to equal the intracellular crosscorrelation coefficient when the model assumptions are satisfied. However, real neurons may violate both the joint normality of the voltage distribution and voltage thresholding as the spikegeneration mechanism. To evaluate the proposed method in a more biologically realistic situation, we conducted simulations using two identical conductancebased leaky integrateandfire neurons. The membrane potential (V) of each neuron evolves according to the equation
Data analysis
Using the (time stationary) simulated data from the leaky integrateandfire neurons, we calculated the entries in a contingency table by counting the total number of times each cell fired simultaneously (N
_{11}), the number of times only one of them fired (N
_{10} andN
_{01}), and the number of times when neither cell fired (N
_{00}). The sum of these numbers, using a time bin of 1 ms, was N = 9 × 10^{5}. Given these counts, we calculated the maximum likelihood estimate of ρ and estimated 95% confidence intervals of ρ by data resampling (Efron and Tibshirani 1993). The extracellular correlation coefficientJ_{N}
was calculated according to
Simulation results
Three simulation series, in which the cells had different mean firing rates, were run on the pair of leaky integrateandfire neurons. The joint normal distributions from which excitatory and inhibitory conductances were drawn had means of
In Fig. 3 A, spike crosscorrelation coefficient (J _{ N }) values are shown against actual membrane potential correlation (ρ_{V} ) for simulation series 1 (thick line) and 3 (thin line). The curves are monotonic but strongly nonlinear as expected from the asymmetry in positive and negative bounds of J_{N} . At any one intracellular correlation value, different firing rates may produce different spike crosscorrelations. Similarly, Fig. 3 A shows that the same extracellular correlation coefficient can be produced by different combinations of intracellular membrane correlations and firing rates. These relationships between J_{N} and ρ_{V} are also predicted by our model (seeappenidx).
Figure 3, B–D, illustrates the relation between the tetrachoric correlation coefficient ρ and the membrane voltage correlation ρ_{V} for each of the three simulation series. The points lie very close to the unity line in all cases. Thus the tetrachoric correlation ρ is a very good estimator of the membrane voltage correlation, ρ_{V}. It can be seen that confidence intervals get larger as the membrane voltage correlations become more negative, implying that estimates of negative correlations are less precise than positive correlations; this occurs when simultaneous counts are very small. Simulations yieldingN _{11} = 0 (obtained for large negative values of ρ) were ignored as the tetrachoric correlation cannot be estimated in this situation.
Simulations shown in Fig. 3 were run with Gaussian whitenoise input conductances, and the correlation coefficients were calculated only at zero time lag. We also performed some simulations where the inputs driving the conductance of the model were filtered in various ways. A typical outcome of these simulations is depicted in Fig.4. In this graph, we compare the spike crosscorrelation (ρ_{V}, dotted line), the tetrachoric correlation coefficient with 95% confidence intervals (ρ, thin line), and the membrane potential correlation (J_{N} , thick line). The tetrachoric correlation slightly underestimates the membrane voltage correlation. We suspect these biases may be due to the fact our model is “static” and does not include any dynamics in it. Nevertheless, the tetrachoric analysis captured the overall shape of the temporally dynamic correlation, especially when compared with the spike crosscorrelation (dotted line in Fig. 4), resulting in a good estimate of membrane voltage correlations through time in these more complex cases.
These simulations suggest that, despite its simplicity, the tetrachoric coefficient provides a good estimate of the intracellular correlation in conditions where the assumptions of the model are violated. One clear violation of the model's assumption is that the joint distributions of membrane potentials in the simulated neurons are not Gaussian (Fig. 5). Figure 5, A–D, shows membrane voltage distributions from simulations run with the temporally dynamic conductance inputs as in Fig. 4 (see Simulation results); firing rates were 21.3 and 11.5 spikes/s. Figure 5 A is a density plot of the joint distribution; B is the corresponding contour plot. Figure 5,C and D, shows probability densities for marginal distributions of membrane voltage for cells 1 and 2, respectively. The departure from normality, in both the joint and marginal distributions, is marked. A lilliefors test rejected the null hypothesis of normality for both marginal distributions (P < 0.01). Figure 5,E–H, shows membrane voltage distributions from simulations run with Gaussian whitenoise conductance inputs with correlation 0.8; firing rates were 13.7 and 13.6 spikes/s. Departure from normality is less marked in these data, as the greatest density of points appears to be distributed normally. However, the marginal probability densities are negatively skewed, and a lilliefors test rejected normality for each (P < 0.01). Figure 5, I–L,shows membrane voltage distributions from simulations run with Gaussian whitenoise conductance inputs with correlation –0.8; firing rates were 13.5 and 13.6. Marginal distributions are again skewed, and normality is rejected (P < 0.01).
DISCUSSION
The question about the “proper” way of normalizing the JPSTH and the crosscorrelogram has a long history going back to the pioneering work of Perkel et al. (1967). As pointed out by Ito and Tsuji (2000), the crux of this problem is defining the most “appropriate” measure of correlation in a 2 × 2 contingency table with the frequenciesp _{00},p _{01},p _{10}, andp _{11}. Interestingly, a similar discussion occurred in statistics, where a heated debate on the “proper” way of defining measures of association in a 2 × 2 contingency table took place between G. Udny Yule and Karl Pearson (Agresti 1990; Pearson 1900;Pearson and Heron 1913; Yule 1912). Yule advocated nonparametric measures of association such as his coefficient of association and coefficient of colligation (Yule 1912). Along these lines, a number of nonparametric methods to “normalize” the crosscorrelogram and the JPSTH have been proposed (Bair et al. 2001; Brosch and Schreiner 1999; Eggermont and Smith 1996; Epping and Eggermont 1987; Ghose et al. 1994;Kruger and Aiple 1988; Palm et al. 1988;Roe and Ts'o 1999; Salinas and Sejnowski 2000; van den Boogaard 1996). The relative advantages and disadvantages of these different measures have been reviewed elsewhere (Aertsen et al. 1989;Brillinger 1976; Hubalek 1982;Palm et al. 1988).
On the other side of the statistics debate, Karl Pearson favored a parametric model of association where the proportions in the contingency table arose by the thresholding of a continuous underlying distribution (Pearson 1900; Pearson and Heron 1913). The elegant work of Ito and Tsuji (2000)is one of the few examples of a parametric framework measuring the association between cell pairs. These authors convincingly argued that it is not possible to derive a modelfree normalization of JPSTH that corrects for the influence of firing modulations. In their work, they adopted a model for the interaction between the spike trains (see alsoAertsen and Gerstein 1985) and devised a procedure to identify the interaction parameter. Our model differs from theirs in that we do not postulate a specific interaction among point processes (the spike trains). Instead, we put forward a simple model for the physical generation of the spikes, where correlations in the spike trains arise by virtue of correlations in the membrane potentials of the neurons. An additional difference is that our model considers both cells on an equal footing (the tetrachoric correlation is “symmetric” on both cells), while spike interaction models require one cell to be identified as “driving” a second one. Finally, the tetrachoric correlation technique has the advantage that it separates the contribution of spike rate modulation and of the membrane voltage correlation to the nJPSTH, allowing not only for comparisons of correlations within the same pair at different times, but allowing the investigator to pool data across a population of cell pairs with different firing rates (see ).
The model presented here also may help clarify basic aspects of the relationship between the extracellular (spike) and intracellular (voltage) correlations, J_{N} and ρ_{V}. Figure 3 A shows that the relationship between these variables is highly nonlinear and accelerating; this relationship is predicted even in our simple model (Fig. 6). In a recent study, Lampl et al. (1999) studied experimentally the relationship between the membrane voltage and spike correlations of neurons in vivo in cat area 17 and found that the extracellular covariograms (based on the spikes) tend to be narrower in time than the intracellular covariograms. Part of their results may be explained by the accelerating nonlinear mapping in Figs. 3 A and 6, which would predict such an effect.
We have shown that a simple model of neural firing in cell pairs reduces to the mathematical framework introduced by Pearson to derive the tetrachoric correlation coefficient, which has since been used as a measure of association in binary data (Agresti 1990;Everitt 1992). The mean firing rates of the neurons are taken into account in the method by the threshold parameters, and the bounds on the tetrachoric correlation coefficient are always –1 and +1. This makes the comparison of correlation coefficients within and between pairs straightforward with no further normalization necessary. Although it is based on an extremely simple model of spike generation, we have shown that the tetrachoric correlation coefficient is a good estimate of intracellular correlation when applied to leaky integrateandfire neurons in situations that violate many of the model's assumptions. The distributions of membrane voltage in our simulations departed from strict normality to different degrees (Fig.5), and we observed a slight underestimation of the intracellular correlations in an extreme situation where normality was clearly violated (Figs. 4 and 5, A–D). Published data suggest that the distribution of membrane voltage during spontaneous activity may well be approximated by a Gaussian distribution at least in simple cells (Anderson et al. 2000). Thus we believe that the tetrachoric correlation can offer a natural framework for the study of neural correlations. We would like to emphasize that the key idea behind our proposal is to use an estimate of the intracellular membrane voltage correlation as the measure of association, and there may be other ways of obtaining such estimates that improve on the tetrachoric correlation. Adding some simple dynamics to the model may be one way to proceed in the future.
We end with a word of caution. Additional work is needed to establish with complete confidence that the tetrachoric correlation is a reasonable method to apply in real neurons. Evaluating the method in more complex models of real neurons is one possibility. However, we feel actual recordings from cells in a slice preparation, where one could study spontaneous correlated activity or inject currents to induce various degrees of correlation, may provide a solid foundation to test the performance of competing methods that estimate intracellular correlations from extracellular spike trains.
Acknowledgments
This work was supported by National Eye Institute Grant EY12816 to D. L. Ringach, National Science Foundation Grant IBN9720305 to D. L. Ringach, and a Research to Prevent Blindness grant to the Jules Stein Eye Institute at UCLA.
Footnotes

Address for reprint requests: D. L. Ringach, Dept. of Neurobiology and Psychology, Franz Hall, Rm 8441B, University of California, Los Angeles, CA 900951563 (Email: dario{at}ucla.edu).
 Copyright © 2003 The American Physiological Society
Appendix
Definition of the problem
Assume we collect the spike trains of two neurons in a number of identical experimental trials. Neural responses are represented as binary sequences of 0s and 1s that result from binning time at a fine scale (Perkel 1967). Following the notation ofBrody (1999), we denote byS
Model description
The main idea of our proposal is to interpret the relative frequencies within the context of a simple model of the joint activity of the neurons. First, we assume that each cell has an associated membrane potential, denoted by V
_{1} andV
_{2}, which are jointly normal with density
If one had exact knowledge of the probabilities in Eq. EA2, a reasonable strategy for estimating the model parameters would be to first identify the thresholds for each cell by inverting Eqs.EA9
and
EA10, which yields
Predicted relationship between intracellular and extracellular correlation
Given two fixed firing rates, the model predicts a particular relationship between the intracellular membrane correlation, ρ, and the extracellular correlation of the spikes as defined byJ_{N} . This relationship helps to clarify some of issues that arise in the interpretation ofJ_{N} . For example, consider an instance where two cells are firing at 10 and 20 spikes/s. The relationship between ρ and J_{N} predicted by the model is given by the thin curve in Fig. 6. It can be seen that the function is monotonic but strongly nonlinear as expected from the asymmetry in positive and negative bounds of J_{N} . The thick line in Fig. 6 illustrates the predicted relationship for a combination of firing frequencies of 2 and 20 spikes/s. These two curves make evident the fact that a particular value ofJ_{N} (for example,J_{N} = 0.1) can be obtained by different combinations of the mean firing rates and intracellular voltage correlation (horizontal dashed line in Fig. 6). Similarly, a fixed intracellular voltage correlation (for example, ρ = 0.5) can lead to different values of J_{N} depending on the mean firing rates of the cells (vertical dashed line in Fig. 6). Thus a modulation ofJ_{N} (t _{1},t _{2}) overt _{1} andt _{2} can be caused purely by changes in the correlation between the membrane potentials, purely by changes in the instantaneous firing rates of the neurons or both.
Our method provides a way to separate the contributions of firing rate modulation and membrane voltage correlation toJ_{N}
(t
_{1},t
_{2}). This can be done by estimating the parameters of the model at each pair of times,t
_{1} andt
_{2}. The result of this procedure are three functions, θ_{1}(t
_{1},t
_{2}), θ_{2}(t
_{1},t_{2}
), and ρ(t
_{1},t
_{2}). The function θ_{1}(t
_{1},t_{2}
) is expected to vary mainly alongt
_{1} and be constant along thet
_{2} axis. Thus we can estimate the mean firing rate of cell 1 by averaging acrosst
_{2}, θ_{1}(t
_{1}) =