Pairwise Models Are Not Sufficient to Describe the Distribution of Activity Patterns of Large Networks Responding to Natural Stimuli.

The most general description of the population activity of n neurons, which uses all possible correlation functions among cells, can be written using the maximum entropy principle (36, 37) as:

where must be found numerically, such that the corresponding averages, such as <x_i>, <x_ix_j>, or <x_ix_jx_k> over , agree with the experimental ones; the partition function Z is a normalization factor. However, has an exponential number of parameters (interactions). Thus, this representation would be useful only if a small subset of parameters would be enough to capture the network behavior.

For small groups of 10–20 neurons, it has been shown that the minimal model of the network that relies only on pairwise interactions, namely, the maximum entropy pairwise model, with only s and s in Eq. 1, captures more than 90% of the correlation structure of the network patterns (19, 20). Learning , which has only parameters, is a known hard computational problem (29). For 100 cells, we cannot even enumerate all network states and must rely on a Monte Carlo algorithm to find the parameters of (Methods). The resulting model gave a seemingly accurate description of the network spatial patterns (Fig. 1D), correcting the orders-of-magnitude mistakes that makes. We obtained similar results for the spatiotemporal activity patterns of 10 neurons over 10 consecutive time bins, which we learned using the temporal correlations between cells (Fig. 1E and Methods).

However, closer inspection reveals that the most common population activity patterns were actually misestimated (Fig. 1 D and E), whereas the rare patterns were predicted well by. This reflects a real failure of the pairwise model in this case and not a numerical artifact of our algorithm (Fig. S1). Moreover, is an accurate model of the population's response to artificial movies without spatial structure, such as spatiotemporal white noise (24) (Fig. 1F and Fig. S1). Thus, the inaccuracy of the pairwise model for the retina responding to natural stimuli is the signature of the contribution of higher-order interactions to the neural code of large populations.

High-order interactions, which reflect the tendency of groups of triplets and quadruplets, for example, to fire synchronously (or to be silent together) beyond what can be explained from the pairwise relations between cells, have been studied before in small networks (4, 26, 27, 36, 37). These synchronous events are believed to be crucial for neural coding and vital for the transmission of information to downstream brain regions (4, 6, 7, 10). However, unlike the case of small groups of neurons, where one can systematically explore all interaction orders, how can we find the high-order interactions that define the neural population code in large networks among the exponential number of potential ones?

Learning the Functional Interactions Underlying the Neural Code from Reliably Sampled Activity Patterns.

We uncovered the structure of the population code of the retina using a network model that gives an extremely accurate approximation to the empirical distribution of population activity patterns. We note that if we were given the true distribution, , over all the network activity patterns, we could then infer the interactions among neurons simply by solving a set of linear equations for the parameters. For a specific activity pattern , substituting the appropriate 0’s and 1’s in Eq. 1, we get that , where the sums are only over terms where all x_is are 1. The interaction parameters that govern the network could then be calculated recursively, starting from the low-order ones and going up (36, 38):

In general, this approach would not be useful for learning large networks, because the number of possible activity patterns grows exponentially with network size. We would therefore rarely see any pattern more than once, regardless of how long we observe the system, and as a result we would not be able to estimate the parameters in Eq. 2. However, the distribution of pattern occurrences in large networks of neurons had many patterns, both spatial and spatiotemporal ones, that appeared many times during the experiment (Fig. 1B). This feature of the pattern count distribution is a direct result of the sparseness of the neural code and the correlation structure of neural activity, and does not hold for other kinds of distributions (Fig. S2).

We can then use the frequently occurring patterns to estimate the interactions governing the joint activity of neural populations accurately and to build the RI model that approximates the distribution of neural population responses:

where are the dependencies of any order among neurons, whose interaction parameters, , are inferred, as in Eq. 2, from patterns μ that appeared more than n_RI times (our learning threshold) in the training data. The order of interactions in the model is therefore determined by the nature of the patterns in the data and not by an arbitrarily imposed finite order; higher-order interactions beyond pairs, if needed, would stand for synchronous events that could not be explained by pairwise statistics.

Sparse Low-Order Interaction Network Model Captures the Code of Large Populations of Neurons.

