Skip to main content
Log in

Information-geometric measure of 3-neuron firing patterns characterizes scale-dependence in cortical networks

  • Published:
Journal of Computational Neuroscience Aims and scope Submit manuscript

Abstract

To understand the functional connectivity of neural networks, it is important to develop simple and incisive descriptors of multineuronal firing patterns. Analysis at the pairwise level has proven to be a powerful approach in the retina, but it may not suffice to understand complex cortical networks. Here we address the problem of describing interactions among triplets of neurons. We consider two approaches: an information-geometric measure (Amari 2001), which we call the “strain,” and the Kullback-Leibler divergence. While both approaches can be used to assess whether firing patterns differ from those predicted by a pairwise maximum-entropy model, the strain provides additional information. Specifically, when the observed firing patterns differ from those predicted by a pairwise model, the strain indicates the nature of this difference—whether there is an excess or a deficit of synchrony—while the Kullback-Leibler divergence only indicates the magnitude of the difference. We show that the strain has technical advantages, including ease of calculation of confidence bounds and bias, and robustness to the kinds of spike-sorting errors associated with tetrode recordings. We demonstrate the biological importance of these points via an analysis of multineuronal firing patterns in primary visual cortex. There is a striking scale-dependent behavior of triplet firing patterns: deviations from the pairwise model are substantial when the neurons are within 300 microns of each other, and negligible when they are at a distance of >600 microns. The strain identifies a consistent pattern to these interactions: when triplet interactions are present, the strain is nearly always negative, indicating that there is less synchrony than would be expected from the pairwise interactions alone.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4

Similar content being viewed by others

