ARTICLE

Estimating Membrane Voltage Correlations From Extracellular Spike Trains

Interdepartmental Program for Neuroscience, Brain Research Institute and

and

Interdepartmental Program for Neuroscience, Brain Research Institute and

Departments of Neurobiology and Psychology, and Jules Stein Eye Institute, University of California, Los Angeles, California 90095

Search for more papers by this author

Published Online:01 Apr 2003https://doi.org/10.1152/jn.000889.2002

Abstract

The cross-correlation coefficient between neural spike trains is a commonly used tool in the study of neural interactions. Two well-known 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 model-based approach to the interpretation of spike train correlations that circumvents these problems. The basic idea of our proposal is to estimate the cross-correlation 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 conductance-based leaky integrate-and-fire 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 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₁ andt₂, 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 )

J_{N}^{max} = \frac{min (p_{1 \cdot}, p_{\cdot 1}) - p_{1 \cdot} p_{\cdot 1}}{\sqrt{p_{1 \cdot} (1 - p_{1 \cdot}) p_{\cdot 1} (1 - p_{\cdot 1})}}

Equation 1

J_{N}^{min} = \frac{- p_{1 \cdot} p_{\cdot 1}}{\sqrt{p_{1 \cdot} (1 - p_{1 \cdot}) p_{\cdot 1} (1 - p_{\cdot 1})}}

Equation 2

where J

_{N}^{max}

is the upper bound on the cross-correlation coefficient,J

_{N}^{min}

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 Table1, which shows the theoretical minimum and maximum correlations for four pairs of cells with different mean firing rates. A direct comparison of J_Nvalues 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.

**Table 1. Examples of upper and lower bounds on the correlation coefficients for different combinations of firing rates**
Rate 1, Spikes/s	Rate 2, Spikes/s	p_1·	p_·1	J $_{N}^{max}$	J $_{N}^{min}$
5	10	0.005	0.010	+0.7053	−0.0071
10	10	0.010	0.010	+1.0000	−0.0101
5	20	0.005	0.020	+0.4962	−0.0101
2	20	0.002	0.020	+0.3134	−0.0064

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 ). Figure1B 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).

**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 , these bounds can be approximated byJ $_{N}^{max}$ = $\sqrt{p_{1 \cdot} / p_{\cdot 1}}$ and J $_{N}^{min}$ = − $\sqrt{p_{1 \cdot} p_{\cdot 1}}$ , which is the reason for the parallel structure when the bounds are plotted in double logarithmic coordinates.

Download figure Download PowerPoint

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.

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 (θ′ $_{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, andp₁₀, the probability that cell 1 fired when cell 2 did not is given byp₁₀ =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 byp₀₁ =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.

**Fig. 2.**The model. The points represent examples of a bivariate normal distribution of the membrane voltages in a pair of cells. θ′₁ is the firing threshold of cell 1; θ′₂ is the firing threshold of cell 2; ρ is the membrane potential correlation coefficient (equal to +0.5 in this example); p₁₁ is the probability that both cells will fire; p₁₀ is the probability that cell 1 fired and cell 2 did not; p₀₁represents the probability that cell 2 fired and cell 1 did not; andp₀₀ is the probability that neither cell fired. In our analysis, the parameters of the model, θ′₁, θ′₂, and ρ are simultaneously estimated using maximum likelihood estimation.

Download figure Download PowerPoint

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₀₀,N₀₁,N₁₀, andN₁₁, 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₁₀, andN₁₁ can be downloaded fromhttp://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

C \frac{d V}{d t} = g_{e} (V_{e} - V) + g_{i} (V_{i} - V) + g_{leak} (V_{leak} - V)

Equation 3

where C is the capacitance,g_e,g_i, andg_leak are the excitatory, inhibitory, and leak conductances, respectively, andV_e,V_i, andV_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, andV_th = −55 mV. We departed fromWielaard et al. (2001) and usedV_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. SeeSimulation 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₁₀ andN₀₁), 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 coefficientJ_N was calculated according to

J_{N} = \frac{p_{11} - p_{1 \cdot} p_{\cdot 1}}{\sqrt{p_{1 \cdot} (1 - p_{1 \cdot}) p_{\cdot 1} (1 - p_{\cdot 1})}}

Equation 4

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

ρ_{V} = \frac{〈 V_{1} (t) V_{2} (t) 〉 - 〈 V_{1} (t) 〉 〈 V_{2} (t) 〉}{\sqrt{〈 V_{1} (t)^{2} 〉 - 〈 V_{1} (t) 〉^{2}} \sqrt{〈 V_{2} (t)^{2} 〉 - 〈 V_{2} (t) 〉^{2}}}

Equation 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 $\bar{g}$ _exc, $\bar{g}$ _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.

**Table 2. Standard deviations of the excitatory and inhibitory conductances in the different simulation series together with the resulting mean firing rates**
	Cell 1			Cell 2
	Excitatory, μS/cm²	Inhibitory, μS/cm²	Firing rate, spikes/s	Excitatory, μS/cm²	Inhibitory, μS/cm²	Firing rate, spikes/s
Series 1	7	4	6.2	6	4	2.8
Series 2	10	4	13.5	6	5	4.2
Series 3	10	10	13.5	4	4	13.5

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 (seeappenidx).

Download figure Download PowerPoint

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₁₁ = 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.

**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.

Download figure Download PowerPoint

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).

