A Tractable Method for Describing Complex Couplings between Neurons and Population Rate

Christophe Gardella; Olivier Marre; Thierry Mora

doi:10.1523/ENEURO.0160-15.2016

Abstract

Neurons within a population are strongly correlated, but how to simply capture these correlations is still a matter of debate. Recent studies have shown that the activity of each cell is influenced by the population rate, defined as the summed activity of all neurons in the population. However, an explicit, tractable model for these interactions is still lacking. Here we build a probabilistic model of population activity that reproduces the firing rate of each cell, the distribution of the population rate, and the linear coupling between them. This model is tractable, meaning that its parameters can be learned in a few seconds on a standard computer even for large population recordings. We inferred our model for a population of 160 neurons in the salamander retina. In this population, single-cell firing rates depended in unexpected ways on the population rate. In particular, some cells had a preferred population rate at which they were most likely to fire. These complex dependencies could not be explained by a linear coupling between the cell and the population rate. We designed a more general, still tractable model that could fully account for these nonlinear dependencies. We thus provide a simple and computationally tractable way to learn models that reproduce the dependence of each neuron on the population rate.

Significance Statement

The description of the correlated activity of large populations of neurons is essential to understand how the brain performs computations and encodes sensory information. These correlations can manifest themselves in the coupling of single cells to the total firing rate of the surrounding population, as was recently demonstrated in the visual cortex, but how to build this dependence into an explicit model of the population activity is an open question. Here we introduce a general and tractable model based on the principle of maximum entropy to describe this population coupling. By applying our approach to multielectrode recordings of retinal ganglion cells, we find complex forms of coupling, with the unexpected tuning of many neurons to a preferred population rate.

Introduction

An important feature of neural population codes is the correlated firing of neurons. Manifestations of collective activity are observed in the correlated firing of individual pairs of neurons (Arnett, 1978), and through the coupling of single neurons to the activity in its surrounding population (Arieli et al., 1996; Tsodyks et al., 1999). These correlations, whether they are evoked by common inputs or result from interactions between neurons, imply that the neural code must be studied through the collective patterns of activity rather than by individual neuron.

As the number of possible firing patterns in a population grows exponentially with its size, they cannot be sampled exhaustively for large populations. Several modeling approaches have been suggested to describe the collective activity patterns of a neural population (Martignon et al., 1995; Schneidman et al., 2003, 2006; Pillow et al., 2008; Cocco et al., 2009; Tkačik et al., 2014). In these approaches, a small number of statistics (e.g., mean firing rate, pairwise correlations) is measured to constrain the parameters of the model. Models are then evaluated on their ability to predict statistics of the population activity that were not fitted to the data. These models are computationally hard to infer, and one must usually have recourse to approximate methods to fit them.

Recently, Okun et al. (2015) investigated how the activity of the whole population influenced the behavior of single neurons in the primary visual cortex of awake mice and monkeys. In particular, they studied the role of the correlation between neurons and the summed activity of the population, called the population rate. To assess whether these couplings between neurons and population activity were sufficient to describe the correlative structure of the code, synthetic spike trains preserving these couplings were generated and compared to data. However, the numerical method used to generate synthetic spike trains is computationally heavy, and is unable to predict the probability of particular patterns of spikes, as most of them are unlikely to ever occur.

Here we introduce a new method, based on the principle of maximum entropy, to exactly account for the coupling between individual neurons and the population rate. This model is tractable, meaning that predictions for the statistics of the activity can be derived analytically. The gradient and Hessian of the likelihood of the model can thus also be computed efficiently, allowing for fast inference using Newton’s method. Compared with previous methods (Okun et al., 2015), our method can fit hours of large-scale recordings of large populations in a few seconds on a standard laptop computer. We tested it on recordings of the salamander retina (160 neurons). We uncovered new ways for individual neurons to be coupled to the population, where a single neuron is tuned to a particular value of the population rate, rather than being monotonically coupled to the population.

Materials and Methods

Recordings from retinal ganglion cells