References

  • Abeles, M. (1991). Corticonics: neural circuits of the cerebral cortex. Cambridge: Cambridge U P.

    Google Scholar 

  • Amari, S. I. (2001). Information geometry on hierarchy of probability distributions. IEEE Transactions on Information Theory, 47, 1701–1711.

    Article  Google Scholar 

  • Arabzadeh, E., Petersen, R. S., & Diamond, M. E. (2003). Encoding of whisker vibration by rat barrel cortex neurons: implications for texture discrimination. The Journal of Neuroscience, 23(27), 9146–9154.

    CAS  PubMed  Google Scholar 

  • Bonifazi, P., Goldin, M., Picardo, M. A., Jorquera, I., Cattani, A., Bianconi, G., et al. (2009). GABAergic hub neurons orchestrate synchrony in developing hippocampal networks. Science, 326(5958), 1419–1424.

    Article  CAS  PubMed  Google Scholar 

  • Buzsaki, G. (2004). Large-scale recording of neuronal ensembles. Nature Neuroscience, 7(5), 446–451.

    Article  CAS  PubMed  Google Scholar 

  • Fee, M. S., Mitra, P. P., & Kleinfeld, D. (1996). Automatic sorting of multiple unit neuronal signals in the presence of anisotropic and non-Gaussian variability. Journal of Neuroscience Methods, 69(2), 175–188.

    Article  CAS  PubMed  Google Scholar 

  • Feinerman, O., Rotem, A., & Moses, E. (2008). Reliable neuronal logic devices from patterned hippocampal cultures. Nature Physics, 4(12), 967–973.

    Article  CAS  Google Scholar 

  • Franke, F., Natora, M., Boucsein, C., Munk, M. H., & Obermayer, K. (2009). An online spike detection and spike classification algorithm capable of instantaneous resolution of overlapping spikes. Journal of Computational Neuroscience.

  • Gutig, R., Aertsen, A., & Rotter, S. (2003). Analysis of higher-order neuronal interactions based on conditional inference. Biological Cybernetics, 88(5), 352–359.

    Article  CAS  PubMed  Google Scholar 

  • Kennel, M. B., Shlens, J., Abarbanel, H. D., & Chichilnisky, E. J. (2005). Estimating entropy rates with Bayesian confidence intervals. Neural Computation, 17(7), 1531–1576.

    Article  PubMed  Google Scholar 

  • Kohn, A., & Smith, M. A. (2005). Stimulus dependence of neuronal correlation in primary visual cortex of the macaque. The Journal of Neuroscience, 25(14), 3661–3673.

    Article  CAS  PubMed  Google Scholar 

  • Li, C. Y., Poo, M. M., & Dan, Y. (2009). Burst spiking of a single cortical neuron modifies global brain state. Science, 324(5927), 643–646.

    Article  CAS  PubMed  Google Scholar 

  • Marre, O., El Boustani, S., Fregnac, Y., & Destexhe, A. (2009). Prediction of spatiotemporal patterns of neural activity from pairwise correlations. Physical Review Letters, 102(13), 138101.

    Article  CAS  PubMed  Google Scholar 

  • Martignon, L., Von Hasseln, H., Grun, S., Aertsen, A., & Palm, G. (1995). Detecting higher-order interactions among the spiking events in a group of neurons. Biological Cybernetics, 73(1), 69–81.

    Article  CAS  PubMed  Google Scholar 

  • Martignon, L., Deco, G., Laskey, K., Diamond, M., Freiwald, W., & Vaadia, E. (2000). Neural coding: higher-order temporal patterns in the neurostatistics of cell assemblies. Neural Computation, 12(11), 2621–2653.

    Article  CAS  PubMed  Google Scholar 

  • Meister, M., Lagnado, L., & Baylor, D. A. (1995). Concerted signaling by retinal ganglion cells. Science, 270(5239), 1207–1210.

    Article  CAS  PubMed  Google Scholar 

  • Montani, F., Ince, R. A., Senatore, R., Arabzadeh, E., Diamond, M. E., & Panzeri, S. (2009). The impact of high-order interactions on the rate of synchronous discharge and information transmission in somatosensory cortex. Philosophical Transactions A Mathematical, Physical and Engineering Sciences (1901), 367, 3297–3310.

    Article  Google Scholar 

  • Nakahara, H., & Amari, S. (2002). Information-geometric measure for neural spikes. Neural Computation, 14(10), 2269–2316.

    Article  PubMed  Google Scholar 

  • Ohiorhenuan, I. E., Mechler, F., Purpura, K. P., Schmid, A. M., Hu, Q., & Victor, J.D. (2010). Sparse coding and high-order correlations in fine-scale cortical networks. Nature. doi:10.1038/nature09178.

  • Palm, G. (1981). Evidence, information, and surprise. Biological Cybernetics, 42(1), 57–68.

    Article  CAS  PubMed  Google Scholar 

  • Palm, G., Aertsen, A. M., & Gerstein, G. L. (1988). On the significance of correlations among neuronal spike trains. Biological Cybernetics, 59(1), 1–11.

    Article  CAS  PubMed  Google Scholar 

  • Pastalkova, E., Itskov, V., Amarasingham, A., & Buzsaki, G. (2008). Internally generated cell assembly sequences in the rat hippocampus. Science, 321(5894), 1322–1327.

    Article  CAS  PubMed  Google Scholar 

  • Prut, Y., Vaadia, E., Bergman, H., Haalman, I., Slovin, H., & Abeles, M. (1998). Spatiotemporal structure of cortical activity: properties and behavioral relevance. Journal of Neurophysiology, 79(6), 2857–2874.

    CAS  PubMed  Google Scholar 

  • Reid, R. C., Victor, J. D., & Shapley, R. M. (1997). The use of m-sequences in the analysis of visual neurons: linear receptive field properties. Visual Neuroscience, 14(6), 1015–1027.

    Article  CAS  PubMed  Google Scholar 

  • Ringach, D. L., Shapley, R. M., & Hawken, M. J. (2002). Orientation selectivity in macaque V1: diversity and laminar dependence. The Journal of Neuroscience, 22(13), 5639–5651.

    CAS  PubMed  Google Scholar 

  • Roudi, Y., Nirenberg, S., & Latham, P. E. (2009). Pairwise maximum entropy models for studying large biological systems: when they can work and when they can‘t. PLoS Computational Biology, 5(5), e1000380.

    Article  PubMed  Google Scholar 

  • Roudi, Y., Tyrcha, J., & Hertz, J. (2009). Ising model for neural data: model quality and approximate methods for extracting functional connectivity. Physical Review. E, Statistical, Nonlinear, snd Soft Matter Physics, 79(5 Pt 1), 051915.

    Article  Google Scholar 

  • Rulkov, N. F., & Bazhenov, M. (2008). Oscillations and synchrony in large-scale cortical network models. Journal of Biological Physics, 34(3–4), 279–299.

    Article  PubMed  Google Scholar 

  • Schmid, A. M., Purpura, K. P., Ohiorhenuan, I. E., Mechler, F. & Victor, J. D. (2009). Subpopulations of neurons in visual area V2 perform differentiation and integration operations in space and time. Frontiers in Systems Neuroscience, 3, 15. doi:10.3389/neuro.06.015.2009.

    Article  PubMed  Google Scholar 

  • Schneidman, E., Berry, M. J., 2nd, Segev, R., & Bialek, W. (2006). Weak pairwise correlations imply strongly correlated network states in a neural population. Nature, 440(7087), 1007–1012.

    Article  CAS  PubMed  Google Scholar 

  • Shlens, J., Field, G. D., Gauthier, J. L., Grivich, M. I., Petrusca, D., Sher, A., et al. (2006). The structure of multi-neuron firing patterns in primate retina. The Journal of Neuroscience, 26(32), 8254–8266.

    Article  CAS  PubMed  Google Scholar 

  • Shlens, J., Field, G. D., Gauthier, J. L., Greschner, M., Sher, A., Litke, A. M., et al. (2009). The structure of large-scale synchronized firing in primate retina. The Journal of Neuroscience, 29(15), 5022–5031.

    Article  CAS  PubMed  Google Scholar 

  • Tang, A., Jackson, D., Hobbs, J., Chen, W., Smith, J. L., Patel, H., et al. (2008). A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. The Journal of Neuroscience, 28(2), 505–518.

    Article  CAS  PubMed  Google Scholar 

  • Vaadia, E., Bergman, H., & Abeles, M. (1989). Neuronal activities related to higher brain functions–theoretical and experimental implications. IEEE Transactions on Biomedical Engineering, 36(1), 25–35.

    Article  CAS  PubMed  Google Scholar 

  • Victor, J. D., Mechler, F., Repucci, M. A., Purpura, K. P., & Sharpee, T. (2006). Responses of V1 neurons to two-dimensional Hermite functions. Journal of Neurophysiology, 95(1), 379–400.

    Article  PubMed  Google Scholar 

  • Vogels, T. P., & Abbott, L. F. (2005). Signal propagation and logic gating in networks of integrate-and-fire neurons. The Journal of Neuroscience, 25(46), 10786–10795.

    Article  CAS  PubMed  Google Scholar 

  • Yu, S., Huang, D., Singer, W., & Nikolic, D. (2008). A small world of neuronal synchrony. Cerebral Cortex, 18(12), 2891–2901.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