**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.

Download figure Download PowerPoint

DISCUSSION

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 frequenciesp₀₀,p₀₁,p₁₀, andp₁₁. 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 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 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.

**Fig. 6.**Relationship between cross-correlation coefficient (J_N) 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:J_N = 0.1

Download figure Download PowerPoint

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.

This work was supported by National Eye Institute Grant EY-12816 to D. L. Ringach, National Science Foundation Grant IBN-9720305 to D. L. Ringach, and a Research to Prevent Blindness grant to the Jules Stein Eye Institute at UCLA.

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 $_{1}^{r}$ (t) andS $_{2}^{r}$ (t) the spike trains of two cells at time t, for trial r. The nJPSTH is defined by (Aertsen et al. 1989)

J_{N} (t_{1}, t_{2}) = \frac{〈 S_{1}^{r} (t_{1}) S_{2}^{r} (t_{2}) 〉 - 〈 S_{1}^{r} (t_{1}) 〉 〈 S_{2}^{r} (t_{2}) 〉}{ς_{1} (t_{1}) ς_{2} (t_{2})}

Equation A1

where 〈·〉 denotes averaging across trials and ς_i(t) =

\sqrt{〈 S_{i}^{r} (t)^{2} 〉 - 〈 S_{i}^{r} (t) 〉^{2}}

is the SD of the response for cell i at time t. J_N(t₁,t₂) is a correlation coefficient that under normal circumstances ranges between –1 (perfect anti-correlation) and +1 (perfect correlation). However, as we discuss in the following text, tighter bounds apply in this case due to the binary nature of the signalsS

_{i}^{i}

(t) (Palm et al. 1988). In the discussion that follows, it helps to fix the values of t₁ andt₂ and consider the activity of the cells at bins t₁ andt₂ at trial r as the realization of two random variables,S₁ andS₂. The joint distribution ofS₁ andS₂ is fully specified by the probabilities

p_{11} = Pr (S_{1} = 1 \land S_{2} = 1)

p_{10} = Pr (S_{1} = 1 \land S_{2} = 0)

p_{01} = Pr (S_{1} = 0 \land S_{2} = 1)

p_{00} = Pr (S_{1} = 0 \land S_{2} = 0)

Equation A2

Note that this distribution has only 3 degrees of freedom as the sum of all probabilities must equal one. Using the joint distribution of S₁ andS₂ in Eq. EA1, we obtain

J_{N} = \frac{p_{11} - p_{1 \cdot} p_{\cdot 1}}{\sqrt{p_{1 \cdot} (1 - p_{1 \cdot}) p_{\cdot 1} (1 - p_{\cdot 1})}}

Equation A3

where p_1· =p₁₁+p₁₀ is the probability of firing for the first cell, and p_·1 =p₁₁+p₀₁ is the probability of firing for the second cell. Then, given two fixed mean firing rates (i.e., firing probabilities), p_1· andp_·1, it is easily seen thatJ_S is bounded above by

