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

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.


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 popu-lation 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., 2003Schneidman et al., , 2006Pillow 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.

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 , 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, i , 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 ϭ ͑ 1 , ѧ, N ͒ of all neuron responses (Fig. 1A). We define the population rate K as the number of neurons spiking in the interval K͑͒ ϭ ͚ iϭ1 N i . We build three models for the probability of responses P͑͒. 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, P͑ i ϭ 1͒ ϭ ͗ i ͘, 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): 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: *OM and TM contributed equally to this work.
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 ͗ i ͘ and P(K) of the data. Z is a normalization factor. Note that ␤ K͑͒ depends on the state through K͑͒. 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 ͗ i ͘ and P(K) as before, as well as the linear correlation ͗K· i ͘ between each neuron response i and the population rate K, for i ϭ 1,. . .,N. It takes the following form (see Mathematical derivations): Analogously to the minimal model, the parameters ␣ i , ␤ K , and ␥ i are inferred so that the model agrees with the mean statistics ͗ i ͘, P(K), and ͗K· i ͘ of the data. Importantly, despite their common notation, the values of the fitted parameters ␣ i and ␤ K are different from the ones

Complete-coupling model
The third model reproduces the joint probability distributions of the response of each neuron and the population rate P͑ i , K͒. It takes the following form (see Mathematical derivations): The parameters ͑h iK ͒ iϭ1, ѧ, N;Kϭ0, ѧ, N are inferred so that the model agrees with the data on P͑ i , K͒ for each (i, K) pair. Note that h iK͑͒ 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 ␣ i ϩ ␤ K . 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, allowing for its rapid computation (e.g. in polynomial time in N). All statistics of the model, such as P( i ,K) or covariances ͗ i j ͘ Ϫ ͗ i ͗͘ j ͘ 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: where Coeff͓Q, X n ͔ denotes the coefficient of polynomial Q of order X n .

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, L ϭ ͟ ␣ϭ1 n P͑ ͑␣͒ ͒, where ͑ ͑1͒ , ѧ, ͑n͒ ͒ are the n activity patterns recorded in the experiment, assumed to be independently drawn.
In practice, we maximized the normalized log-likelihood ᏸ ϭ ͑1/n͒log L 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 ͑␣ i ͒ iϭ1, ѧ, N and ͑␤ K ͒ Kϭ0, ѧ, N . For the linear-coupling model, the optimization was performed over these two sets of parameters, as well as ͑␥ i ͒ iϭ1, ѧ, N . For the complete-coupling model, the optimization was performed over all elements of the matrix ͑h iK ͒ iϭ1, ѧ, N;Kϭ0, ѧ, N . We stopped the algorithm when the fitting error was smaller than 10 Ϫ6 (see Mathematical derivations).