We analyzed previously published ex vivo recordings from retinal ganglion cells of the tiger salamander (Ambystoma tigrinum; Tkačik et al., 2014). In brief, animals were killed according to institutional animal care standards. The retina was extracted from the animal, maintained in an oxygenated Ringer’s solution, and recorded on the ganglion cell side with a 252-electrode array. Spike sorting was performed with custom software (Marre et al., 2012), and N = 160 neurons were selected for the stability of their spike waveforms and firing rates, and the lack of refractory period violation.

Maximum entropy models

We are interested in modeling the probability distribution of population responses in the retina. The responses are first binned into 20 ms time intervals. The response of neuron i in a given interval is represented by a binary variable, , which takes value 1 if the neuron spikes in this interval, and 0 if it is silent. The population response in this interval is represented by the vector of all neuron responses (Fig. 1A). We define the population rate K as the number of neurons spiking in the interval .

Figure 1.

A maximum entropy model for population coupling. A, Spikes trains are recorded with a multielectrode array and binned in 20 ms time windows. We study the dependence between the binned response of each neuron and the population rate K, defined as the summed activity of all neurons. B–D, The linear-coupling model fits the following three observables with high accuracy: the population rate distribution (B); the cells firing rates (C); and population couplings (D). For each observable, the model fit is plotted against empirical values.

We build three models for the probability of responses . These models reproduce some chosen statistics, meaning that these statistics have same value in the model and in empirical data. The first model reproduces the firing rate of each neuron and the distribution of the population rate. The second model also reproduces the correlation between each neuron and the population rate. The third model reproduces the whole joint probability of single neurons with the population rate. It is a hierarchy of models, because the statistics of each model are also captured by the next one.

Minimal model

We build a first model that reproduces the firing rate of each neuron, , and the distribution of the population rate, P(K). We also want the model to have no additional constraints, and thus be as random as possible. In statistical physics and information theory, the randomness of a distribution P is measured by its entropy S(P): (1)

where the sum runs over all possible states. The maximum entropy model is the distribution that maximizes this entropy while reproducing the constrained statistics. Using the technique of Lagrange multipliers (see Mathematical derivations), one shows that the model must take the following form: (2)

where the parameters α_i , i = 1,…,N and β_K , K = 0,…,N must be fitted so that the distribution of Equation 2 matches the statistics and P(K) of the data. Z is a normalization factor. Note that depends on the state through . We refer to this distribution as the minimal model, as no explicit dependency between the activity of individual neurons and the population rate is constrained.

Linear-coupling model

The second model reproduces and P(K) as before, as well as the linear correlation between each neuron response and the population rate K, for i = 1,…,N. It takes the following form (see Mathematical derivations): (3)

Analogously to the minimal model, the parameters α_i , β_K , and γ_i are inferred so that the model agrees with the mean statistics , P(K), and of the data. Importantly, despite their common notation, the values of the fitted parameters α_i and β_K are different from the ones fitted in the minimal model (see Mathematical derivations).

Complete-coupling model

The third model reproduces the joint probability distributions of the response of each neuron and the population rate . It takes the following form (see Mathematical derivations): (4)

The parameters are inferred so that the model agrees with the data on for each (i, K) pair. Note that depends on the state . We refer to this model as the complete-coupling model since it reproduces exactly the joint probability between each neuron and the population rate.

Model solution

The minimal and linear-coupling models can be written in the same form as the complete-coupling model (Eq. 4), but with constraints on the form of h_iK . In the minimal model, the matrix h_iK is constrained to have the form . In the linear-coupling model, it is constrained to have the form α_i + β_K + γ_i K. In the complete-coupling model, the matrix h_iK has no imposed structure, and all its elements must be learned from the data.

Since all the considered models can be viewed as subcases of the complete-coupling model (Eq. 4), we only describe the mathematical solution to this general case. First, we describe how to solve the direct problem, i.e., how to compute statistics of interest, such as P(σ_i ,K), from the parameters h_iK . In the next section, we explain how to solve the inverse problem—the reverse task of inferring the model parameters from the statistics—which relies on the solution to the direct problem.