J_{N}^{max} = \frac{min (p_{1 \cdot}, p_{\cdot 1}) - p_{1 \cdot} p_{\cdot 1}}{\sqrt{p_{1 \cdot} (1 - p_{1 \cdot}) p_{\cdot 1} (1 - p_{\cdot 1})}}

Equation A4

and bounded below by

J_{N}^{min} = \frac{- p_{1 \cdot} p_{\cdot 1}}{\sqrt{p_{1 \cdot} (1 - p_{1 \cdot}) p_{\cdot 1} (1 - p_{\cdot 1})}}

Equation A5

Simplified expressions for Eqs. EA4 and EA5can be obtained assuming the firing probabilities are relatively low (p_1·,p_·1 ≪ 1) as we usually find in the cortex, and if we consider, without loss of generality, thatp_1· ≤p₁. Then, these bounds can be approximated by J

_{N}^{max}

≈

\sqrt{p_{1 \cdot} / p_{\cdot 1}}

and J

_{N}^{min}

≈ −

\sqrt{p_{1 \cdot} p_{\cdot 1}}

, which combined lead to the expression ‖J

_{N}^{max}

_{N}^{min}

‖ ≈ 1/p_·1. This demonstrates a large asymmetry between the maximum and minimum absolute values attainable byJ_N, as already noted in the work ofPalm et al. (1988).

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₁ andV₂, which are jointly normal with density

f (V_{1}, V_{2}) = \frac{1}{2 π \sqrt{1 - ρ^{2}} ν_{1} ν_{2}} exp \{- \frac{1}{2 (1 - ρ^{2})}

Equation A6

\times (\frac{(V_{1} - μ_{1})^{2}}{ν_{1}^{2}} + \frac{(V_{2} - μ_{2})^{2}}{ν_{2}^{2}} - \frac{2 ρ (V_{1} - μ_{1}) (V_{2} - μ_{2})}{ν_{1} ν_{2}})\}

Here, μ₁ are the mean membrane voltages, ν₁ their SDs, and ρ the correlation coefficient between the membrane voltages (i = 1, 2). Second, we assume that each cell fires if the voltage exceeds a threshold, V_i > θ_i. Otherwise, the cell does not fire. Under the assumptions of the model, the probability that both cells fire simultaneously is given by

p_{11} = \int_{θ_{1}}^{\infty} \int_{θ_{2}}^{\infty} f (V_{1}, V_{2}) d V_{1} d V_{2}

= \frac{1}{2 π \sqrt{1 - ρ^{2}} ν_{1} ν_{2}} \int_{θ_{1}}^{\infty} \int_{θ_{2}}^{\infty} exp \{- \frac{1}{2 (1 - ρ^{2})} (\frac{(V_{1} - μ_{1})^{2}}{ν_{1}^{2}} + \frac{(V_{2} - μ_{2})^{2}}{ν_{2}^{2}} - \frac{2 ρ (V_{1} - μ_{1}) (V_{2} - μ_{2})}{ν_{1} ν_{2}})\} d V_{1} d V_{2}

Equation A7

Substituting V′_i = (V_i − μ_i)/ς_i we obtain

p_{11} = \frac{1}{2 π \sqrt{1 - ρ^{2}}} \int_{θ_{1}^{'}}^{\infty} \int_{θ_{2}^{'}}^{\infty} exp \{- \frac{1}{2 (1 - ρ^{2})} (V_{1}^{' 2} + V_{2}^{' 2} - 2 ρ V_{1}^{'} V_{2}^{'})\} d V_{1}^{'} d V_{2}^{'}

Equation A8

where θ′_i = (θ_i − μ_i)/ν_i. Similarly, the mean probability of firing for each cell is given by

p_{1 \cdot} = \frac{1}{\sqrt{2 π}} \int_{θ_{1}^{'}}^{\infty} exp \{- \frac{V_{1}^{' 2}}{2}\} d V_{1}^{'} = \frac{1}{2} (1 - erf (\frac{θ_{1}^{'}}{\sqrt{2}}))

Equation A9

and