This work was supported by NIH F31-EY019454 and NIH GM7739 (IEO), and NIH R01-EY9314 (JDV). The portion of this work related to the Kullback-Leibler divergence has been presented at the 2007 meetings of CoSyNe and the Society for Neuroscience.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Ifije E. Ohiorhenuan.

Additional information

Action Editor: David Golomb

Electronic Supplementary Materials

Below is the link to the electronic supplementary material.

Figure S1

(PPT 171 kb)

Appendices

Appendix

This Appendix expands on a number of technical points mentioned in the main text. First, we demonstrate the basic formula for the strain, eq. (9); this is well-known ((Amari 2001) (Nakahara and Amari 2002)) and we demonstrate it here for the reader’s convenience and to illustrate the notation used in the remainder of the Appendix. We then develop asymptotic expressions for the bias and variance of the “plugin” estimator for the strain (eqs. (11) and (12) of the main text), and for the effect of a lockout correction. Finally, we comment on ways that the strain can be applied to recordings from more than three neurons.

Our notational setup is the same as in the main text: we consider a network of M neurons and population state vectors \( \vec{\sigma } = \left( {{\sigma_1},{\sigma_2}, \ldots, {\sigma_M}} \right) \), with the convention that σ k  = +1 means that neuron k is firing, and σ k  = −1 means that it is not. We also use \( s\left( {\vec{\sigma }} \right) \) to represent the parity of \( \overrightarrow \sigma \),

$$ s\left( {\vec{\sigma }} \right) = \prod\limits_j {{\sigma_j}} . $$
(A1)

\( s\left( {\overrightarrow \sigma } \right) \) is +1 if all neurons fire, or if an even number of neurons do not fire, and −1 otherwise.

For a three-neuron network (M = 3), the firing probabilities can always be expressed in the form

$$ p\left( {\overrightarrow \sigma } \right) = \frac{{Z\left( {\overrightarrow \sigma } \right)}}{Z}, $$
(A2)

where

$$ Z\left( {\overrightarrow \sigma } \right) = \exp \left( {\sum\limits_{j = 1}^M {{\alpha_j}{\sigma_j} + \sum\limits_{j < k}^M {{\beta_{jk}}{\sigma_j}{\sigma_k}} } + \gamma {\sigma_1}{\sigma_2}{\sigma_3}} \right) $$
(A3)

and

$$ Z = \sum\limits_{\vec{\sigma }} {Z\left( {\overrightarrow \sigma } \right)} . $$
(A4)

See, for example, (Amari 2001; Nakahara and Amari 2002). For notational simplicity, we focus on the three-neuron case, but the derivations immediately extend to M > 3 neurons, as we indicate.

Calculating the strain from the firing probabilities

Our first objective is to derive eq. (9) of the main text. With the above notation, it can be written

$$ \gamma = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log p\left( {\overrightarrow \sigma } \right)} . $$
(A5)

To demonstrate eq. (A5), we proceed as follows. First, we note that exactly half of the \( \overrightarrow \sigma \)’s have even parities (\( s\left( {\overrightarrow \sigma } \right) = + 1 \)), and half have odd parities (\( s\left( {\vec{\sigma }} \right) = - 1 \)). Thus, \( \sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)} = 0 \), and therefore (using eq. (A2))

$$ \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log p\left( {\overrightarrow \sigma } \right)} = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\left( {\log p\left( {\overrightarrow \sigma } \right) + \log Z} \right) = } \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log Z\left( {\overrightarrow \sigma } \right)} . $$
(A6)