A model is considered tractable if there exists an analytical expression for the normalization factor, (5)

allowing for its rapid computation (e.g. in polynomial time in N). All statistics of the model, such as P(σ_i ,K) or covariances between pairs of neurons, can then be calculated efficiently through derivatives of Z (see Mathematical derivations). In general, maximum entropy models are not tractable, because sums of the kind in Equation 5 involve a sum over an exponential number of terms (2 ^N ). Fortunately, in our case, the technique of probability-generating functions provides an expression for Z that is amenable to fast computation using polynomial algebra (see Mathematical derivations), as follows: (6)

where denotes the coefficient of polynomial Q of order Xⁿ .

Model inference

We now describe how to fit the models to experimental data. The inference of the model parameters is equivalent to a problem of likelihood maximization (Ackley et al., 1988). The model reproduces the empirical statistics exactly when the parameters maximize the likelihood of experimental data measured by the model, , where are the n activity patterns recorded in the experiment, assumed to be independently drawn.

In practice, we maximized the normalized log-likelihood instead of L, which is equivalent theoretically but is more convenient for computation. We used Newton’s method to perform the maximization. This method requires computation of the first and second derivatives of the normalized log-likelihood. These derivatives can be expressed as functions of mean statistics of the model and can be calculated using the solution to the direct problem sketched in the previous section, and detailed in the Mathematical derivations section. Because the model is tractable, these mean statistics can be computed quickly, and the model can be inferred rapidly.

For the minimal model, the optimization was performed over the parameters and . For the linear-coupling model, the optimization was performed over these two sets of parameters, as well as . For the complete-coupling model, the optimization was performed over all elements of the matrix . We stopped the algorithm when the fitting error was smaller than 10⁻⁶ (see Mathematical derivations).

Regularization

Prior to learning the model, we regularized the empirical population rate distribution P(K) and conditional firing rates to mitigate the effects of low sampling noise. This regularization allowed us to remove zeros from the mean statistics, avoiding issues with the fitting procedure. We performed this regularization using pseudocounts (see Mathematical derivations).

Tuning curves in the population rate

We define the tuning curve of neuron i in the population rate as the conditional probability of neuron i to spike given the summed activity of all neurons but i, . It is equal to the following: (7)

where we can use and . Each of these quantities can be computed using the solution to the direct problem (see Mathematical derivations).

We then tested for each neuron whether its tuning curve had significant local maxima. We first identified the set of for which was significantly larger than points below and above . To assess significance, we measured the SD of the difference across 100 training sets consisting of random halves of the dataset. The difference was said to be significant when it was 5 SDs above 0.

For the cells for which the presence of a maximum was determined, we then evaluated the location of the maximum, , by taking the median of the maxima determined for each training set. We inferred the presence and position of minima in a similar way.

To estimate the quality of the model prediction for the tuning curve, we quantified how the model differed from the data. We trained the model on 100 random training sets and computed , the difference between in the testing data and predicted by the model, measured by the Kullback–Leibler (KL) divergence. The KL divergence between two distributions P and Q of a random variable x is: . We regularized in the testing set before computing the KL divergence. To measure sampling noise, we computed the difference between the testing and the training sets, D_KL (test train), where was regularized in both sets. The normalized KL divergence, z, is defined as the difference between D_KL (test model) and D_KL (test train), divided by the SD, as follows: (8)

In other words, it measures by how many SDs the data differ from the model.

Quality of the model

Pairwise correlations

In order to measure the quality of the predictions of correlations between pairs of neurons and , we used cross-validation. We randomly divided the dataset into 100 training and testing sets half the size of the data, and learned the model on the training sets. The correlations of each testing set were then predicted with the model . The quality of the model prediction was measured by a goodness-of-fit index quantifying the amount of correlations predicted by the model. We define it as follows: (9)