Regularization
Prior to learning the model, we regularized the empirical population rate distribution P(K) and conditional firing rates P͑ iԽ K͒ 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, K \i ϭ ͚ j i j . It is equal to the following: 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 K \i for which P͑ i ϭ 1 Խ K \i ͒ was significantly larger than points below and above K \i . 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, K \i ‫ء‬ , 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 D KL ͑testmodel͒, the difference between P͑ i , K͒ 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: D KL ͑PQ͒ ϭ ͚ x P͑x͒log ͓P͑x͒/Q͑x͔͒. We regularized P͑ i , K͒ 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 P͑ i , K͒ 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: 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 i and j , 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 c test,ij were then predicted with the model c model,ij . 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: where c train, ij 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 ͗͗log P complete ͑͒ Ϫ log P minimal ͑͒͘ /͗log P linear ͑͒ Ϫ log P minimal ͑͒͘ ͘ test , where ͗·͘ test 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: the entropy of the spike patterns measured by their frequencies, P data ͑͒, 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, S model Ͼ S 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, 0 Ͻ I model /I Ͻ 1, where I model ϭ S indep Ϫ S model , 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 (N Յ 20), because P data requires evaluation of 2 N pattern frequencies, which is prohibitive for large networks.

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): where i equals 1 when neuron i spikes within the time window, and 0 otherwise, K ϭ ͚ i i is the population rate, and Z is a normalization constant. The parameters ͑␣ i ͒ iϭ1, ѧ, N , ͑␥ i ͒ iϭ1, ѧ, N and ͑␤ K ͒ Kϭ0, ѧ, N must be fitted to empirical data. We refer to this model as the linear-coupling model, because of the linear term ␥ i K i 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 Ϫ6 (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 singleneuron 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 K \i ϭ ͚ j i j , denoted by P͑ i ϭ 1ԽK \i ͒. 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.
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 i and K \i . A majority of cells (85 of 160) displayed a local maximum in their empirical tuning curves, at some preferred value K \i ‫ء‬ 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, K \i ‫ء‬ , 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, K \i ‫ء‬ Ͼ 0. The model performance can be quantified by computing the KL divergence between the data and the model for the joint probability P͑ i , K͒ 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, P͑ i , K͒. This model takes the following form (see Materials and Methods): 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 h iK ϭ ␣ i ϩ ␤ K ϩ ␥ i K. 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, ͗ i j ͘ Ϫ ͗ i ͗͘ j ͘, can be calculated analytically from the model parameters (see Materials and Methods) and directly compared with the data (Fig. 3A,B).
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 This simpler maximum entropy model reads as follows: 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 linearcoupling 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 A B C 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.
continued 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 P͑ i , K͒ 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.
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 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.

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 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 Shlens et al., 2006Shlens et al., , 2009Ganmor 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 3 , 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.

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
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 2 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 , 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
Derivation of the model form Maximum entropy models A maximum entropy model is defined by a distribution that maximizes its entropy, as follows: while reproducing a set of chosen statistics. In the case where these statistics are the means of some observables ᏻ 1 ͑͒, ѧ, ᏻ M ͑͒, the form of the model is given by the following: where Z is a normalization factor. Equation 14 is obtained by maximizing the entropy while constraining the chosen statistics using the method of Lagrange multipliers. The Lagrange multipliers a are model parameters that must be adjusted so that the mean observables predicted by the model, ͗ᏻ a ͘ agree with those of the data, ͗ᏻ a ͘ data ϭ ͑1/n͒ ͚ ␣ϭ1 n ᏻ͑ ͑␣͒ ͒, where ͑ ͑1͒ , ѧ, ͑n͒ ͒ are the n activity patterns recorded in the experiment. This fitting procedure is equivalent to maximizing the likelihood of the data under the model L ϭ ͟ ␣ϭ1 n P͑ ͑␣͒ ͒, assuming that the patterns are independently drawn. The likelihood maximization problem is convex, and the distribution P͑͒ maximizing the likelihood is always unique. However, if the constrained observables are linearly related, the optimal set of a is not unique (even though the resulting distribution is), and must be set by choosing a convention.

Minimal model
In the minimal model, the statistics we constrain are P͑ i ϭ 1͒ for each neuron i, and P͑K ϭ k͒ for each k ϭ 0,...,N. They correspond to the means of the following observables: where ␦ x, y is Kronecker's delta, which is equal to 1 if x ϭ y, and 0 otherwise. Note that while in the main text we use K both as a short-hand for ͚ i i and its realization as a random variable, here we distinguish the two by using K and k, respectively. Applying Equation 14 to this choice of observables ͑ i , ␦ K, k ͒ yields the following: where each k is associated with the constraint on ͗␦ K, k ͘ and each ␣ i is associated with the constraint on ͗ i ͘. In the second line, we have used the fact that in the first sum, the only term which is nonzero is the one for which k ϭ K. For convenience, we rescale the parameters K , which will give a common form to our three models. We first set 0 ϭ 0, which is possible because the model is invariant when adding a constant to all k parameters (this changes only the normalization factor Z). We then introduce the rescaled parameters ␤ K , which are defined as ␤ 0 ϭ 0 and the model takes the following form: This model has 2N Ϫ 1 parameters: there are N coefficients ͑␣ i ͒ iϭ1 N and N ϩ 1 coefficients ͑␤ k ͒ kϭ0 N , but ␤ 0 is not used, and the model is invariant under a change in parameters, ␣ i Ј ϭ ␣ i ϩ c, ␤ k Ј ϭ ␤ k Ϫ c, for any number c.

Linear-coupling model
The linear-coupling model reproduces P͑ i ͒ and P(K), and also the linear correlation between the neuron response i and the population rate K, ͗ i K͘. The three sets of constrained observables are thus ͑ i ͒ iϭ1, ѧ, N , ͑␦ K, k ͒ kϭ0, ѧ, N , and ͑ i K͒ iϭ1, ѧ, N . With this choice of observables, Equation 14 reads as follows: where, in addition to the ␣ i and k parameters, each ␥ i parameter is associated with the constraint on ͗ i K͘. Note that, in general, the inferred parameters ␣ i and k will be different from the ones inferred in the minimal model. This is due to the fact that the set of observables i , ␦ K, k and i K are not independent. Therefore, the parameters ␥ i cannot be learned independently from ␣ i and k .
As for the minimal model, we rescale the parameters K with ␤ 0 ϭ 0 and ␤ K ϭ K /K for K Ͼ 0, as follows: This model has 3N Ϫ 2 parameters: there are 2N coefficients, ͑␣ i ͒ iϭ1 N and ͑␥ i ͒ iϭ1 N , and N ϩ 1 coefficients ͑␤ k ͒ kϭ0 N , but ␤ 0 is not used, and the model is invariant under changes in parameters ␣ i Ј ϭ ␣ i ϩ c, ␤ k Ј ϭ ␤ k Ϫ c ϩ dK, ␥ i Ј ϭ ␥ i Ϫ d for any numbers c and d.

Complete-coupling model
The third maximum entropy model reproduces the joint probability distributions between the response of each neuron and the population rate, P͑ i , K͒. The problem reduces to matching P͑ i ϭ 1, K͒ for all i ϭ 1,. . .,N and K ϭ 0,. . .,N, since P͑ i ϭ 0, K͒ can be determined through the following: where the distribution P(K) is set by the following: This holds true because K is the number of neurons spiking, so: (25) where we can then multiply both sides by P(K). Here ͗ . Խ K͘ stands for the mean conditioned by K.
Therefore, we impose that the model reproduces only the statistics P͑ i ϭ 1, K͒, which are the means of the observables i ␦ K, k . Using Equation 14 with this set of observables yields the following:

Regularization
We regularized the empirical population rate distribution P(K) and conditional firing rates P͑ iԽ K͒ using pseudocounts. If we denote by n ϭ 2.8 10 5 the total number of responses recorded during the experiment and by n K the number of responses with K spikes in the population, the distribution of population rates K was computed as follows: where P indep ͑K͒ is the distribution of K in a model of independent neurons reproducing the empirical firing rates ͗ i ͘. Similarly, if we denote by n iK the number of responses in which neuron i spiked and in which the population rate was K, the conditional firing rates were estimated as follows: where again P indep ͑ i ϭ 1 Խ K͒ is the estimate of the conditional firing rate according to the independent model. The terms scaling as play the role of pseudocounts. These pseudocounts are not taken to be uniform, but rather follow the prediction of a model of independent neurons. We used ϭ 1 so that the total weight of pseudocounts is equivalent to a single observed pattern.
Calculating statistics from the model We start by providing an analytical expression for the normalization factor, which is defined as follows: All useful statistics predicted by the model can be derived from the expression of Z, as we shall see below.
To calculate Z, we decompose it as a sum over groups of patterns with the same population activity K: Z ϭ ͚ kϭ0 N Z k with: We introduce the polynomial Q͑X͒ ϭ ͟ iϭ1 N ͑1 ϩ e h ik X͒.
Expanding Q, we can calculate its coefficient of order X k , denoted by Coeff[Q,X k ]. This coefficient is the sum of all the terms having exactly k factors e h ik , as follows: It is obtained by recursively computing the coefficients of ͟ iϭ1 n ͑1 ϩ e h ik X͒, of order up to X k , for n ϭ 1 to N, using the following relation: for any number b, polynomial F, and order X l . Z k can thus be computed in time linear in kN, and Z ϭ ͚ k Z k can be computed rapidly.
Many statistics of the model can then be calculated by deriving Z. For example, the mean observables according to the model in Equation 14 are given by the following: This formula gives the following expression for the joint distribution of i and K: (1 ϩ Xe h j,K ), X K ͬ . (39) Similarly, pairwise correlations are computed using the following formula:

Model inference
To learn the model parameters from the data, we maximized the normalized log-likelihood ᏸ ϭ ͑1/n͒logL using Newton's method. The update equation for the parameter values in Newton's method read as follows: where a is an adjustable step size taken typically between 0.1 and 1. ͑t͒ ϭ ͑ a ͑t͒ ͒ aϭ1, ѧ, M is the vector of the parameters at iteration t; and ٌᏸ and Ᏼ are the gradient and Hessian of ᏸ with respect to the parameters a . In the general context of maximum entropy models (Eq. 14), one can show that the gradient and Hessian read as follows: where we used Equation 37 for Equation 43 and a similar formula for Equation 45. Both quantities can readily be computed as derivatives of the normalization factor Z.
For time efficiency, we only updated the Hessian every 100 iterations of the algorithm. We stopped the algorithm when the fitting error reached 10 Ϫ6 . The fitting error was defined as the maximum error on P(K) and P͑ i ͒ for the minimal model; on P(K), P͑ i ͒ , and ͗K i ͘ for the linearcoupling model; and on P(K) and P͑ iԽ K͒ for the completecoupling model.
The code for the models inference is available at https:// github.com/ChrisGll/MaxEnt_Model_Population_Coupling.