New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods

Communicated by Roderick MacKinnon, The Rockefeller University, New York, NY, June 18, 2009 (received for review February 15, 2009)
Abstract
Complexity of neural systems often makes impracticable explicit measurements of all interactions between their constituents. Inverse statistical physics approaches, which infer effective couplings between neurons from their spiking activity, have been so far hindered by their computational complexity. Here, we present 2 complementary, computationally efficient inverse algorithms based on the Ising and “leaky integrateandfire” models. We apply those algorithms to reanalyze multielectrode recordings in the salamander retina in darkness and under random visual stimulus. We find strong positive couplings between nearby ganglion cells common to both stimuli, whereas longrange couplings appear under random stimulus only. The uncertainty on the inferred couplings due to limitations in the recordings (duration, small area covered on the retina) is discussed. Our methods will allow realtime evaluation of couplings for large assemblies of neurons.
A vertebrate retina is a structured, complex network of interacting neurons that process visual input stimuli at the photoreceptors into an output pattern of action potentials of the retinal ganglion cells (1–2). It is now a wellestablished fact that retinal cells process information in a collective fashion: The firing of one ganglion cell is correlated with the firing pattern of other cells (3–4). Multielectrode recordings have made accessible hourslong, simultaneous spiking activity of tens of retinal ganglion cells and thus have become a powerful tool to investigate the information processing performed by a vertebrate retina (5–7). The analysis of pairwise correlations in the activity has revealed different patterns of synchrony between 2 cells that have been related to different retina circuits (7).
Analyzing the concerted activity of all of the recorded cells is, however, a very challenging task. Recently Schneidman et al. (8) and Shlens et al. (9) pointed out that correlations in the firing activity of cell populations can be reconstructed from the average firing rates, f_{i}, and 2cell correlations, c_{ij}, alone. The theoretical model, which has been used to generate the frequencies of all possible 2^{N} spiking configurations for a system of N neurons, is the wellknown Ising model. It is characterized by a reduced set of ≈N^{2} parameters: N “fields,” h_{i}, experienced by individual cells, and N(N − 1)/2 “couplings,” J_{ij}, between pairs of cells. Computing the parameters h_{i} and J_{ij} from the firing patterns can be viewed as an example of the inverse statistical physics method.
The existence of a lowdimensional parameterization of the retinal activity (8–11) is an interesting and encouraging result. Naively, one may be tempted to assign to the inferred parameters a simple interpretation: The fields could represent the external stimuli, and the couplings could reflect the physiological interactions between the cells. However, because most of the neural circuitry (all cells in intermediate layers and most of the cells in the ganglion layer) is not recorded, the inferred fields and couplings are only “effective” quantities, and many questions remain unanswered. First, how do the effective couplings depend on the visual stimulus? Second, to what extent are the inferred couplings affected by the incomplete sampling of the activity, both from temporal (finite duration of the recordings) and spatial (small area of the retina covered by the electrode array) points of views? Third, do the couplings strongly depend on the model used for the inference? The main difficulty in addressing these questions in a systematic manner is connected with the computational complexity of the applied inverse statistical algorithm.
In the present work, we propose 2 efficient algorithms to calculate the effective couplings from the spiking activity of a population of neurons. The first method is based on the abovementioned Ising model, the second one on the wellknown “leaky integrateandfire” (I&F) model. We then use these algorithms to characterize the effective couplings between ganglion cells from previously published recordings of the salamander retinal activity under different visual stimuli.
Results
Ising and I&F Models.
Existing algorithms for inverse Ising problems are based on timeconsuming learning schemes (11–12). It is possible, however, to drastically improve inverse Ising methods if one takes advantage of the fact that the neurons are more often silent than active. This situation corresponds to the behavior of the Ising model at high fields. Based on a highfield expansion of the Ising thermodynamic potential at fixed f_{i} and c_{ij}, we have developed an efficient algorithm, which calculates the couplings, as well as their relative accuracy, in a time polynomial in N [Methods and supporting information (SI) Appendix, Section 1]. This means in practice that for a typical recording of the activity of N = 32 ganglion cells in the salamander retina, N(N − 1) ≈ 10^{3} couplings can be inferred (on a personal computer) in a couple of minutes rather than in many hours.
The Ising approach takes into account correlations in the “same time bin” (of the order of 20 ms) but completely ignores timedelayed correlations. We have thus developed a second approach based on the I&F model (13–15) (Methods and SI Appendix, Section 3). I&F is a dynamical model describing spike generation by a neuron in the presence of the synaptic inputs from other neurons; couplings are pairwise but are not a priori symmetric. Using techniques from statistical field theory, we have devised an efficient inference procedure for calculating the couplings, G_{ij}, from all of the S recorded spikes, in a time polynomial in both N and S. For the recording mentioned above (S ≈ 10^{5} spikes) the running algorithm time is of the order of 20 s on a personal computer.
As an example of potential application, we reanalyze 3 recordings from salamander ganglion cells: (i) a 2,000slong recording of the spontaneous activity of 32 cells in total darkness (5); (ii) a 4,450slong recording of the same 32 cells in the same retina illuminated with randomly flickering bright squares (5) (the locations of the receptive field centers of those cells are known); and (iii) a recording from 40 cells in another retina presented with a 120slong natural movie repeated 20 times (8).
Both the Ising and I&F models are approximations to the real neural activity. We have checked that multicell correlations, calculated in the Ising model within the same time bin, are rather faithfully reproduced, although slightly overestimated (SI Appendix, Section 2). The I&F model, on the other hand, despite its simplicity, is capable of reproducing some important features of the longertime correlations between spiking events (SI Appendix, Section 4).
Correlations Versus Interactions.
The task of inferring couplings from the dataset is more complex than the longestablished process of analyzing correlations in the firing activity (3–7). The coupling between 2 cells is indeed not only a function of the activity of the pair itself but depends on the activity of all of the recorded cells. It is, however, interesting to compare the couplings inferred through the Ising model and the correlation indices (Methods) previously defined (6) in the analysis of multielectrode array data (see SI Appendix, Section 4). The logarithm of the correlation index is the first contribution to the Ising coupling in our largefield expansion. This 2cell approximation to the coupling corresponds to calculating the coupling only from the average firing rates and correlation of the 2 cells in the pair (see Methods). As shown in Fig. 1, some Ising couplings J_{ij} can be accurately approximated from their 2cell approximation. In general, however, the coupling between 2 cells cannot be deduced directly from the correlation of their firing activities, because this correlation may result from indirect couplings via other neurons. By adding the consecutive expansion terms to this 2cell approximation, one can systematically generalize the correlation index analysis by including correlations within the larger clusters of cells (Methods). We have found that pairs of cells with positive and large correlation indices have, in general, large couplings with values close to their 2cell approximation, whereas pairs with positive and small correlation indices may have negative couplings (SI Appendix, figure 5). The presence of negative couplings is an interesting finding that cannot be deduced from the analysis of the data in terms of correlations (see histogram of couplings, and of their 2cell approximations in SI Appendix, figures 3 and 4).
Note that the difference between couplings and correlations indices and the size of the clusters of cells necessary for an accurate inference depend, in general, on the structure of the neural circuit.
Comparison of Couplings Obtained in Dark and Flicker Conditions.
Although the cell activities are different in flicker and dark conditions, many pairs have similar Ising couplings J_{ij} in both conditions (Fig. 2); this is the case both for strong positive couplings (such as pair “a,” whose cells often spike together within a time window Δt = 20 ms), and for several negative couplings (e.g., pair “b,” showing no correlation). Note, however, that few large and positive couplings under flicker stimulus have small or even negative values in dark (e.g., pair “‘c,” with a strong correlation in flicker and anticorrelation in dark), and few positive couplings in the dark have negative values under flicker conditions (e.g., pair “d,” whose cells are positively correlated but with a long delay >20 ms). The difference of behaviors of pairs a, b, c, and d between the 2 stimuli presumably correspond to different classes of physiological interactions, as appears from the fulltime histogram of crosscorrelations.
Comparison Between I&F and Ising Couplings.
The presence of numerous strong positive couplings is corroborated by the I&F analysis (Fig. 3). In the I&F model, G_{ij} is rarely equal to G_{ji}, but the couplings are grouped along the symmetry axis. Overall asymmetry of couplings is larger in dark than in flicker conditions (Fig. 3 and SI Appendix) The symmetrized I&F interactions (G_{ij} + G_{ji})/2 are proportional to positive Ising couplings J_{ij}, but the agreement between the Ising and I&F models is poorer for negative couplings, especially under flicker stimuli, because the Ising model ignores correlations with delays larger than Δt (e.g., pair d of Fig. 2). We find this discrepancy to be particularly important for structured stimuli, such as natural movie, which generally lead to larger delays in the spiking activity (Fig. 3).
Dependence of Couplings on Removal of Cells from the Recordings.
Multicellular recordings of neurons interrogate only a small portion of retinal ganglion cells. To study how restricting the number of measured neurons affects the inferred couplings, we “removed” a portion of the recorded cells and recalculated J_{ij}. Within a linear response approximation, the removal of a cell does not change couplings beyond 600 μm (the data are limited to distances <1.4 mm). However, couplings are often affected by removal, or inclusion, of nearby cells (SI Appendix, Sections 6 and 7). We conclude that in order for the couplings between central cells to be reliable, the electrode arrays have to be dense enough (with spacing typically <40 μm) and have to probe most ganglion cells in a surrounding region >1 mm.
Dependence of Couplings on Distance Between Cells.
Fig. 4 shows the Ising and I&F couplings between pairs of ganglion cells as a function of the distance between their receptive field centers for dark and flicker conditions. The analyzed recordings do not lead to any clear dependence of the couplings on the cell type, probably because most recorded cells (80%) are of the OFF type. We therefore group all cells together, independently of their type. The results are compatible with the conclusions of Meister et al. (6) based on the study of pairwise correlations. The dependences of the Ising and I&F couplings on distance are very similar, although, as explained in Fig. 3, the Ising negative couplings have larger amplitudes than the I&F couplings. We observe that whereas strong positive couplings are found at small distances (<500 μm) both in dark and flicker conditions, the positive couplings at larger distances (up to 1,000 μm) exist mainly in flicker but not in dark conditions and are hence stimulus induced. At the same time, negative couplings, including the few conserved interactions identified above, are acting at distances >200 μm. The absence of negative interactions at short distances in dark and in flicker conditions could come from the homogeneity in the recorded cell types.
Spatial Networks of Couplings.
Spatial information about the couplings can be depicted in form of 2D maps in the retinal plane obtained by connecting the receptive field centers with lines representing the couplings. Remarkably, the maps of the strongest couplings obtained with the I&F and the Ising inverse statistical approaches include largely the same pairs of cells (Fig. 5). The map of the strongest couplings in dark conditions (Fig. 5 A and C) has a simple connectivity: with shortrange couplings between neighboring cells. In contrast, the strong interactions present in flicker and absent in dark (e.g., pair c) are long range (Fig. 5 B and D). The appearance of these couplings seems to be at odds with recent studies in primate retina, which have found that the quasitotality of the informational entropy of spikes can be accounted for by assuming only interactions between adjacent cells (9).
Discussion
We have presented 2 computationally efficient inverse statistical approaches, thereby establishing a solid basis for the rapid, lowdimensional parameterization of neural activity. We have classified the inferred couplings with respect to their sign, dependence on the stimulus, on the inference model, and their spatial range. Although most pairwise correlations are positive, both the Ising and I&F models predict the existence of negative couplings, which can be long range and stimulus dependent. In addition to the stimulusdependent couplings, we have located strong, positive, shortrange couplings, which are similar both in the dark and under flicker stimuli.
A crucial issue connected with this lowdimensional parameterization of neural activity remains: Can the effective couplings J_{ij} be given a simple physiological interpretation (2)? For instance, it will be important to explore the physiological meaning of the longrange couplings, which appear under flicker stimuli, and could be, e.g., mediated by amacrine cells. It would be also interesting to see whether symmetric I&F couplings correspond to direct interactions between neurons (e.g., via gap junctions) (7). Additional physiological measurements could be performed to explore the underlying nature of the interneuronal couplings.
Most of the inferred couplings are expected to reflect the effective interactions through the intermediate retina layers, which are not easily accessible to measurement and the circuitry of which cannot be reconstructed (1–2). Yet, effective couplings between ganglion cells could give important insights on the way the retina encodes and transmits visual information. The computational efficiency of our inverse statistical approaches should make possible experiments in which the couplings are calculated in real time, while the firing activity is being recorded. Visual stimuli could be adjusted during experiments to explore particular aspects of retinal function. This possibility is not restricted only to multielectrode recordings in vertebrate retina: The same approach could be also used to analyze the activity of other neuronal systems, such as the recordings from slices of vertebrate cortex (12).
Methods
Algorithm for the Inverse Ising Problem.
Data are encoded into bit configurations: s_{i}^{τ} = 1 if cell i is active in bin τ, 0 otherwise, where time bins are of width Δt; τ = 1,2,…,B is in the index of the time bin in the recording of duration B × Δt. In the Ising model, the loglikelihood of a Nbit configuration (in the same bin) is where the parameters h_{i} and J_{ij} are called, respectively, effective fields and couplings; the value of F[{h_{i}},{J_{ij}}], called free energy in statistical physics, is such that the probabilities of the 2^{N} configurations sum up to 1. The effective fields and couplings must be inferred to reproduce the experimentally measured average activities, p_{i} = 〈s_{i}〉, where p_{i} = f_{i} × Δt is the probability that cell i fires in a bin, and cell–cell correlations, p_{ij} = 〈s_{i}s_{j}〉 here 〈·〉 denotes the average, with the probability P^{like}. Those equations are ill defined because of imperfect sampling, e.g., J_{ij} would be infinite if cells i and j never spike together in the set of B configurations. To avoid this problem, we introduce a Gaussian prior over the parameters and infer the effective fields and couplings through the maximization of the Bayesian logposterior probability where Γ is expected to be of the order of 1/B (SI Appendix, Section 1).
Because the calculation of the free energy F, Eq. 1, for a given set of fields and couplings requires a computational effort growing exponentially with the number N of cells, P^{post} cannot be directly maximized from the definition (Eq. 2) when N exceeds, say, 15. Fortunately, we need to know F for large fields only: The firing probability p_{i} is small because the bin width Δt (≈20 ms) is much smaller than the average interspike interval (>1 s for dark and flicker recordings). We have extended largefield expansions existing in statistical physics to: (i) tackle the case of nonuniform (celldependent) couplings and fields, and (ii) perform the minimization over the (large) fields and the couplings in Eq. 2 to express F in terms of the firing probabilities p_{i} and correlations p_{ij}.
The outcome of the procedure is that the coupling J_{ab} can be decomposed into a sum of contributions coming from all clusters of k + 2 distinct cells including cells a and b, The first term in the expansion is J_{ab}^{(2)} = \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\frac{1}{4}$$ \end{document} logCI_{ab}, where CI_{ab} = p_{ab}/(p_{a}p_{b}) is the correlation index of cells a and b (SI Appendix, Section 1) and is called here the 2cell approximation to the coupling. Terms with k ≥ 1 correspond to network contributions to J_{ab}: Informally speaking, J_{ab}^{(k+2)}[a,b,i_{1},i_{2},…,i_{k}] is the contribution to the coupling between cells a and b that cannot be obtained from the knowledge of the activities of smaller subsets of those k + 2 cells.
Our expression for J_{ab}^{(k+2)}[a,b,i_{1},i_{2},…,i_{k}] is an infinite series in powers of the connected correlation c_{ij} − p_{i}p_{j}. Although we are not able to obtain a closed analytical expression we can show the following fact. Assume that all firing probabilities are small, that is, p_{i} < ε, for some small positive ε. Then the contribution from a cluster of k + 2 cells to the effective coupling J_{ab} is of the order of ε^{k} and decreases very quickly with k. This property allows us to truncate the sum in Eq. 2. In practice, we numerically sum up the series and calculate the contributions to the couplings coming from all clusters with k + 2 ≤ 7 cells. A similar procedure gives the expansion of the fields h_{i}.
Notice that the above algorithm not only provides us with the most likely values for the couplings and fields but also with the uncertainty on those values because of the finite number of recorded data. The error estimates for the couplings are obtained from the inverse of the Hessian matrix of log P^{post} (Fisher information matrix), see SI Appendix, Section 1.
Algorithm for the Inverse I&F Problem.
In the I&F model, the membrane potential of each cell obeys the following equation, where g is the leak conductance of the membrane, C is the capacitance, I^{syn} is the current because of synaptic inputs from the other registered cells, where t_{jl} is the time at which cell j emits its lth spike. The integration kernel K(t) vanishes for t < 0 (causality) and for t > t_{s} where t_{s} is the synaptic integration time; the integral of K over the [0;t_{s}] interval is normalized to unit. The synaptic strength G_{ij}, which may be positive or negative, measures the overall change in the potential of cell i resulting from a spike emitted by cell j. The current coming from the other cells is modeled as a white noise process, with average value I_{i,} and Gaussian fluctuations uncorrelated from cell to cell. When the potential V_{i} exceeds the threshold value V_{th}, cell i emits a spike and the potential is reset to its rest value (equal to 0). Up to a rescaling of the couplings and currents we can set V_{th} = C = 1 without loss of generality; g can then be interpreted as the inverse of the leaking time of the membrane.
Eqs. 4–6 implicitly define the likelihood P[{t_{jl}}∣{G_{ij},I_{i}}] that our population of N neurons emit spikes at times {t_{jl}} given the couplings G_{ij} and the average currents I_{i}. In the absence of a priori probability on the couplings and currents, the most likely values for the latter is found through the maximization of P. In principle, P can be calculated through the resolution of many Fokker–Planck equations, one for each spike in the dataset, associated to 1dimensional Ornstein–Uhlenbeck processes with moving boundaries. In practice, this approach and related numerical approximations (15) are inadequate to process datasets with hundreds of thousands of spikes.
We have, therefore, resorted to an approximation for P, asymptotically exact when the amplitude σ of the synaptic noise Eq. 6 tends to zero: P is then dominated by the contribution coming from particular trajectories V*_{i}(t) of the potentials, called classical paths in physics (16) and optimal paths in largedeviation theory (17). The difficulty in finding the optimal potentials V*_{i}(t) consists in determining whether and when they touch the threshold (without crossing it) at intermediate times t ≠ t_{ik} (18). We have devised a fast and analytical procedure to calculate the optimal potentials, and thus the probability P, for a population of coupled cells, see SI Appendix, Section 3. As the likelihood P is a convex function of G_{ij} and I_{i} the couplings and currents with maximal probability are easily found by using a local ascent procedure, e.g., the Newton–Raphson algorithm, see SI Appendix, Section 3.
The procedure is able to determine N × (N − 1) couplings and N currents from S spikes in a computational time scaling as S × N^{2}. The algorithm has been tested on artificially generated data from networks with known couplings and currents with up to N = 160 cells and S = 10^{7} spikes. The inferred couplings are largely independent on the initial conditions over the potentials used in the I&F simulation. For moderate noise levels, the interactions are correctly inferred, with error bars of the order of 1/ \batchmode \documentclass[fleqn,10pt,legalpaper]{article} \usepackage{amssymb} \usepackage{amsfonts} \usepackage{amsmath} \pagestyle{empty} \begin{document} $$\sqrt{S}$$ \end{document} .
Acknowledgments
M. Meister and M. Berry kindly made their multielectrode recordings available to us. We thank them for discussions and encouragement. We benefited from discussions with C. Bargmann, D. Butts, T. Holy, S. Kuehn, B. Knight, and R. Yuste and from their comments on our manuscript.
Footnotes
 ^{1}To whom correspondence should be addressed. Email: cocco{at}lps.ens.fr

Author contributions: S.C., S.L., and R.M. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/cgi/content/full/0906705106/DCSupplemental.
References
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Meister M,
 Lagnado L,
 Baylor DA
 ↵
 ↵
 ↵
 Shlens J,
 et al.
 ↵
 ↵
 Tkacik G,
 Schneidman E,
 Berry MJ, II,
 Bialek W
 ↵
 Tang A,
 et al.
 ↵
 ↵
 ↵
 Paninski L,
 Pillow JW,
 Simoncelli EP
 ↵
 ZinnJustin J
 ↵
 Dembo A,
 Zeitouni O
 ↵
Citation Manager Formats
Sign up for Article Alerts
Article Classifications
 Physical Sciences
 Physics
 Biological Sciences
 Neuroscience