We now use eq. (A3) for \( Z\left( {\overrightarrow \sigma } \right) \):

$$ \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log p\left( {\overrightarrow \sigma } \right)} = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {{\sigma_1}{\sigma_2}{\sigma_3}} \left( {\sum\limits_{j = 1}^M {{\alpha_j}{\sigma_j} + \sum\limits_{j < k}^M {{\beta_{jk}}{\sigma_j}{\sigma_k}} } + \gamma {\sigma_1}{\sigma_2}{\sigma_3}} \right). $$
(A7)

Expanding the right-hand-side leads to a simplification, because \( {\sigma_k} = \pm 1 \), and therefore \( \sigma_k^2 = 1 \):

$$ \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log p\left( {\overrightarrow \sigma } \right)} = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {\left( {{\alpha_1}{\sigma_2}{\sigma_3} + {\alpha_2}{\sigma_1}{\sigma_3} + {\alpha_3}{\sigma_2}{\sigma_3} + {\beta_{12}}{\sigma_3} + {\beta_{13}}{\sigma_2} + {\beta_{23}}{\sigma_1} + \gamma } \right)} . $$
(A8)

In calculating the sum over all 2M network states \( \overrightarrow \sigma \), the α-terms vanish because each of the alternatives \( {\sigma_j}{\sigma_k} = \pm 1 \) occurs half of time, and the β-terms vanish because each of the alternatives \( {\sigma_k} = \pm 1 \) occur half the time. Thus,

$$ \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\vec{\sigma }} \right)\log p\left( {\vec{\sigma }} \right)} = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} \gamma, $$
(A9)

which is equivalent to eq. (A5) and eq. (9) of the main text.

Note that the above derivation immediately generalizes to the highest coefficient of a complete exponential expansion, e.g., the coefficient of σ 1 σ 2 σ 3 σ 4 for M = 4, etc.

Asymptotic estimates of bias and variance

In this section, we determine the asymptotic bias and variance of estimates of γ, based on the “plug-in” estimator corresponding to eq. (A5).

To formalize the problem, we assume that our data consist of N observations (“bins”), each of which are independent, and (as in the main text) that the number of observations that are in the state \( \overrightarrow \sigma \) is \( N\left( {\overrightarrow \sigma } \right) \), so \( N = \sum {N\left( {\overrightarrow \sigma } \right)} \). Our estimate of \( p\left( {\overrightarrow \sigma } \right) \) is the naïve one, namely

$$ {p_{est}}\left( {\vec{\sigma }} \right) = \frac{{N\left( {\vec{\sigma }} \right)}}{N} $$
(A10)

With this estimate for \( p\left( {\vec{\sigma }} \right) \), the “plug-in” estimator corresponding to eq. (A5) is

$$ {\gamma_{est}} = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log \frac{{N\left( {\overrightarrow \sigma } \right)}}{N}} . $$
(A11)

The bias of the estimator γ est is, by definition, the expected departure of the estimator from its true value γ:

$$ B \equiv \left\langle {{\gamma_{est}}} \right\rangle - \gamma, $$
(A12)

where <> denotes an average over all experiments with N observations, and true state probabilities \( p\left( {\overrightarrow \sigma } \right) \). Its variance is, by definition,

$$ V \equiv \left\langle {{{\left( {{\gamma_{est}} - \left\langle {{\gamma_{est}}} \right\rangle } \right)}^2}} \right\rangle, $$
(A13)

Thus, our goal is to determine the mean and variance of γ est , as defined by eq. (A11).

We begin by noting that eq. (A11) simplifies, because half of the \( \overrightarrow \sigma \)’s have even parities, and half have odd parities. Thus,

$$ {\gamma_{est}} = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log N\left( {\overrightarrow \sigma } \right)} . $$
(A14)

The right-hand-side of eq. (A14) is a linear combination of logarithms of values drawn from a multinomial process of 2M symbols (the states \( \overrightarrow \sigma \)), each with probability \( p\left( {\overrightarrow \sigma } \right) \). Therefore, to determine the mean and variance of γ est , we determine the mean and variance of the logarithms of values drawn from a multinomial distribution.

We determine these quantities from the well-known expressions for the mean, variance, and covariance of a multinomial process. In N observations, the expected number of times that a symbol with probability \( p\left( {\overrightarrow \sigma } \right) \) is encountered is

$$ \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle = Np\left( {\overrightarrow \sigma } \right). $$
(A15)

With \( \delta \left( {\overrightarrow \sigma } \right) = N\left( {\overrightarrow \sigma } \right) - \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle \), the variance of each count is given by

$$ \left\langle {\delta {{\left( {\overrightarrow \sigma } \right)}^2}} \right\rangle = Np\left( {\overrightarrow \sigma } \right)\left( {1 - p\left( {\overrightarrow \sigma } \right)} \right) $$
(A16)

