Box-scaling as a proxy of finite-size correlations
Daniel A. Martin, Tiago L. Ribeiro, Sergio A. Cannas, Tomas S. Grigera, Dietmar Plenz, Dante R. Chialvo
BBox-scaling as a proxy of finite-size correlations
Daniel A. Martin , , Tiago L. Ribeiro , Sergio A. Cannas , ,Tomas S. Grigera , , , Dietmar Plenz , and Dante R. Chialvo , Center for Complex Systems & Brain Sciences (CEMSC ),Universidad Nacional de San Mart´ın, San Mart´ın, (1650) Buenos Aires, Argentina Section on Critical Brain Dynamics, National Institute of Mental Health,National Institutes of Health, Bethesda, MD, 20892, USA Facultad de Matem´atica Astronom´ıa F´ısica y Computaci´on,Universidad Nacional de C´ordoba, Instituto de F´ısica Enrique Gaviola (IFEG-CONICET),Ciudad Universitaria. (5000) C´ordoba, Argentina. Instituto de F´ısica de L´ıquidos y Sistemas Biol´ogicos (IFLySiB),CCT CONICET La Plata and Universidad Nacional de La Plata,Calle 59 no 789, B1900BTE La Plata, Argentina Departamento de F´ısica, Facultad de Ciencias Exactas,Universidad Nacional de La Plata, Argentina and Consejo Nacional de Investigaciones Cient´ıficas y T´ecnicas (CONICET), Godoy Cruz 2290, Buenos Aires, Argentina (Dated: July 17, 2020)The scaling of correlations as a function of system size provides important hints to understandcritical phenomena on a variety of systems. Its study in biological systems offers two challenges:usually they are not of infinite size, and in the majority of cases sizes can not be varied at will.Here we discuss how finite-size scaling can be approximated in an experimental system of fixed andrelatively small size by computing correlations inside of a reduced field of view of various sizes (i.e.,“box-scaling”). Numerical simulations of a neuronal network are used to verify such approximation,as well as the ferromagnetic 2D Ising model. The numerical results support the validity of theheuristic approach, which should be useful to characterize relevant aspects of critical phenomena inbiological systems.
Complex biological phenomena at all levels includ-ing macroevolution, neuroscience at different scales, andmolecular biology are of increasing interest. In somesystems, the origin of such complexity has been tracedto critical phenomena via models and theory [1–11].Nonetheless, the connection between complexity and crit-icality still needs to be established carefully in each case.Among others, a very distinctive indicator of the presenceof critical phenomena is the observation of an increase inthe correlation length as a function of the size of the sys-tem under study [12–15]. Such observation exposes oneof the hallmarks of criticality: a complex dynamic whichlacks a characteristic scale. Less evident but equally rel-evant is the fact that at criticality the only scales are theones “imposed” from outside, i.e., the finite size of thesystem and the limited time of system observation.In most physical systems, criticality can be studiedthrough the variation of system properties as some ex-ternal parameter (say temperature) is changed. Howeverthis kind of tuning is usually off-limits in biological sys-tems. Alternatively, one can in principle establish thelack of an intrinsic scale by demonstrating the size-scalingdirectly, i.e. by studying the correlation function of sys-tems of increasing size. While this can be done withrelative ease in numerical studies, it is much harder toachieve in experiments [16]. This is especially true inbiological systems, which in most cases can neither becut in small pieces, nor can they easily be enlarged. Thebrain is a prototypical biological system for which criticaldynamics has been suggested to hold the key to its corefunctions. [3, 4, 6, 10, 11]. The correlations reported for ongoing and evoked brain activity have been foundto depend in a unique way on the size of the observationwindow [4, 11, 17]. However, there has been no system-atic analysis whether and how system-size scaling can beapproximated by varying the size of an observation win-dow (without changing a system’s size). We will refer tothis latter approach as “box-scaling”, since it resemblesthe fractal box-counting algorithm [18].The letter is organized as follows. First the connectedcorrelation function (CCF) is defined. Then the CCF isstudied for a neuronal network model under two scenar-ios: in the first one we proceed in the standard manner,increasing the system size and determining its correlationbehaviour. In the second setting, the CCF is examinedusing a fixed system size (relatively large) while varyingthe size of the field of view (changing the observable win-dow size). After that similar calculations are described inthe ferromagnetic 2D Ising model. The letter concludeswith a summary of the main results.
Correlation analysis:
The connected correlation func-tion measures how a local quantity loses spatial correla-tion as distance is increased [19]. Here we compute it as[8, 15] C ( r ) = 1 c (cid:80) i,j u i u j δ ( r − r ij ) (cid:80) i,j δ ( r − r ij ) (1)where δ ( r − r ij ) is a smoothed Dirac δ function selectingpairs of neurons (or region of interest) at mutual distance r ; r ij is the Euclidean distance from the site i to site j ; u i is the value of the signal s of site i at time t , after a r X i v : . [ c ond - m a t . d i s - nn ] J u l subtracting the overall mean of signals s at that time t , u i ( t ) = s i ( t ) − s ( t ); c is a normalization factor to ensurethat C ( r = 0) = 1.The definition Eq. (1) of the CCF differs from the sta-tistical mechanic definition of CCF (or its normalizedcounterpart, the Pearson correlation coefficient) in an im-portant aspect: the u i are the fluctuations with respect tothe instantaneous space average , s ( t ) = (1 /N ) (cid:80) Ni s i ( t ),as opposed to the ensemble, or phase, average (cid:104) s i (cid:105) , de-fined through an appropriate probability distribution.We also note that the effect of the denominator is tocompensate for the density fluctuations, i.e. the latticespatial structure, thus disentangling the intrinsic fluctu-ations of s i from those due to the discrete distribution ofparticles in space. It is particularly convenient when theeffects of open boundaries are expected to be relevant [8].The correlation length ξ measures the scale at whichtwo points start to become uncorrelated. In a critical sys-tem the correlation length is infinite, meaning that thedecay of the correlation lacks a characteristic scale. Inthis case, the position of the zero of C ( r ) in Eq. 1 is a use-ful length scale, because then it increases proportionallyto the system size L (the CCF defined with instanta-neous space averages substracted must have a zero [8]).This functional dependence, attesting scale invariance,suggests the presence of critical dynamics. However, itcan not be used when the system size is fixed, as in thecase of brain networks. Instead, we can define a charac-teristic scale r defined by the first zero of C ( r ), but withthe correlation function computed for partial regions ofthe entire system. For that the system is subdivided inboxes of side W , (with W < L ) and r determined by C W ( r ) = 0, where C W ( r ) is the CCF Eq. 1 restricted toa box of size W . Thus, the implementation follows thesame logic and limitations than the box-counting algo-rithm commonly used to compute the fractal dimensionof a data set, image or object [18].The hypothesis tested here is that the behavior of r when W is varied with L fixed is the same as wouldbe obtained with W = L and varying L , at least when W (cid:28) L , namely r ∼ (cid:40) ξ log( W/ξ ) , ξ (cid:28) W < L,W, ξ (cid:29) L (cid:29) W. (2)This behavior can be justified for physical systems inequilibrium by extending the arguments of [8, Sec. 2.3.3]to the box-scaling case (see Appendix).In the following we will study the scaling behaviourof the characteristic length r in two models: a 2D neu-ronal network as well as in the 2D ferromagnetic Isingmodel. In all simulations, C ( r ) was measured for all in-teger values of r . Subsequently, the smallest value of r for which C ( r ) is negative, r m , was computed, and r was found as the zero crossing of the linear fit of C ( r )between r = r m − r = r m .
2D Neuronal network model:
We study a neuronal net-work described previously [20, 21]. In short, the model is a cellular automata network on a square lattice, in whicheach neuron can be in one of three states at each timestep: 0 for resting, 1 for active (lasting one time step),and 2 for refractory (lasting two time steps). Each neuronconnects to K other neurons (here always K = 16), inwhich Euclidean nearest neighbour neurons are favouredby an exponential decaying function of the distance r between them (fixed here to r d = 5). To ease the inter-pretation we imposed a cutoff in the interaction probabil-ity by preventing neurons to be connected at distances r greater than a given value (called here interaction length I , fixed here at r = 20).The network overall rate of activity is set by a verysmall ( h = 10 − step − ) independent Poisson perturba-tion to each neuron. We have verified that the resultsare robust over a wide range of h values (e.g., h = 10 − to 10 − step − ). The model includes a control parameter σ = K × T , in which T is the probability that an activeneuron (i.e. in state 1) can excite one of the K neighborsthat are connected to it. Therefore, as shown previously[20], for any given value of K , the model can be madecritical by changing the transmission probability T suchthat σ = 1. We study two different scenarios: in thefirst we compute the CCF of the neuronal activities col-lected from systems of increasing sizes, from L = 40 upto L = 640. In the second case, a system of (fixed) largesize ( L = 1000, i.e., L × L neurons) was simulated, fromwhich the activity of each neuron inside of windows ofsizes smaller than L (from W = 40 up to W = 640) wereextracted for the correlation analysis. In all cases we con-sidered periodic boundary conditions. Results presentedbelow correspond to averages of twenty realizations last-ing 2 . × steps for each parameter value.There are four scale lengths to consider in this model:The first is the interaction length I , it is the scale atwhich neurons interact via direct connections. The sec-ond is the system size L . The third is the size of the ob-servation window ( W < L ) which determines the subsetof neurons selected to compute the CCF. The last scaleis the characteristic length r which will be determinedfrom the CCF.Fig. 1 shows the connected correlation functions com-puted for various system and window sizes, considering s = 1 in Eq. 1, if neuron is active or refractory, and0 otherwise, and three values of the control parameter σ . The results correspond to representative values of thecontrol parameter: sub-critical σ = 0 .
64, super-critical( σ = 1 .
6) and critical ( σ ∼
1) as indicated in the legend.Top panels (A, B and C) correspond to computationsfrom increasing system sizes. Bottom panels (D, E andF) correspond to CCF computed using various windowlength sizes from a system of size L = 1000. It can beseen that the functions obtained when changing systemsize L or changing window size W are qualitatively verysimilar. In particular we note that, as expected fromEq. 2, for both sub-critical and super-critical values of σ (panels A, D and C, F respectively) correlations do notgrow much beyond the model interaction length I = 20. C (r) C (r) L=40L=80L=160L=225L=320L=450L=640 C (r) C (r) C (r) W=40W=80W=160W=225w=320W=450W=640 C (r) A B C r D E F
FIG. 1.
Connected correlation function of the neuronal network model . Curves in panels A, B and C for differentsystem sizes L and those in D, E and F for different window sizes W computed on a system of size L = 1000. Results are forthree control parameter values corresponding to sub-critical ( σ =0.64, panels A and D), critical ( σ =1.024, panels B and E) andsuper-critical ( σ = 1 .
6, panels C and F) regimes of the model. Arrow in panel B illustrates the value of r for L = 640.
10 100 1000System Length (L)0100200 10 100 100010100 σ = 0.64 σ = 0.96 σ = 1.024 σ = 1.12 σ = 1.6 C ha r a c t e r i s t i c Leng t h (r ) C ha r a c t e r i s t i c Leng t h (r )
10 100 10001010010 100 1000Window Length (W)0100200 B CAD E F
FIG. 2.
Characteristic Length r of the neuronal network model: The zero crossings of the CCF shown in Fig.1 areplotted in linear-linear (left), log-linear(middle) and log-log (right) axis. Top three panels correspond to different system size L and bottom three panels to different box length W . Different symbols correspond to the values of the control parameter σ denoted in the legends. Dashed lines are visual aids to emphasize the predicted logarithmic behaviour for both sub-critical andsuper-critical regimes (open circles) and the linear dependence expected for the critical regime (filled circles). Open trianglesare used to denote results obtained for intermediate values of σ . At criticality, on the other hand r increases when eitherthe system itself (panel B) or windows increase in size(panel E).The values of r extracted from the curves in Fig. 1 (aswell as two other σ values) are plotted in Fig. 2. Herethe same data is presented using different axis formatsto visualize the different functional dependency near andaway the critical point. Dashed lines in panels B and E are a guide to the eye illustrating the expected loga-rithmic behaviour of r for sub-critical and super-criticalregimes (open circles). The dashed lines in panels A andD denote the linear dependence expected for r (filledcircles) in the critical regime. Finally, the same data isplotted in log-log axis in panels C and F to reveal thecrossover behaviour for W values close to and smallerthan the interaction length ( I = 20), denoted by the de- C (r) C (r) L=16L=32L=64L=96L=128L=192L=256 C (r) C (r) C (r) W=16W=32W=64W=96W=128W=192W=256 C (r) A B CFED
FIG. 3.
Ferromagnetic 2D Ising model. Connected correlation function:
Typical results for three temperatures T = 2 .
00 (panels A, D); T = 2 .
27 (panels B, E) and T = 3 . L and those in the bottom panels computed on a system of size L = 600 usingthe window sizes W indicated in the legend.
10 100 1000System Length (L)0100
T=2.0T=2.2T=2.27T=2.4T=3.0
10 100 1000101000 300 6000100 C ha r a c t e r i s t i c Leng t h (r ) C ha r a c t e r i s t i c Leng t h (r )
10 100 10001010010 100 1000Window Length (W)0100200AD B CFE
FIG. 4.
Ferromagnetic 2D Ising model. Characteristic Length r . The zero crossings of the computed CCF areplotted in linear-linear (left), log-linear (middle) and log-log (right) axis. Top three panels correspond to different system sizes L and bottom three panels to different window lengths W . Different symbols correspond to the values of the temperature T denoted in the legends. Dashed lines are visual aids to emphasize the predicted logarithmic behaviour for both sub-critical andsuper-critical regimes (open circles), and the linear dependence expected for the critical regime (filled circles). Open trianglesare used to denote results obtained for intermediate values of temperature indicated in the legend. viation from the asymptotic linear dependency for large W . Also it is worth to notice the small deviation of thelinear scaling observed at criticality when the value of W approaches the system size L ( W = 640 in panel D ofFig 2). This deviation is expected from the theory, asdiscussed in the Appendix.Overall, these results show that the scaling of the char-acteristic length r follows a similar functional depen- dence with either the box-scaling or the system-size.
2D ferromagnetic Ising model:
The results obtainedfrom the neuronal model were replicated in numericalsimulations of the ferromagnetic 2D Ising model. Simi-lar to the previous model, the simulations used two sce-narios: in the first the CCF was computed in the stan-dard way from a model running on square lattices of in-creasing sizes from L = 16 up to L = 512. In the sec-ond setup, a relatively large L = 600 square lattice (i.e.600 ×
600 spins) was simulated, and the CCF was com-puted from square window of smaller sizes from W = 16up to W = 512. Results correspond to averages of fiverealizations each one lasting at least 5 × Monte Carlosteps. In all cases we considered periodic boundary con-ditions.Fig. 3 shows representative results for the two scenar-ios and three different temperatures: sub-critical ( T =2 .
0, Panels A and D), critical ( T = 2 .
27, Panels B andE) and super-critical ( T = 3 .
0, Panels C and F). Thetop panels represent results computed for increasing sys-tem sizes and the bottom panels for a fixed lattice sizeand various window sizes. Note that, as already seen inthe simulations of the neuronal model, the computationof the CCF by changing system size L or by changingwindow size W produces very similar results.The dependence of r with system and window sizeis shown in Fig. 4 using the same format as in Fig. 2for the neuronal model. It is clear that the results ob-tained from varying the system or the window size show astriking similarity, suggesting that for this system the ap-proximation is also valid. Notice the small deviation fromlinearity observed at criticality for W = 512 in panel D ofFig. 4 which is similar to that exhibited by the neuronalmodel for W sizes near the value of L .In summary, the results obtained from a neuronal net-work model and from the ferromagnetic 2D Ising modelshow that the finite-size scaling of the correlation length ξ can be approximated —near the small-size limit— bythe dependence of the characteristic length r on windowsize. The results are particularly relevant at the experi-mental level in neuroscience, in which techniques to mapdifferent areas of the brain cortex are now available [22],while changing system size is not feasible. In that di-rection, the present analysis is fully consistent with theexperimental observations being reported in [17]. Acknowledgements:
Work supported by1U19NS107464-01 BRAIN Initiative (USA), ZIAMH00297 of the DIRP, NIMH (USA) and CONICET(Argentina). DAM acknowledges additional supportfrom ANPCyT Grant No. PICT-2016-3874 (AR).
APPENDIX: Box-size scaling in equilibrium.
The relationship between r and W observed in our nu-merical experiments can be understood in the case of aphysical system in equilibrium, where one can relate theCCF computed with space averages (as we do here) to theusual phase-average connected correlation. We show herehow to adapt the argument of [8, Sec. 2.3.3] to the caseof W < L . We start with the relationship between C ( r )and the correlation C ph ( r ), defined as C ( r ) but using thefluctuations with respect to the phase average (cid:104) s i (cid:105) , whichfor W (cid:29) a (where a represents the microscopic lengthssuch as lattice spacing or interaction range) is [8] C ( r ) = C ph ( r ) − (cid:10) [ s − (cid:104) s (cid:105) ] (cid:11) . (3)The variance of s computed over a volume W d can be written in terms of C ph ( r ), (cid:10) [ s − (cid:104) s (cid:105) ] (cid:11) = 1 W d (cid:90) W d d d r C ph ( r ) g ( r ) , (4)where g ( r ) is the radial distribution function describingthe density correlations of the lattice, and arises herebecause the definition of the CCF includes a denomina-tor which is essentially r d − g ( r ). The definition of r is C ( r ) = 0, so that C ph ( r ) = 1 W d (cid:90) W d d d r C ph ( r ) g ( r ) . (5)This equation is useful because, for equilibrium physicalsystems near a critical point, we know the scaling formof C ph ( r ), which we can use to obtain the relationship weseek.We must distinguish two cases: (i) ξ (cid:28) W < L : In this case we can write for r (cid:29) a , C ph ( r ) = r − d +2 − η e r/ξ . At this scale the system ishomogeneous, and we can approximate g ( r ) ≈
1. Clearly r will depend on ξ and W but not on L . Due to the shortrange of C ph ( r ) the integral in Eq. 5 can be extended toinfinity, so that r − d +2 − η e − r /ξ = 1 W d ξ − η (cid:90) ∞ dx x − η e − x , (6)which gives to leading order r ∼ ξ log W. (7) (ii) ξ (cid:29) L (cid:29) W . This is the critical case, where C ph ( r ) = r − d +2 − η h ( r/L ) for r (cid:29) a [16]. For L → ∞ , thescaling function h ( x ) goes to a constant and the decay isa pure power law, but for finite L the decay is modulatedby the scaling function. Plugging into Eq. 5 and usingagain g ( r ) ≈ r − d +2 − η h ( r /L ) = W − d L − η (cid:90) W/L h ( u ) u − η du. (8)If W = L the integral reduces to some constant, and wesee that r ∼ L is a solution, justifying the claim that thezero of C ( r ) is proportional to L when the correlation iscomputed over the whole sample. If W < L (i.e. if C ( r )is computed over a box smaller than the whole system),then in general r will depend on both L and W . How-ever if W (cid:28) L we are in a regime where C ( r ) shoulddecay almost as a pure power law, because the modulat-ing effects of h ( x ) will be noticeable only for r ≈ L . Thismeans that we can replace h ( u ) with a constant insidethe integral, so that r − d +2 − η h ( r /L ) ∼ W − d L − η (cid:90) W/L u − η du ∼ W − d +2 − η , (9)which gives r ∼ W . [1] T. Mora & W. Bialek, Are biological systems poised atcriticality? J. Stat. Phys.
268 (2011).[2] P. Bak,
How nature works: The science of self-organizedcriticality.
Springer Science, New York (1996).[3] D.R. Chialvo, Emergent complex neural dynamics.
Na-ture Physics
744 (2010).[4] A. Haimovici, E. Tagliazucchi, P. Balenzuela, D.R.Chialvo, Brain organization into resting state networksemerges at criticality on a model of the human connec-tome.
Phys. Rev. Lett.
Biochim. Biophys. Acta , 53 (2009).[6] J.M. Beggs & D. Plenz, Neuronal avalanches in neocor-tical circuits.
Journal of Neuroscience , 11167 (2003).[7] Q.Y. Tang, Y.Y. Zhang, J. Wang, W. Wang, D.R.Chialvo, Critical fluctuations in the native state of pro-teins. Phys. Rev. Lett.
Physics Reports
Rev. Mod. Phys.
Frontiers inPhysiology , 15.[11] D. Fraiman & D.R. Chialvo (2012) What kind of noiseis brain noise: anomalous scaling behavior of the restingbrain activity fluctuations. Frontiers in physiology , 307.[12] B. Schmittmann & R.K.P. Zia, Statistical mechanics ofdriven diffusive systems, in: C. Domb and J.L.Lebowitz (eds.), Phase transitions and critical phenomena , Vol. 17,Academic Press, London (1995)[13] U.C. Tuber,
Critical Dynamics,
Cambridge UniversityPress, Cambridge (2014).[14] M.N. Barber, Finite-size scaling, in: C. Domb and J.L.Lebowitz (eds.),
Phase transitions and critical phenom-ena , Academic Press, London (1983).[15] A. Attanasi, A. Cavagna, L. Del Castello, I. Giardina,S. Melillo, L. Parisi, O. Pohl, B. Rossaro, E. Shen, E.Silvestri, M. Viale, Finite-size scaling as a way to probenear criticality in natural swarms.
Phys. Rev. Lett.
Finite-size Scaling,
North Holland,Amsterdam (1988).[17] T.L. Ribeiro, S. Yu, D.A. Martin, D. Winkowski, P.Kanold, D.R. Chialvo, D. Plenz, Trial-by-trial variabil-ity in cortical responses exhibits scaling in spatial cor-relations predicted from critical dynamics. submittedbioRxiv, doi: https://doi.org/10.1101/2020.07.01.182014[18] B. Mandelbrot
The Fractal Geometry of Nature , WHFreeman, San Francisco (1983).[19] N. Goldenfeld,
Lectures on Phase Transitions and theRenormalization Group,
Addison-Wesley (1992).[20] O. Kinouchi & M. Copelli. Optimal dynamical range ofexcitable networks at criticality.
Nat. Phys. (5) 348351(2006).[21] T.L. Ribeiro, S. Ribeiro, H. Belchior, F. Caixeta, M.Copelli. Undersampled critical branching processes onsmall-world and random networks fail to reproduce thestatistics of spike avalanches. PloS one (4) e94992(2014).[22] T. Bellay, A. Klaus, S. Seshadri, D. Plenz, Irregularspiking of pyramidal neurons organizes as scale-invariantneuronal avalanches in the awake state. eLife4