where is the correlation in the corresponding training set. The numerator of Equation 9 is the part of the correlations in the testing set that is predicted by the model, and the lower one is a normalization correcting for sampling noise. We have C = 1 when the model perfectly accounts for the correlations of the training set. When the model completely ignores correlations, as in a model of independent neurons, c_ij = 0, then C = 0.

Likelihood

Using the models learned on the same 100 training sets, we computed the likelihood of responses in the testing sets for the minimal, linear-coupling, and complete-coupling models. In this article, the log-likelihood is expressed in bits, using binary logarithms. We then computed the improvement in mean log-likelihood compared to the minimal model, for complete versus linear models as the ratio , where is the mean over training sets and is the mean over responses in each testing set.

Multi-information

The multi-information (Cover and Thomas 1991; Schneidman et al., 2003) quantifies the amount of correlative structure captured by a model. It is defined as the difference between the entropy of a model of independent neurons reproducing firing rates and the empirical data: I = S_indep − S_data. Here is the entropy of the spike patterns measured by their frequencies, , in the data, and S_indep is the entropy if all neurons were independent.

The entropy of a maximum entropy model is by construction higher than that of the real data, , because the model has maximum entropy given the statistics it reproduces. Its entropy is also smaller than S_indep, provided that constraints include the spike rates, because the model has more structure and reproduces more statistics than if neurons were independent. Thus, the fraction of correlations that is accounted by the maximum entropy model, , where , can be viewed as a measure of how well the model captures the correlative structure of responses. The true multi-information I can only be calculated for small groups of neurons ( ), because P_data requires evaluation of 2 ^N pattern frequencies, which is prohibitive for large networks.

Results

Tractable maximum entropy model for coupling neuron firing to population activity

The principle of maximum entropy (Jaynes, 1957a,b) provides a powerful tool to explicitly construct probability distributions that reproduce key statistics of the data, but are otherwise as random as possible. We introduce a novel family of maximum entropy models of spike patterns that preserve the firing rate of each neuron, the distribution of the population rate, and the correlation between each neuron and the population rate, with no additional assumptions (Fig. 1A). Under these constraints, the maximum entropy distribution over spike patterns in a fixed 20 ms time window is given by the following (see Materials and Methods): (10)

where equals 1 when neuron i spikes within the time window, and 0 otherwise, is the population rate, and Z is a normalization constant. The parameters and must be fitted to empirical data. We refer to this model as the linear-coupling model, because of the linear term in the exponential.

Unlike maximum entropy models in general, this model is tractable, meaning that its prediction for the statistics of spike patterns has an analytical expression that can be computed efficiently using polynomial algebra. This allows us to infer the model parameters rapidly for large populations on a standard computer, using Newton’s method (see Materials and Methods). We learned this model in the case of a population of N = 160 salamander retinal ganglion cells, stimulated by a natural movie. It took our algorithm 14 s to fit the 3N–2 model parameters (see Mathematical derivations) so that the maximum discrepancy between the model and the data was <10⁻⁶ (Fig. 1B–D).

The linear-coupling model provides a rigorous mathematical formulation for the hypotheses underlying the modeling approach of Okun et al. (2015) applied to cortical populations. In that work, synthetic spike trains were generated by shuffling spikes from the original data so as to match the three constraints listed above on the single-neuron spike rates, the distribution of population rates, and their linear correlation. Shuffling data (i.e., increasing randomness and hence entropy while constraining mean statistics) have previously been shown to be equivalent to the principle of maximum entropy in the context of pairwise correlations (Bialek and Ranganathan, 2007). Our formulation provides a fast way to learn the model and to make predictions from it, as we shall see below. In addition, it allows us to calculate the probability of individual spike patterns (Eq. 10), which a generative procedure such as the one in Okun et al. (2015) cannot.

Tuning curves of single neurons to the population activity