p_{\cdot 1} = \frac{1}{\sqrt{2 π}} \int_{θ_{2}^{'}}^{\infty} exp \{- \frac{V_{2}^{' 2}}{2}\} d V_{2}^{'} = \frac{1}{2} (1 - erf (\frac{θ_{2}^{'}}{\sqrt{2}}))

Equation A10

where

erf (x) = (2 / \sqrt{π}) \int_{0}^{t} e^{- t^{2}} dt

is the error function. Equations EA8-EA10 establish the relationship between the parameters of our model θ′₁, θ′₂, and ρ and the probabilities of firing p_1·,p_·1, andp₁₁.

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

θ_{1}^{'} = \sqrt{2} {erf}^{- 1} (1 - 2 p_{1 \cdot})

θ_{2}^{'} = \sqrt{2} {erf}^{- 1} (1 - 2 p_{\cdot 1})

Equation A11

where erf⁻¹(x) is the inverse of the error function. Once the parameters θ′₁ and θ′₂ have been identified, Eq. EA8 can be solved numerically to obtain the correlation coefficient ρ. However, we do not normally have access to the actual probabilities. In this case, it is appropriate to estimate all the parameters of the model simultaneously using maximum likelihood estimation (Tallis 1962).

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₁,t₂) overt₁ andt₂ 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₁,t₂). This can be done by estimating the parameters of the model at each pair of times,t₁ andt₂. The result of this procedure are three functions, θ₁(t₁,t₂), θ₂(t₁,t₂), and ρ(t₁,t₂). The function θ₁(t₁,t₂) is expected to vary mainly alongt₁ and be constant along thet₂ axis. Thus we can estimate the mean firing rate of cell 1 by averaging acrosst₂, θ₁(t₁) = $\frac{1}{T} \int_{0}^{T}$ θ₁ (t₁,t₂) dt₂. Similarly, we can estimate θ₂(t₂) = $\frac{1}{T} \int_{0}^{T}$ θ₂ (t₁,t₂) dt₁. The functions θ₁(t) and θ₂(t) represent the changes in the firing rates of the cells, whereas ρ(t₁,t₂) represents the change in their intracellular voltage correlation. We refer to ρ(t₁,t₂) as the joint tetrachoric histogram (JTH).

FOOTNOTES

Address for reprint requests: D. L. Ringach, Dept. of Neurobiology and Psychology, Franz Hall, Rm 8441B, University of California, Los Angeles, CA 90095-1563 (E-mail: [email protected]edu).

REFERENCES

