The activity of a neural network is defined by patterns of spiking and silence from the individual neurons. Because spikes are (relatively) sparse, patterns of activity with increasing numbers of spikes are less probable, but, with more spikes, the number of possible patterns increases. This tradeoff between probability and numerosity is mathematically equivalent to the relationship between entropy and energy in statistical physics. We construct this relationship for populations of up to N = 160 neurons in a small patch of the vertebrate retina, using a combination of direct and model-based analyses of experiments on the response of this network to naturalistic movies. We see signs of a thermodynamic limit, where the entropy per neuron approaches a smooth function of the energy per neuron as N increases. The form of this function corresponds to the distribution of activity being poised near an unusual kind of critical point. We suggest further tests of criticality, and give a brief discussion of its functional significance. entropy | information | neural networks | Monte Carlo | correlation O ur perception of the world seems a coherent whole, yet it is built out of the activities of thousands or even millions of neurons, and the same is true for our memories, thoughts, and actions. It is difficult to understand the emergence of behavioral and phenomenal coherence unless the underlying neural activity also is coherent. Put simply, the activity of a brain-or even a small region of a brain devoted to a particular task-cannot be just the summed activity of many independent neurons. How do we describe this collective activity? Statistical mechanics provides a language for connecting the interactions among microscopic degrees of freedom to the macroscopic behavior of matter. It provides a quantitative theory of how a rigid solid emerges from the interactions between atoms, how a magnet emerges from the interactions between electron spins, and so on (1, 2). These are all collective phenomena: There is no sense in which a small cluster of molecules is solid or liquid; rather, solid and liquid are statements about the joint behaviors of many, many molecules. At the core of equilibrium statistical mechanics is the Boltzmann distribution, which describes the probability of finding a system in any one of its possible microscopic states. As we consider systems with larger and larger numbers of degrees of freedom, this probabilistic description converges onto a deterministic, thermodynamic description. In the emergence of thermodynamics from statistical mechanics, many microscopic details are lost, and many systems that differ in their microscopic constituents nonetheless exhibit quantitatively similar thermodynamic behavior. Perhaps the oldest example of this idea is the "law of corresponding states" (3). The power of statistical mechanics to describe collective, emergent phenomena in the inanimate world led many people to hope that it might also provide a natural language for describing networks of neurons (4-6). However, if one takes the language of statistical mechanics seriously, then as we consider networks with larger and larger numbers of neurons, we should see the emergence of something like thermodynamics. At first sight, the notion of a thermodynamics for neural networks seems hopeless. Thermodynamics is about temperature and heat, both of which are irrelevant to the dynamics of these complex, nonequilibrium systems. However, all of the thermodynamic variables that we can measure in an equilibrium system can be calculated from the Boltzmann distribution, and hence statements about thermodynamics are equivalent to statements about this underlying probability distribution. It is then only a small jump to realize that all probability distributions over N variables can have an associated thermodynamics in the N → ∞ limit. This link between probability and thermodynamics is well-studied by mathematical physicists (7), and has been a useful guide to the analysis of experiments on dynamical systems (8, 9). To be concrete, consider a system with N elements; each element is described by a state σ i , and the state of the entire system is σ ≡ fσ 1 , σ 2 , ⋯, σ N g. We are interested in the probability PðσÞ that we will find the system in any one of its possible states. It is natural to think not about the probability itself but about its logarithm, EðσÞ = −ln PðσÞ. [1] In an equilibrium system, this is precisely the energy of each state (in units of k B T), but we can define this energy for any probability distribution. As discussed in detail in Supporting Information, all The activity of a brain-or even a small region of a brain devoted to a particular task-cannot be just the summed activity of many independent neurons. Here we use methods from statistical physics to describe the collective activity in the retina as it responds to complex inputs such as those encountered in the natural environment. We find that the distribution of messages that the retina sends to the brain is very special, mathematically equivalent to the behavior of a material near a critical point in its phase diagram. If we imagine a family of systems in which the number of degrees of freedom N varies, then a thermodynamic limit will exist provided that both the entropy and the energy are proportional to N at large N. The existence of this limit is by no means guaranteed. In most systems, including the networks that we study here, there are few states with high probability, and many more states with low probability. At large N, the competition between decreasing probability and increasing numerosity picks out a special value of E = E * , which is the energy of the typical states that we actually see; E * is the solution to For most systems, the energy EðσÞ has only small fluctuations around E * , hðδEÞ 2 i=ðE * Þ 2 ≈ 1=N, and, in this sense, most of the states that we see have the same value of log probability per degree of freedom. However, hidden in the function SðEÞ are all of the parameters describing the interactions among the N degrees of freedom in the system. At special values of these parameters, ½d 2 SðEÞ=dE 2 E=E * → 0, and the variance of E diverges as N becomes large. This is a critical point, and it is mathematically equivalent to the divergence of the specific heat in an equilibrium system (10). These observations focus our attention on the "density of states" N ðEÞ. Rather than asking how often we see specific combinations of spiking and silence in the network, we ask how many states there are with a particular probability. The vertebrate retina offers a unique system in which the activity of most of the neurons comprising a local circuit can be monitored simultaneously using multielectrode array recordings. As described more fully in ref. 11, we stimulated salamander retina with naturalistic grayscale movies of fish swimming in a tank (Fig. 1A) , while recording from 100 to 200 retinal ganglion cells (RGCs); additional experiments used artificial stimulus ensembles, as described in Supporting Information. Sorting the raw data (12), we identified spikes from 160 neurons whose activity passed our quality checks and was stable for the whole ∼ 2 h duration of the experiment; a segment of the data is shown in Fig. 1B . These experiments monitored a substantial fraction of the RGCs in the area of the retina from which we record, capturing the behavior of an almost complete local population responsible for encoding a small patch of the visual world. The experiment collected a total of ∼ 2 × 10 6 spikes, and time was discretized in bins of duration Δτ = 20 ms; all of the results discussed below are substantially the same at Δτ = 10 ms and Δτ = 40 ms (Fig. S1 ). For each neuron i, σ i = 1 in a bin denotes that the neuron emitted at least one spike, and σ i = 0 denotes that it was silent. Conceptually, estimating the function N ðEÞ and hence the entropy vs. energy is easy: We count how often each state occurs, thus estimating its probability, and then count how many states have (log) probabilities in a given range. In Fig. 1 C and D, we show the first steps in this process. We identify the unique patterns of activitycombinations of spiking and silence across all 160 neurons-that occur in the experiment, and then count how many times each of these patterns occurs. Even without trying to compute SðEÞ, the results of Fig. 1D are surprising. With N neurons that can either spike or remain silent, there are 2 N possible states. Not all these states can be visited equally often, because spikes are less common than silences, but even taking account of this bias, and trying to capture the correlations among neurons, our best estimate of the entropy for the patterns of activity we observe is s ≈ 0.15 bits=neuron (see below). With N = 160 cells, this means that the patterns of activity are spread over 2 Ns ≈ 1.67 × 10 7 possibilities, 100 times larger than the number of samples that we collected during our experiment. Indeed, most of the states that we saw in the full population occurred only once. However, roughly one thousand states occurred with sufficient frequency that we can make a reasonable estimate of their probability just by counting across ∼2 h. Thus, the probability distribution PðσÞ is extremely inhomogeneous. To probe more deeply into the tail of low-probability events, we can construct models of the distribution of states, and we have done this using the maximum entropy method (11, 13): We take from experiment certain average behaviors of the network, and then search for models that match these data but otherwise have as little structure as possible. This works if matching a relatively small number of features produces a model that predicts many other aspects of the data. The maximum entropy approach to networks of neurons has been explored, in several different systems, for nearly a decade (14-23), and there have been parallel efforts to use this approach in other biological contexts (24) (25) (26) (27) (28) (29) (30) (31)(32)(33)(34)(35). Recently, we have used the maximum entropy method to build models for the activity of up to N = 120 neurons in the experiments described above (11); see Fig. 2 . We take from experiment the mean probability of each neuron generating a spike (hσ i i), the correlations between spiking in pairs of neurons (hσ i σ j i), and the probability that K out of the N neurons spike in the same small window of time [PðKÞ]. Mathematically, the maximum entropy models consistent with these data have the form where K = P N i=1 σ i counts the number of neurons that spike simultaneously, and Z is set to ensure normalization. All of the parameters fh i , J ij , V ðKÞg are determined by the measured averages This model accurately predicts the correlations among triplets of neurons (figure 7 in ref. 11), and how the probability of spiking in individual neurons depends on activity in the rest of the population (figure 9 in ref. 11). One can even predict the time-dependent response of single cells from the behavior of the population, without reference to the visual stimulus ( figure 15 in ref. 11 ). Most important for our present discussion, the distribution of the energy Eðfσ i gÞ across the observed patterns of activity agrees with the distribution predicted by the model, deep into the tail of patterns that occur only once in the 2-h-long experiment ( figure 8 in ref. 11 ). This distribution is closely related to the plot of entropy vs. energy that we would like to construct, and so the agreement with experiment gives us confidence. The direct counting of states ( Fig. 1 ) and the maximum entropy models (Fig. 2) give us two complementary ways of estimating the function N ðEÞ and hence the entropy vs. energy in the same data set. Results are in Fig. 3 ; see also Supporting Information. As emphasized above, the plot of entropy vs. energy contains all of the thermodynamic behavior of a system, and this has a meaning for any probability distribution, even if we are not considering a system at thermal equilibrium. Thus, Fig. 3 is as close as we can get to constructing the thermodynamics of this network. With the direct counting of states, we see less and less of the plot at larger N, but the part we can see is approaching a limit as N → ∞, and this is confirmed by the results from the maximum entropy models. This, by itself, is a significant result. If we write down a model like Eq. 4, then, in a purely theoretical discussion, we can scale the couplings between neurons J ij with N to guarantee the existence of a thermodynamic limit (5), but with J ij constructed from real data, we can't impose this scaling ourselves-either it emerges from the data or it doesn't. We can make the emergence of the thermodynamic limit more precise by noting that, at a fixed value of S=N, the value of E=N extrapolates to a well-defined limit in a plot vs. 1=N, as in Fig. 3A , Inset. The results of this extrapolation are strikingly simple: The entropy is equal to the energy, within (small) error bars. If the plot of entropy vs. energy is a straight line with unit slope, then Eq. 2 is solved not by one single value of E but by a whole range. Not only do we have d 2 S=dE 2 = 0, as at an ordinary critical point, but all higher-order derivatives also are zero. Thus, the results in Fig. 3 suggest that the joint distribution of activity across neurons in this network is poised at a very unusual critical point. We expect that states of lower probability (e.g., those in which more cells spike) are more numerous (because there are more ways to arrange K spikes among N cells as K increases from very low values). However, the usual result is that this trade-offwhich is precisely the trade-off between energy and entropy in thermodynamics-selects typical states that all have roughly the same probability. The statement that SðEÞ = E, as suggested in Fig. 3 , is the statement that states which are ten times less probable are exactly 10 times more numerous, and so there is no typical value of the probability. The vanishing of d 2 S=dE 2 corresponds, in an equilibrium system, to the divergence of the specific heat. Although the neurons obviously are not an equilibrium system, the model in Eqs. 3 and 4 is mathematically identical to the Boltzmann distribution. Thus, we can take this model seriously as a statistical mechanics problem, and compute the specific heat in the usual way. Further, we can change the effective temperature by considering a one-parameter family of models, with Eðfσ i gÞ as before (Eq. 4). Changing T is just a way of probing one direction in the parameter space of possible models, and is not a physical temperature; the goal is to see whether there is anything special about the model (at T = 1) that describes the real system. Results for the heat capacity of our model vs. T are shown in Fig. 4 . There is a dramatic peak, and, as we look at larger groups of neurons, the peak grows and moves closer to T = 1, which is the model of the actual network. Importantly, the heat capacity grows even when we normalize by N, so that the specific heat, or heat capacity per neuron, is growing with N, as expected at a critical point, and these signatures are clearer in models that provide a more accurate description of the population activity; for details, see Supporting Information. The temperature is only one axis in parameter space, and, along this direction, there are variations in both the correlations among neurons and their mean spike rates. As an alternative, we consider a family of models in which the strength of correlations changes but spike rates are fixed. To do this, we replace the energy function in Eq. 4 with where α controls the strength of correlations, and we adjust all of the h ′ i ðαÞ to hold mean spike rates to their observed values. At α = 0, the model describes a population of independent neurons, so that the correlation coefficients are all C = 0 (Fig. 5A, Left) . For α > 1, the distribution of correlations broadens such that at α = 2, some pairs are very strongly correlated (Fig. 5A, Right) . This is reflected in a distribution of states that cluster around a small number of prototypical patterns, much as in the Hopfield network (4, 5). The entropy vs. energy plot, shown in Fig. 5B , singles out the ensemble at α = 1 (Fig. 5A , Middle): Going toward independence (smaller α) gives rise to a concave bump at low energies, whereas α > 1 ensembles deviate away from the equality line more at high energies. Correspondingly, we see in Fig. 5D that there is a peak in the specific heat of the model ensemble near α = 1. As we look at larger and larger networks, this peak rises and moves toward α = 1, which describes the real system. The evidence for criticality that we find here is consistent with extrapolations from the analysis of smaller populations (15, 18). Those predictions were based on the assumption that spike probabilities and pairwise correlations in the smaller populations are drawn from the same distribution as in the full system, and that these distributions are sufficient to determine the thermodynamic behavior (36). Signs of criticality also are observable in simpler models, which match only the distribution of summed activity in the network (22), but less accurate models have weaker signatures of criticality (Fig. S2) . In equilibrium thermodynamics, poising a system at a critical point involves careful adjustment of temperature and other parameters. Finding that the retina seems to be poised near criticality should thus be treated with some skepticism. Here we consider some ways in which we could be misled into thinking that the system is critical when it is not (see also Supporting Information). Part of our analysis is based on the use of maximum entropy models, and one could worry that the inference of these models is unreliable for finite data sets (37, 38). Expanding on the discussion of this problem in ref. 11, we find a clear peak in the specific heat when we learn models for N = 100 neurons from even one-tenth of our data, and the variance across fractions of the data are only a few percent (Fig. S3) . Although the inference of maximum entropy models is accurate, less interesting models might mimic the signatures of criticality. In particular, it has been suggested that independent neurons with a broad distribution of spike rates could generate a distribution of N neuron activity patterns fσ i g that mimics some aspects of critical behavior (39). However, in an independent model built from the actual spike rates of the neurons, the probability of seeing the same state twice would be less than one part in a billion, dramatically inconsistent with the measured P c ≈ 0.04. Such independent models also cannot account for the faster than linear growth of the heat capacity with N (Fig. 4) , which is an essential feature of the data and its support for criticality. In maximum entropy models, the probability distribution over patterns of neural activity is described in terms of interactions between neurons, such as the terms J ij in Eq. 4; an alternative view is that the correlations result from the response of the neurons to fluctuating external signals. Testing this idea has a difficulty that has nothing to do with neurons: In equilibrium statistical mechanics, models in which spins (or other degrees of freedom) interact with one another are mathematically equivalent to a collection of spins responding independently to fluctuating fields (see Supporting Information for details). Thus, correlations always are interpretable as independent responses to unmeasured fluctuations, and, for neurons, there are many possibilities, including sensory inputs. However, the behavior that we see cannot be simply inherited from correlations in the visual stimulus, because we find signatures of criticality in response to movies with very different correlation structures (Fig. S4) . Further, the pattern of correlations among neurons is not simply explained in terms of overlaps among receptive fields (Fig. S5) , and, at fixed moments in the stimulus movie, neurons with nonzero spike probabilities have correlations across stimulus repetitions that can be even stronger than across the experiment as a whole (Fig. S6) . When we rewrite a model of interacting spins as independent spins responding to fluctuating fields, criticality usually requires that the distribution of fluctuations be very special, e.g., with the variance tuned to a particular value. In this sense, saying that correlations result from fluctuating inputs doesn't explain our observations. Recently, it has been suggested that sufficiently broad distributions of fluctuations lead generically to critical phenomenology (40). As explained in Supporting Information, mean field models have the property that the variance of the effective fields becomes large at the critical point, but more general models do not, and the correlations we observe do not have the form expected from a mean field model. The fact that quantitative changes in the strength of correlations would drive the system away from criticality ( Fig. 5D) suggests that the distribution of equivalent fluctuating fields must be tuned, rather than merely having sufficiently large fluctuations. The traditional formulation of the neural coding problem makes an analogy to a dictionary, asking for the meaning of each neural response in terms of events in the outside world (41). However, before we can build a dictionary, we need to know the lexicon, and, for large populations of neurons, this already is a difficult problem: With 160 neurons, the number of possible responses is larger than the number of words in the vocabulary of a well-educated English speaker, and is more comparable to the number of possible short phrases or sentences. In the same way that the distribution of letters in words embodies spelling rules (28), and the distribution of words in sentences encodes aspects of grammar (42) and semantic categories (43), we expect the distribution of activity across neurons to reveal structures of biological significance. In the small patch of the retina that we consider, no two cells have truly identical input/output characteristics (44). Nonetheless, if we count how many combinations of spiking and silence have a given probability in groups of N > 20 cells, this relationship is reproducible from group to group, and simplifies at larger N. This relationship between probability and numerosity of states is mathematically identical to the relationship between energy and entropy in statistical physics, and the simplification with increasing N suggests that we are seeing signs of a thermodynamic limit. If we can identify the thermodynamic limit, we can try to place the network in a phase diagram of possible networks. Critical surfaces that separate different phases often are associated with a balance between probability and numerosity: States that are a factor F times less probable also are a factor F times more numerous. At conventional critical points, this balance occurs only in a small neighborhood of the typical probability, but, in the network of RGCs, it extends across a wide range of probabilities (Fig.  3) . In model networks with slightly stronger or weaker correlations among pairs of neurons, this balance breaks down (Fig. 5) , and less accurate models have weaker signatures of criticality (Fig. S2) . The strength of correlations depends on the structure of visual inputs, on the connectivity of the neural circuit, and on the state of adaptation in the system. The fact that we see signatures of criticality in response to very different movies, but not in model networks with stronger or weaker correlations, suggests that adaptation is tuning the system toward criticality. A sudden change of visual input statistics should thus drive the network to a noncritical state, and, during the course of adaptation, the distribution of activity should relax back to the critical surface. This can be tested directly. Is criticality functional? The extreme inhomogeneity of the probability distribution over states makes it possible to have an instantaneously readable code for events that have a large dynamic range of likelihoods or surprise, and this may be wellsuited to the the natural environment; it is not, however, an efficient code in the usual sense. Systems near critical points are maximally responsive to certain external signals, and this sensitivity also may be functionally useful. Most of the systems that exhibit criticality in the thermodynamic sense also exhibit a wide range of time scales in their dynamics, so that criticality may provide a general strategy for neural systems to bridge the gap between the microscopic time scale of spikes and the macroscopic time scales of behavior. Critical states are extremal in all these different senses, and more; it may be difficult to decide which is relevant for the organism. Related signatures of criticality have been detected in ensembles of amino acid sequences for protein families (29), in flocks of birds (33) and swarms of insects (45), and in the network of genes controlling morphogenesis in the early fly embryo (46); there is also evidence that cell membranes have lipid compositions tuned to a true thermodynamic critical point (47). Different, dynamical notions of criticality have been explored in neural (48, 49) and genetic (50, 51) networks, and in the active mechanics of the inner ear (52-54); recent work connects dynamical and statistical criticality, with the retina as an example (55). These results hint at a general principle, but there is room for skepticism. A new generation of experiments should provide decisive tests of these ideas. Experiments were performed on the larval tiger salamander, Ambystoma tigrinum tigrinum, in accordance with institutional animal care standards at Princeton University. ACKNOWLEDGMENTS. We thank A. Cavagna, D. S. Fisher, I. Giardina, M. Ioffe,
REF could be very relevant.
Thermodynamics and signatures of criticality in a network of neurons
Proceedings of the National Academy of Sciences of the United States of America
Proceedings of the National Academy of Sciences of the United States of America
Mathematics Medicine