We wondered whether the linear-coupling model could explain how the response of single neurons depended on the population rate. We examined the firing probability of neuron i as a function of the summed activity of the other neurons , denoted by . This quantity can be viewed as the tuning curve of neuron i in response to the rest of the population. It can be calculated analytically from the parameters of the model (see Materials and Methods) and compared with empirical values. The tuning curves of four representative cells are shown in Figure 2A–D.

Figure 2.

Tuning curves of single neurons as a function of the population rate. A–D, Spiking probability of neuron i conditioned on the summed activity of all other neurons , as observed in the data (black curves; SE is shaded in gray) and predicted by linear-coupling model (red curves). Each subfigure corresponds to a different representative cell. E, Histogram of the difference between the preferred population rate—at which the tuning curve is maximal—observed in the data, , and predicted by the linear-coupling model, . Data are shown for the 38 cells that had at least one local maximum both in the linear model and in the data. When the empirical tuning curve had two local maxima, the closest one to the model prediction was chosen. F, Histogram of the normalized KL divergence between the observed joint distributions and its prediction by linear-coupling model. The arrows indicate the value for the four example cells A–D. The vertical line shows a normalized divergence of 2, meaning that cells sitting on its right deviate from the linear-coupling model by >2 SDs.

The linear-coupling model predicts a variety of tuning curves (in red), from sublinear to superlinear. Although its prediction was qualitatively close to the empirical value for some cells (Fig. 2A), the model generally did not account well for the coupling between and . A majority of cells (85 of 160) displayed a local maximum in their empirical tuning curves, at some preferred value of the population activity to which the neuron is tuned. The model did not predict the existence of this maximum in 47 of these 85 cells (Fig. 2C). Even when it did, the location of the maximum, , was poorly predicted, as can be seen by the distribution of the difference between the model and the data (Fig. 2E). In six cases, the tuning curve had two local maxima, while the model only predicted one. Another 27 cells had a minimum in their empirical tuning curve, which was never reproduced by the model (Fig. 2D). Interestingly, no cells were tuned to fire when the rest of the population is silent; even cells whose spiking activity was anticorrelated with the rest of the population had a nonzero preferred population rate, .

The model performance can be quantified by computing the KL divergence between the data and the model for the joint probability of the neuron and population activity (Fig. 2F). The KL divergence is a measure of the dissimilarity between two distributions, P and Q (Cover and Thomas, 1991), which quantifies the amount of information that is lost if we use Q to approximate P. We calculated a normalized KL divergence (see Materials and Methods) measuring by how many SDs the KL divergence between the linear-coupling model and the data deviated from what would be expected from sampling noise (Fig. 2F). A majority of cells (143 of 160) deviated by >2 SDs, meaning that their tuning curve was not well accounted for by the linear-coupling model. This observation is consistent with the failure of the model to account for the qualitative properties of their tuning curves.

Together, these results indicate that the full dependency between single cells and the population rate cannot be explained by their linear correlation only.

A refined maximum entropy model

To overcome the limitations of the linear-coupling model, and to fully account for the variety of tuning curves found in data, we built a maximum entropy model constrained to match all joint probabilities of the population rate with each single neuron response, . This model takes the following form (see Materials and Methods): (11)

where the parameters h_iK for i = 1,…,N and K = 0,…,N are fitted to empirical data, and Z is a normalization constant. Note that the linear-coupling model can be viewed as a particular case of this model, with parameters h_iK constrained to take the form . By construction, this model exactly reproduces the tuning curves of Figure 2A–D.

Although this model has many more parameters than the simpler linear-coupling model, it is still tractable, and we could infer its N(N – 1) + 1 parameters (see Mathematical derivations) in 7 s for the whole population of 160 neurons. Hereafter, we refer to this model as the complete-coupling model.

Pairwise correlations

The models introduced thus far are only constrained to reproduce the firing rate of each neuron, the distribution of population rates, and the coupling of each neuron with the population rate. We asked whether these simple models could account for correlations between individual pairs of cells, which were not fitted to the data. The correlation between two neurons, , can be calculated analytically from the model parameters (see Materials and Methods) and directly compared with the data (Fig. 3A,B).