and the covariances are given by

$$ \left\langle {\delta \left( {\overrightarrow \sigma } \right)\delta \left( {\overrightarrow \sigma } \prime \right) } \right\rangle = - Np\left( {\overrightarrow \sigma } \right)p\left( {\overrightarrow \sigma \prime } \right),\overrightarrow \sigma \ne \overrightarrow \sigma \prime . $$
(A17)

We now calculate the corresponding statistics of the logarithms asymptotically:

$$ \left\langle {\log N\left( {\overrightarrow \sigma } \right)} \right\rangle = \left\langle {\log \left( {\left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle + \delta \left( {\overrightarrow \sigma } \right)} \right)} \right\rangle = \log \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle + \left\langle {\log \left( {1 + \frac{{\delta \left( {\overrightarrow \sigma } \right)}}{{\left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle }}} \right)} \right\rangle $$

Expanding the logarithm via a Taylor series yields

$$ \left\langle {\log N\left( {\overrightarrow \sigma } \right)} \right\rangle = \log \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle + \left\langle {\frac{{\delta \left( {\overrightarrow \sigma } \right)}}{{\left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle }}} \right\rangle - \frac{1}{2}\left\langle {{{\left( {\frac{{\delta \left( {\overrightarrow \sigma } \right)}}{{\left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle }}} \right)}^2}} \right\rangle + O\left( {N{{\left( {\overrightarrow \sigma } \right)}^{ - 3/2}}} \right), $$
(A18)

or

$$ \left\langle {\log N\left( {\overrightarrow \sigma } \right)} \right\rangle = \log \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle - \frac{{Np\left( {\overrightarrow \sigma } \right)\left( {1 - p\left( {\overrightarrow \sigma } \right)} \right)}}{{2{{\left( {Np\left( {\overrightarrow \sigma } \right)} \right)}^2}}} + O\left( {{N^{ - 3/2}}} \right) = \log Np\left( {\overrightarrow \sigma } \right) - \frac{{1 - p\left( {\overrightarrow \sigma } \right)}}{{2Np\left( {\overrightarrow \sigma } \right)}} + O\left( {{N^{ - 3/2}}} \right), $$
(A19)

where we have used \( \left\langle {\delta \left( {\overrightarrow \sigma } \right)} \right\rangle = 0 \), eq. (A15), and eq. (A16).

Thus, the bias in the logarithm of the counts of the state \( \overrightarrow \sigma \) is

$$ B\left( {\overrightarrow \sigma } \right) \equiv \left\langle {\log N\left( {\overrightarrow \sigma } \right)} \right\rangle - \log \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle = - \frac{{1 - p\left( {\overrightarrow \sigma } \right)}}{{2Np\left( {\overrightarrow \sigma } \right)}} + O\left( {{N^{ - 3/2}}} \right). $$
(A20)

The bias estimate (text eq. (11)) now follows from the definition (A12) by combining eq. (A20) with eq. (A14):

$$ B \equiv \left\langle {{\gamma_{est}}} \right\rangle - \gamma = \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)B\left( {\overrightarrow \sigma } \right)} = - \frac{1}{{{2^{M + 1}}N}}\sum\limits_{\vec{\sigma }} {\frac{{s\left( {\overrightarrow \sigma } \right)}}{{p\left( {\overrightarrow \sigma } \right)}}} + O\left( {{N^{ - 3/2}}} \right) $$
(A21)

with a simplification resulting from the fact that the parities sum to zero (\( \sum {s(\vec{\sigma })} = 0 \)).

Similar algebra provides asymptotic expressions for the variances of the logarithms of the counts,

$$ \begin{array}{*{20}{c}} {V\left( {\overrightarrow \sigma ,\overrightarrow \sigma } \right) \equiv \left\langle {{{\left( {\log N\left( {\overrightarrow \sigma } \right) - \left\langle {\log N\left( {\overrightarrow \sigma } \right)} \right\rangle } \right)}^2}} \right\rangle = } \hfill \\ {\left\langle {{{\left( {\left( {\log \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle + \log \left( {1 + \frac{{\delta \left( {\overrightarrow \sigma } \right)}}{{\left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle }}} \right)} \right) - \left( {\log \left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle - \frac{{1 - p\left( {\overrightarrow \sigma } \right)}}{{2Np\left( {\overrightarrow \sigma } \right)}}} \right)} \right)}^2}} \right\rangle + O\left( {{N^{ - 3/2}}} \right) = } \hfill \\ {\left\langle {{{\left( {\log \left( {1 + \frac{{\delta \left( {\overrightarrow \sigma } \right)}}{{\left\langle {N\left( {\overrightarrow \sigma } \right)} \right\rangle }}} \right) - \frac{{1 - p\left( {\overrightarrow \sigma } \right)}}{{2Np\left( {\overrightarrow \sigma } \right)}}} \right)}^2}} \right\rangle + O\left( {{N^{ - 3/2}}} \right) = } \hfill \\ {\frac{1}{{Np\left( \sigma \right)}} - \frac{1}{N} + O\left( {{N^{ - 3/2}}} \right)} \hfill \\ \end{array} $$
(A22)