We learned an RI model for the data from Fig. 1, and the resulting model had far fewer parameters than the pairwise maximum entropy model but included a small set of higher-order interactions. This RI model predicted the probability of the spatial population activity patterns in a 1-h test data extremely well (Fig. 2A). Its accuracy was actually within the experimental sampling error of the data itself, which we quantified using the log-likelihood ratio of the model and the empirical data, . The RI model was equally accurate in predicting the spatiotemporal patterns in the population (Fig. 2B). For comparison, we also present the poor performance of the independent model and the systematic failure of the pairwise model (Fig. 2 A and B). We found similar behavior of the models when we used only the subpopulation of “Off” cells in the salamander retina recorded in our experiment (Fig. S3). Moreover, the interactions among the Off cells network were similar to the ones we found using the full network of all ganglion cells recorded (Fig. S3). The RI model remained a highly accurate model for different time bin values used to analyze neural activity (Fig. S4), whereas became even less accurate for larger time bins of 40 ms.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

RI model is a highly accurate and easy-to-learn model of large neural population activity. (A) For all spatial activity patterns of the cells from Fig. 1A that occur with probability P (abscissa), we plot the mean and SD (error bar) of the log-likelihood ratio of the model and data (ordinate): . The blue funnel marks the estimated 95% confidence interval of the empirical measurement. Independent model (Left), pairwise maximum entropy model (Center), and RI model (Right). (B) Model predictions for population spatiotemporal activity patterns of 10 neurons over 10 time steps. Other details are as in A.

Because of the sparseness of the neural code, the number of interaction parameters that the RI model relies on is very small, much smaller than the number of all possible pairs, and these parameters are low-order ones. The functional architecture of the retina responding to a natural movie, inferred from the RI model, is shown in Fig. 3A, overlaid on the physical layout of the receptive field centers of the cells. We found that beyond a selective group of single-cell biases (α_i) and pairwise interactions (β_ij), there were significant higher-order interactions (e.g., γ_ijk, δ_ijkl). However, importantly, the model was very sparse, using only hundreds of pairwise and triplet interactions, tens of quadruplets, and a single quintuplet (with a learning threshold of ). The pairwise interactions were mostly positive; that is, cells tended to spike together or to be silent together, and their strength decayed with distance (Fig. 3D). These higher-order interactions had a similar decay with distance, but the interactions among nearby cells were mostly negative, indicating that neurons were more silent together than would be expected from pairwise relations alone (Fig. 3D). We note that the decay of interactions with distance is stronger than the decay of spatial correlations in the natural movies we showed the retina (Fig. 3D, Inset). We found a similar decay of pairwise and higher-order temporal interactions among cells with the time difference between them as well as a similar polarity of the interaction sign (Fig. 3E).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Functional interaction maps of retinal ganglion cells responding to a natural movie, based on the RI model. (A) Pairwise interactions, β_ij, learned by the RI model for 99 retinal ganglion cells responding to a natural movie. (Upper) The center of the receptive field of each of the neurons is marked by a black dot (mapped using flickering checkerboard stimuli; Methods). The saturation level of an edge between vertices is set according to the absolute value of the interaction. (Lower) Histogram of the pairwise interaction terms. (B) Same as in A, except for the triple-wise interaction terms, γ_ijk, of the RI model. Note that each of the interactions is shown as a three-line star that connects each of the participating cells and the center of mass of their locations. (C) Same as in A, except for the quadruple-wise interaction terms, δ_ijkl, of the RI model; interactions are shown as a four-line star that connects the cells and the center of mass of their locations. (D) RI model pairwise interaction values are plotted as a function of distance between the receptive field centers of the interacting neurons (blue). Higher-order interaction values are plotted as a function of their average pairwise distance (yellow). (Inset) Decay of pairwise interaction terms (blue; line is a scaled exponential fit to the blue dots in D) and of the correlations between pixels in the movie (black) as a function of the distance on the retina. corr., correlation. (E) RI model pairwise interaction values are plotted as a function of the temporal delay between the interacting neurons (blue dots; solid line corresponds to exponential fit). Higher-order interaction values are plotted as a function of the maximal lag between pairs (yellow dots; solid line corresponds to exponential fit).

Higher-Order Interactions in the Neural Code of the Retina Responding to Natural Scenes.