1 Abbott LF, Sejnowski TJNeural Codes and Distributed Representations: Foundations of Neural Computation.1999MIT PressCambridge, MA
Google Scholar
2 Agresti ACategorical Data Analysis.1990WileyNew York
Google Scholar
3 Aertsen AMHJ, Gerstein GLEvaluation of neuronal connectivity: sensitivity of cross-correlations.Brain Res3401985341354
Crossref | PubMed | ISI | Google Scholar
4 Aertsen AMHJ, Gerstein GL, Habib MK, Palm GDynamics of neuronal firing correlation: modulation of “effective connectivity.”J Neurophysiol611989900917
Link | ISI | Google Scholar
5 Anderson JS, Lampl I, Reichova I, Carandini M, Ferster DStimulus dependence of two-state fluctuations of membrane potential in cat visual cortex.Nat Neurosci32000617621
Crossref | PubMed | ISI | Google Scholar
6 Bair W, Zohary E, Newsome WTCorrelated firing in macaque visual area MT: time scales and relationship to behavior.J Neurosci21200116761697
Crossref | PubMed | ISI | Google Scholar
7 Brillinger DRMeasuring the association of point processes: a case history.Am Math Monthly8319761622
Crossref | ISI | Google Scholar
8 Brody CDCorrelations without synchrony.Neural Comput11199915371551
Crossref | PubMed | ISI | Google Scholar
9 Brosch M, Schreiner CECorrelations between neural discharges are related to receptive field properties in cat primary auditory cortex.Exp J Neurosci11199935173530
ISI | Google Scholar
10 Efron B, Tibshirani RJAn Introduction to the Bootstrap.1993Chapman and HallNew York
Google Scholar
11 Eggermont JJ, Smith GMNeural connectivity only accounts for a small part of neural correlation in auditory cortex.Exp Brain Res1101996379391
Crossref | PubMed | ISI | Google Scholar
12 Epping WJM, Eggermont JJCoherent neural activity in the auditory midbrain of the grassfrog.J Neurophysiol57198714641483
Link | ISI | Google Scholar
13 Everitt BSAnalysis of Contingency Tables2nd ed.1992Chapman and HallNew York
Google Scholar
14 Ghose GM, Freeman RD, Ohzawa ILocal intracortical connections in the cat's visual cortex: postnatal development and plasticity.J Neurophysiol72199412901303
Link | ISI | Google Scholar
15 Hubalek ZCoefficients of association and similarity, based on binary (presence-absence) data: an evaluation.Biol Rev571982669689
Crossref | ISI | Google Scholar
16 Ito H, Tsuji SModel dependence in quantification of spike interdependence by joint peri-stimulus time histogram.Neural Comput122000195217
Crossref | PubMed | ISI | Google Scholar
17 Kruger J, Aiple FMultimicroelectrode investigation of monkey striate cortex: spike train correlations in the infragranular layers.J Neurophysiol601988798828
Link | ISI | Google Scholar
18 Lampl I, Reichova I, Ferster DSynchronous membrane potential fluctuations in neurons of the cat visual cortex.Neuron221999361374
Crossref | PubMed | ISI | Google Scholar
19 Palm G, Aertsen AMHJ, Gerstein GLOn the significance of correlations among neuronal spike trains.Biol Cybern591988111
Crossref | PubMed | ISI | Google Scholar
20 Pearson KMathematical contributions to the theory of evolution. VII. On the correlation of characters not quantitatively measurable.Philos Trans R Soc Lond A Math Phys1951900147
Crossref | Google Scholar
21 Pearson K, Heron DOn theories of association.Biometrika91913159315
Crossref | Google Scholar
22 Perkel DH, Gerstein GL, Moore GPNeuronal spike trains and stochastic point processes. II. Simultaneous spike trains.Biophys J71967419440
Crossref | PubMed | ISI | Google Scholar
23 Ringach DLSpatial structure and symmetry of simple-cell receptive fields in macaque primary visual cortex.J Neurophysiol882002455463
Link | ISI | Google Scholar
24 Roe AW, Ts'o DYSpecificity of color connectivity between primate V1 and V2.J Neurophysiol82199927192730
Link | ISI | Google Scholar
25 Salinas E, Sejnowski TJImpact of correlated synaptic input on output firing rate and variability in simple neuronal models.J Neurosci20200061936200
Crossref | PubMed | ISI | Google Scholar
26 Tallis GMThe maximum likelihood estimation of correlation from contingency tables.Biometrics181962342353
Crossref | ISI | Google Scholar
27 van den Boogaard H, Hesselmans G, Johannesma PSystem identification based on point processes and correlation densities. I. The nonrefractory neuron model.Math Biosci801986143171
Crossref | ISI | Google Scholar
28 Wielaard DJ, Shelley M, McLaughlin D, Shapley RHow simple cells are made in a nonlinear network model of the visual cortex.J Neurosci21200152035211
Crossref | PubMed | ISI | Google Scholar
29 Yule GUOn the methods of measuring association between two attributes.J R Stat Soc751912579652
Crossref | Google Scholar

Journal of Neurophysiology 89 4 cover image

Volume 89Issue 4April 2003Pages 2271-2278

Copyright & Permissions

History

Received 4 October 2002

Accepted 5 December 2002

Published online 1 April 2003

Published in print 1 April 2003

PDF download

Metrics

Downloaded 270 times

Estimating Membrane Voltage Correlations From Extracellular Spike Trains

Abstract

INTRODUCTION

RESULTS

Tetrachoric correlation—underlying model assumptions

Testing the method on leaky integrate-and-fire neurons

Data analysis

Simulation results

DISCUSSION

Definition of the problem

Model description

Predicted relationship between intracellular and extracellular correlation

FOOTNOTES

REFERENCES

More from this issue >

Metrics