Figure 3.

Maximum entropy models of population coupling partly account for pairwise correlations. A, B, The observed correlation coefficient between all pairs of neurons is compared with its prediction according to the linear-coupling (A) and complete-coupling (B) models. C, Distribution of pairwise correlation coefficients, as observed in the data and predicted by minimal, linear-coupling, and complete-coupling models.

In order to understand the importance of the population rate coupling for the prediction of pairwise correlations, we built a null model constrained only by the firing rates of each neuron and the distribution of the population rate. This simpler maximum entropy model reads as follows: (12)

We call it the minimal model. Interactions between neurons derive only from the fluctuations of the population activity, rather than from an explicit coupling. This model has 2N − 1 parameters (see Mathematical derivations), which are inferred using the same techniques as before.

To quantify the performance of the different models, we calculated a goodness-of-fit index ranging from 0, when the correlations were not predicted at all, to 1, when they were predicted perfectly (see Materials and Methods). This index was 0.380 ± 0.001 for the minimal model, 0.526 ± 0.002 for the linear-coupling model, and 0.544 ± 0.002 for the complete-coupling model. Thus, a substantial amount of pairwise correlations could be explained from the coupling of neurons to the population. By this measure, the complete model performed slightly (but significantly) better than the linear-coupling model.

Figure 3C shows the distribution of the pairwise correlations in the data and as predicted by the three models. The minimal model fails to reproduce the long tail of large correlations, and predicts no negative correlations, while 28% of empirical correlations are negative. By contrast, the linear-coupling and complete-coupling models predict 7.6% and 7.7% of negative interactions, respectively, and a longer tail of large correlation coefficients. Thus, the coupling to the population rate is important to reproduce both large correlations and the strong asymmetry of the distribution.

Prediction of probabilities of spike patterns

We quantified the capacity of models to describe population responses by computing the probability of responses predicted by each model. The mean log-likelihood of responses was −33.10 ± 0.06 bits for the minimal model, −30.12 ± 0.06 bits for the linear-coupling model, and −29.49 ± 0.06 bits for the complete-coupling model. The improvement in mean log-likelihood compared with the minimal model was 51.3 ± 0.5% higher for the complete-coupling model than for the linear-coupling model, meaning that nonlinear couplings to the population are important to model the probability of responses.

The multi-information I (Cover and Thomas, 1991) quantifies, in bits, the amount of correlations in the response, whether they are pairwise or of higher order (see Materials and Methods). To assess the performance of the models in capturing the collective behavior of the networks, we calculated the ratio of the multi-information explained by the model to that estimated directly from the data, I_model /I_data. This ratio gives a measure of how well the probability of particular spike patterns is predicted by the model: it is 1 when the model is a perfect description of the data, and 0 when the model assumes independent neurons with no correlation between them. Because it requires the estimation of the probability of all possible spike patterns of the populations, the multi-information can be only calculated for small populations of, at most, 20 cells.

With this measure, the linear coupling model could account for 65% of the multi-information for groups of 10 neurons, and 53% for groups of 20 neurons. The complete model slightly improved these ratios to 68% and 56%, respectively (Table 1). Thus, more than half of the correlative structure in the spike patterns could be explained by the coupling to the population rate alone.

View this table:

Table 1:

Multi-information I estimated either directly from the data or from maximum entropy models, for random subpopulations of 10 and 20 neurons (100 subpopulations each), as well as the ratio of the multi-information between model and data

Discussion

In this study, we have introduced a general computational model for coupling individual neurons to the population rate. This model formalizes and simplifies the generative procedure proposed by Okun et al. (2015) to study population coupling, and overcomes its computational difficulties. In addition, it allows for nonlinear coupling to the population rate.