We found that the RI model was much better than the pairwise model in describing the retinal population's response to different natural movies (Fig. 4 A and B). The benefit of the high-order interactions that the RI models identify was even more pronounced when the whole retina patch was presented with the same stimulus, as we found for a “natural pixel” (NatPix) movie in which the stimulus shown was the luminance of one pixel from a natural movie, retaining natural temporal statistics, and for Gaussian full-field movies (Methods). The RI model's accuracy in describing the population activity patterns was similar to that of the model (but with far fewer parameters) for spatiotemporal white noise stimuli that had no temporal or spatial correlations (Fig. 4 A and B). The RI model was also similar to the model in capturing the retinal response to a manipulated natural movie, which retained the pairwise correlations in the movie but with randomized phases (Methods). Fig. 4B summarizes the average gain in accuracy of the RI model compared with for many groups of neurons in salamander and archer fish retinae responding to artificial, naturalistic, and natural stimuli. Our results indicate that as the stimulus becomes more correlated, the network response becomes more correlated with a growing contribution of high-order interactions; thus, the improvement introduced by the RI model becomes more pronounced. We conclude that the higher-order interactions in the neural response to natural movies are driven by the higher-order statistics in natural scenes. Furthermore, most of the high-order interactions that we found for widely correlated stimuli were gone when we constructed the conditionally independent response of the cells (Fig. 4C). This is done by shuffling the responses of each neuron to a repetitive movie, thereby retaining only stimulus-induced correlations among neurons. Thus, beyond the stimulus statistics, the retinal circuitry does play an important role in generating the functional higher-order interactions in the response.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Higher-order interactions in neural response are driven by long-range stimulus correlations and give more accurate, easier to learn, and hierarchical models of the neural population activity. (A) Comparison of the accuracy of the second-order maximum entropy model and the RI model, measured by the absolute log-likelihood ratio , where is summed only over patterns x, which appear more than once in the experiment. is plotted against , such that dots below the equality line (black) demonstrate better performance of the RI model. Each dot of the same color represents one network of 50 neurons taken from one experiment. Different natural and artificial visual stimuli were used. (B) Overall summary of average and STD of for many groups of 50 neurons from different retinae responding to many different kinds of movies. Datasets used for the salamander retinas had about 100 cells recorded simultaneously (Methods). In the archer fish retina (purple), 38 neurons were recorded; averages are over randomly selected groups of 30 neurons. (C) Histograms of the number of interactions of different orders observed in the real data (green bars) compared with the number observed in conditionally independent surrogate data (gray bars) for a NatPix movie (Upper) and a natural movie (Lower). Cond. indep., conditionally independent. (D) Accuracy of the RI model as a function of the number of parameters in the model. is plotted against the number of parameters in the model for the spatial activity patterns of the 99-neuron network. Results were averaged over 10 randomly selected train and test data (STD error bars are smaller than marker size). The maximum entropy pairwise model (red dot) and the independent model (gray dot) are also shown for comparison. The number of parameters was varied by changing the learning threshold, n_RI. (E) RI model parameters averaged over all 30-neuron subnetworks containing the 29 nearest neighbors of each neuron are plotted against their value estimated from the full 99-neuron network. Error bars represent SDs over different randomly selected subsets of data (x axis) and SDs across different subnetworks (y axis). sub-nets, subnetworks.

We emphasize that f_RI was much more accurate than using a considerably smaller number of parameters. In particular, although the full-pairwise model requires almost 5,000 parameters for 99 neurons, the RI model was already more accurate with just 450 parameters (Fig. 4D). The model continued to improve as we added more parameters or lowered the learning threshold n_RI, or when we added training samples (Fig. S5). The cross-validation performance continued to improve even when we used patterns that appeared only twice in the data to train the model, indicating that the model was not overfit (other measures of model performance and validation are shown in Fig. S6). We found similar behavior of the RI and models for spatiotemporal patterns (Fig. S5).

We note that the RI model is not guaranteed to be a normalized probability distribution over all the population activity patterns. The reason is that by using only the reliable patterns to learn the model, we are inevitably missing nonzero interactions, corresponding to less common patterns that are not used in the training. Unlike other such pseudolikelihood models, which are common in machine learning and statistics (29, 38, 39), the RI model does give an extremely accurate estimate of the actual probability for almost all patterns observed in the test data. The patterns for which it would fail are typically so rare that we would not see them even in a very long experiment.

Hierarchical Modular Organization of the Interaction Network Underlying the Code of Large Neural Populations.

For networks much larger than 100 neurons, we expect that no activity pattern will appear more than once even in a very long experiment because of neural noise; thus, the RI approach could not be simply applied. However, because the interactions of all orders decayed with distance between cells, we hypothesized that the RI model might be applied in a hierarchical manner. Indeed, we found that the interactions inferred from smaller overlapping subnetworks of spatially adjacent neurons (based on their receptive fields) were very similar to those estimated using the full network (Fig. 4E). This surprising result reflects an overlapping modular structure of neuronal interactions (40, 41); namely, neurons directly interact only with a relatively small set of adjacent neurons, known as their “Markov blanket” (28), and the interaction neighborhood of each neuron overlaps with that of its neighbors. This is not an ad hoc approximation but rather a mathematical equality based on the local connectivity structure of the network (42). This suggests that the RI approach may be scalable and useful to model and analyze much larger networks, also beyond the retina.

