Pairwise Maximum Entropy Models for Studying Large Biological Systems: When They Can Work and When They Can't

One of the most critical problems we face in the study of biological systems is building accurate statistical descriptions of them. This problem has been particularly challenging because biological systems typically contain large numbers of interacting elements, which precludes the use of standard brute force approaches. Recently, though, several groups have reported that there may be an alternate strategy. The reports show that reliable statistical models can be built without knowledge of all the interactions in a system; instead, pairwise interactions can suffice. These findings, however, are based on the analysis of small subsystems. Here, we ask whether the observations will generalize to systems of realistic size, that is, whether pairwise models will provide reliable descriptions of true biological systems. Our results show that, in most cases, they will not. The reason is that there is a crossover in the predictive power of pairwise models: If the size of the subsystem is below the crossover point, then the results have no predictive power for large systems. If the size is above the crossover point, then the results may have predictive power. This work thus provides a general framework for determining the extent to which pairwise models can be used to predict the behavior of large biological systems. Applied to neural data, the size of most systems studied so far is below the crossover point.


Introduction
Many fundamental questions in biology are naturally treated in a probabilistic setting. For instance, deciphering the neural code requires knowledge of the probability of observing patterns of activity in response to stimuli [1]; determining which features of a protein are important for correct folding requires knowledge of the probability that a particular sequence of amino acids folds naturally [2,3]; and determining the patterns of foraging of animals and their social and individual behavior requires knowledge of the distribution of food and species over both space and time [4][5][6].
Constructing these probability distributions is, however, hard. There are several reasons for this: i) biological systems are composed of large numbers of elements, and so can exhibit a huge number of configurations-in fact, an exponentially large number, ii) the elements typically interact with each other, making it impossible to view the system as a collection of independent entities, and iii) because of technological considerations, the descriptions of biological systems have to be built from very little data. For example, with current technology in neuroscience, we can record simultaneously from only about 100 neurons out of approximately 100 billion in the human brain. So, not only are we faced with the problem of estimating probability distributions in high dimensional spaces, we must do this based on a small fraction of the neurons in the network.
Despite these apparent difficulties, recent work has suggested that the situation may be less bleak than it seems, and that an accurate statistical description of systems can be achieved without having to examine all possible configurations [2,3,[7][8][9][10][11]. One merely has to measure the probability distribution over pairs of elements and use those to build the full distribution. These ''pairwise models'' potentially offer a fundamental simplification, as the number of pairs is quadratic in the number of elements, not exponential. However, support for the efficacy of pairwise models has, necessarily, come from relatively small subsystems-small enough that the true probability distribution could be measured experimentally [7][8][9]11]. While these studies have provided a key first step, a critical question remains: will the results from the analysis of these small subsystems extrapolate to large ones? That is, if a pairwise model predicts the probability distribution for a subset of the elements in a system, will it also predict the probability distribution for the whole system? Here we find that, for a biologically relevant class of systems, this question can be answered quantitatively and, importantly, generically-independent of many of the details of the biological system under consideration. And the answer is, generally, ''no.'' In this paper, we explain, both analytically and with simulations, why this is the case.