We have used our model to investigate population coupling in large recordings of N = 160 retinal ganglion cells. We found that most cells had a nonlinear coupling to the population rate. In particular, a large fraction of cells were tuned to a preferred value of the population rate. Even more strikingly, a few cells had a least preferred population rate (i.e., they were more likely to spike at lower or higher populations rates). We found no cell that was maximally active when all other neurons were silent, even among cells that were anticorrelated with the population rate. These results emphasize the need for the nonlinear coupling afforded by our model, as they uncover new dependencies that do not fit within the proposed division between soloists and choristers (Okun et al., 2015), such as the tuning to a specific population rate. It would be interesting to test whether these nonlinear couplings can also be found at the cortical level.

Overall, our model reaches a similar predictive performance than what was found in the cortex. The coupling to the population rate accounted for more than half of the correlations between pairs of neurons. In Okun et al. (2015), a custom measure of the fraction of explained pairwise correlations (different from the one used in the present work) gave 0.34. Applying the same measure to our case yields a similar value (0.33). However, this similarity in performance can be due to different underlying mechanisms. In the retina, most correlations are due to common input from previous layers (Trong and Rieke, 2008), while ganglion cells do not make synaptic connections to each other. In contrast, at the cortical level, a larger part of the variability in the activity should be due to internal dynamics generated by recurrent connections (Arieli et al., 1996; van Vreeswijk and Sompolinsky, 1996; Tsodyks et al., 1999). It would be interesting to test our model on cortical data to see whether these differences result in different types of nonlinear population coupling.

Our maximum entropy model of population coupling is complementary to maximum entropy models reproducing correlations between all pairs of neurons. Pairwise models have been shown to accurately describe the collective activity of retinal ganglion cells (Schneidman et al., 2006; Shlens et al., 2006, 2009; Ganmor et al., 2011a,b; Tkačik et al., 2014), and in cortical networks in vitro (Tang et al., 2008) and in vivo (Yu et al., 2008), but they are not tractable, requiring summation over all 2^N possible spiking states in order to implement Boltzmann machine learning (Ackley et al., 1988). Alternative methods based on mean field approximations (Cocco et al., 2009; Cocco and Monasson, 2011) or Monte-Carlo simulations (Broderick et al., 2007) have been proposed. However, Monte-Carlo methods require hours of computations, although recent efforts have tried to lower these computation times for moderately large populations (Ferrari, 2015).

By contrast, the models of population couplings introduced here are much easier to solve. They are tractable, so their predictions can be computed analytically in time N ³, and their parameters can be inferred in a few seconds on a personal computer from large-scale, hour-long recordings of spike trains for a population of N = 160 neurons. These models can then be used to generate synthetic spike trains to calculate analytically response statistics, such as pairwise correlations, or to estimate the probability of particular spike trains. Compared with the shuffling method described by Okun et al. (2015), which is equivalent to the linear-coupling model, our method is simpler and computationally less intensive.

The procedure is general and can be applied to any multineuron recording of individual spikes. The speed of model inference could prove to be an important advantage when studying very large populations, which can now reach 1000 cells (Schwarz et al., 2014). In the case of the linear-coupling model, the number of parameters is also smaller, scaling with the population size N rather than N ² for the pairwise correlations model.

Note that the population coupling models introduced here belong to a different class than the pairwise models. Each class captures the features of neural responses that the other cannot: models of population coupling should be sufficient for studying the global properties of collective activity, while pairwise models are still needed to account for the detailed structure of the response statistics. Pairwise models have been reported to capture 90% of the correlations, as measured by the multi-information for populations of size N = 10 (Schneidman et al., 2006), while our model captures, at most, 70% (Table 1). Yet, pairwise models can also miss important aspects of the collective activity, such as the probability of large population rates (Tkačik et al., 2014), which is captured by our population model.

Both classes of models consider same-time spike patterns, with no regard for the dynamics of spike trains and their temporal correlations. Generalizations of pairwise maximum entropy models to temporal statistics are even harder to solve computationally (Vasquez et al., 2012; Nasser et al., 2013). By contrast, our models of population coupling are fully compatible with any model describing the dynamics of the population rate, such as that by Mora et al. (2015).

Mathematical derivations