and their covariances,

$$ V\left( {\overrightarrow \sigma, \overrightarrow \sigma \prime } \right) \equiv \left\langle {\left( {\log N\left( {\overrightarrow \sigma } \right) - \left\langle {\log N\left( {\overrightarrow \sigma } \right)} \right\rangle } \right)\left( {\log N\left( {\overrightarrow \sigma \prime } \right) - \left\langle {\log N\left( {\overrightarrow \sigma \prime } \right)} \right\rangle } \right)} \right\rangle = - \frac{1}{N} + O\left( {{N^{ - 3/2}}} \right),\overrightarrow \sigma \ne \overrightarrow \sigma \prime . $$
(A23)

The variance estimate (text eq. (12)) follows from the definition (A13) by combining eqs. (A22) and (A23) with eq. (A14):

$$ V \equiv \left\langle {{{\left( {{\gamma _{est}} - \left\langle {{\gamma _{est}}} \right\rangle } \right)}^2}} \right\rangle = \frac{1}{{{2^{2M}}}}\sum\limits_{\overrightarrow \sigma ,\overrightarrow \sigma \prime } {s\left( {\overrightarrow \sigma } \right)s\left( {\overrightarrow \sigma } \right)V\left( {\overrightarrow \sigma ,\overrightarrow \sigma \prime } \right)} = \frac{1}{{{2^{2M}}N}}\sum\limits_{\overrightarrow \sigma } {\frac{1}{{p\left( {\overrightarrow \sigma } \right)}}} + O\left( {{N^{ - 3/2}}} \right). $$
(A24)

In practice, γ est is approximately Gaussian-distributed, and the above asymptotic expressions for bias (eq. (A21)) and variance (eq. (A24)) are accurate provided that the minimum number of counts, \( N{\min_{\overrightarrow \sigma }}p\left( {\overrightarrow \sigma } \right) \), is at least ~10. Simulations that support this approximation are shown in Fig. 5.

Fig. 5
figure 5

The asymptotic bias and variance estimates of the plug-in estimator for the strain are highly accurate, provided that there are at least 10 counts of each firing pattern. In the upper portion of the Figure, the asymptotic Gaussian distribution (smooth curve) is compared to the empirical distribution of the plug-in estimator, eq. (A11). The Gaussian has a mean given by the asymptotic debiaser (A21) and variance corresponding to the asymptotic estimate (A24). The histograms show the empirical distribution of the plug-in estimator (A11) for 1000 simulated experiments. In each set of simulations, spike pattern counts are distributed in a Poisson fashion according to the probabilities given by eq. (7), with the parameters for the right-hand-side indicated at the top of each column. From top to bottom, the number of bins in the experiment increases from 102 to 105; the latter value is comparable to the number of bins used in the experiments reported here. Note that when the number of bins is small (upper rows), some of the spike patterns are not encountered, so the plug-in estimator diverges; these counts are not included in the histograms. The numbers inset into each plot indicate the expected number of counts of the least likely firing pattern. Note that when there are at least 10 counts of each firing pattern, the asymptotic Gaussian provides an excellent approximation to the empirical distribution. The lower portion of the Figure shows the mean (points) and standard deviation (error bars) of the empirical distribution of the plug-in estimator, as a function of the number of bins in the simulated experiment. Points are only plotted when all of the 1000 simulated experiments had at least one count of all patterns. The horizontal line indicates the true value of the strain. The three smooth curves indicate the asymptotic estimates of the mean and standard deviation, according to eqs. (A21) and (A24). The two rows show the same data plotted on different vertical scales, so that the biases can be clearly seen

The effect of the lockout correction

As mentioned in the main text, spike-sorting introduces a potential confound into the analysis of multineuronal data. The source of the problem is that when two spikes overlap in time, their waveforms sum, and neither can be reliably identified. When this happens, a bin that actually contains two firing events will be registered as containing none; this is known as a “lockout.” (We note that recently, spike-sorting algorithms that may be able to resolve these overlaps have been developed (Franke et al. 2009); the applicability of such methods is likely to depend on the details of the recording situation, such as firing rates and noise. Thus, even with the availability of advanced spike-sorting algorithms, it is desirable to have analytic approaches that are insensitive to spike-sorting errors.)

As described in Methods, an approximate correction for lockout can be constructed, by making use of the fact that in a typical recording, the analysis bin (the interval of time used for tallying spiking states) is many times larger than the spike width. Figure 1 shows that the strain is nearly independent of the lockout correction, indicating that it is relatively insensitive to the effects of the lockout. Our goal here is to show why this is the case.