The extrapolation problem
To gain intuition into the extrapolation problem, let us consider a specific example: neuronal spike trains. Fig. 1A shows a typical spike train for a small population of neurons. Although the raw spike times provide a complete description, they are not a useful representation, as they are too high-dimensional. Therefore, we divide time into bins and re-represent the spike train as 0 s and 1 s: 0 if there is no spike in a bin; 1 otherwise ( Fig. 1B) [7][8][9]11]. For now we assume that the bins are independent (an assumption whose validity we discuss below, and in more detail in the section ''Is there anything wrong with using small time bins?''). The problem, then, is to find p true r ð Þ:p true r 1 ,r 2 , . . . ,r N ð Þ where r i is a binary variable indicating no spike (r i~0 ) or one or more spikes (r i~1 ) on neuron i. Since this, too, is a high dimensional problem (though less so than the original spike time representation), suppose that we instead construct a pairwise approximation to p true , which we denote p pair , for a population of size N. (The pairwise model derives its name from the fact that it has the same mean and pairwise correlations as the true model; see Eq. (15).) Our question, then, is: if p pair is close to p true for small N, what can we say about how close the two distributions are for large N?
To answer this question quantitatively, we need a measure of distance. The measure we use, denoted D N , is defined in Eq. (3) below, but all we need to know about it for now is that if D N~0 then p pair~ptrue , and if D N is near one then p pair is far from p true . In terms of D N , our main results are as follows. First, for small N, in what we call the perturbative regime, D N is proportional to N{2. In other words, as the population size increases, the pairwise model becomes a worse and worse approximation to the true distribution. Second, this behavior is entirely generic: for small N, D N increases linearly, no matter what the true distribution is. This is illustrated schematically in Fig. 2, which shows the generic behavior of D N . The solid red part of the curve is the perturbative regime, where D N is a linearly increasing function of N; the dashed curves show possible behavior beyond the perturbative regime.
These results have an important corollary: if one does an experiment and finds that D N is increasing linearly with N, then one has no information at all about the true distribution. The flip side of this is more encouraging: if one can measure the true distribution for sufficiently large N that D N saturates, as for the dashed blue line in Fig. 2, then there is a chance that extrapolation to large N is meaningful. The implications for the interpretation of experiments is, therefore, that one can gain information about large N behavior only if one can analyze data past the perturbative regime.
Under what conditions is a subsystem in the perturbative regime? The answer turns out to be simple: the size of the system, (B) Spike count in each bin. In this example the bins are small enough that there is at most one spike per bin, but this is not necessary-one could use bigger bins and have larger spike counts. doi:10.1371/journal.pcbi.1000380.g001

Author Summary
Biological systems are exceedingly complicated: They consist of a large number of elements, those elements interact in nonlinear and highly unpredictable ways, and collective interactions typically play a critical role. It would seem surprising, then, that one could build a quantitative description of biological systems based only on knowledge of how pairs of elements interact. Yet, that is what a number of studies have found. Those studies, however, focused on relatively small systems. Here, we ask the question: Do their conclusions extend to large systems? We show that the answer depends on the size of the system relative to a crossover point: Below the crossover point the results on the small system have no predictive power for large systems; above the crossover point the results on the small system may have predictive power. Moreover, the crossover point can be computed analytically. This work thus provides a general framework for determining the extent to which pairwise models can be used to predict the behavior of large biological systems. It also provides a useful heuristic for designing experiments: If one is interested in understanding truly large systems via pairwise interactions, then make sure that the system one studies is above the crossover point.
N, times the average probability of observing a spike in a bin, must be small compared to 1. For example, if the average probability is 1/100, then a system will be in the perturbative regime if the number of neurons is small compared to 100. This last observation would seem to be good news: if we divide the spike trains into sufficiently small time bins and ignore temporal correlations, then we can model the data very well with a pairwise distribution. The problem with this, though, is that temporal correlations can be ignored only when time bins are large compared to the autocorrelation time. This leads to a kind of catch-22: pairwise models are guaranteed to work well (in the sense that they describe spike trains in which temporal correlations are ignored) if one uses small time bins, but small time bins is the one regime where ignoring temporal correlations is not a valid approximation.
In the next several sections we quantify the qualitative picture presented above: we write down an explicit expression for D N , explain why it increases linearly with N when N is small, and provide additional tests, besides assessing the linearity of D N , to determine whether or not one is in the perturbative regime.
Quantifying how well the pairwise model explains the data A natural measure of the distance between p pair and p true is the Kullback-Leibler (KL) divergence [12], denoted D KL p true kp pair À Á and defined as The KL divergence is zero if the two distributions are equal; otherwise it is nonzero.
Although the KL divergence is a very natural measure, it is not easy to interpret (except, of course, when it is exactly zero). That is because a nonzero KL divergence tells us that p pair =p true , but it does not give us any real handle on how good, or bad, the pairwise model really is. To make sense of the KL divergence, we need something to compare it to. A reasonable reference quantity, used by a number of authors [7][8][9], is the KL divergence between the true distribution and the independent one, the latter denoted p ind . The independent distribution, as its name suggests, is a distribution in which the variables are taken to be independent, where p i r i ð Þ is the distribution of the response of the i th neuron, r i . With this choice for a comparison, we define a normalized distance measure-a measure of how well the pairwise model explains the data-as Note that the denominator in this expression, D KL p true kp ind ð Þ , is usually referred to as the multi-information [7,13,14].
The quantity D N lies between 0 and 1, and measures how well a pairwise model does relative to an independent model. If it is 0, the pairwise model is equal to the true model (p pair r ð Þ~p true r ð Þ); if it is near 1, the pairwise model offers little improvement over the independent model; and if it is exactly 1, the pairwise model is equal to the independent model (p pair r ð Þ~p ind r ð Þ), and so offers no improvement.
How do we attach intuitive meaning to the two divergences D KL p true kp pair À Á and D KL p true kp ind ð Þ ? For the latter, we use the fact that, as is easy to show, where S ind and S true are the entropies [15,16] of p ind and p true , respectively, defined, as usual, to be S p ½ ~{ P r p r ð Þ log 2 p r ð Þ. For the former, we use the definition of the KL divergence to write The quantityS S pair has the flavor of an entropy, although it is a true entropy only when p pair is maximum entropy as well as pairwise (see Eq. (6) below). For other pairwise distributions, all we need to know is thatS S pair lies between S true and S ind . A plot illustrating the relationship between D N , the two entropies S ind and S true , and the entropy-like quantityS S pair , is shown in Fig. 3.
Note that for pairwise maximum entropy models (or maximum entropy models for short), D N has a particularly simple interpretation, since in this caseS S pair really is an entropy. Using S maxent to denote the pairwise entropy of a maximum entropy model, for this case we have as is easy to see by inserting Eqs. (4) and (5) into (3). This expression has been used previously by a number of authors [7,9].
D N in the perturbative regime The extrapolation problem discussed above is the problem of determining D N in the large N limit. This is hard to do in general, but there is a perturbative regime in which it is possible. The small parameter that defines this regime is the average number of spikes produced by the whole population of neurons in each time bin. It is given quantitatively by Nndt where dt is the bin size and n the average firing rate, with n i the firing rate of neuron i. The first step in the perturbation expansion is to compute the two quantities that make up D N : D KL p true kp ind ð Þ and D KL p true kp pair À Á . As we show in the section ''Perturbative Expansion'' (Methods), these are given by where Here and in what follows we use O Nndt ð Þ n ð Þto denote terms that are proportional to Nndt ð Þ n in the limit Nndt?0. The N-dependence in Eq. (9a) has been noted previously [7], although the authors did not compute the prefactor, g ind .
The prefactors g ind and g pair , which are given explicitly in Eqs. (42) and (44), depend on the low order statistics of the spike trains: g ind depends on the second order normalized correlation coefficients, g pair depends on the second and third order normalized correlation coefficients (the normalized correlation coefficients are defined in Eq. (16) below), and both depend on the firing rates of the individual cells. The details of that dependence, however, are not important for now; what is important is that g ind and g pair are independent of N and ndt (at least on average; see next section).
Inserting Eq. (8) into Eq. (3) (into the definition of D N ) and using Eq. (9), we arrive at our main result, Note that in the regime Nndt%1, D N is necessarily small. This explains why, in an analytic study of non-pairwise model in which Nndt was small, Shlens et al. found that D N was rarely greater than 0.1 [8].
We refer to quantities with a superscript zero as ''zeroth order.'' Note that, via Eqs. (4) and (5), we can also define zeroth order entropies, These quantities are important primarily because differences between them and the actual entropies indicate a breakdown of the perturbation expansion (see in particular Fig. 4 below).
Assuming, as discussed in the next section, that g ind and g pair are approximately independent of N, n, and dt, Eq. (10) tells us that D N scales linearly with N in the perturbative regime-the regime in which Nndt%1. The key observation about this scaling is that it is independent of the details of the true distribution, p true . This has a very important consequence, one that has major implications for experimental data: if one does an experiment with small ndt and finds that D N is proportional to N{2, then the system is, with very high probability, in the perturbative regime, and one does not know whether p pair will remain close to p true as N increases. What this means in practical terms is that if one wants to know whether a particular pairwise model is a good one for large systems, it is necessary to consider values of N that are significantly greater than N c , where We interpret N c as the value at which there is a crossover in the behavior of the pairwise model. Specifically, if N%N c , the system is in the perturbative regime and the pairwise model is not informative about the large N behavior, whereas if N&N c , the The better the pairwise model, the closerS S pair is to S true . This is reflected in the normalized distance measure, D N , which is the distance between the red and cyan lines divided by the distance between the red and black lines. doi:10.1371/journal.pcbi.1000380.g003 system is in a regime in which it may be possible to make inferences about the behavior of the full system.
The prefactors, g ind and g pair As we show in Methods (see in particular Eqs. (42) and (44)), the prefactors g ind and g pair depend on which neurons out of the full population are used. Consequently, these quantities fluctuate around their true values (in the sense that different subpopulations produce different values of g ind and g pair ), where ''true'' refers to an average over all possible N-neuron sub-populations. Here we assume that the N neurons are chosen randomly from the full population, so any set of N neurons provides unbiased estimates of g ind and g pair . In our simulations, the fluctuations were small, as indicated by the small error bars on the blue points in Fig. 5. However, in general the size of the fluctuations is determined by the range of firing rates and correlation coefficients, with larger ranges producing larger fluctuations.
Because N does not affect the mean values of g ind and g pair , it is reasonable to think of these quantities-or at least their true values-as being independent of N. They are also independent of n, again modulo fluctuations. Finally, as we show in the section ''Bin size and the correlation coefficients'' (Methods), g ind and g pair are independent of dt in the limit that dt is small compared to the width of the temporal correlations among neurons. We will assume this limit applies here. In sum, then, to first approximation, g ind and g pair are independent of our three important quantities: N, n, and dt. Thus, we treat them as effectively constant throughout our analysis.

The dangers of extrapolation
Although the behavior of D N in the perturbative regime does not tell us much about its behavior at large N, it is possible that other quantities that can be calculated in the perturbative regime, g ind , g pair , and S ind (the last one exactly), are informative, as others have suggested [7]. Here we show that this is not the case-they also are uninformative.
The easiest way to relate the perturbative regime to the large N regime is to ignore the corrections in Eqs. (8a) and (8b), extrapolate the expressions for the zeroth order terms, and ask what their large N behavior tells us. Generic versions of these extrapolations, plotted on a log-log scale, are shown in Fig. 4A, along with a plot of the independent entropy, S ind (which is necessarily linear in N).
The first thing we notice about the extrapolations is that they do not, technically, have a large N behavior: one terminates at the point labeled N ind , which is where D 0 KL p true kp ind ð ÞS ind (and thus, via Eq. (0a), S 0 true~0 ; continuing the extrapolation implies negative true zeroth order entropy); the other at the point labeled N pair , which is where D 0 KL p true kp pair À Á S ind (and thus, via Eq. (5) and the fact thatS S 0 pair ƒS ind , S 0 true ƒ0). Despite the fact that the extrapolations end abruptly, they still might provide information about the large N regime. For example, N pair and/or N ind might be values of N at which something interesting happens. To see if this is the case, in Fig. 4B we plot the naive extrapolations ofS S pair and S true (that is, the zeroth order quantities given in Eq. (11),S S 0 pair and S 0 true ), on a linear-linear plot, along with S ind . This plot contains no new information compared to Fig. 4A, but it does elucidate the meaning of the extrapolations. Perhaps its most striking feature is that the naive extrapolation of S true has a decreasing portion. As is easy to show mathematically, entropy cannot decrease with N (intuitively, that is because observing one additional neuron cannot decrease the entropy of previously observed neurons). Thus, N ind , which occurs well beyond the point where the naive extrapolation of S true is decreasing, has essentially no meaning, something that has been pointed out previously by Bethge and Berens [10]. The other potentially important value of N is N pair . This, though, suffers from a similar problem: when N~N pair , S 0 true is negative. How do the naively extrapolated entropies-the solid lines in Fig. 4B-compare to the actual entropies? To answer this, in Fig. 4B we show the true behavior of S true andS S pair versus N (dashed lines). Note that S true is asymptotically linear in N, even though the neurons are correlated, a fact that forcesS S pair to be linear in N, as it is sandwiched between S true and S ind . (The asymptotically linear behavior of S true is typical, even in highly correlated systems. Although this is not always appreciated, it is easy to show; see the section ''The behavior of the true entropy in the large N limit,'' Methods.) Comparing the dashed and solid lines, we see that the naively extrapolated and true entropies, and thus the naively extrapolated and true values of D N , have extremely different behavior. This further suggests that there is very little connection between the perturbative and large N regimes.
In fact, these observations follow directly from the fact that g ind and g pair depend only on correlation coefficients up to third order (see Eqs. (42) and (44)) whereas the large N behavior depends on correlations at all orders. Thus, since g ind and g pair tell us very little, if anything, about higher order correlations, it is not surprising that they tell us very little about the behavior of D N in the large N limit.

Numerical simulations
To check that our perturbation expansions, Eqs. (8)(9)(10), are correct, and to investigate the regime in which they are valid, we performed numerical simulations. We generated, from synthetic data, a set of true distributions, computed the true distance measures, D KL p true kp ind ð Þ , D KL p true kp pair À Á , and D N numerically, and compared them to the zeroth order ones, D 0 KL p true kp ind ð Þ , D 0 KL p true kp pair À Á , and D 0 N . If the perturbation expansion is valid, then the true values should be close to the zeroth order values whenever Nndt is small. The results are shown in Fig. 5, and that is, indeed, what we observed. Before discussing that figure, though, we explain our procedure for constructing true distributions. The set of true distributions we used were generated from a third order model (so named because it includes up to third order interactions). This model has the form where Z true is a normalization constant, chosen to ensure that the probability distribution sums to 1, and the sums over i, j and k run from 1 to N Ã . The parameters h true i ,J true ij and K true ijk were chosen by sampling from distributions (see the section ''Generating synthetic data,'' Methods), which allowed us to generate many different true distributions. In all of our simulations we calculate the relevant quantities directly from Eq. (13) . Consequently, we do not have to worry about issues of finite data, as would be the case in realistic experiments.
For a particular simulation (corresponding to a column in Fig. 5), we generated a true distribution with N Ã~1 5, randomly chose 5 neurons, and marginalized over them. This gave us a 10-neuron true distribution. True distributions with Nv10 were constructed by marginalizing over additional neurons within our 10-neuron population. To achieve a representative sample, we considered all possible marginalizations (of which there are 10 choose N, or ). The results in Fig. 5 are averages over these marginalizations.
For neural data, the most commonly used pairwise model is the maximum entropy model. Therefore, we use that one here. To emphasize the maximum entropy nature of this model, we replace the label ''pair'' that we have been using so far with ''maxent.'' The maximum entropy distribution has the form Fitting this distribution requires that we choose the h i and J ij so that the first and second moments match those of the true distribution. Quantitatively, these conditions are where the angle brackets, S . . .T maxent and S . . .T true , represent averages with respect to p maxent and p true , respectively. Once we have h i and J ij that satisfy Eq. (15), we calculate the KL divergences, Eqs. (1) and (4), and use those to compute D N . The results are shown in Fig. 5. The rows correspond to our three quantities of interest: D KL p true kp ind ð Þ , D KL p true kp pair À Á , and D N (top to bottom). The columns correspond to different values of ndt, with the smallest ndt on the left and the largest on the right. Red circles are the true values of these quantities; blue ones are the zeroth order predictions from Eqs. (9) and (10b).
As suggested by our perturbation analysis, the smaller the value of ndt, the larger the value of N for which agreement between the true and zeroth order values is good. Our simulations corroborate this: the left column of Fig. 5 has ndt~0:024, and agreement is almost perfect out to N~10; the middle column has ndt~0:029, and agreement is almost perfect out to N~7; and the right column has ndt~0:037, and agreement is not good for any value of N. Note that the perturbation expansion breaks down for values of N well below N C (defined in Eq. (12)): in the middle column of Fig. 5 it breaks down when N=N c &0:23, and in the right column it breaks down when N=N c &0:15. This is not, however, especially surprising, as the perturbation expansion is guaranteed to be valid only if N=N c %1.
These results validate the perturbation expansions in Eqs. (8) and (10), and show that those expansions provide sensible predictions-at least for some parameters. They also suggest a natural way to assess the significance of one's data: plot D KL p true kp ind ð Þ , D KL p true kp pair À Á , and D N versus N, and look for agreement with the predictions of the perturbation expansion. If agreement is good, as in the left column of Fig. 5, then one is in the perturbative regime, and one knows very little about the true distribution. If, on the other hand, agreement is bad, as in the right column, then one is out of the perturbative regime, and it may be possible to extract meaningful information about the relationship between the true and pairwise models.
That said, the qualifier ''at least for some parameters'' is an important one. This is because the perturbation expansion is essentially an expansion that depends on the normalized correlation coefficients, and there is an underlying assumption that they don't exhibit pathological behavior. The k th order normalized correlation coefficient for the distribution p r ð Þ, denoted r p i1i2...i k , is written A potentially problematic feature of the correlation coefficients is that the denominator is a product over mean activities. If the mean activities are small, the denominator can become very small, leading to very large correlation coefficients. Although our perturbation expansion is always valid for sufficiently small time bins (because the correlation coefficients eventually becomes independent of bin size; see the section ''Bin size and the correlation coeffcients,'' Methods), ''sufficiently small'' can depend in detail on the parameters. For instance, at the maximum population size tested (N~10) and for the true distributions that had ndtv0:03, the absolute error of the prediction had a median of approximately 16%. However, about 11% of the runs had errors larger than 60%. Thus, the exact size of the small parameter at which the perturbative expansion breaks down can depend on the details of the true distribution.
External fields and pairwise couplings have a simple dependence on firing rates and correlation coefficients in the perturbative regime Estimation of the KL divergences and D N from real data can be hard, in the sense that it takes a large amount of data for them to converge to their true values. In addition, as discussed above, in the section ''The prefactors g ind and g pair '', there are fluctuations in D N associated with finite subsampling of the full population of neurons. Those fluctuations tend to keep D N from being purely linear, as can seen, for example, in the blue points in Fig. 5F and 5I. We therefore provide a second set of relationships that can be used to determine whether or not a particular data set is in the perturbative regime. These relationships are between the parameters of the maximum entropy model, the h i and J ij , and the mean activity and normalized second order correlation coefficient (the latter defined in Eq. (19) below).
Since the quantity ndt plays a central role in our analysis, we replace it with a single parameter, which we denote d, In terms of this parameter, we find (using the same perturbative approach that led us to Eqs. (8)(9)(10); see the section ''External fields, pairwise couplings and moments,'' Methods), that where r ij , the normalized second order correlation coefficient, is defined in Eq. (16) with k~2; it is given explicitly by (We don't need a superscript on r or a subscript on the angle brackets because the first and second moments are the same under the true and pairwise distributions.) Equation (18a) can be reconstructed from the low firing rate limit of analysis carried out by Sessak and Monasson [17], as can the first three terms in the expansion of the log in Eq. (18b). Equation (18) tells us that the N-dependence of the h i and J ij , the external fields and pairwise couplings, is very weak. In Fig. 6 we confirm this through numerical simulations. Equation (18b) also provides additional information-it gives us a functional relationship between the pairwise couplings and the normalized pairwise correlations function, r ij . In Fig. 7A-C we plot the pairwise couplings, J ij , versus the normalized pairwise correlation coefficient, r ij (blue dots), along with the prediction from Eq. (18b) (black line). Consistent with our predictions, the data in Fig. 7A-C essentially follows a line-the line given by Eq. (18b).
A relationship between the pairwise couplings and the correlations coefficients has been sought previously, but for the more standard Pearson correlation coefficient [7,9,11]. Our analysis explains why it was not found. The Pearson correlation coefficient, denoted c ij , is given by In the small Sr i T limit-the limit of interest-the right hand side of Eq. (20) is approximately equal to Sr i TSr j T À Á 1=2 r ij . Because and there is a one-to-one relationship between r ij and J ij (Eq. (18b)), there can't be a one-to-one relationship between c ij and J ij . We verify the lack of a relationship in Fig. 7D and 7E, where we again plot J ij , but this time versus the standard correlation coefficient, c ij . As predicted, the data in Fig. 7D and 7E is scattered over two dimensions. This suggests that r ij , not c ij , is the natural measure of the correlation between two neurons when they have a binary representation, something that has also been suggested by Amari based on information-geometric arguments [18]. Note that the lack of a simple relationship between the pairwise couplings and the standard correlation coefficient has been a major motivation in building maximum entropy models [7,11]. This is for good reason: if there is a simple relationship, knowing the J ij 0 s adds essentially nothing. Thus, plotting J ij versus r ij (but not c ij ) is an important test of one's data, and if the two quantities fall on the curve predicted by Eq. (18b), the maximum entropy model is adding very little information, if any.
As an aside, we should point out that the N-dependence is a function of the variables used to represent the firing patterns. Here we use 0 for no spike and 1 for one or more spikes, but another, possibly more common, representation, derived from the Ising model and used in a number of studies [7,9,11], is to use 21 and +1 rather than 0 and 1. This amounts to making the change of variables s i~2 r i {1. In terms of s i , the maximum entropy model has the form p r The second term on the right side of Eq. (21a) is proportional to N{1, which means the external fields in the Ising representation acquire a linear N-dependence that was not present in our 0/1 representation. The two studies that reported the N-dependence of the external fields [7,9] used this representation, and, as predicted by our analysis, the external fields in those studies had a component that was linear in N.

Is there anything wrong with using small time bins?
An outcome of our perturbative approach is that our normalized distance measure, D N , is linear in bin size (see Eq. (10b)). This suggests that one could make the pairwise model look better and better simply by making the bin size smaller and smaller. Is there anything wrong with this? The answer is yes, for reasons discussed above (see the the section ''The extrapolation problem''); here we emphasize and expand on this issue, as it is an important one for making sense of experimental results.
The problem arises because what we have been calling the ''true'' distribution is not really the true distribution of spike trains. It is the distribution assuming independent time bins, an assumption that becomes worse and worse as we make the bins smaller and smaller. (We use this potentially confusing nomenclature primarily because all studies of neuronal data carried out so far have assumed temporal independence, and compared the pairwise distribution to the temporally independent-but still correlated across neurons-distribution [7][8][9]11]. In addition, the correct name ''true under the assumption of temporal independence,'' is unwieldy.) Here we quantify how much worse. In particular, we show that if one uses time bins that are small compared to the characteristic correlation time in the spike trains, the pairwise model will not provide a good description of the data. Essentially, we show that, when the time bins are too small, the error one makes in ignoring temporal correlations is larger than the error one makes in ignoring correlations across neurons.
As usual, we divide time into bins of size dt. However, because we are dropping the independence assumption, we use r t , rather than r, to denote the response in bin t. The full probability distribution over all time bins is denoted 2 r 1 , . . . ,r M À Á . Here M is the number of bins; it is equal to T=dt for spike trains of length T. If time bins are approximately independent and the distribution of r t is the same for all t (an assumption we make for convenience only, but do not need; see the section ''Extending the normalized distance measure to the time domain,'' Methods), we can write Furthermore, if the pairwise model is a good one, we have Combining Eqs. (22) and Eq. (23) then gives us an especially simple expression for the full probability distribution: 2 r 1 , . . . ,r M À Á & P t p pair r t ð Þ. The problem with small time bins lies in Eq. (22): the right hand side is a good approximation to the true distribution when the time bins are large compared to the spike train correlation time, but it is a bad approximation when the time bins are small (because adjacent time bins become highly correlated). To quantify how bad, we compare the error one makes assuming independence across time to the error one makes assuming independence across neurons. The ratio of those two errors, denoted c, is given by It is relatively easy to compute c in the limit of small time bins (see the section ''Extending the normalized distance measure to the time domain,'' Methods), and we find that As expected, this reduces to our old result, D N , when there is only one time bin (M~1). When M is larger than 1, however, the second term is always at least one, and for small bin size, the third term is much larger than one. Consequently, if we use bins that are small compared to the temporal correlation time of the spike trains, the pairwise model will do a very bad job describing the full, temporally correlated spike trains.

Discussion
Probability distributions over the configurations of biological systems are extremely important quantities. However, because of the large number of interacting elements comprising such systems, these distributions can almost never be determined directly from experimental data. Using parametric models to approximate the true distribution is the only existing alternative. While such models are promising, they are typically applied only to small subsystems, not the full system. This raises the question: are they good models of the full system?
We answered this question for a class of parametric models known as pairwise models. We focused on a particular application, neuronal spike trains, and our main result is as follows: if one were to record spikes from multiple neurons, use sufficiently small time bins and a sufficiently small number of cells, and assume temporal independence, then a pairwise model will almost always succeed in matching the true (but temporally independent) distribution-whether or not it would match the true (but still temporally independent) distribution for large time bins or a large number of cells. In other words, pairwise models in the ''sufficiently small'' regime, what we refer to as the perturbative regime, have almost no predictive value for what will happen with large populations. This makes extrapolation from small to large systems dangerous.
This observation is important because pairwise models, and in particular pairwise maximum entropy models, have recently attracted a great deal of attention: they have been applied to salamander and guinea pig retinas [7], primate retina [8], primate cortex [9], cultured cortical networks [7], and cat visual cortex [11]. These studies have mainly operated close to the perturbative regime. For example, Schneidman et al. [7] had Nndt&0:35 (for the data set described in their Fig. 5), Tang et al. [9] had Nndt&0:06 to 0.4 (depending on the preparation), and Yu et al. [11] had Nndt&0:2. For these studies, then, it would be hard to justify extrapolating to large populations.
The study by Shlens et al. [8], on the other hand, might be more amenable to extrapolation. This is because spatially localized visual patterns were used to stimulate retinal ganglion cells, for which a nearest neighbor maximum entropy models provided a good fit to their data. (Nearest neighbor means J ij is zero unless neuron i and neuron j are adjacent.) Our analysis still applies, but, since all but the nearest neighbor correlations are zero, many of the terms that make up g ind and g pair vanish (see Eqs. (42) and (44)). Consequently, the small parameter in the perturbative expansion becomes Kndt (rather than Nndt), where K is the number of nearest neighbors. Since K is fixed, independent of the population size, the small parameter will not change as the population size increases. Thus, Shlens et al.may have tapped into the large population behavior even though they considered only a few cells at a time in their analysis. Indeed, they have recently extended their analysis to more than 100 neurons, and they still find that nearest neighbor maximum entropy models provide very good fits to the data [19].

Time bins and population size
That the pairwise model is always good if Nndt is sufficiently small has strong implications: if we want to build a good model for a particular N, we can simply choose a bin size that is small compared to 1=Nn. However, one of the assumptions in all pairwise models used on neural data is that bins at different times are independent. This produces a tension between small time bins and temporal independence: small time bins essentially ensure that a pairwise model will provide a close approximation to a model with independent bins, but they make adjacent bins highly correlated. Large time bins come with no such assurance, but they make adjacent bins independent. Unfortunately, this tension is often unresolvable in large populations, in the sense that pairwise models are assured to work only up to populations of size 1= nt corr ð Þwhere t corr is the typical correlation time. Given that n is at least several Hz, for experimental paradigms in which the correlation time is more than a few hundred ms, 1= nt corr ð Þ is about one, and this assurance does not apply to even moderately sized populations of neurons.
These observations are especially relevant for studies that use small time bins to model spike trains driven by natural stimuli. This is because the long correlation times inherent in natural stimuli are passed on to the spike trains, so the assumption of independence across time (which is required for the independence assumption to be valid) breaks badly. Knowing that these models are successful in describing spike trains under the independence assumption, then, does not tell us whether they will be successful in describing full, temporally correlated, spike trains. Note that for studies that use stimuli with short correlation times (e.g., nonnatural stimuli such as white noise), the temporal correlations in the spike trains are likely to be short, and using small time bins may be perfectly valid.
The only study that has investigated the issue of temporal correlations in maximum entropy models does indeed support the above picture [9]: for the parameters used in that study (Nndt~0:06 to 0.4), the pairwise maximum entropy model provided a good fit to the data (D N was typically smaller than 0.1), but it did not do a good job modeling the temporal structure of the spike trains.

Other systems-Protein folding
As mentioned in the Introduction, in addition to the studies on neuronal data, studies on protein folding have also emphasized the role of pairwise interactions [2,3]. Briefly, proteins consist of strings of amino acids, and a major question in structural biology is: what is the probability distribution of amino acid strings in naturally folding proteins? One way to answer this is to approximate the full probability distribution of naturally folding proteins from knowledge of single-site and pairwise distributions. One can show that there is a perturbative regime for proteins as well. This can be readily seen using the celebrated HP protein model [20], where a protein is composed of only two types of amino acids. If, at each site, one amino acid type is preferred and occurs with high probability, say 1{d with d%1, then a protein of length shorter than 1=d will be in the perturbative regime, and, therefore, a good match between the true distribution and the pairwise distribution for such a protein is virtually guaranteed.
Fortunately, the properties of real proteins generally prevent this from happening: at the majority of sites in a protein, the distribution of amino acids is not sharply peaked around one amino acid. Even for those sites that are sharply peaked (the evolutionarily-conserved sites), the probability of the most likely amino acid, 1{d, rarely exceeds 90% [21,22]. This puts proteins consisting of only a few amino acids out of the perturbative regime, and puts longer proteins-the ones usually studied using pairwise models-well out of it.
This difference is fundamental: because many of the studies that have been carried out on neural data were in the perturbative regime, the conclusions of those studies-specifically, the conclusion that pairwise models provide accurate descriptions of large populations of neurons-is not yet supported. This is not the case for the protein studies, because they are not in the perturbative regime. Thus, the evidence that pairwise models provide accurate descriptions of protein folding remain strong and exceedingly promising.

Open questions
In our analysis, we sidestepped two issues of practical importance: finite sampling and alternative measures for assessing the quality of the pairwise model. These issues are beyond the scope of this paper, but in our view, they are natural next steps in the analysis of pairwise models. Below we briefly expand on them.
Finite sampling refers to the fact that in any real experiment, one has access to only a finite amount of data, and so does not know the true probability distribution of the spike trains. In our analysis, however, we assumed that one did have full knowledge of the true probability distribution. Since a good estimate of the probability distribution is crucial for assessing whether the pairwise model can be extrapolated to large populations, it is important to study how such estimates are affected by finite data. Future work is needed to address this issue, and to find ways to overcome data limitation-for example, by finding efficient methods for removing the finite data bias that affects information theoretic quantities such as the Kullback-Leibler divergence.
There are always many possible ways to assess the quality of a model. Our choice of D N was motivated by two considerations: it is based on the Kullback-Leibler divergence, which is a standard measure of ''distance'' between probability distributions, and it is a widely used measure in the field [7][8][9][10]. It suffers, however, from a number of shortcomings. In particular, D N can be small even when the pairwise model assigns very different probabilities to many of the configurations of the system. It would, therefore, be important to study the quality of pairwise models using other measures.

Methods
The behavior of the true entropy in the large N limit To understand how the true entropy behaves in the large N limit, it is useful to express the difference of the entropies as a mutual information. Using S N to denote the true entropy of N neurons and I 1; N ð Þ to denote the mutual information between one neuron and the other N neurons in a population of size Nz1, we have If knowing the activity of N neurons does not fully constrain the firing of neuron Nz1, then the single neuron entropy, S 1 , will exceed the mutual information, I 1; N ð Þ, and the entropy will be an increasing function of N. For the entropy to be linear in N, all we need is that the mutual information saturates with N. Because of synaptic noise, this is a reasonable assumption for networks of neurons: even if we observed all the spikes from all the neurons, there would still be residual noise associated with synaptic failures, jitter in release time, variability in the amount of neurotransmitter released, stochastic channel dynamics, etc. Consequently, in the large N limit, we may replace I 1; N ð Þ by its average, denoted SI 1; ? ð ÞT. Also replacing S 1 by its average, denoted SS 1 T, we see that for large N, the difference between S Nz1 and S N approaches a constant. Specifically, where this expression is valid in the large N limit and the corrections are sublinear in N.

Perturbative expansion
Our main quantitative result, given in Eqs. (8)(9)(10), is that the KL divergence between the true distribution and both the independent and pairwise distributions can be computed perturbatively as an expansion in powers of Nd in the limit Nd%1. Here we carry out this expansion, and derive explicit expressions for the quantities g ind and g pair .
To simplify our notation, here we use p r ð Þ for the true distribution. The critical step in computing the KL divergences perturbatively is to use the Sarmanov-Lancaster expansion [23][24][25][26][27][28] for p r ð Þ, where dr i~ri {r i ð29cÞ This expansion has a number of important, but not immediately obvious, properties. First, as can be shown by a direct calculation, where the subscripts p and ind indicate an average with respect to p r ð Þ and p ind r ð Þ, respectively. This has an immediate corollary, This last relationship is important, because it tells us that if a product of dr 0 s contains any terms linear in one of the dr i , the whole product averages to zero under the independent distribution. This can be used to show that from which it follows that Thus, p r ð Þ is properly normalized. Finally, a slightly more involved calculations provides us with a relationship between the parameters of the model and the moments: for i=j=k, Sdr Virtually identical expressions hold for higher order moments. It is this last set of relationships that make the Sarmanov-Lancaster expansion so useful.
Note that Eqs. (32a) and (32b), along with the expression for the normalized correlation coefficients given in Eq. (16), imply that These identities will be extremely useful for simplifying expressions later on.
Because the moments are so closely related to the parameters of the distribution, moment matching is especially convenient: to construct a distribution whose moments match those of p r ð Þ up to some order, one simply needs to ensure that the parameters of that distribution, H i , J ij , K ijk , etc., are identical to those of the true distributions up to the order of interest. In particular, let us write down a new distribution, q r ð Þ, We can recover the independent distribution by letting j q r ð Þ~0, and we can recover the pairwise distribution by letting J q ij~J p ij . This allows us to compute D KL pkq ð Þ in the general case, and then either set j q to zero or set J q ij to J p ij . Expressions analogous to those in Eqs. (31)(32)(33) exist for averages with respect to q r ð Þ; the only difference is that p is replaced by q. Defining and inserting this expression into Eq. (43), we recover Eq. (8b).

External fields, pairwise couplings and moments
In this section we derive, to leading order in Nd, expressions relating the external fields and pairwise couplings of the maximum entropy model, h i and J ij , to the first and second moments; these are the expressions reported in Eq. (18). The calculation proceeds along the same lines as in the previous section. There is, though, one extra step associated with the fact that the quadratic term in the maximum entropy distribution given in Eq. (14) is proportional to r i r j , not dr i dr j . However, to lowest order in Nd, this doesn't matter. That's because X ivj J ij r i r j~X ivj J ij dr i dr j zr i X j=i J ij r j zconstants: where r i is defined as in Eq. (29d) except with H p i replaced by h i , and we used the fact that J ij~Jji . The second term introduces a correction to the external fields, h i . However, the correction is O Nd ð Þ, so we drop it. We should keep in mind, though, that our final expression for h i will have corrections of this order.
Using Eq. (14), but with r i replaced by dr i where it appears with J ij , we may write (after a small amount of algebra) where p ind r ð Þ is the same as the function p ind r ð Þ defined in Eq. (29a) except that H p i is replaced by h i , the subscript ''ind'' indicates, as usual, an average with respect to p ind r ð Þ, and the two functions j 2 r ð Þ and y x ð Þ are defined by and y x ð Þ:e x {1{x: Given this setup, we can use Eqs. (55) and (56) below to compute the moments under the maximum entropy model. That's because both y x ð Þ and its first derivative vanish at x~0, which are the two conditions required for these equations to be valid. Using also the fact that Sdr i T ind~0 , Eqs. (55) and (56) imply that where the first term in Eq. (48b) came from Eq. (29d) with H p i replaced by h i , the term ''zJ ij '' in Eq. (48c) came from Sdr i dr j j 2 r ð ÞT ind , and for the second two expressions we used the fact that, to lowest order in Nd, the denominator in Eq. (45) is equal to 1.
Finally, using Eq. (19) for the normalized correlation coefficient, dropping the subscript ''maxent'' (since the first and second moments are the same under the maxent and true distributions), and inverting Eqs. (48b) and (48c) to express the external fields and coupling coefficients in terms of the first and second moments, we arrive at Eq (18).

Averages of powers of j p and j q
Here we compute Sj m p j n q T ind , which, as can be seen in Eq. (37), is the key quantity in our perturbation expansion. Our starting point is to (formally) expand the sums that make up j p and j q (see Eqs. (29b) and (34b)), which yields The sum over m 1 , . . . ,m l f g is a sum over all possible configurations of the m i . The coefficient y l ð Þ m1,...,m l are complicated functions of the J p ij ,J q ij ,K p ijk,K q ijk , etc. Computing these functions is straightforward, although somewhat tedious, especially when l is large; below we compute them only for l~2 and 3. The reason l starts at 2 is that mzn §2; see Eq. (37).
We first show that all terms with superscript l ð Þ are O d l À Á . To do this, we note that, because the right hand side of Eq. (49) is an average with respect to the independent distribution, the average of the product is the product of the averages, Then, using the fact that dr i~1 {r i ð Þwith probability r i and {r i with probability 1{r i ð Þ(see Eq. (29c)), we have The significance of this expression is that, for mw1, Consequently, if all the m i in Eq. (50) are greater than 1, then the right hand side is O d l À Á . This shows that, as promised above, the superscript l ð Þ labels the order of the terms.
As we saw in the section ''The KL divergence in the Sarmanov-Lancaster representation'', we need to go to third order in d, which means we need to compute the terms on the right hand side of Eq. (49) with l~2 and 3. Let us start with l~2, which picks out only those terms with two unique indices. Examining the expressions for j p and j q given in Eqs. (29b) and (34b), we see that we must keep only terms involving J ij , since K ijk has three indices, and higher order terms have more. Thus, the l~2 contribution to the average in Eq. (49), which we denote Sj p r ð Þj q r ð ÞT 2 ð Þ ind , is given by Pulling J p ij and J q ij out of the averages, using Eq. (33a) to eliminate J p ij and J q ij in favor of r p ij and r q ij , and applying Eq. (51) (while throwing away some of the terms in that equation that are higher than second order in d), the above expression may be written In the perturbative regime, a distribution generated with these values of the external fields has firing rates approximately equal to the r Ã i ; see Eq. (18a) and Fig. 6. To generate the J true ij and K true ijk , we drew them from Gaussian distributions with means equal to 0.05 and 0.02 and standard deviations of 0.8 and 0.5, respectively. Using non-zero values for K ijk means that the true distribution is not pairwise.

Bin size and the correlation coefficients
One of our main claims is that D N is linear in bin size, dt. This is true, however, only if g ind and g pair are independent of dt, as can be seen from Eq. (10b). In this section we show that independence is satisfied if dt is smaller than the typical correlation time of the responses. For dt larger than such correlation times, g ind and g pair do depend on dt, and D N is no longer linear in dt. Note, though, that the correlation time is always finite, so there will always be a bin size below which the linear relationship, D N *dt, is guaranteed.
Examining Eqs. (42) and (44), we see that g ind and g pair depend on the normalized correlation coefficients, r ij andr r ijk (we drop superscripts, since our discussion will be generic). Thus, to understand how g ind and g pair depend on bin size, we need to understand how the normalized correlation coefficients depend on bin size. To do that, we express them in terms of standard crosscorrelograms, as the cross-correlograms contain, in a very natural way, information about the temporal timescales in the spike train.
We start with the second order correlation coefficient, since it is simplest. The corresponding cross-correlogram, which we denote C ij t ð Þ, is given by where t k i is the time of the k th spike on neuron i (and similarly for t l j ), and d : ð Þ is the Dirac d-function. The normalization in Eq. (57) is slightly non-standard-more typical is to divide by something with units of firing rate (n i , n j or n i n j À Á 1=2 ), to give units of spikes/s.
The normalization we use is convenient, however, because in the limit of large t, C ij t ð Þ approaches one. It is slightly tedious, but otherwise straightforward, to show that when dt is sufficiently small that only one spike can occur in a time bin, r ij is related to C ij t ð Þ via The (unimportant) factor 1{ t j j=dt ð Þcomes from the fact that the spikes occur at random locations within a bin. Equation (58) has a simple interpretation: r ij is the average height of the central peak of the cross-correlogram relative to baseline. How strongly r ij depends on dt is thus determined by the shape of the cross-correlogram. If it is smooth, then r ij approaches a constant as dt becomes small. If, on the other hand, there is a sharp peak at t~0, then r ij *1=ndt~1=d for small dt, so long as dt is larger than the width of the peak. (The factor of n included in the scaling is approximate; it is a placeholder for an effective firing rate that depends on the indices i and j. It is, however, sufficiently accurate for our purposes.) A similar relationship exists between the third order correlogram and the correlation coefficient. Thus, r r ijk is also independent of dt in the small dt limit, whereas if the central peak is sharp it scales as 1 d 2 . The upshot of this analysis is that the shape of the crosscorrelogram has a profound effect on the correlation coefficients and, therefore, on D N . What is the shape in real networks? The answer typically depends on the physical distance between cells. If two neurons are close, they are likely to receive common input and thus exhibit a narrow central peak in their cross-correlogram. Just how narrow depends on the area. Early in the sensory pathways, such as retina [29][30][31] and LGN [32], peaks can be very narrowon the order of milliseconds. Deeper into cortex, however, peaks tend to broaden, to at least tens of milliseconds [33,34]. If, on the other hand, the neurons are far apart, they are less likely to receive common input. In this case, the correlations come from external stimuli, so the central peak tends to have a characteristic width given by the temporal correlation time of the stimulus, typically 100 s of milliseconds.
Although clearly both kinds of cross-correlograms exist in any single population of neurons, it is convenient to analyze them separately. We have already considered networks in which the cross-correlograms were broad and perfectly flat, so that the correlation coefficients were strictly independent of bin size. We can also consider the opposite extreme: networks in which the the cross-correlograms (both second and higher order) among nearby neurons exhibit sharp peaks while those among distant neurons are uniformly equal to 1. In this regime, the correlation coefficients depend on dt: as discussed above, the second order ones scale as 1=d and the third as 1=d 2 . This means that the arguments of Consequently, f r ij ,0 À Á *1=d and fr r true ijk ,r r pair ijk *1 d 2 . This implies that g ind and g pair scale as Nd and N 2 d, respectively, and so D N *N, independent of d. Thus, if the bin size is large compared to the correlation time, D N will be approximately independent of bin size.
Extending the normalized distance measure to the time domain In this section we derive the expression for c given in Eq. (25). Our starting point is its definition, Eq. (24). It is convenient to define R to be a concatenation of the responses in M time bins, R: r 1 ,r 2 , . . . ,r M À Á ð59Þ where, as in the section ''Is there anything wrong with using small time bins?'', the superscript labels time, so 2 R ð Þ is the full, temporally correlated, distribution.
With this definition, we may write the numerator in Eq. (24) as where S M true is the entropy of 2 R ð Þ, the last sum follows from a marginalization over all but one element of 2 R ð Þ, and p t true r ð Þ is the true distribution at time r (unlike in the section ''Is there