INTRODUCTION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION APPENDIX REFERENCES

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 cross-correlogram (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 shuffled-corrected cross-correlogram, also called the normalized cross-covariogram (Brody 1999). The correlation coefficient in the nJPSTH is a function of two time variables, t₁ and t₂, 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 cross-covariogram is a function of a single time variable, , representing a relative time lag between the activity of the neurons.

The interpretation of the cross-correlation coefficients obtained from either the nJPSTH or the normalized cross-covariogram is complicated by two well-known facts (Palm et al. 1988). First, due to the binary nature of the spike train representation, the cross-correlation coefficients are not bounded by -1 and +1 (representing perfect anti-correlation 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 APPENDIX)

(1)

(2)

where J is the upper bound on the cross-correlation coefficient, J represents the lower bound, p_1· is the probability of cell 1 firing in a small time interval (1-ms bins are used throughout this paper), and p_·1 is the probability of cell 2 firing in a second small interval.

Some of the consequences of the dependence of the correlation bounds on the firing rates of the cells are exemplified by Table 1, 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 and J_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 cross-covariogram, 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 cross-correlation coefficient as a function of the firing rates of the cells is shown in Fig. 1. In Fig. 1A, the upper bounds, calculated from Eq. 1, are represented as a pseudo-color 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 APPENDIX). Figure 1B 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 APPENDIX).

View larger version (45K): [in this window] [in a new window]   Fig. 1. Dependence on mean firing rate of cross-correlation bounds. A: theoretical upper bounds (shown in color code) for pairs with firing rates ranging from 1 to 100 spikes/s (shown on a double logarithmic plot). B: theoretical lower bounds for the same mean firing rates. As shown in the APPENDIX, these bounds can be approximated by J =  and J = , which is the reason for the parallel structure when the bounds are plotted in double logarithmic coordinates.

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 cross-covariogram. 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 conductance-based, leaky integrate-and-fire neurons. Details of the mathematical derivation of the method are provided in the APPENDIX.

	RESULTS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION APPENDIX REFERENCES