Contribution of Higher-Order Interactions to Stimulus Coding.

The average stimulus preceding the joint spiking events of pairs and of triplets of ganglion cells (pattern-triggered average) was significantly different from what would be expected from conditionally independent cells (Fig. 5 A–C), reflecting that the population codebook carried information differently than what could be read from the single cells or pairs (4, 7). To demonstrate the role of higher-order interactions in stimulus coding under natural conditions, we presented the retina repetitively with a 50-s natural movie and trained an RI model on the responses of groups of 10 neurons to each 1-s clip of the movie using half of the repeats. We found that these models could discriminate between the clips very accurately (Fig. 5D). For larger networks, more samples are required; thus, we presented the same retina with two long and different full-field stimuli: one with natural temporal statistics (NatPix) and one in which each frame was sampled from a Gaussian distribution [full field flicker (FFF); Methods]. For each of the stimuli, an RI model was learned from the population responses f_NatPix and f_FFF, respectively. New stimuli (not used for the training) were sampled from both stimulus classes and presented to the retina, and they were then classified according to the log-likelihood ratio, . Using the RI models, we found that we could accurately classify new stimulus segments within a few hundred milliseconds (Fig. 5E) far better than (almost threefold faster) or the independent model (almost 10-fold faster). Finally, we also show that a limited memory-based estimation of f_RI, which is learned using only frequent patterns that appeared at least once every few minutes in the training data (Fig. 5F), proves to be highly accurate. This limited memory approach can be readily converted into an online learning algorithm: As new samples are presented, only the parameters corresponding to patterns that appear frequently are retained, whereas parameters learned from patterns that have not appeared in the recent past are completely forgotten.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Role of high-order interactions in encoding and decoding of natural stimuli. (A) Examples of the average stimulus preceding a specific triplet (Left, blue) or quadruplet (Right, blue) of cells spiking within 20 ms for FFF stimulus; stimuli are shown in contrast-normalized units. The STAs of the single neurons comprising the triplet/quadruplet are shown in gray. The vertical scale bar is expressed in units of stimulus STD. STA, spike-triggered average. (B) Absolute value of the triplet STA peak is plotted against the absolute value of the peaks of the average STAs of the subsets comprising the triplet (single neurons, yellow; pairs, purple; normalized contrast units). As in A, triplet STAs correspond to significantly greater changes in luminance (P < 1e-9, n = 70, two-sided paired Wilcoxon signed rank test). (C) Norm of the difference between pair STAs (in normalized contrast units) and the average STAs of the single neurons is plotted for real neural data (x axis) vs. conditionally independent surrogate data (y axis). Conditionally independent pair STAs were significantly more similar to the single-neuron STAs (P < 1e-3, n = 57, two-sided paired Wilcoxon signed rank test). (D) Accurate classification of short natural movie clips by the RI model. A 50-s natural movie was presented repetitively to the retina 101 times. We divided the movie into 50 consecutive 1-s clips, and an RI model was constructed for each clip, based on training data, for 10 randomly chosen groups of 10 neurons. Later, the likelihood of each trained model was estimated on test data taken from the different clips. The (i, j)th entry in the matrix is the log ratio between the likelihood of the data taken from clip j according to the model trained on clip i and the likelihood of the same data according to the model trained on clip j. Negative values indicate that the model trained on the tested clip gave a higher likelihood. (E) For large groups of neurons, the RI model rapidly outperforms the independent and pairwise models. The same group of neurons was presented with naturalistic (NatPix) and artificial (FFF) stimuli (Methods). RI models were trained on each stimulus. The average log-likelihood ratio, , is plotted for test data taken from the NatPix (red) or FFF (green) stimuli as a function of time (shaded area represents SEM over 100 repeats). Classification by and models is shown for comparison. Note that a difference of 10 in the ordinate means a likelihood ratio of over 1,000. Indep., independent. (F) Performance of a limited memory (“online”) version of the RI model, where each interaction is estimated from the corresponding activity pattern only if it has appeared at least once every min in the experiment (Methods). (Inset) of the limited-memory RI model and the data as a function of the time window used for training the model. Data from the main panel are marked by the orange dot.