Testing the critical brain hypothesis using a phenomelogical renormalization group

We present a systematic study to test a recently introduced phenomenological renormalization group, proposed to coarse-grain data of neural activity from their correlation matrix. The approach allows, at least in principle, to establish whether the collective behavior of the network of spiking neurons is described by a non-Gaussian critical fixed point. We test this renormalization procedure in a variety of models focusing in particular on the contact process, which displays an absorbing phase transition at λ = λc between a silent and an active state. We find that the results of the coarse-graining do not depend on the presence of long-range interactions, but some scaling features persist in the super-critical system up to a distance of 10% from λc. Our results provide insights on the possible subtleties that one needs to consider when applying such phenomenological approaches directly to data to infer signatures of criticality.

The possibility that living systems may be poised at criticality is a fascinating hypothesis [1][2][3], and in recent years it has been explored in a vast variety of areas [4][5][6][7].
Tools from statistical mechanics, such as the renormalization group [8][9][10], teach us that at criticality the macroscopic, collective behavior of the system is described by a few relevant attributes, such as the embedding dimension of the system and its symmetries, while most the microscopic details of the system become irrelevant. At the critical point the physical properties are determined by a non-trivial fixed point in the space of the possible models. However, in the broad landscape of natural systems one often has to deal directly with data without an explicit model, and the systems are typically finite, so that most of the time it is hard to come up with a definitive answer about whether they are poised near a critical point [3].
Recently, a phenomenological coarse-graining procedure was introduced in [11,12] to deal with the longrange interactions that one reasonably expects in a network of neurons, but of which the full interaction network is not necessarily known. Data from single-neuron recordings, from the hippocampus of a mouse running along a virtual track, were directly analyzed by the authors in order to understand if this coarse-graining procedure (which we recall in Section 1) drives the system towards a non-trivial fixed point in the renormalization group sense, hence if the neural dynamics is critical and details independent. Indeed, the brain is probably one of the most impressively complex system we are able to study and the idea that the collective behavior of neurons might emerge from a self-organized critical state has been widely studied in the last year [2,[13][14][15].
One of the first evidences that suggested this hypothesis is the presence of neuronal avalanches that spontaneously occurs in the brain, i.e. during spontaneous activity, that show a spatio-temporal power law distribution with exponents compatible with those of a mean field branching process [6]. However, this conclusion is highly debated, in fact, and the same exponents may stem, for instance, from an underlying non-critical neutral [16] or random [17] dynamics and, in general, the subject is far from being settled.
In this paper, we aim to test this phenomenological renormalization group (PRG) method by applying it to a well known non-equilibrium statistical model, the contact process [18]. This model belongs to the universality class of directed percolation and displays an absorbing phase transition, which has been widely studied [19]. Its critical behavior is well understood and the exponents are known from numerical studies, so we shall regard it as a "control case" to investigate the ability of this procedure to extract the relevant information and infer signatures of a critical state in out-of-equilibrium systems.
On a d-dimensional hyper-cubic lattice with nearest neighbor interactions it is sufficient to introduce longrange connections to change the topology, for instance, to that of a small-world network. Hence by simulating the contact process we are able in particular to probe the impact of short and long range interactions on the coarse grained system behaviour, and in particular we are able to shed some light on the possible outcomes and interpretations of the emergent fixed point describing the system collective behavior. Along the road, we also test the PRG procedure in other models to better characterize its results.

I. THE COARSE-GRAINING PROCEDURE
In this section we briefly describe the coarse-graining procedure introduced in [11,12] that we aim to test. The authors propose to build clusters of variables by grouping together neurons that are most correlated, so that the overall correlation structure tends do be preserved.
Let us consider a system (e.g. neural circuit) of N vari-arXiv:2001.04353v1 [cond-mat.stat-mech] 13 Jan 2020 ables (e.g. neurons) connected by a given, but typically unknown, interaction network. Denoting state variables of the neurons as σ (1) i for i = 1, ..., N , where the superscript 1 denotes that we are at the first step of the renormalization procedure, we search for the maximal non-diagonal element of the normalized correlation matrix where · represents the average over the time-series of neural activity. The pair (i, j * (i)) of maximally correlated variables is removed and we search again, ending up with a set of pairs {i, j * (i)}. The coarse-grained variables are defined as where i = 1, . . . , N/2. We iterate this process, producing clusters of K = 1, 2, 4, . . . , 2 k−1 variables. Each one defines a new variable σ (k) i as the summed activity of cluster i.
Under this coarse-graining procedure, in [11,12] the behavior of various quantities is analyzed in order to make some parallels with the behavior of critical systems. In particular, the following observables are studied: the mean variance of the neural activity; the distribution of the individual coarse-grained variables; the spectrum of the covariance matrix; and the mean autocorrelation function.
The mean variance of the activity is defined as where N k is the number of variables after k steps of the coarse-graining procedure. If the variables are independent one would obtain a variance scaling as M 2 (K) ∝ Kα withα = 1.
More generally, we can study the full distribution of the individual coarse-grained variables. Since a coarsegrained variable σ is the probability distribution of the normalized activity. Thus we look at an effective (reduced) free energy (this formula is based on the assumption that the energy of the system is zero when no activity is present) and at its possible scaling F ∼ −Kβ. For independent variables, we expectβ = 1. A scaling behavior of the ranked spectrum of the covariance matrix at the critical point is expected. In fact, the correlation function decays algebraically as G(x) ∼ |x| −(d−2+η) , and one can show (see Appendix A) that in translational invariant systems the eigenvalues of the covariance matrix scale as where r is the rank of λ r , ordered from the highest to the smallest. If we consider the variables inside the clusters that we build along the coarse-graining, the highest possible rank r is given by the number of variables K that make up each cluster. Hence at criticality we should find with µ = (2 − η)/d, and this is a direct consequence of the power law decay of the correlation function in space. Finally, the mean autocorrelation function is obtained by where C (k) Since we are grouping correlated variables, the decay of the autocorrelation is slower in clusters of bigger size. However, in a critical system we might expect dynamical scaling, which would imply a power law scaling of the autocorrelation times τ c ∝ Kz.
A different test can be performed by exploiting the fact that, in systems with translational invariance, the eigenvalues λ k of the covariance matrix in momentum space are the Fourier transform of the correlation function G(k). Since coarse-graining in momentum space amounts to average over the Fourier modes with small wavelength, we expect that averaging over low variance contributions in the covariance matrix should lead to an equivalent result. Hence, we consider the set of eigenvectors of the covariance matrix {u r }, ordered according to the value of the corresponding eigenvalue, from the the highest to the smallest one, and we introduce the projectors where P ij (N ) is the identity, hence the eigenvectors are orthonormalized. The authors of [11,12] propose to consider a cutoffK < N , in analogy to the cutoff in momentum space, in such a way that the low variance contributions do not enter the projector (6. Then the coarsegrained variables are defined as where z i (K) assures that the coarse-grained variables have unitary variance, i.e., φ 2 i (K) = 1. By means of the Young-Eckart theorem [21], the above procedure allows one to find the best decomposition with rankK of the original data matrix. In this setting it is interesting to look at the distribution as we change the cutoffK. In fact, a renormalization group transformation typically drives the joint probability towards a fixed point, and if the variables are weakly correlated such fixed point is the one obtained from the central limit theorem [22]. Hence, the authors of [11,12] propose to use this PRG approach to test whether the joint distribution converges towards a non-Gaussian critical fixed point.

II. THE MODEL
The contact process [18,19,23] is possibly the simplest non-equilibrium model used to describe the propagation of neural activity on a network [15,16]. Hence, it is a proper modelling framework to test the coarse-graining procedure just described. Consider a collection of N nodes of a given network. Each node can be either active (occupied) or inactive (empty), and we identify its state by means of a binary variable σ i (t) = 1, 0 respectively. The activity spreads via a nearest neighbors interaction, and it depends on the number of active neighbors n i (t) = j∈ i σ j (t), whereas each active site is emptied at a unitary rate. The rates w[σ i (t) → σ i (t + dt)|n i (t)] that define the process for a node with k i neighbors are given by where λ is the spreading rate.
The configuration with all empty sites is an absorbing state, since the system cannot escape from it. In particular, if λ > λ c the stationary state is an active fluctuating phase, whereas if λ < λ c the system eventually gets trapped in the absorbing configuration. Hence, we choose the density ρ of the active sites as an order parameter, while λ is the control parameter. In fact, exactly at λ = λ c the density of active sites undergoes large fluctuations and the system is often close to the absorbing state [24].
The nature of the phase diagram can be readily understood from the mean field approximatioṅ This equation has two stationary solutions: ρ v st = 0 and the active state ρ a st = (λ − 1)/λ. The former is stable if λ < 1, and the latter if λ > 1. Hence the mean field critical point is λ MF c = 1. The contact process is not exactly solvable even in one dimension, therefore we need to rely on numerical studies. We implement the usual scheme [26]: an occupied site i is randomly chosen, and with probability 1 − p λ = 1/(1 + λ) the site is emptied. With probability p λ = λ/(1 + λ) one of the neighbors is picked at random and, if empty, is occupied. The time is increased by 1/N occ , where N occ is the number of occupied sites.
We are interested in two different types of interaction network topology to test their effect on the coarsegraining procedure: a 2D lattice with periodic boundary conditions and a small-world network. The former is a rather standard choice, and the estimated critical point is λ 2D c ≈ 1.6488 [19]. On the other hand, the latter setting is more realistic, given the existence of long synaptic connections occurring in a network of neurons.
In general, it would be ideal to have a coarse-graining procedure that works both for short-range and long-range interactions, especially if one needs to deal directly with neural activity data and the specific network architecture is not accessible. In the small-world case the critical point λ SW c depends on the rewiring probability, and it has been studied numerically in [27].
We perform all the simulation with N = 40 2 sites and analyze clusters of size K = 2, . . . , 256. In momentum space, we keep up to N/128 ≈ 12 eigenvalues, which is less than 1% of the original modes.  (3). Both are fitted with the corresponding power laws shown in Section 1. Notice that both the critical and the super-critical regime show power law behaviour but with different exponents: in particular, the silence probability is smaller and decays much faster in the super-critical case due to the proliferation of the activity. However, we do not find full compatibility between the super-critical contact process and the independent case.

III. RESULTS
We now consider the contact process in a 2D lattice, both at λ = λ 2D c ≈ 1.6488 and in the super-critical phase at λ = 3. We will show that the PRG is indeed able to distinguish between these two different phases as the coarse-graining procedure gives indeed results that are rather different. However, we find the existence of some caveats that is important to take into account.
In Figure 2 we see that in the super-critical regime the exponents of the variance and of the free energy are not exactly compatible to the independent caseα = 1 = β. Nevertheless, the profile of the free energy clearly shows that the underlying dynamics is different, as for clusters of size K > 32 the silence probability vanishes in the active phase. Notably, if we compare the critical exponentβ = 0.65 ± 0.02 with the one obtained for real neurons in [11,12],β neurons = 0.893 ± 0.003 we see that in the contact process the decay of the silence probability with the cluster size is slightly slower, meaning that the sites tend to be more active than real neurons. Figure 3 shows the correlation structure of the system's quasi-stationary state. The change in the spectrum of the covariance matrix is more evident, since in the supercritical case the eigenvalues span a smaller set of values. In the critical case, instead, we find a power law decay with an exponent µ = 0.63 ± 0.02, with µ = (2 − η)/d. In real neurons, [11,12] report µ neurons == 0.71 ± 0.06. We note that since one of the hyper-scaling relations of the contact process yields [28] we expect µ ≈ 0.6, which is compatible with what we find using the PRG procedure. The time-autocorrelation function shows an evident change as well: in the supercritical regime the autocorrelation decays exponentially, whereas at criticality we find a power scaling with an exponentz = 0.50±0.06. We note also that in the supercritical regime a power-law seems to be present, but the small exponent is compatible with the absence of scaling [29].
The evolution of the joint probability distribution of the coarse-grained variables in Figure 4 shows once more the differences in the underlying dynamics. The most notable result is the convergence in momentum space: for λ > λ 2D c the fixed point is Gaussian in accord with the central limit theorem, whereas at λ = λ 2D c we see distinct non-Gaussian tails. The last coarse-graining step in momentum space only keeps N/128 modes, so the fact that we still see non-trivial tails is significant.
These results prove to be very stable when we change the underlying topology and we introduce long-range connections by choosing a small-world network. In particular, we implement a Watts-Strogatz model [30] with a rewiring probability p = 0.01. The critical point is λ SW c ≈ 1.7961, as given by [27]. We find that all the considerations we made so far hold in the small-world topology as well, and the presence of long-range interactions does not affect the results of the PRG coarse-graining procedure.
As a sanity check, we also implement a synchronous update algorithm. The results show no difference with respect to the asynchronous one used insofar [31].

A. Persistence of the scaling near a critical point
A natural question one may ask is how sensible this PRG approach is, i.e., how easy it is to distinguish a truly critical system from a super-critical one. We test this in the contact process by moving the control parameter from the critical point λ c to λ nc ≈ 1.1λ c , that is a 10% increase. Notice that, although it is not trivial to define a finite-size critical point [32] for the transition in the contact process [33], at λ nc we do see distinctive features of a super-critical dynamics: in fact, we do not see  (2). Right panel: evolution of the probability distribution of the coarse-grained variables in momentum space, as in Equation (8). We show the results for K = 32, . . . , 256 and for N/8, . . . , N/128 modes (from brighter to darker color). The distribution in direct space is very different due to the critical contact process being typically close to the absorbing state. In momentum space, the super-critical contact process converges to a Gaussian fixed point in agreement with to the central limit theorem [22], whereas at criticality we see the presence of non-Gaussian tails.
considerable fluctuations in the density of sites, nor the system constantly approaches the absorbing state as at λ = λ c . Hence, at λ nc the dynamical evolution is significantly super-critical, and we shall refer to this as a near-critical case to distinguish it from the super-critical regime we described before.
In Figure 5 we see non-trivial scaling behaviors of both the variance and the free energy. If we compare them to Figure 2, they are arguably more similar to the critical regime rather than the super-critical one. Noticeably, the exponentsα andβ are in between the two cases (i.e., critical and super-critical ones), suggesting that as λ smoothly changes from λ = λ c to λ = +∞, the exponents smoothly approach 1. This result calls for carefulness as the scaling inferred from the PRG of the variance and of the free energy are not related only to the system critical state [34]. Indeed, the difficulty to distinguish between critical or quasi critical states is confirmed also from other studies (using different approaches) [35].
On the other hand, in Figure 6 the eigenvalues of the covariance matrix do not display an evident power law scaling as we change the cluster size, and the scaling of the autocorrelation time function is not significant, especially for larger clusters.
The most convincing results to discriminate between critical and quasi-critical state are the joint probability distributions (2,8) that we show in Figure 7, in particular the one in momentum space. We do not see the non-Gaussian tails that we previously found at the critical point, which is expected since away from criticality the variables are much less correlated with one another and they are eventually dominated by the central limit theorem.
Nevertheless, following [17] in Appendix B we introduce a simple model of conditionally independent neurons that shows a non-trivial form of the joint probabil-  (3) at λ = λnc. Even if this is a super-critical contact process, the exponents are far from being comparable with the independent casẽ α = 1 =β. Overall, the scaling is more similar to the real critical case, see Figure 2. Notice in particular how the silence probability is non vanishing even for large clusters.
ity distribution along the coarse-graining in momentum space, although the system is clearly non in a critical state. At the same time, for this model we find that the clustering of maximally correlated pairs of variable seems to work in identifying the model as not critical, because there is no interaction between the neurons themselves.
All these results therefore suggest that in order to infer the system state is crucial not to focus on a single observable, but analyze all the different coarse-grained quantities.

IV. DISCUSSION
The phenomenological approach introduced in [11,12] has two considerable advantages: it is model independent so it can be applied directly to the data, and it is stable with respect to the presence of long-range interactions. In this work we tested the PRG both in equilibrium models (see Appendix C for the Ising model), where we expect that it is be able to distinguish between critical and noncritical phases, but also in non-equilibrium models, such as the contact process. We have found that the super and sub-critical regimes can be easily recognized, even though the nature of the phase transition is qualitatively different from the one of the Ising model.
At the same time we have highlighted that quasicritical states are difficult to infer, especially if only a subset of physical quantities is analyzed. In non-trivial dynamical models, such as the contact process, the strategy that works best seems to be the one related to the correlation structure. For instance, the presence of non-Gaussian tails in the joint probability distribution of the coarse-grained variables in momentum space, in Figure  4 and Figure 7, might be a good signature of a possible underlying criticality, but at the same time clustering maximally correlated variables fails as we approach the critical point, Figure 2 and Figure 5. Interestingly, in considerably simpler models (such as the one of Appendix B) the situation is reversed. Hence, in principle,  6. Scaling of the eigenvalues of the covariance matrix for K = 32, 64, 128 (from brighter to darker color) and of the autocorrelation times during the coarse-graining for λ λc. The spectrum of the covariance matrix does not show a power law decay, and the scaling of the autocorrelation times is not as convincing as the critical case shown in Figure 3.
one needs to study both the approach via maximally correlated pairs and in momentum space. Notably, this is the case of [11,12], hence the claim of the authors that the neuronal dynamics is critical should still hold.
However, we believe that this approach gives in general a set of necessary conditions for criticality rather than sufficient ones. For instance, the presence of a non-trivial distribution of the coarse-grained variables in momentum space is a necessary condition for criticality because it implies that the underlying variables are strongly correlated, but the convergence in the critical case is hardly clear and calls for particular attention when dealing with experimental data. One might also wonder, inspired by simple non-critical models as the one studied in [17], if there are distinct contributions to the dynamical evolution of the system -an intrinsic one, given by the interaction between the microscopic degrees of freedom, and an extrinsic one, given by some external driver [36]. The binomial model we propose in Appendix B is an archetypal example of the latter, and further research is needed if we want to distinguish the two contributions and, eventually, understand which of them contributes to poising the system at criticality.
Overall we believe that, as it is, this PRG might be considered as a better method to infer the presence of a critical state with respect to typical inference methods based on the identification of avalanches in both size and duration with particular exponents [3,6,[37][38][39].
The existence in the data of the signatures of criticality we have highlighted -such as non-Gaussian tails in the distribution of coarse-grained variables in momentum space -are possibly a powerful and stable indicator to characterize the state of a neuronal network. Extending these methods and test them systematically as in the present work might provide further insights in the understanding of the role of criticality in living systems. For example, recent works suggest that the actual transition in the brain dynamics is not between low and high neural activity states, but rather between an asynchronous and synchronous states [39,40]. An interesting future Let us briefly show that a power law spectrum of the covariance matrix is the consequence of the algebraic decay of the spatial correlation function at the critical point. Consider the covariance matrix In a system with translational invariance each element of the covariance matrix is given by C ij = C(x i − x j ) for some function C, whose Fourier transform is given by Hence the covariance matrix has entries given by which means that in Fourier space the covariance matrix is diagonal. In fact, it is easy to show that the eigenvalues are given by the Fourier transform of the correlation hence e ikx is a eigenfunction of eigenvalue G(k). This has a non trivial implication for the eigenvalue spectrum of the covariance matrix in a critical system, where we expect the algebraic decay G(r) ∼ r −(d−2+η) . Since the eigenvalues are the Fourier transform of the correlation function, we shall write If this is a decreasing function of |k|, that is if η < 2, then we consider a ranking of eigenvalues from small momentum to large momentum. Hence the highest eigenvalue has rank r = 1, which implies This implies that the eigenvalues of the covariance matrix decay as a power law of their rank, namely where λ 1 ≥ λ 2 ≥ . . . λ N and µ = (2 − η)/d.

APPENDIX B: MODEL OF CONDITIONALLY INDEPENDENT NEURONS
The convergence of the joint probability distribution in Equation (8) is related to the spectrum of the covariance matrix, but one should be careful when considering its relation with criticality. Consider N random variables (σ t 1 , . . . , σ t N ). At each time t the distribution of the i-th variable, which we can think of as a neuron that can be either active or inactive, is a simple binomial distribution with parameter ξ i (t). However, we consider the case in which also ξ(t) = (ξ 1 (t), . . . , ξ N (t)) is itself a random variable distributed according to some distribution p(ξ), so that the N neurons are conditionally independent. We can think of this case as the one of neurons that follow a common external dynamics, represented by p(ξ), but otherwise show no intrinsic dynamical features.
The probability that a neuron is either active or inactive is then a binomial distribution conditioned to the value of ξ i (t), that is , and the each neuron is described by the joint probability p(σ t i , ξ) = p(σ t i | ξ = ξ(t)) p(ξ). Since the coarse-graining procedure depends on the equal-time covariance of the neurons cov(σ i , σ j ) = C ij , we need the marginal probabilities It is then trivial to see that, even if the neurons are not correlated, their covariance is not vanishing but depends on the covariance of p(ξ), Let us consider the simple case of ξ i = ξ j ∀i, j, so that at each time all the neurons fire with the same probability ξ, and take p(ξ i ) to be a uniform distribution. In this case the covariance matrix is simply with a = 1/4 and b = 1/12. The eigenvalues of this matrix are given by where m is the corresponding multiplicity. Therefore, there are N − 1 eigenvalues with the same value. The eigenvector associated to the highest eigenvalue is 1/ √ N (1, . . . , 1) T , but there is no obvious choice for the other eigenvectors in Equation (6) because the ranking is ill-defined. However, from a numerical standpoint the spectrum of the covariance matrix will not be degenerate, so if we simulate the model we can try to apply the procedure regardless.
As we can see from Figure 8 the joint probability does not converge to a Gaussian, even though there is nothing critical about the underlying dynamics. Hence the proposed coarse-graining in momentum space fails for a simple set of conditionally independent binomial variables, albeit it seemed to be the most promising procedure for the supercritical contact process in the vicinity of the critical point. On the other hand, and perhaps not surprisingly, in this model the proposed coarse-graining procedure via maximally correlated variables does work: in fact, since the off-diagonal elements of the covariance matrix are all equal we are randomly pairing neurons together and no scaling property emerges. This simple model shows once more how careful one should be when employing these kind of procedures. All in all, the two approaches combined seem to work well, in the sense that a system might really be critical if both of them indicate the presence of an underlying scale invariance. However, when we take them individually they might point in the wrong direction.

APPENDIX C: ISING MODEL
We also test the coarse-graining procedure in the 2D Ising model. We do not show the results explicitly for the sake of brevity, but they are in line with what one might expect.
In the disordered phase at T > T c the coarse-graining drives the system towards a behavior that is comparable with the one of independent random variables, very much like one would expect from a usual block-spin transformation in real space.
For T < T c , instead, the behavior of coarse-grained variables resembles the one of a perfectly ordered system.
However, as we lower the temperature we see a nontrivial effect due to the spontaneous symmetry breaking that occurs at the transition. In fact, the Ising model in its ordered phase is essentially low dimensional [41] in the sense that one single eigenvalue eventually dominates the spectrum of the covariance matrix. Once we take this into account, we find that at criticality its spectrum scales with an exponent µ = 0.88 ± 0.03 which is perfectly compatible with the value µ = 7/8 one gets from the exact solution. Hence, in this case of an equilibrium phase transition with spontaneous symmetry breaking, this procedure does identify two distinct phases [42].