Tetrachoric correlationunderlying 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 ('_i, where i = 1, 2); otherwise, the cell does not fire. This simple model of spike generation in a neural pair should be considered more of a conceptual construct rather than a faithful representation of the underlying biophysical processes. Clearly real cells are much more complicated than mere thresholding devices. The advantage of this oversimplification is that within the context of this model, one can estimate the intracellular membrane correlation from the spike data and that, perhaps surprisingly, the technique works remarkably well in more complex situations when the assumptions of the model are violated.

Figure 2 shows an example of a joint normal distribution of the membrane voltages in a pair of neurons. In this diagram, p₁₁ 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 '_i is exceeded, and p₁₀, the probability that cell 1 fired when cell 2 did not is given by p₁₀ = p_1·  p₁₁. Similarly, p_·1 represents the probability of cell 2 firing and is obtained by integration of the distribution over the area in which '₂ is exceeded. p₀₁ represents the probability that cell 2 fired when cell 1 did not and its value is given by p₀₁ = p_·1  p₁₁. The three parameters of the model are: '₁ (the threshold for cell 1), '₂ (the threshold for cell 2), and  (the correlation coefficient between the cells' membrane potentials). Our task is to estimate these parameters given the extracellular spike trains.

View larger version (16K): [in this window] [in a new window]   Fig. 2. The model. The points represent examples of a bivariate normal distribution of the membrane voltages in a pair of cells. '1 is the firing threshold of cell 1; '2 is the firing threshold of cell 2;  is the membrane potential correlation coefficient (equal to +0.5 in this example); p11 is the probability that both cells will fire; p10 is the probability that cell 1 fired and cell 2 did not; p01 represents the probability that cell 2 fired and cell 1 did not; and p00 is the probability that neither cell fired. In our analysis, the parameters of the model, '1, '2, and  are simultaneously estimated using maximum likelihood estimation.

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 N instances (assumed independent) in the time bins under consideration. These numbers are denoted by N₀₀, N₀₁, N₁₀, and N₁₁, 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 cross-correlation 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₀₀, N₀₁, N₁₀, and N₁₁ can be downloaded from http://manuelita.psych.ucla.edu/~dario/.

Testing the method on leaky integrate-and-fire neurons

The tetrachoric correlation coefficient is only expected to equal the intracellular cross-correlation 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 spike-generation mechanism. To evaluate the proposed method in a more biologically realistic situation, we conducted simulations using two identical conductance-based leaky integrate-and-fire neurons. The membrane potential (V) of each neuron evolves according to the equation

(3)

where C is the capacitance, g_e, g_i, and g_leak are the excitatory, inhibitory, and leak conductances, respectively, and V_e, V_i, and V_leak are the excitatory, inhibitory, and leak reversal potentials. We followed the work of Wielaard et al. (2001) and used the following values for the parameters: C = 10⁶ F cm², g_leak = 50 × 10⁶ S cm², V_leak = 70 mV, V_e = 0 mV, V_i = 80 mV, and V_th = 55 mV. We departed from Wielaard et al. (2001) and used V_reset = 66.3 mV. This was done to minimize the large transient observed in the membrane potential traces with a reset equal to V_leak that is not observed in actual intracellular recordings. We simulated a stationary process where the excitatory and inhibitory conductances were Gaussian white noise. The SDs of the conductances were adjusted between simulation runs to vary the mean firing rates of the model neurons. To induce correlated activity of the model neurons, the excitatory conductances were drawn from a joint normal distribution with specified correlation; the inhibitory conductances were drawn from a separate joint normal distribution with the same correlation. The degree of correlation between the inputs was varied in different runs to yield different degrees of membrane voltage correlation. See Simulation results for actual values used in the simulations. A total of 900 s of data was simulated to match the amount of data we collect in our reverse correlation experiments (Ringach 2002). A Matlab Simulink model implementing the integrate-and-fire model used in these simulations can be downloaded from http://manuelita.psych.ucla.edu/~dario.

Data analysis

Using the (time stationary) simulated data from the leaky integrate-and-fire neurons, we calculated the entries in a contingency table by counting the total number of times each cell fired simultaneously (N₁₁), the number of times only one of them fired (N₁₀ and N₀₁), and the number of times when neither cell fired (N₀₀). The sum of these numbers, using a time bin of 1 ms, was N = 9 × 10⁵. 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 coefficient J_N was calculated according to

(4)

where the probabilities were estimated as p_ij = N_ij/N, for i, j  {0,1} (see also the APPENDIX). The intracellular correlation between the membrane voltages, V₁(t) and V₂(t), was calculated by the usual formula

(5)

where  ·  indicates averaging across time.

Simulation results

Three simulation series, in which the cells had different mean firing rates, were run on the pair of leaky integrate-and-fire neurons. The joint normal distributions from which excitatory and inhibitory conductances were drawn had means of _exc, _inh = 20 × 10⁶ S cm². Firing rates were varied between the series by changing the SD of the excitatory and inhibitory conductances (Table 2). Seventeen individual runs in each series were simulated at excitatory and inhibitory conductance correlations from 0.8 to +0.8, resulting in membrane voltage correlations from 0.25 to 0.6. The membrane voltage correlation coefficient (_V), tetrachoric correlation coefficient (), and the correlation coefficient at zero time lag from the cross-covariogram (J_N) were calculated in each case.

In Fig. 3A, spike cross-correlation 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 cross-correlations. Similarly, Fig. 3A 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 (see APPENIDX).

View larger version (8K): [in this window] [in a new window]   Fig. 3. Simulation results for white-noise input at 0 time lag ( = 0). A: spike train cross-correlation coefficient (JN) against actual membrane voltage correlation coefficient (V) for simulation series 1 and 3. Series 1: thin line, mean spike rates 6.2 and 3.8. Series 3: thick line, mean spike rates 13.5 and 13.5. B-D: plot of the estimated tetrachoric correlation coefficient () against membrane voltage correlation coefficient for series 1-3, respectively. The diagonal line represents the identity line. Error bars represent 95% confidence intervals obtained by data resampling.

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 yielding N₁₁ = 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 white-noise 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 cross-correlation (_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 cross-correlation (dotted line in Fig. 4), resulting in a good estimate of membrane voltage correlations through time in these more complex cases.

View larger version (17K): [in this window] [in a new window]   Fig. 4. Performance of tetrachoric correlation when the conductances are not white. Firing rates were 21.3 and 11.5 spikes/s. Dotted line, spike cross-correlation; thin line, tetrachoric correlation coefficient with 95% confidence intervals; thick line, membrane potential correlation coefficient.

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 5A 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 white-noise 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 white-noise 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).

View larger version (16K): [in this window] [in a new window]   Fig. 5. Joint membrane voltage distributions from leaky integrate-and-fire model neurons. A-D: membrane voltage distributions for the simulation depicted in Fig. 4. A: density plot of joint distribution of membrane potentials. B: contour plot of the joint distribution. C: marginal probability density for the membrane voltage of cell 1. D: marginal probability density for the membrane voltage of cell 2. E-H: membrane voltage distributions of simulation run with Gaussian white noise with values as in series 3 (see Table 2) and conductance correlation = 0.8; firing rates were 13.7 and 13.6 spikes/s. E: density plot of the joint distribution of membrane potentials. F: contour plot of joint distribution. G: marginal probability density for cell 1. H: marginal probability density for cell 2. I-L: membrane voltage distributions of simulation run with Gaussian white noise with values as in series 3 (see Table 2) and conductance correlation = 0.8; firing rates were 13.5 and 13.6 spikes/s. I: density plot of joint distribution. J: contour plot of joint distribution. K: marginal probability density for cell 1. L: marginal probability density for cell 2.

	DISCUSSION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION APPENDIX REFERENCES

The question about the "proper" way of normalizing the JPSTH and the cross-correlogram 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 frequencies p₀₀, p₀₁, p₁₀, and p₁₁. 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 cross-correlogram 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 model-free 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 also Aertsen 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 APPENDIX).

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 3A 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. 3A and 6, which would predict such an effect.

View larger version (11K): [in this window] [in a new window]   Fig. 6. Relationship between cross-correlation coefficient (JN) and membrane voltage correlation () predicted by the model. Thin line: firing rates 10 and 20 spikes/s. Thick line: firing rates 2 and 20 spikes/s. Vertical dashed line:  = 0.5. Horizontal dashed line: JN = 0.1

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 integrate-and-fire 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.