We do this by first describing the lockout correction strategy in detail, and then deriving a formula that relates the lockout-corrected value of the strain to the observed probability of each spiking state.

The lockout correction is based on an estimate of how many events are miscounted due to overlapping spikes. To obtain this estimate, we assume that firing rates are constant for the duration of the analysis bin, that spikes occur in a Poisson fashion, and that no bin contains more than one spike from each neuron. We then partition each analysis bin into W sub-intervals, where each sub-interval corresponds to the spike width. (In the data presented here, W = 8 for the datasets analyzed with 10 ms bins and W = 12 for the datasets analyzed with 14.8 ms bins.) We estimate the miscounts by relating the total frequency of multineuronal events to the observed frequency of the multineuronal events. In order for a multineuronal event to be observed, each sub-interval can contain no more than one spike.

Consider, for example, an event in which neurons 1 and 2 fire. The spikes from neurons 1 and 2 can fall into the W sub-intervals in W 2 equally-likely ways, since each of the W sub-intervals is equally likely to contain either spike. Usually, the spikes from the two neurons fall into different sub-intervals, so both are registered; these events properly contribute to the estimate of p 110. This happens in W(W-1) of the W 2 possible configurations in which the spikes can fall into sub-intervals. But in the remaining W configurations, both spikes fall into the same sub-interval, so neither spike is registered, and the event is an erroneous contribution to p 000. As a result of these overlap events, p 110 is under-estimated by a factor of \( \frac{{{W^2}}}{{W\left( {W - 1} \right)}} \approx 1 + 1/W \); the resulting undercounts contribute to p 000. (Here, we write p 110 as the probability of a state in which σ 1 = +1, σ 2 = +1, and σ 3 = −1, etc.).

Similarly, triplet firing events can populate the W sub-intervals in W 3 ways, but only \( W\left( {W - 1} \right)\left( {W - 2} \right) \) of these are observed configurations, in which no sub-interval has 2 or more spikes. Thus, p 111 is undercounted by a factor of \( \frac{{{W^3}}}{{W\left( {W - 1} \right)\left( {W - 2} \right)}} \approx 1 + 3/W \). Most of the undercounts occur when two spikes occur in one sub-interval, and a third occurs in another sub-interval. The two spikes that occur in the same sub-interval occlude each other. The third spike, which occurs in its own sub-interval, is correctly detected. Thus, the three-spike event appears to contain only one spike, and contributes to one of the singleton event probabilities p 100, p 010, and p 001. A small number of undercounts are events in which all three spikes occur within the same sub-interval. The three spikes occlude each other and thus the event contributes to p 000. We neglect these events, since their proportion is O(W −2).

Putting the above facts together leads to formulas that relate the raw estimates of the state probabilities to those that would have been encountered in the absence of a lockout. We add the undercounts to the triplet firing events

$$ {p_{111}}^{corrected} = {p_{111}}^{raw}\left( {1 + \frac{3}{W}} \right), $$
(A25)

and to the doublet firing events

$$ {p_{110}}^{corrected} = {p_{110}}^{raw}\left( {1 + \frac{1}{W}} \right),{p_{101}}^{corrected} = {p_{101}}^{raw}\left( {1 + \frac{1}{W}} \right)\,and\,{p_{011}}^{corrected} = {p_{011}}^{raw}\left( {1 + \frac{1}{W}} \right). $$
(A26)

The counts that are added to the triplet firing events (eq. (A25)) are removed from the singlet firing events

$$ {p_{100}}^{corrected} = {p_{100}}^{raw} - \frac{{{p_{111}}^{raw}}}{W},{p_{010}}^{corrected} = {p_{010}}^{raw} - \frac{{{p_{111}}^{raw}}}{W},\,and\,{p_{001}}^{corrected} = {p_{001}}^{raw} - \frac{{{p_{111}}^{raw}}}{W}, $$
(A27)

and the counts that are added to the doublet firing events (eq. (A26)) are removed from the null events

$$ {p_{000}}^{corrected} = {p_{000}}^{raw} - \frac{{{p_{011}}^{raw} + {p_{101}}^{raw} + {p_{110}}^{raw}}}{W}. $$
(A28)

Thus, (as mentioned in the main text), we anticipate that the lockout correction will be small: it shifts events within the numerator (eqs. (A25) and (A27)) or within the denominator (eqs. (A26) and (A28)), but not between numerator and denominator.

To formalize this observation, we substitute the above formulae (eqs. (A25) to (A28)) into

$$ {\gamma^{corrected}} = \frac{1}{8}\log \frac{{{p_{111}}^{corrected}{p_{100}}^{corrected}{p_{010}}^{corrected}{p_{001}}^{corrected}}}{{{p_{011}}^{corrected}{p_{101}}^{corrected}{p_{011}}^{corrected}{p_{000}}^{corrected}}} $$
(A29)

