Statistical models of neural activity, criticality, and Zipf's law

TL;DR

This paper explores the relationship between criticality in neural dynamics, statistical models like maximum entropy, and Zipf's law, highlighting differences between statistical and dynamical criticality and their implications for neural activity analysis.

Contribution

It clarifies the distinction between statistical and dynamical criticality in neural systems and examines the emergence of Zipf's law and Fisher information as indicators of criticality.

Findings

Zipf's law observed in retinal activity under various conditions

Divergence of specific heat as a signature of statistical criticality

Examples of systems critical in one sense but not the other

Abstract

In this overview, we discuss the connections between the observations of critical dynamics in neuronal networks and the maximum entropy models that are often used as statistical models of neural activity, focusing in particular on the relation between "statistical" and "dynamical" criticality. We present examples of systems that are critical in one way, but not in the other, exemplifying thus the difference of the two concepts. We then discuss the emergence of Zipf laws in neural activity, verifying their presence in retinal activity under a number of different conditions. In the second part of the chapter we review connections between statistical criticality and the structure of the parameter space, as described by Fisher information. We note that the model-based signature of criticality, namely the divergence of specific heat, emerges independently of the dataset studied; we suggest this is compatible with previous theoretical findings.

Figures

Equations

$${\sigma }_n(t)=\{\begin{array}{ll}1& \text{ if neuron n spikes between t and t+δt }\\ 0& \text{ otherwise.}\end{array}$$

$$\mathbf{\sigma }(t)=({\sigma }_1(t),\dots ,{\sigma }_N(t))$$

$$P(\mathbf{\sigma })=\frac{1}{Z}{e}^{-E(\mathbf{\sigma })},$$

$$P(\mathbf{\sigma })=\frac{1}{Z(h,J)}\exp \left(\sum _{i=1}^N{h}_i{\sigma }_i+\sum _{i\ne j}{J}_{ij}{\sigma }_i{\sigma }_j\right).$$

$$P(\mathbf{\sigma })=\frac{1}{Z(h,J)}\exp \left(\sum _{i=1}^N{h}_i{\sigma }_i+\sum _{i,j}{J}_{ij}^{(2)}{\sigma }_i{\sigma }_j+\sum _{i,j,k}{J}_{ijk}^{(3)}{\sigma }_i{\sigma }_j{\sigma }_k+\dots \right)$$

$$P_T(\mathbf{\sigma })=\frac{1}{Z(h/T,J/T)}\exp \left(\frac{1}{T}\left[\sum _{i=1}^N{h}_i{\sigma }_i+\sum _{i\ne j}{J}_{ij}{\sigma }_i{\sigma }_j\right]\right)$$

$$P_{T\to \infty }(\mathbf{\sigma })=\frac{1}{2^N},\forall \mathbf{\sigma }.$$

$$E_k=-\log p_k+\text{const.}$$

$$\log p_k$$

$$E_k$$

$$P_T(x)=\frac{1}{Z(T)}{e}^{-E(x)/T}.$$

$$E(\mathbf{v},\mathbf{h})=-{\mathbf{a}}^{⊺}\mathbf{v}-{\mathbf{b}}^{⊺}\mathbf{h}-{\mathbf{v}}^{⊺}J\mathbf{h}.$$

$$P_T(\mathbf{v})=\frac{1}{Z}\sum _{h_{\mathrm{1...}N}=0,1}{e}^{\frac{1}{T}({\mathbf{a}}^{⊺}\mathbf{v}+{\mathbf{b}}^{⊺}\mathbf{h}+{\mathbf{v}}^{⊺}J\mathbf{h})}=\frac{e^{{\mathbf{a}}^{⊺}\mathbf{v}/T}}{Z}\prod _{j=1}^{N_h}\left(1+e^{b_j/T+({\mathbf{v}}^{⊺}J{)}_j/T}\right)$$

$$E_T(\mathbf{v})=\frac{{\mathbf{a}}^{⊺}\mathbf{v}}{T}+\sum _{j=1}^{N_h}\log \left(1+e^{b_j/T+({\mathbf{v}}^{⊺}J{)}_j/T}\right)$$

$$c(T)=\frac{\mathrm{Var}(E_T)}{N_v{T}^2}.$$

$$D_{KL}(p;q)={\int }_{S}p(x)\log \frac{p(x)}{q(x)}dx.$$

$$D_{KL}(P_{\theta +\delta \theta };P_{\theta })\approx \frac{1}{2}\delta {\theta }^{T}F(\theta )\delta \theta ,$$

$$F_{ij}(\theta )=-\int P_{\theta }(x)\frac{\partial \log P_{\theta }(x)}{\partial {\theta }_i}\frac{\partial \log P_{\theta }(x)}{\partial {\theta }_j}dx.$$

$$P_{\theta }(x)=\frac{e^{-H_{\theta }(x)}}{Z_{\theta }},H_{\theta }(x)=\sum _{k=1}^n{\theta }_k{f}_k(x),$$

$$F_{ij}(\theta )=\mathrm{Cov}[f_i(x),f_j(x)].$$

$$\frac{\partial }{\partial \beta }=\frac{1}{n}\sum _{k=1}^n\frac{\partial }{\partial {\theta }_k}$$

$$c(\beta )=\frac{{\beta }^2}{N}\mathrm{Var}[E].$$

$$D_{KL}(P_{{\theta }_0+\alpha v_k};P_{{\theta }_0})\approx \frac{{\alpha }^2}{2}{v}_k^{T}F({\theta }_0)v_k=\frac{{\alpha }^2}{2}{\lambda }_k,$$

Full text

2018

1 Introduction

The debate about criticality in neural systems began with the observation of power laws in a number of experimentally measured variables related to neural activity. The first experimental observation of neuronal avalanches beggs2003neuronal found that their size distribution follows a power law with exponent of about $-3/2$, and their duration distribution follows one of exponent near $-2$ in cortical slices. These values are compatible with the exponents expected in critical branching processes — a well-studied topic in the field of complex systems physics Athreyab . Similar observations have been consistently reported in literature; moreover, the presence of power-law avalanche statistics was found to be theoretically justified by functional arguments on numerous occasions beggs2008criticality ; shew2009neuronal ; shew2013functional ; gautam2015maximizing ; vazquez2017stochastic , and was shown to differ in different brain states priesemann2013neuronal ; hahn2017spontaneous . For an excellent high-level discussion of the topic, see beggs2012being .

An equally interesting instance of a power law is the finding that, in the population statistics of a neural network’s activity, the rank of a state (the first being the most frequently observed, and so on) and its frequency are inversely proportional. This phenomenon, known as Zipf’s law, was first observed by Auerbach in 1913: “If one sorts individuals by a given property in a descending fashion and stops doing so at rank $n_{1}$, or at $n_{2}$, or generally at rank $n_{x}$, where the property has gone down to values $p_{1}$, $p_{2}$, $p_{x}$, then a certain law exists between $n_{x}$ and $p_{x}$. In our case, this law is especially simple, it is expressed by the formula: $n_{x}\cdot p_{x}=\rm{constant}$” Auerbach1913 . This author already alluded to the possibility of more complex forms of the same law (e.g. for the distribution of wealth), but did not speculate why it assumes its simplest form in the studied data. In the mid-1930s, the American linguist George Kingsley Zipf discovered that the frequency of occurrence of words in Joyce’s Ulysses and American newspapers follows the same law Zipf1949 , which is today called Zipf’s law: the frequency of each word decays as a power law of its frequency rank. After Auerbach’s original example, city sizes gabaix1999zipf ; jiang2011zipf , Zipf’s law was confirmed in a variety of fields, including citation counts in scientific literature redner1998popular , earthquake magnitudes, wealth, solar flare size, number of emails and phone calls received, and many others newman2005power .

Over the years several attempt have been made to understand Zipf’s law. Zipf himself explains it by the principle of least effort: If words are stored in a linear array, then the low-frequency items are optimally located in a more distant place than more often used ones. The product of distance and frequency can be considered as a measure of the effort necessary to retrieve the word which he claims to be a constant. However, the assumption of an array, where the effort needed for retrieving an item is linear, which is necessary in order to obtain an inverse relationship rather than a general power law, seems unnatural when considering how items are stored in a neural network.

Another potential cause for the law can be seen in the idea of preferential attachment Barabasi1999 . If, for example, the probability to move to or away from a city is assumed to be independent of its size, then Zipf’s law for city sizes emerges. Other assumptions have been discussed and been shown to provide a better match for the distribution of city sizes Vitanov2015 , but again this may not easily carry over to states in a neural network.

Li Li1992 demonstrated that the words of an artificial language that simply consists of randomly chosen letters including a space sign tend to obey Zipf’s law. However, the ‘space’ sign, which separates words in Li’s approach, plays no such role in the analysis of neural data.

The authors of Ref. aitchison2016zipf aim at an explanation of Zipf’s law by the existence of latent variables. Differently from the above attempts, this study is directly relevant for the analysis of neural activity. It also subsumes the scheme proposed by Li Li1992 .

Although all of these attempts have their interest, there is some agreement that a deeper understanding is still lacking. In addition, there seems to be no clear justification on why criticality in the statistical sense and Zipf’s law have been observed in neural data, or what brain function might benefit from it. It is interesting in this context that Zipf’s law is a system property, i.e. it depends on the number of elements in the system and does not automatically apply to subset or unions of Zipfian sets. It can not be reduced to the mere presence of a particular probability distribution (such as $P(x)\sim x^{-2}$), but requires a conditional sampling procedure to be reproduced in a simulation Cristelli2012 . The observations in neural data as well as a number of unsolved problems with this subject make it a very interesting subject of further investigation.

In what follows, we will discuss the connections between the observations of critical dynamics and maximum entropy models that are often used as statistical models of neural activity, reviewing the recent literature on the matter, and debate the possible relationship between this statistical concept of criticality and the dynamical criticality related to avalanche statistics. First, we will illustrate the concept of Zipf’s distribution, its origin, and its applicability to neural data. In section 2, we will introduce maximum entropy models, and show their connection to criticality and Zipf’s law. In section 3, we will make three observations that emphasise the difference between statistical and dynamical criticality: (3.1) a system that shows dynamical, but not statistical criticality, (3.2) the process of fitting an energy-based model, and (3.3) the application of the theory of a large-scale corpus of biological data, where Zipf’s law appears to hold, although the system is not dynamically critical. In section 4, we will show connections between statistical criticality and the structure of the parameter space, as described by Fisher information. Finally, in Section 5, we will return to the question debate whether there is a relationship between statistical and dynamical criticality and conclude with an outlook on the problem.

2 Statistical description of spike trains

2.1 Zipf’s law in neural data

For the specific case of neural activity, Zipf’s law refers to the rank-probability law for the occurrence of each possible pattern of activity, which has been observed to follow a power law in the same sense as for words in the English language tkacik2015thermodynamics . To understand what we mean by pattern or state, we need to adopt a simplified way of representing spike trains that we can call digital: discretising time in bins of equal size $\delta t$, we can define a Boolean variable

[TABLE]

At any given time, then, the population activity is described by a codeword

[TABLE]

which describes, up to a precision of $\delta t$, the spiking state of the $N$ neurons considered (Figure 1). Modelling the system statistically, in this framework, means giving a full account of the probability of each possible codeword to appear. Note that, typically, we are not concerned with the dynamics of the system, and we disregard temporal correlations on scales larger than $\delta t$: this approach is suited to describe short-time correlations across space or properties of the encoding. Needless to say, the choice of $\delta t$ can have important consequences on the results: in the limit of very large bin size, the pattern where all neurons fire simultaneously will be the only one to be observed; in the opposite limit of small $\delta t$, the silent pattern will be the most common, patterns with a single active neuron arbitrarily rare, and multi-neuron patterns absent. The results we discuss hold for bin sizes of the order of 5–20 ms, i.e. of the same order of magnitude as the typical correlation length between neurons; this value is commonly adopted in the literature schneidman2006weak .

To understand why Zipf laws are considered a signature of criticality, we will now illustrate the relationship, exposed by recent literature, between them and the critical points of models that have been used to describe neural activity, and are well known in physics.

2.2 Statistical modelling

The activity patterns of individual neurons and neural networks invariably display stochastic characteristics. A common approach, which we can call top-down, of modelling the nature of this activity is to make (simplifying) assumptions on the actual workings of neurons, synapses, and networks, in order to set up a computational model the results of which can then be compared with experimental observations. A large part, perhaps the largest, of computational neuroscience is based on this paradigm, predominantly by simulations of spiking neural networks.

Here, on the contrary, we are concerned with what we call bottom-up modelling, which seeks to infer properties of neural activity in an entirely data-driven way. Understanding the correlation structure, the distribution of firing rates, or the repetition of identical patterns from experimental data are examples of this approach. In other words, the data is described in terms of probabilities and other statistical descriptors, instead of parameters directly implied by the biological or physical theory.

In the bottom-up approach, a very broad family of models is available. We will restrict ourselves to energy-based statistical models, a number of models developed in the last decade which adopt a log-linear relation between probability and state variables. In an energy-based model probabilities are expressed in terms of an energy function $E$, in analogy with statistical physics:

[TABLE]

where $Z$ is the relevant normalisation factor. Many energy-based models used in neuroscience adopt the aforementioned digital description of spike trains in terms of binary variables: we will focus on these. In this case, for $N$ neurons, $\boldsymbol{\sigma}$ can take $2^{N}$ different values, and determining the full population probability distribution requires specifying $2^{N}$ probabilities, which is an unrealistic task even for modest population sizes. Assumptions on the analytical form of the distribution are therefore required in order to infer a complete distribution from a relatively small number of samples.

2.3 Maximum entropy models

The first, and perhaps more elegant, strategy developed to this end is to adopt a maximum entropy approach jaynes1957information , in which one first selects what features of the data should be exactly reproduced, and determines then the highest-entropy probability distribution consistent with those constraints. Schneidman et al. schneidman2006weak and Shlens et al. shlens2006structure first applied this approach to neural data, using a Pairwise Maximum Entropy (PME) model, which exactly fits all $\langle\sigma_{i}\rangle$ and $\langle\sigma_{i}\sigma_{j}\rangle$, i.e. firing rates and pairwise correlations. Indeed, the question behind that research was primarily related to the importance of correlations in the vertebrate retina, including the study of higher-order interactions.

The PME probability distribution over all codewords has the following form:

[TABLE]

The expression above is mathematically identical, in statistical physics, to that of the canonical ensemble for the Ising model with arbitrary couplings, a generalisation of the model originally used to describe ferromagnetism in solids ising1925beitrag .

By definition, a successful fit of a PME model correctly reproduces all firing rates and pairwise correlations present in the data from the considered neural population. Fitting based solely on second-order statistics does not imply that third-order correlations and other statistical measures are correctly reproduced. Reports that higher-order correlations are largely irrelevant were thus very surprising schneidman2006weak ; shlens2006structure ; Tang2008 , although these observation may be restricted to low activity and high pairwise correlations Yu2011 . Assessing whether a maximum entropy model can capture additional statistics of the data provides a source of interpretability: If, say, a PME model can account for third-order correlations, then the latter are not constrained further by the data. If, conversely, third-order correlations diverge from the PME prediction, we learn that the neural activity uses higher-order statistics to encode information. Whether this is the case depends on the system and on the distance between the neurons considered ohiorhenuan2010sparse .

Several attempts have been made at improving the quality of the fit of statistical models, using different features as known statistics. The generalisation of equation (2)

[TABLE]

can describe any probability distribution of binary variables exactly. However, finding the values of $J^{(n)}$ is computationally expensive and the benefits typically do not outweigh the costs for $n\geq 3$.

As a different way to assess at least some aspects of the higher-order statistics, we can consider, for instance, the probability distribution of the number of neurons firing in a time bin, $p(K)$, where $K(t)=\sum_{i=1}^{N}\sigma_{i}$, was used as a target. This can be introduced as a further constraint in a maximum entropy model in combination with firing rates and pairwise correlations, leading to the $K$-pairwise model tkacik2014searching ; mora2015dynamical . It typically produces significantly better fits than a pure PME, and is much less computationally expensive than attempting to fit higher order cumulants. Another related approach, the population tracking model, fits $p(K)$ together with the conditional probabilities $P(\sigma_{i}=1|K)$ of each neuron firing, given the current population firing rate, providing a lightweight and interpretable model odonnell2016population .

An example of an energy-based model which does not rely on the maximum entropy principle, finally, is to use a restricted or semi-restricted Boltzmann machine (RBM/sRBM). Despite not directly aiming at fitting correlations, $p(K)$, and cumulants, as a maximum-entropy model would, RBMs were shown to perform at least comparably well in fitting all these aspects koster2014modeling . An advantage is that their complexity can be tuned, offering a choice of various degrees of accuracy and the corresponding computational costs. Additionally, contrastive divergence, the algorithm used for fitting, is an approximate but relatively fast and reliable algorithm, which lets one fit the simultaneous activity of a large number of units (up to several hundreds, whereas the exact learning algorithm for a PME model is not usable in practice over $N\approx 40$, although more efficient methods have been studied). Finally, RBMs can be interpretable models, specifically by studying the roles taken by hidden units.

Although a detailed discussion is beyond the scope of this chapter, we should at least mention the efforts to reproduce the time dynamics of the system, so that the statistical model fits both the distribution of single-time bin patterns and the conditional distribution of the pattern given the pattern in the previous time bin marre2009prediction ; nasser2013spatio ; gardella2018blindfold .

2.4 Phase transitions in models

Although this may initially seem less relevant from the point of view of research in neuroscience, we should remind ourselves that the Ising model is one of the earliest and most commonly studied paradigms of a phase transition. To understand its behaviour, let us make the temperature dependence of equation (2) explicit:

[TABLE]

Note that, in the high temperature limit, this converges to a uniform probability:

[TABLE]

Conversely, when $T\to 0$, only a small number of states with non-zero probability survive, the others becoming infinitely rare. If a system that obeys $P_{T\to 0}(\boldsymbol{\sigma})$ is perturbed in any way, it will eventually converge to this stable set, under any reasonable dynamics. In other words, the distribution becomes, in the zero temperature limit, a finite set of stable attractors, the same as the stationary distribution of a Hopfield network Hopfield1982 , where the attractors play the role of memory patterns.

Clearly, neither of the two limiting cases can be a realistic description of neural statistics, and the truth lays in between them, in a regime where the model is much more informative. The physics literature shows that there is sharp phase transition between a disordered phase and a spin glass phase, with the exact location of the critical point depending on the statistics of $h$ and $J$ nishimori2001statistical . It is then a natural question to ask whether the Ising model that results from a fit to neural activity is in one of the two phases, or poised near the critical point, and whether this relates to other concepts of criticality in neural systems.

The divergence of specific heat, also called heat capacity, in a macroscopic system is a classic signature of discontinuity in the properties of the system upon variation of a single parameter, typically temperature (generalisations of this idea will be discussed in the next section). The most classic example is the case of a change in the state of matter — solid to liquid, liquid to gas, etc. — where an infinitesimal change in temperature through the critical point requires a finite amount of energy (the latent heat). This is equally true for spin systems of the form we examined above. Tkačik et al. tkacik2015thermodynamics fitted a model of the form (2) to binned spike trains from recordings of the salamander retina subjected to movies of naturalistic stimuli. They varied the temperature of the model around $T=1$ (this value corresponding to the fit to neural data), and studied the specific heat as a function of $T$ for an increasing number of neurons.

Their result clearly showed a peak in the specific heat of their models, with the peak temperature approaching $T{=}1$ as $N$ is increased. This is evidence that $T{=}1$ coincides with the critical point, and therefore, the model is poised at criticality for parameter values exactly corresponding to those that fit the neural data. Similar observations were independently repeated, e.g. in mora2015dynamical ; nonnenmacher2017signatures ; hahn2017spontaneous , generating a debate on the nature of this observation and its biological interpretation, as will be discussed in later sections.

2.5 Model criticality and Zipf’s law

Zipf laws can be related to statistical criticality in the sense of models, as shown in tkacik2015thermodynamics (supplementary information), as follows. Call $p_{1},...,p_{k},...,p_{2^{N}}$ the probability of occurrence for each of the $2^{N}$ possible codewords. In statistical physics, microcanonical entropy can be defined as $S=\log\Omega$, where $\Omega=\Omega(E)$ is the number of states with energy lower than $E$. On the other hand, the energy level associated with a pattern is a function of its probability:

[TABLE]

Now, Zipf’s law states that, for every pattern, its rank $r_{k}\propto 1/p_{k}$. In the notation used above, note that $r_{k}=\Omega(E_{k})$. Therefore, Zipf’s law implies

[TABLE]

If the above linear relation holds, then $d^{2}S/dE^{2}{=}0$. Since both specific heat and the variance of energy are inversely proportional to $d^{2}S/dE^{2}$, these thermodynamic quantities diverge. This is the classic signature of a second order phase transition.

The rank-probability relation defined by Zipf’s law, therefore, is a model-independent way of showing criticality in this statistical sense. Its appearance guarantees the divergence of the specific heat of a PME model fit to the same data, but does not require complex and computationally expensive fitting procedures, and relies only on the statistical properties of the data.

3 Statistical and dynamical criticality

As we have mentioned, most energy-based models do not account for dynamics, as they are concerned only with fitting a single-time bin distribution. The formulation of the Ising model in physics describes a stationary distribution and does not include any dynamics. Transition probabilities from a state to another can be added through additional assumptions about the dynamics of the system. For instance, Glauber dynamics Glauber1963 generates a Markov chain whose stationary distribution coincides with the distribution (2). This is useful, for example, when sampling states from that probability distribution. Avalanches can be observed in high-dimensional Ising systems when they are driven out of stationarity by a change in temperature or applied magnetic field. Therefore, it is not expected that maximum entropy models reproduce any aspect related to avalanche dynamics.

It is not clear a priori, then, whether the observation of Zipf laws and diverging specific heats should be related to power-laws in the dynamics. In fact, finding a connection between the two concepts seems challenging. In the next sections, we will provide examples of how the two might be entirely distinct, which prompts questions on the nature, meaning, and relevance of statistical criticality.

3.1 The Eurich model is dynamically, but not statistically critical

For a discussion of the relationship between the two concepts of criticality, it is interesting to consider the Eurich model for neural avalanches as a “testbed” eurich2002finite . It is mathematically well-understood and can be conveniently tuned, because the parameter values for (quasi-)critical as well as sub- or super-critical behaviour are known analytically (even for finite systems). The very definition of the model is such that all neurons have identical properties, and the same for pairs, triplets, etc., of neurons. As a consequence, all $N$ patterns with exactly one active neuron appear with the same frequency; all $N(N-1)/2$ patterns with exactly two active neurons appear with the same frequency, and so on, giving the rank-probability plot a step-like appearance, which cannot follow a Zipf law (Figure 2).

It is tempting to consider the tail of the rank distribution. Although the number of states increases with the activity (for neurally plausible activity levels), their probability decreases strongly if the firing rate (per time bin) is low. Therefore, steps will disappear for higher ranks, which may or may not produce a power-law-like behaviour. However, in the statistical approach, typically small values of $N\delta t$ (see equation 1) are used, such that the potentially Zipf-like tail (figure 2, right) will be statistically irrelevant. It has also been claimed that the statistical approach is most likely restricted to low-activity patterns Yu2011 .

Note that the rank curve would be less step-like if some form of heterogeneity is introduced, as opposed to the complete symmetry between neurons that characterises the Eurich model. However, this would amount to an additional assumption that is, as the Eurich model shows, not necessary for criticality. In particular, we would need to assume that the patterns with a single active neuron already follow Zipf’s law. This is a particular, but not unreasonable assumption, as this can be expected, for example, in a scale-free neural network. In fact, we cannot rule out that a Zipf profile, or at least an approximation, could be found just by tuning the distribution of firing rates, even in the absence of correlations. Such a finding would entirely rule out any relation to dynamical criticality, which appears exclusively as a consequence of emergent phenomena deriving from complex interactions. However, it would still require a specific distribution of firing rates among the neurons, an assumption that in itself would prompt questions about its functional reasons. Firing rate distributions have been studied extensively mizuseki2013preconfigured , and found to be highly skewed, with a small fraction of neurons responsible for the majority of emitted spikes. A clear theory on why this is an advantage for the encoding is still missing. In section 3.3, we will consider a case in which the Zipf relation holds even when correlations are destroyed, which suggests the long-tailed firing rate distribution is sufficient for it to hold.

3.2 Fitting energy-based models to critical activity

A natural way of checking if dynamical and statistical criticality are related could involve fitting a statistical model to neural models that exhibit various kinds of dynamics, and can be tuned to a supercritical (noisy), subcritical or critical regime. This was one of the goals of in Ref. hahn2017spontaneous . The authors identified five different dynamical states of the cat and monkey cortex, studied their avalanche statistics, and evaluated the temperatures corresponding to peak specific heats of Ising models fitted to each dataset. The results did show a small but significant relationship between avalanche dynamics and specific heat peak location. However, it should be noted that some of the results in this work were obtained with small datasets of six neurons only, which may not offer insight on what happens in the thermodynamic limit.

We attempted a similar task fitting Bernoulli RBMs to the activity generated by a tunable model of a neural network, similar to the binary, non-leaky, integrate and fire model used by gautam2015maximizing . In this model, the strength of inhibition can be tuned, leading to a network with low, random-looking activity and a short-tailed avalanche distribution (high inhibition); or a network generating activity in large bursts (low or no inhibition). The critical regime lies in between the two. Details about the implementation of the model are given below.

We found that although the absolute value of the peak does depend on the correlations of the data, its location is always at a temperature near $T=1$ (which is the value corresponding to the original fit), and further approaches this temperature as the number of units increases, i.e. in the thermodynamic limit. These results are compatible with what was shown by nonnenmacher2017signatures , using a different dataset, for the $K$-pairwise model.

It seems, then, not only that the statistical model that we fitted does not accurately detect criticality in the dynamical sense, but it also exhibits statistical criticality no matter the dataset it was fitted to. This implies, on the one hand, that the dynamical criticality of a dataset and the statistical criticality of a model fitted to it are unrelated, and, on the other hand, that a model fitted to datasets of very different nature all tend to exhibit statistical criticality. This is compatible with an argument that was put forward by theoreticians mastromatteo2011criticality , as will be discussed in section 4.2.

Methods: network model

The binary neuron model used for the simulations is similar to the one presented by gautam2015maximizing . In this model, each neuron has a probability of firing given by a weighted sum of its inputs, divided by a factor dependent on its own firing history. A fifth of all neurons are inhibitory, while the rest are excitatory; the value of inhibition was tuned to 0.0, 1.0, or 2.0, to enforce different regimes, corresponding to different correlations and avalanche statistics. While in the original work the connectivity was all-to-all, with weights drawn from a uniform distribution, we modified it to identical couplings, but set on a network with scale-free degree distribution. This enforced larger variability of firing rates between neurons; both experimental evidence and the theory in larremore2011predicting suggest this choice does not affect the location of the critical point. We simulated 1000 neurons, from which we took subsets of the sizes required for analysis, for 1 million time steps (conventionally taken to equal 1 ms). The resulting activity was re-binned in 5 ms bins, to reduce sparseness.

Methods: RBM specific heats

As we will argue in section 4.2, the direction that best indicates the critical point coincides with the fist eigenvector (the one corresponding to the largest eigenvalue) of the Fisher information tensor. However, in practice, this is never orthogonal to the direction of increasing/decreasing temperature: thus, varying temperature is an acceptable way to look for a phase transition.

In statistical physics, the general expression for the probability of a pattern in an energy-based model at temperature $T$ is

[TABLE]

The expression for the energy in the case of RBMs is

[TABLE]

Where $\mathbf{v}$ and $\mathbf{h}$ are vectors of visible and hidden binary variables respectively. Since this expression is linear in the parameters $a_{i},b_{j}$ and $J_{ij}$ for all $i,j$, changing the temperature of a model coincides with rescaling these parameters by a linear factor $\beta=1/T$. In the following, we have adopted the standard strategy of fitting an RBM to neural data, obtaining values for its parameters, and then rescaling them — this means $T=1$ (no rescaling) coincides with the parameters as they were fitted, the values corresponding to a model that correctly reproduces the given data. Fits were obtained by 1-step persistent contrastive divergence.

We can then compute the specific heats at different temperatures. The marginal probability of $\mathbf{v}$ is

[TABLE]

Disregarding an additive constant, the energy of a visible pattern can be expressed as the logarithm:

[TABLE]

In accordance with statistical physics, we can define

[TABLE]

This quantity can be computed from a sample. For each temperature value, in each dataset, we extracted 20 chains of 2000 samples, taken every 10 steps in order to reduce spurious correlations, and computed $c(T)$ by the expression above.

3.3 Retinal activity obeys Zipf’s law, but is not dynamically critical

The mammalian retina, a system that is often chosen when studying the statistics of neural activity, and whose encoding and dynamical properties are well known, is an example of the opposite case: It was the first system in which statistical criticality was observed, but it does not exhibit dynamical criticality.

Avalanches arise in the mammalian retina only during the period of development: for mice, in the first few days after birth, before eye opening, when the retina does not respond to light and the network activates spontaneously. During this stage, the activity of the retina consists of the so-called retinal waves, which are effectively power-law distributed avalanches. Direct comparison with a computational model showed that these are indeed the signature of a critical state between locally and globally connected activity hennig2009early . However, these disappear in a functional retina: Figure 4 shows the statistics of a 20-minute recording of an untreated, adult mouse retina under an uncorrelated black-and-white checkerboard stimulation. It is evident that the avalanche statistics is short-tailed, and, at the same time, the probability-rank plot of pattern frequencies is well compatible with a Zipf law. Note that correlations between the activities of retinal ganglion cells change significantly with the statistics of the stimulus, and the avalanche statistics will consequently appear different. The example of adult retinas is complementary to Section 3.1 in the sense that, here, a system that does not show dynamical criticality can well obey Zipf’s law.

It is worth mentioning that the observation of Zipf law in retinas is very robust to a number of external factors. We found no significant differences in the rank-frequency plots of patterns observed when the retina was treated with bicuculline (a GABA blocker) compared to a control; analogously for retinas under stimuli characterised by very different level of spatial correlations.

If Zipf’s laws have a functional role, there is no expectation this phenomenon would survive in a non functioning neural system, such as a retina that has been pharmacologically treated in a way that breaks its normal operative mode. Here, we took data from the same mouse retina, before and after treatment with a 20 $\mu$M solution of bicuculline, which is a gabaA antagonist. The results are shown in figure 6 (left): as it is evident, there is no clear difference between the two rank-probability plots. Of course, the only strong argument against the functional role of Zipf laws would be finding a functional retina in which this law is broken, which is not the case here. However, we can notice that even an intervention that significantly disrupts the retina’s activity, by blocking inhibitory interactions, doesn’t prevent this phenomenon to arise. This is despite the large change in the correlation between neurons induced by bicuculline.

Likewise, one may expect a dependence of pattern frequency-rank statistics on stimulus statistics. The retina, after all, is a neural system design to encode a stimulus — and the correlation structure of its neurons’ activity strongly depends on the correlations in the stimulus. However, we found no significant difference in Zipf laws under different stimulations. Figure 6 (centre) shows a single group of 100 neurons selected in a retina that was stimulated with light patterns of different kinds. All stimulus presentations consisted of black-and-white random checkerboards, which are binarised versions of random noise of given frequency spectrum $f^{-a}$ with $a=0.5,1.0,1.5,2.0$ in space and time: from near-white noise to the statistics of natural images (figure 5).

The independence from correlations is evident in the right panel of figure 6: here, the “control” curve is the same as in the left panel, and is compared with the rank-probability plot for a “shuffled” version of the same data, where the firing rates were kept the same, but spikes were moved in time in order to cancel neuron-neuron correlations. The difference between the two curves is clearly not significant. This demonstrates how a firing rate distribution which is long-tailed (approximately log-normal) can in itself produce a Zipf-like plot. More research is needed to show whether this holds in general.

Methods

The handling of the retinas, experimental apparatus, and the first part of the data analysis pipeline were performed as illustrated in hilgen2017unsupervised . Starting from detected and sorted spikes, we removed those with very low amplitudes, by selecting a threshold corresponding roughly to the lowest 10%. This was to ensure only good-quality events were left. Then, we selected, for each Zipf plot, $N=100$ clusters all pertaining to the same area of the retina (figure 7).

At this stage, spikes were binarised into a $N\times T$ matrix $S$ of boolean variables, with $S(n,t)=1$ if neuron $n$ spiked between times $t$ and $t+\delta t$ and $S(n,t)=0$ otherwise. When multiple spikes from the same neuron occurred in a single time bin, the extra spikes were disregarded. For recordings shown in this chapter, $T=120000$ or more, and $\delta t=10$ ms, implying at least 20 minutes of neural activity were recorded.

4 Parametric sensitivity

The basic fitting procedure of a maximum entropy model minimises the quadratic difference between the data moments and the moments predicted by the model. During fitting, any model is updated by exploring the parameter space, following a direction given by the loss function. When a model admits a phase transition, the parameter space is characterised by (at least) two regions, corresponding to the phases, separated by a critical surface. From a theoretical point of view, asking why a model is poised at criticality coincides with asking why the fitting process tends to lead towards the critical surface in the parameter space. This has been discussed by mastromatteo2011criticality ; before introducing their argument, we provide some theoretical background.

4.1 Model distance

Intuitively, a phase transition occurs at a location (the critical point or critical surface) where an arbitrarily small change in the parameters yields a sharp, qualitative change in the behaviour of the model. In this section, we will formalise this idea, and link the notion of model distance to the statistical physics framework that we have introduced above.

A common measure of the distance (in model space) of a probability distribution $p$ from a given one $q$, both defined on a set $S$, is the Kullback-Leibler divergence

[TABLE]

It measures the amount of information that is lost when approximating $p$ by $q$. The name divergence stresses that this quantity does not have the mathematical properties of a distance, namely not being symmetric. If the model space is parametrised by ${\theta}$, and $p$ and $q$ are close to each other in this space, so that $q=P_{\theta}$ and $p=P_{\theta+\delta\theta}$, at second order in $\delta\theta$, the divergence can be approximated as

[TABLE]

where $F$ is called Fisher information tensor (FIT) which is given by

[TABLE]

Fisher information is here expressed as a statistical quantity, but it has an important relation to the physics of statistical models. Consider a Hamiltonian model, where the probability distribution is given by

[TABLE]

which is an obvious generalisation of maximum entropy models. Calculating the Fisher information for this form of $P$ (6), we retrieve the direct connection between the covariance (with respect to $P_{\theta}$) of the physical quantities $f$ and the FIT that is characteristic for probability distributions of the exponential type:

[TABLE]

This means that the FIT characterises the variances and correlations of the functions $f$ which are now considered as stochastic variables and depend on the state $x$ of the system.

Additionally, note that changing the temperature in the traditional canonical ensemble corresponds to scaling the Hamiltonian by a factor $\beta=1/T$, similarly to equation (3). In the formulation above (6), this is equivalent to scaling all the $\theta_{i}$ by $\beta$. Given a point $\theta$ on the parameter manifold, the direction $\partial/\partial\beta$ can be expressed as

[TABLE]

which is just a linear combination. The specific heat is given by

[TABLE]

Thus, we can arrange for the specific heat to be one of the entries of the Fisher information matrix, with a change of basis, which includes the $\beta$ direction as a base vector together with other $n-1$ linearly independent ones. Analogous considerations can be made for the magnetic field and magnetic susceptibility. In this sense, the Fisher information tensor is a generalisation of specific heats and susceptibilities.

4.2 Fisher information and criticality

We can now look at the relationship between statistical criticality and the model’s parameter space. Suppose any generalised susceptibility (i.e. a component of the Fisher tensor) diverges at a point $\theta_{0}$. Then an eigenvalue of the Fisher information, say $\lambda_{k}$, diverges at $\theta_{0}$. Call $v_{k}$ the corresponding normalised eigenvector. For small $\alpha$,

[TABLE]

and the r.h.s. diverges. This means that, moving from $\theta_{0}$ in the $v_{k}$ direction by an arbitrarily small step yields a model $P_{\theta_{0}+\alpha v_{k}}$ that is completely different from $P_{\theta_{0}}$, as indicated by an infinite KL divergence.

We introduced this description in terms of Fisher information in order to give an interpretation of criticality from the point of view of modelling. A model is at a critical point whenever there is a direction in parameter space that leads to an infinitely different model by a finite change in parameters. This, incidentally, shows that the best way of measuring the distance from a critical point is not to vary temperature, but to use the first eigenvalue of the FIT and move in the direction of the corresponding eigenvector. Temperature is not always the most relevant control parameter.

Mastromatteo and Marsili mastromatteo2011criticality have argued that, because of this special property critical points have in the parameter space, they are particularly favoured by model fitting. In particular, they show that distinguishable models accumulate near critical points, whereas models farther away from these are largely indistinguishable. Their argument, in brief, goes a follows. Two models are considered indistinguishable if their Kullback-Leibler divergence is less than a given value $\epsilon$. For small $\epsilon$, $D_{KL}$ is approximated by Fisher information, and the volume of parameter space occupied by models indistinguishable from $\theta_{0}$ turns out to be proportional to $(\det F(\theta_{0}))^{-\frac{1}{2}}$. This quantity diverges at critical points due to the first eigenvalue diverging as explained above. Thus, most models actually are poised near a critical point, according to this metric. They conclude that criticality may be a feature induced by the inference process, rather than one intrinsic to the real system being studied by the model. This may be the reason why statistical models seem to be poised at a critical point, for a variety of training datasets, as we showed in section 3.2. However, it does not affect Zipf laws, which are directly observed in the data.

4.3 Criticality and parameter ‘sloppiness’

It is well known that the parameter spaces of many models often show only a small number of directions (linear combinations of parameter changes) along which the overall properties of the model strongly change (“stiff”), and a large number of directions which have little influence on the model (“sloppy”). This phenomenon, termed “model sloppiness”, has been observed in a wide number of cases in systems science gutenkunst2007universally ; machta2013parameter .

For the specific case of neuronal networks, in Ref. panas2015sloppiness , although for small numbers of neurons, “stiff” dimensions corresponding to large FIT eigenvalues were identified. The remaining “sloppy” dimensions, on the other hand, can change without much effect on the goodness of fit of the model. A further development of this approach has been reported in Ref. herzog2018dimensionality , where it was shown that about half of the dimensions in the data manifold are irrelevant for the modelling. As shown in paragraph 4.2, near a critical point, the direction pointing towards the critical surface has a diverging FIT eigenvalue, while the others are smaller. This hints there may be a connection between sloppiness and criticality, which, at the moment, we can only leave at the level of speculation.

Additionally, however, sloppiness indicates that a fitting algorithm for the data may be improved if different dimensions are differently weighted during the optimisation process. We can then ask whether using a natural gradient in the fitting procedure would lead to a different result while evaluating model criticality. In natural gradient optimisation, the components of the gradient are compensated by the inverse Fisher information, i.e. the divergence near a critical point of the model would disappear, at least theoretically when the Fisher information is exactly known. As a result, the fitting procedure is not homogeneous with respect to the set of the parameters, but with respect to the space of the parameters, taking into account its geometrical properties, and parameters can be identified equally well in all regions. In this way, the problem discussed in Ref. mastromatteo2011criticality may disappear — more research will be needed in order to verify this.

5 Discussion

Neuronal avalanches are an experimentally well-studied phenomenon, that can be explained as a consequence of the optimisation of information processing in the brain. It should be noted that an understanding of how the potential functional benefits of this “dynamical” criticality are realised is missing shew2013functional — however, it has been shown that the maximisation of the dynamical range happens at criticality Kinouchi2006 .

Statistical criticality is an equally complex phenomenon to explain theoretically. Like dynamical criticality, it can be taken to indicate the complexity of the neural data and the relevance of higher-order correlations or latent variables, but its functional implications are less clear. In this chapter, we have reviewed the concept, both in the context of fitted statistical models, and as a direct observation of Zipf laws in neural population data. Through experiments on restricted Boltzmann machines, we suggested that the divergence of model specific heat is not a reliable way to infer properties of the data. We mentioned how Fisher information provides the correct description of the parameter space and the critical surfaces, and reviewed a possible explanation of why statistical models tend to poise themselves at a critical point. Then, we tried to describe the connection between statistical and dynamical criticality, and argued there is no clear connection, by showing examples where one of the two was present without the other. Further insight on this matter might come from models that are capable of both, provided they can reproduce not only the equilibrium distribution of the data, but also the dynamics. A multi-time maximum entropy model might provide a starting point for this work.

Of course, it may well be that the observation of Zipf laws is simply a consequence of problems related to how we describe the data — these include the typically small sets of observables, the choice of binning size, failure to account for the real dynamics, and biases introduced by sampling. However, the ubiquity of Zipf laws in complex systems means its emergence in biological neural networks should not surprise us, and it could be explained in terms of mechanisms such as the one described by aitchison2016zipf , or perhaps with preferential attachment. Conversely, an important open problem is an explanation on whether statistical criticality is something that is actively sought by the system because of some functional relevance. On this matter, we tried to analyse the Zipf profile of retinal activity under various conditions (various stimulus statistics, pharmacological treatment), but we found no significant differences in the cases examined. Interestingly, it seems to be possible to generate a Zipf profile simply by enforcing a long-tailed firing rate distribution, despite the absence of correlations. Even if this observation were confirmed, the question would simply shift towards finding a reason for such a skewed distribution of firing rates, which has not yet found a justification in terms of function.

Notably, recent research has started showing how Zipf laws appear in different kinds of parametric models, including “deep” ones, as soon as learning occurs. It has been shown that the Zipf property arises to different degrees in different layers of a deep network, and is maximal in the layers that attain an optimal trade-off between resolution and accuracy in generating samples song2017emergence . This is a starting point in linking statistical criticality to function. It is not known whether similar principles are relevant in the case of biological neural networks, and finding such a link could be an interesting direction of future research.