to obtain (after some algebra)

$$ {\gamma^{corrected}} = {\gamma^{raw}} + \frac{1}{{8W}}\left( {\left( {\frac{{{p_{011}}^{raw}}}{{{p_{000}}^{raw}}} - \frac{{{p_{111}}^{raw}}}{{{p_{100}}^{raw}}}} \right) + \left( {\frac{{{p_{101}}^{raw}}}{{{p_{000}}^{raw}}} - \frac{{{p_{111}}^{raw}}}{{{p_{010}}^{raw}}}} \right) + \left( {\frac{{{p_{110}}^{raw}}}{{{p_{000}}^{raw}}} - \frac{{{p_{111}}^{raw}}}{{{p_{001}}^{raw}}}} \right)} \right) + O\left( {{W^{ - 2}}} \right). $$
(A30)

Equation (A30) provides further reasons that the lockout correction is small: γ corrected-γ raw is a sum of three terms, each of which is a difference of probability ratios. These differences are nonzero only to the extent to which triplet interactions are present. For example, the first term is nonzero only if the probability of joint firing neuron 2 and 3 depends on whether neuron 1 fires. Moreover, each of these ratios has a quadratic dependence on the analysis bin width, and their total contribution is attenuated by the multiplier 1/W.

More than 3 neurons

In the final section, we mention ways in which the strain can be extended to analyses of M > 3 neurons. Above we have taken the simple approach of analyzing larger datasets by considering only three neurons at a time, but other approaches are possible. Perhaps the most natural alternative is to enlarge the pairwise maximum-entropy model (text eq.(3)) to include all third-order terms:

$$ {p_{model}}\left( {\overrightarrow \sigma; {{\mathbf{\alpha }}^{model}},{{\mathbf{\beta }}^{model}},{{\mathbf{\gamma }}^{model}}} \right) = \frac{1}{Z}\exp \left( {\sum\limits_{j = 1}^M {\alpha_j^{model}{\sigma_j} + \sum\limits_{j < k}^M {\beta_{jk}^{model}{\sigma_j}{\sigma_k}} + \sum\limits_{j < k < l}^M {\gamma_{jkl}^{model}{\sigma_j}{\sigma_k}{\sigma_l}} } } \right). $$
(A31)

The resulting model is again maximum-entropy model, but now subject to the constraints of matching all triplet firing probabilities p(j, k, l) to those observed in the data, as well as matching the singlet firing probabilities p(j) and the pairwise firing probabilities p(j, k). The \( M\left( {M - 1} \right)\left( {M - 2} \right)/6 \) third-order parameters γ jkl are a comprehensive description of how triplet firing patterns depart from the predictions of the pairwise model. Note that for M > 3, there is no explicit formula comparable to eq. (9) for calculating the parameters γ jkl from the state probabilities \( p\left( {\overrightarrow \sigma } \right) \); they must be determined via fitting the third-order model.

In order to preserve the feature that γ can be calculated directly from the state probabilities \( p\left( {\overrightarrow \sigma } \right) \), a different generalization is necessary. This generalization is based on including terms up to Mth order in a maximum-entropy model. (At this point it is no longer a “model”; it is guaranteed to fit the set of state probabilities \( p\left( {\overrightarrow \sigma } \right) \) exactly (Amari 2001)). Once terms of all orders are included, the highest-order term (i.e., the coefficient of σ 1 σ 2 σ 3...σ M ) can be calculated explicitly, as

$$ {\gamma_{\left[ M \right]}} \equiv \frac{1}{{{2^M}}}\sum\limits_{\vec{\sigma }} {s\left( {\overrightarrow \sigma } \right)\log p\left( {\overrightarrow \sigma } \right)} . $$
(A32)

γ[M] indicates how the network’s firing patterns differ from those that could be predicted from the interactions up to order M-1. Like γ = γ[3], its sign indicates whether there is an excess or a deficit of synchronous events. As is the case for γ, estimation of γ[M] requires sufficient data so that even the rarest event occurs often enough (e.g., 50 times) for its probability can be accurately estimated. This is likely to be on the edge of practicality for M = 4, and impractical for larger M. We do not advocate estimating γ[M] unless all firing patterns have been observed many times, and thus, do not address the situation in which the plug-in estimator corresponding to eq. (A32) diverges.

Rights and permissions

Reprints and permissions

About this article

Cite this article

Ohiorhenuan, I.E., Victor, J.D. Information-geometric measure of 3-neuron firing patterns characterizes scale-dependence in cortical networks. J Comput Neurosci 30, 125–141 (2011). https://doi.org/10.1007/s10827-010-0257-0

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s10827-010-0257-0

Keywords

Navigation