Heterogeneous excitable systems exhibit Griffiths phases below hybrid phase transitions
aa r X i v : . [ c ond - m a t . d i s - nn ] J a n Heterogeneous excitable systems exhibit Griffiths phases below hybrid phasetransitions
G´eza ´Odor (1) and Beatriz de Simoni (2) (1) Institute of Technical Physics and Materials Science,Center for Energy Research, P. O. Box 49, H-1525 Budapest, Hungary(2) Budapest University of Technology and Economics, M˝uegyetem rkp. 3., H-1111 Budapest, Hungary (Dated: January 25, 2021)In d > /t power-law activity decay and other power-laws, thus it is calledmixed-order or hybrid type. It has recently been shown that the introduction of quenched disorderrounds the discontinuity and second order phase transition and Griffiths phases appear. Herewe provide numerical evidence, that even in case of high graph dimensional hierarchical modularnetworks a Griffiths phase in the K = 2 threshold model is present below the hybrid phase transition.This is due to the fragmentation of the activity propagation by modules, which are connected viasingle links. This provides a widespread mechanism in case of threshold type of heterogeneoussystems, modeling the brain or epidemics for the occurrence of dynamical criticality in extendedGriffiths phase parameter spaces. We investigate this in synthetic modular networks with andwithout inhibitory links as well as in the presence of refractory states. PACS numbers: 05.70.Ln 89.75.Hc 89.75.Fb
I. INTRODUCTION
Phase transitions in genuine nonequilibrium systemshave often been investigated among reaction-diffusion(RD) type of models exhibiting absorbing states [1, 2].In many cases mapping to surface growth, spin systemsor stochastic cellular automata have been used. Critical-ity allows us to define universality classes, defined by thescaling exponents, which have been explored in homoge-neous systems [3, 4]. In heterogeneous network modelsthe situation is less clear. Hybrid phase transition (HPT)means that at the transition point the order parameterexhibits a jump, in conjunction with critical phenomenarelated to it. It can mean avalanches of activity at thetransition point with power-law (PL) size distribution forexample. Such type of transitions have been known for along time [5], for example at tricriticality [6, 7], but hadnot been the focus of research and the term appearedlater. HPT-s have been found in network science in caseof k -cores [8], interdependent networks [9] and multi-plexes [10].The ”mixed-order” naming for the same phenomenain statistical physics arouse some years ago [11] bythe exactly soluble one-dimensional Ising model withlong range interactions. It is also known to appear innonequilibrium models, exhibiting transition to absorb-ing states [12]. Further examples include critical mod-els at extended surface defects [13, 14] and synchroniza-tion [15–17].Criticality is an ubiquitous phenomenon in nature assystems can benefit many ways from it. As correlationsand fluctuations diverge [18] in neural systems workingmemory and long-connections can be generated spon-taneously [19] and the sensitivity to external signals ismaximal. Furthermore, it has also been shown thatinformation-processing capabilities are optimal near the critical point. Therefore, systems tune themselves closeto criticality via self-organization (SOC) [20, 21], pre-sumably slightly below to avoid blowing over excitation.Besides, if quenched heterogeneity (that is called disor-der compared to homogeneous system) is present, rare-region (RR) effects [22] and an extended semi-critical re-gion, known as Griffiths Phase (GP) [23] can emerge.RR-s are very slowly relaxing domains, remaining in theopposite phase than the whole system for a long time,causing slow evolution of the order parameter. In theentire GP, which is an extended control parameter re-gion around the critical point, susceptibility diverges andauto-correlations exhibit fat tailed, power-law behavior,resulting in bursty behavior [24], frequently observed innature [25]. Even in infinite dimensional systems, wheremean-field behavior is expected, Griffiths effects [26] canoccur in finite time windows.It is known that strong disorder can round or smearphase transitions [22]. According to the arguments byImry-Ma [27] and Aizenman-Wehr [28], first-order tran-sitions do not exist in low-dimensional disordered sys-tems. It has recently been shown [29, 30] that this istrue in genuinely nonequilibrium systems [1, 4].Experimental and theoretical research provide evi-dence that the brain operates in a critical state betweensustained activity and an inactive phase [18, 31–36].Criticality in general occurs at continuous, second orderphase transitions. On the other hand, meta-stability andhysteresis are also common in the brain behavior. Theyare related to the ability to sustain stimulus-selectivepersistent activity for working memory [37]. The brainrapidly switches from one state to another in response tostimulus, and it may remain in the same state for a longtime after the end of the stimulus. It suggests the exis-tence of a repertoire of meta-stable states. There havebeen several model describing this [38, 39]. It introducesan apparent contradiction, because meta-stability andhysteresis occur in general at first order, discontinuousphase transitions. But the brain can operate at differentregimes close to the critical point which can provide thedesired advantages for biological systems. Another pos-sible resolution for the above controversy is the operationat a transition of hybrid type. It has also been suggestedin a recent theoretical work [40].Threshold type of systems, like the integrate and firemodels of the brain [41], are also suggested to describeother phenomena, like power-grids [42–44], crack andfracture formation [45], contagion [46], etc. In these mod-els HPT can emerge naturally, thus the present resultscan also be relevant.Heterogeneity effects are very common in nature andresult in dynamical criticality in extended GP-s, in caseof quasi static quenched disorder approximation [47].This leads to avalanche size and time distributions, withnon-universal PL tails. It has been shown within theframework of modular networks [47–49] and a large hu-man connectome graph [50]. In this study we re-use thehierarchical modular network of [48] and provide nu-merical evidence that above the GP a HPT emerges.Meta-stable states and hysteresis behavior can also befound, thus this system can oscillate between up anddown states, depending on the level of oscillations, with-out the need of oscillators at the nodes, as in case ofthe Ginzburg-Landau theory, suggested to model corti-cal dynamics [51]. By extending our model we will showthat the proposed mechanism is very general, providingan explanation for the observed wide range of scale-freebehavior below the transition point. II. MEAN-FIELD APPROXIMATION
Discrete threshold models can be defined as two-statesystems: x i = 0 , i , with a con-ditional activation rule, depending on the sum of activityof neighbors compared to the threshold value K : X j x j w i,j ≥ K , (1)where w i,j is the weight of the link connecting site j to i . In interacting homogeneous systems w i,j is justthe adjacency matrix element: A i,j , which is 1 if nodesare connected or 0 otherwise. To describe stochasticitythis activity creation can be accepted with probability λ , competing with an activity removal process of prob-ability ν . The mean-field description of the thresholdmodel of N nodes can be obtained in a similar way asin case of RD systems [52]. That work is defined on thelattice, but we can apply it for other graphs. In [52] itwas shown that discontinuous jump occurs in mean-fieldmodels of the n > m RD systems, in which n neighbor-ing particles are needed for creation and m neighbors forremoval. Here we don’t have diffusion and particles, butthe activity can be considered as site occupancy and we can map the threshold model with K = 2 to an RD sys-tem with n = 2 active neighbors necessary for creationand m = 1 for spontaneous removal. In the presence ofinhibitions n > ρ , and two active neighboring sites can occurin a ( N − N − / N − N − ρ (1 − ρ ) . (2)Let us call : Λ( N − N − / λ . For a full graph of N nodes we can setup the rate equation dρdt = λρ (1 − ρ ) − νρ , (3)which in the N → ∞ limit provides λ c = 0, but for finitegraphs λ c >
0. In the steady state we have λρ (1 − ρ ) − νρ = 0 . (4)By imposing the condition ν = 1 − λ , (5)we obtain λρ (1 − ρ ) − (1 − λ ) = 0 , (6)which can be solved as ρ = λ ± p λ − λ (1 − λ )2 λ . (7)The solution is real and positive if λ > λ c = 4 / , (8)providing a threshold within a system of size N Λ c = 85( N − N −
2) (9)and an order parameter for Λ → Λ + c ρ c = 1 / . (10)It is important to realize that in the N → ∞ limit Λ c →
0, thus there is no inactive phase in the thermodynamiclimit. But as real systems are always finite sized, we canobserve this hybrid phase transition in them even in themean-field limit. Similar phenomenon has recently beenreported in contagion models [53]. Furthermore, in caseof the presence of inhibitory couplings the HPT at finitetransition rate may survive the N → ∞ limit in highdimensional systems.At this transition point we can determine how the den-sity approaches ρ c : ρ ( t ) − ρ c ∼ t − . (11)Thus here we find PL dynamical behavior even thoughthe transition is discontinuous, as in other known hybridor mixed order transitions. For λ < λ c we have ρ = 0 sta-ble solution and exponentially decaying activity. Rightabove the transition the steady state density vanisheswith a square-root singularity as in case of k-core [8] ormultiplex percolation hybrid transitions [10]:( ρ − ρ c ) ∝ ( λ − λ c ) / , (12)but unlike the contact process [54] near multiple junc-tions [12], or the Kuramoto model with uniform fre-quency distribution [15], which thus belong to anotherhybrid universality class. In the following sections weinvestigate what happens to this HPT if we implementthe threshold model on a hierarchical modular network(HMN). III. THRESHOLD MODEL ON HIERARCHICALMODULAR NETWORKS
In this section we describe the HMN models we usefor the simulations. It is important to note that we be-lieve that hierarchy is not relevant, but that modularity iswhat enhances RR effects as in case of the study [49]. Themodels are motivated by brain networks originated fromKaiser and Hilgetag, who performed numerical studiesto investigate the effects of different topologies on theactivity spreading [55]. Their hierarchical model reflectsgeneral features of brain connectivity at large and meso-scopic scales, where the nodes were intended to representcortical columns instead of individual neurons. The con-nections between them were modeled excitatory, sincethere appears to be no long-distance, inhibitory connec-tions within the cerebral cortex [56].The network was generated beginning with the high-est level and adding modules to the next lower level withrandom connectivity within modules. Kaiser and Hilge-tag explored hierarchical networks with different num-bers of hierarchical levels and numbers of sub-modulesat each level. The total, average node degree was setto a fixed value, motivated by comparative experimentalstudies. However, they investigated different topologiesby varying the edge density across the levels. All thetested HMN-s were small-world type, i.e. exhibited infi-nite topological dimensions.The spreading model they investigated was a two-statethreshold model, in which nodes became activated withprobability λ , when at least K of their directly connectednode neighbors were active at the same time or sponta-neously deactivated, with probability ν . Note that thismodel is very similar to RD models known in statisti-cal physics [3, 4], with a synchronous cellular automaton(SCA) update. Without loss of generality this algorithmproduces faster dynamical scaling results for thresholdmodels than those with random sequential updates.In this paper we investigate versions of HMN-s, whichpossess increasing edge density from top to bottom lev- els. Clearly, such topologies can be expected to be moresuitable for activity localization and for RR effects. i j FIG. 1: Plot of the adjacency matrix of an N = 1024 sizedsample of the HMN2d graph used. Black dots denote connec-tions between nodes i and j . The 4-level structure is clearlyvisible by the blocks along the diagonal, additional long-rangeedges are scattered points around it. One can also make a correspondence with the spatiallyembedded networks [57]. These networks have long links,with algebraically decaying probabilities in the Euclideandistance R as p ( R ) ∼ R − s . (13)We added random long links by level-to-level from topto bottom, similar to in [48]. The levels l = 0 , , ..., l max are numbered from bottom to top. The size of domains,i.e the number of nodes in a level, grows as N l = 4 l +1 inthe case of the 4-module construction, related to a tilingof the 2d base lattice, due to the rough distance levelrelation: p l = b ( 12 s ) l . (14)Here b is related to the average degree h k i of nodes , whichwas prescribed to be h k i = 12 for this construction.We connected nodes in a hierarchical modular way as ifthey were embedded in a regular two-dimensional lattice(HMN2d) as shown by the adjacency matrix on Fig. 1,similarly as in [48]. The 4 nodes of the level l = 0 werefully connected. The single connectedness of the net-works is guaranteed by additional linking of these 4-nodemodules, by two edges to the subsequent ones: the firstand the last nodes of module (i) to the first node of mod-ule (i+1). Accidental multiple connections were removedand self-connections were not allowed. Note that the sin-gle connectedness at low level does not result in stablesteady states in case of the threshold value K = 2.The in-degree distribution of 4 randomly selectedgraphs with N = 4096 nodes can be seen on Fig. 2. The k in c oun t FIG. 2: In-degree distribution of 4 randomly selected l = 5HMN2d graphs. lowest in-degree is always k ini ≥
5. The modularity quo-tient of the networks is high:
Q > .
9, defined by Q = 1 N h k i X ij (cid:18) A ij − k i k j N h k i (cid:19) δ ( g i , g j ) , (15)where A ij is the adjacency matrix and δ ( i, j ) is the Kro-necker delta function. The Watts-Strogatz clustering co-efficient [58] of a network of N nodes is C = 1 N X i n i /k i ( k i − , (16)where n i denotes the number of direct edges intercon-necting the k i nearest neighbors of node i . C = 0 . C r = 0 . C r = h k i /N . Theaverage shortest path length is defined as L = 1 N ( N − X j = i d ( i, j ) , (17)where d ( i, j ) is the graph distance between vertices i and j . In case of this typical network L = 6 .
74, about twicelarger than that of the random network of same size: L r = 3 . L r = ln( N ) − . h k i + 1 / . (18)So this is a small-world network, according to the defini-tion of the coefficient [60]: σ = C/C r L/L r , (19)because σ = 5 .
363 is much larger than unity.We estimated the effective topological dimension usingthe breadth-first search (BFS) algorithm: d = 4 . N ( r ) ∼ r d , where we counted the numberof nodes N ( r ) with chemical distance r or less from the seeds and calculated averages over the trials. Note, thatthis is just an estimate for the finite sized graph, becausewe know that d → ∞ is expected for s = 3. It renders thismodel into the mean-field region, because for thresholdmodels the upper-critical dimension is d c ≤
4. Still, dueto the heterogeneous structure, we find very non-trivialdynamical GP scaling behavior as will be shown in thefollowing sections.
IV. DYNAMICAL SIMULATIONS
Time dependent simulations were performed for singleactive seed initial conditions. It means that a pair ofnodes is activated at neighboring sites: x i = x i +1 = 1,in an otherwise fully inactive system. It can triggeran avalanche, a standard technique in statistical physicsto investigate critical initial slip [2]. We measured thespatio-temporal size s = P Ni =1 P Tt =1 x i and the durationof the avalanches ( T ) for tens of thousands of randominitial conditions: both initial sites and initial graph con-figurations. The graphs we investigated had l = 4 , , , N = 1024 , , , h k i ≃
12, afterthe low level linking and the removal of accidental mul-tiple edges. The ratio of short and long links was ≃ . ν = 1 − λ and updated the sites at discretetime steps, i.e. set the state variables x ′ ( i ) = 1 if it wasinactive x ( i ) = 0 and the sum of active neighbors P j x ( j )exceeded K = 2 with probability λ or to x ′ ( i ) = 0 withprobability 1 − λ if it was active x ( i ) = 1. Following a fullsweep of sites we wrote x ( i ) = x ′ ( i ) for all nodes, corre-sponding to one Monte Carlo step (MCs), throughout thestudy we measure time in MCs units. We have measuredthe density of active nodes ρ ( t ) = 1 /N P Ni =1 x i . A. Excitatory model
The simulations were run for T = 10 MCs, or untilthe system goes to a fully inactive state, corresponding tothe end of the avalanche. We computed the probabilitydensity functions of avalanche sizes p ( s ) and final survivaltime distributions p ( t ). We repeated these simulationsfor different λ branching rates, by increasing its value. AsFig. 3 shows we don’t see exponential decays as shouldbe in the inactive phase of a mean-field model. Instead,there are PL-like tails for λ > .
31, modulated slightlyby oscillations, which is a well known phenomenon whendiscrete spatial periodicity is present, here the size of themodules. The slopes of the PL tails vary from τ = 2 . τ = 1 .
39 as we increase λ from 0 .
315 to 0 . λ = 0 .
315 with δ = 1 . λ = 0 .
33 with δ = 0 . s −10 −8 −6 −4 −2 p ( s ) −1.39(1) s −1.56(1) s −2.02(1) FIG. 3: Avalanche size distributions at different λ branchingrates, denoted by the symbols, in the presence of excitatorylinks in the HMN2d with l = 5 , λ = 0.33, 0.325, 0.322, 0.32 ( l = 5 cyan and l = 6green), 0.315, 0.31. Dashed lines show PL fits for the tails: s > λ = 0.315, 0.322, 0.33. integral, thus δ is related to the duration exponent of P ( t ) ∝ t − τ t as τ t = 1 + δ . (20)These non-universal PL-s suggest that Griffiths effectsare present, as reported in [48] for this model at differ-ent parameters. By repeating the simulations at differ-ent sizes: l = 5 , t −5 −4 −3 −2 −1 p ( t ) −0.172(1) t −0.41(1) t −0.73(1) t −1.80(1) FIG. 4: Survival probability of the activity at differentbranching rates in the K = 2 threshold model with excita-tory links. From top to bottom curves: λ = 0.33, 0.325,0.322, 0.32 ( l = 5 and l = 6), 0.315. Dashed lines show PLfits for the tails: s > at λ = 0.315, 0.322, 0.33. The seminal experiments by Beggs and Plenz [31] re-ported neuronal avalanches with size ( s ) dependence,defined as either the number of electrodes with supra-threshold activity or the sum of the potentials, accordingto a power law, p ( s ) ∝ s − . . For the duration distribu-tion of such events P ( t ) ∝ t − , PL tails were observed.These exponents are in agreement with the mean-field(MF) exponents of the Directed Percolation (DP) crit-icality: τ = 3 / τ t = 2 see [3]. Mean-field exponentsare expected to occur if the fluctuation effects are weak,when the system dimension is above the upper criticaldimension d c .On the other hand, Palva et al. [61] have found thatsource-reconstructed M/EEG data exhibit robust power-law, long-range time correlations and scale-invariantavalanches with a broad range of exponents: 1 < τ < . . < τ t < .
4. These broad range exponent resultshave also been found in a recent cortical electrode ex-periment study on rodents [62]. An obvious explanationfor this wide spread of critical exponents can be hetero-geneity, which in the GP causes non-universal dynamicalexponents [48, 63–65].
B. Inhibitory model
Although inhibitory links are not expected at longrange links of the brain [56], we believe that our syn-thetic model may describe smaller cortical scales as well.Besides, inhibitory mechanisms can occur in other phe-nomena with threshold dynamics. In case of power-grids,for example, feedback is applied to prevent catastrophicblackout avalanches, or in models of social/epidemic con-tagion, nodes with inhibitory properties may also exist.For simplicity we modeled the inhibitions by the intro-duction of links with negative weight contribution ( w i,j )in the threshold comparison rule given by Eq. 1, althoughwe think our results are easily transferred to the case ofinhibitory nodes. As in [64, 65], we randomly flippedthe sign of 20% of the edges after the generation of thenetwork.The same analysis resulted in similar behavior as forthe excitatory case. One can see p ( s ) distributionswith non-universal PL-s ranging from λ = 0 . τ = 1 . λ = 0 .
55 with τ = 1 . N = 4096 to N = 16384.Usually it is believed that overlapping avalanches dis-tort the scaling behavior. From the point of view of sta-tistical physics, this would contradict universal asymp-totic scaling behavior. Indeed we can see the same cumu-lative p ( s ) distribution tails even in case of starting thesystem from half filled active state as shown on the insetof Fig. 5. The only difference is that the tails are shiftedto larger s values following an initial growth, which mightnot be observed in case of short time measurements.The p ( t ) decays show GP behavior from λ = 0 . δ = 1 . λ = 0 .
52 with δ = 0 . s −10 −8 −6 −4 −2 p ( s ) −1.168(1) s −1.651(5) s −4 −2 p ( s ) FIG. 5: Avalanche size distributions at different λ branchingrates, denoted by the symbols, in the presence of inhibitorylinks in HMN2d with l = 5 , λ =0.55, 0.54, 0.53 ( l = 5 green and l = 6 cyan),0.52, 0.51, 0.50 ( l = 5 triangle and l = 6 diamond). Dashedlines show power-law fits for the tails of λ = 0 . , . t > λ = 0 . , . , . , .
525 (bottomto top symbols).
These values do not correspond to the ends of the GP,we did not aim to determine them precisely as the expo-nents are non-universal. Furthermore, as we will show inSect. V the determination of the upper limit of the GP,corresponding to the critical decay is almost impossibleby numerical simulations. t −5 −4 −3 −2 −1 p (t) −1.10(3) t −0.70(3) FIG. 6: Survival probability of the activity at differentbranching rates and ν = 1 − λ for the K = 2 threshold modelwith levels: l = 5 , λ = 0.5, 0.505, 0.510, 0.515,0.520 ( l = 5 purple cross and l = 6 blue circle), 0.525 ( l = 5brown cross and l = 6 brown circle). Again the τ and the τ t = 1 + δ values lie within therange obtained by experiments. C. Inhibitory-refractory model
Finally, we extended the inhibitory case study withthe possibility of refractory states. Refractorieness meansthat, following an activation, nodes cannot fall back im-mediately to inactive state on the next update, insteadthey stay for time ∆ t in a refractory state. Thus theycannot be reactivated by the neighbors they excited. Thisrefractoriness is generic in excitable systems and has beenused in most of the neural studies [18, 66, 67]. One of theconsequences of refractoriness is to induce oscillatory dy-namics if ∆ t is large enough and the spreading propertiesresemble to annular growth, corresponding to Dynami-cal Isotropic Percolation (DIP) [4]. However, real DIPoccurs if re-activation is not possible, i.e. in the limit∆ t → ∞ , and in high dimension the avalanche scalingexponents of DIP are the same as those of DP [4, 68].Thus analytic studies or simulations in high dimensionsdo not show differences. In the extensive GP simulationswe used ∆ t = 1, but on the inset of Fig. 7 we show os-cillatory activity behavior of a single run for ∆ t = 10, λ = 0 . l = 6.In [65] the GP behavior of inhibitory-refractory thresh-old model was investigated on a large human connectomegraph numerically. Non-universal p ( s ) decays were re-ported with 1 . < τ < .
91. Here we can see this inthe range λ = 0 .
39 with τ = 1 . λ = 0 .
43 with τ = 1 . s −10 −8 −6 −4 −2 p ( s ) −1.96(2) s −1.69(2) s −1.555(4) s −1.30(1) t ρ FIG. 7: Avalanche size distributions at different λ branchingrates, denoted by the symbols, in case of the refractory model,in the presence of inhibitory links in HMN2d-s with l = 5 , λ = 0.39, 0.40 ( l = 5 lefttriangle and l = 6 up triangle), 0.41, 0.42, 0.43. Dashed linesare PL fits for the tails of λ = 0 . , . t > ρ ( t ) of a single run for∆ t = 10. The avalanche survival probabilities (see Fig. 8), ex-hibit PL decay from λ = 0 .
40 with δ = 0 . λ = 0 .
43 with δ = 0 . . ≤ τ t ≤ . t −5 −4 −3 −2 −1 p ( t ) −0.92(1) t −0.726(3) t −0.39(1) t −5 −3 −1 ρ FIG. 8: Survival probability of the activity at differentbranching rates λ for the levels: l = 5 ,
6, in case of theinhibitory-refractory model. From bottom to top symbols: λ =0.40, 0.41 ( l = 5 and l = 6), 0.42 ( l = 5 light greenand l = 6 dark green), 0.43. Dashed lines show PL fits for t > λ = 0.4, 0.41, 0.43 cases. Inset: ρ ( t ) at λ = 1, l = 7 averaged over 10 realizations. Blue boxes: excitatory,red diamonds: inhibitory. Black bullets: BFS ρ ( r ) results.Dashed lines are PL fits for the initial regions: 1 ≤ t < d eff = 1 . d eff = 1 . d = 4 . < r < for similar models in [66, 67] complex phase diagramsand non-universal PL-s have also been found and thepossibility of GP has been pointed out. V. STEADY STATE SIMULATIONS
In order to determine the steady state behavior wefirst performed long runs, up to T = 10 MCs, by start-ing the system from fully active state or from randomlyhalf filled activity: ρ (0) = 0 .
5. Fig. 9 shows the resultsfor the excitatory model. At λ = 0 . . ≤ λ < .
33, in agreement with the seedsimulations. At λ = 0 .
33 the density does not saturateto a constant value. Examining it on log.-lin. scale andperforming average over thousands of independent sam-ples it turns out that even the λ = 0 .
34 curve decaysvery slowly. Only for λ ≥ .
35 we can see saturation,corresponding to active steady state, thus, we estimate λ c = 0 . t −3 −2 −1 ρ ( t ) FIG. 9: Evolution of ρ ( t ) for different λ -s in case of startingfrom fully active state in the excitatory model with levels: l = 5 ,
6. From bottom to top symbols: λ = 0.30, 0.32, 0.321( l = 6), 0.322, 0.322 ( l = 6), 0.325, 0.33, 0.34 (l=6), 0.35, 0.4,0.5, 0.6, 0.7. of λ c to higher values as the consequence of the modelmodifications. We show the results for the inhibitorynetwork on Fig. 11. Again, slow activity decays were ob-served, ending up with visible PL tails for 0 . ≤ λ ≤ .
54, while saturation starts from λ c ≥ . ρ c = 0 . . < λ < .
80 the curves do notsaturate up to T = 10 MCs, neither reach a scaling re-gion. They belong to the inactive phase, but it is verydifficult to distinguish them from other (logarithmic) de-cay forms. λ ρ inhibitory−refractoryinhibitoryexcitatory t00.20.40.60.81 ρ FIG. 10: Steady state behavior for the excitatory, inhibitoryand refractory-inhibitory cases. Inset: evolution of ρ in aninhibitory HMN2d with N = 4096 for different initial activitydensities: ρ (0) = 0.0005, 0.001, 0.01, 0.1, 1 (bottom to topcurves). We can see large jumps at the transition points in allcases, suggesting a discontinuous transition above theGP. It is very hard to locate the exact location of thetransition points as stability disappears very slowly. Thissuggests that at the critical point logarithmic decay oc-curs like in case of the disordered DP [22]. t −4 −3 −2 −1 ρ ( t ) FIG. 11: Evolution of ρ ( t ) for different λ -s, shown by legendsin case of starting from active states in the inhibitory model.Thick, normal, thin lines correspond to l = 7, l = 6, l = 5,accordingly. From bottom to top curves: λ = 0.51, 0.52, 0.53,0.54, 0.55, 0.56, 0.6, 0.78, 0.8, 0.85, 0.9, 0.999. The graphshows results with initial condition ρ (0) = 0 . λ ≤ . λ = 0 . l = 7. In all the other cases ρ (0) = 1 isapplied. We have also tried to start from other initial condi-tions than the full and single seed ones. As the inset ofFig. 10 shows we can see different saturation values for ρ (0) < .
1, that means we have multi-stability in case oflow initial densities. This is the consequence of the factthat for low ρ (0)-s only parts of the graph can be acti-vated. Even though the networks are simple connectedand the lowest in-degree is k ini = 5, not all nodes have 2incoming links from the same neighbor, which is neces-sary for the activation. These nodes cannot be activatedby a neighbor if they are on the ”border” of an active ter-ritory, thus the graphs are practically fragmented fromthe activity point of view. This provides a mechanism forthe emergence of GP even in high dimensions, withoutbreaking conjecture provided in [47], according to whichGPs and similar RR effects do not exist in networks withan infinite topological dimension.Furthermore, we can see the emergence of discontinu-ous transition with multi-stable states, which can be con-sidered bi-stable, for initial excitation with node fraction ρ (0) > . VI. CONCLUSIONS
In conclusion we provided numerical evidence thatstrong heterogeneity effects in networks, coming from themodular structure can result in GP even if the topolog-ical dimension is high, where mean-field scaling wouldbe expected. This is the consequence of fragmented ac-tivity propagation caused by the modular topology andthe threshold. We can define effective dimensions of thethese graphs by running seed simulations with λ = 1, ν = 0 and measuring ρ ( t ) ∼ t d eff . For this compactgrowth ρ ( t ) ∼ N ( r ), so d eff provides an estimate for thedimensionality, similarly to the BFS algorithm. This isreminiscent of similar methods, for instance computingthe spectral dimension of a network from random walksimulations [69–72]. While the topological dimension isa purely structural measure, d eff , as well as the spec-tral dimension are observables of processes operating onnetworks, providing insights to dynamical signatures oflocalization, slowing down and dynamical fragmentation.However, it has recently been shown that in models ofHMN-s the spectral dimension is not defined [73], thusour d eff can be a candidate to clarify relation of struc-ture and slow dynamics. The inset of Fig. 8 shows that aninitial scaling can be fitted for the excitatory case with: ρ ( t ) ∼ t . , while for inhibitory: ρ ( t ) ∼ t . . Thusthese effective, activity dimensions are less than d c , muchsmaller than the topological dimension obtained by theBFS, which is also shown on the graph as the function of r . Furthermore, the threshold type models allow for thepossibility to observe hybrid phase transitions, where or-der parameter discontinuity and multi-stability can co-exist with dynamical scaling in a GP, thus they can modelbrain criticality as well as up/down states. External ac-tivation can then push the model among the multi-stablestates if it is poised near the transition point.The investigated K = 2 discrete threshold model re-sults can obviously be extended for higher K values andwe expect to find similar behavior in continuous, inte-grate and fire type models on modular networks. Con-versely, by duplicating the links we get back effectivelythe contact process [54] without RR effects and GP. Forneural systems our results imply that the functional andstructural connectivity can be different. The effects of in-hibition and refractive states have also been studied andemergence of oscillatory states have been shown. Ourmodel results are applicable to a wide range of phenom-ena, like power-grids, crack and fracture dynamics andcontagion. Acknowledgments
We thank R´obert Juh´asz for useful comments and dis-cussions, Miguel Mu˜noz and Dante Chialvo for their valu-able feedbacks. Support from the Hungarian NationalResearch, Development and Innovation Office NKFIH(K128989) is acknowledged. We thank access to the Hun- garian National Supercomputer Network. [1] J. Marro and R. Dickman,
Nonequilibrium Phase Tran-sitions in Lattice Models , Al´ea-Saclay (Cambridge Uni-versity Press, 2005), ISBN 9780521019460.[2] M. Henkel, H. Hinrichsen, and S. L¨ubeck,
Non-equilibrium phase transition: Absorbing Phase Transi-tions (Springer Verlag, Netherlands, 2008).[3] G. ´Odor, Reviews of Modern Physics , 663 (2004).[4] G. ´Odor, Universality in nonequilibrium lattice systems:Theoretical foundations (World Scientific, 2008).[5] J. L. Cardy, Journal of Physics A: Mathematical andGeneral , 1407 (1981).[6] H.-K. Janssen, M. M¨uller, and O. Stenull, Physical Re-view E , 026114 (2004).[7] W. Chan, F. Ghanbarnejad, and P. Grassberger, NaturePhysics , 936– (2015).[8] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes,Physical Review Letters , 040601 (2006).[9] D. Lee, S. Choi, M. Stippinger, J. Kert´esz, and B. Kahng,Physical Review E , 042109 (2016).[10] G. J. Baxter, S. N. Dorogovtsev, J. F. F. Mendes, andD. Cellai, Physical Review E , 042801 (2014).[11] A. Bar and D. Mukamel, Physical Review Letters ,015701 (2014).[12] R. Juh´asz and F. Igl´oi, Physical Review E , 022109(2017).[13] P. Grassberger, Physical Review E , 010102(R) (2017).[14] R. Juh´asz and F. Igl´oi, Physical Review E , 012111(2018).[15] D. Paz´o, Physical Review E , 046211 (2005).[16] J. G´omez-Garde˜nes, S. G´omez, A. Arenas, andY. Moreno, Physical Review Letters , 128701 (2011).[17] B. C. Coutinho, A. V. Goltsev, S. N. Dorogovtsev, andJ. F. F. Mendes, Physical Review E , 032106 (2013).[18] A. Haimovici, E. Tagliazucchi, P. Balenzuela, and D. R.Chialvo, Physical Review Letters , 178101 (2013).[19] J. J. Johnson S., Torres and J. Marro, PLoS ONE ,e50276 (2013).[20] P. Bak, C. Tang, and K. Wiesenfeld, Physical ReviewLetters , 381 (1987).[21] D. R. Chialvo, Nature Physics , 744 (2010).[22] T. Vojta, Journal of Physics A: Mathematical and Gen-eral , R143 (2006).[23] R. B. Griffiths, Physical Review Letters , 17 (1969),ISSN 0031-9007.[24] G. ´Odor, Physical Review E , 042102 (2014).[25] M. Karsai, H.-H. Jo, and K. Kaski, SpringerBriefs inComplexity (2018), ISSN 2191-5334.[26] W. Cota, S. C. Ferreira, and G. ´Odor, Physical ReviewE , 032322 (2016).[27] Y. Imry and S.-k. Ma, Physical Review Letters , 1399(1975).[28] M. Aizenman and J. Wehr, Physical Review Letters ,2503 (1989).[29] P. Villa Mart´ın, J. A. Bonachela, and M. A. Mu˜noz,Physical Review E , 012145 (2014).[30] P. Villa Mart´ın, M. Moretti, and M. A. Mu˜noz, Journalof Statistical Mechanics p. P01003 (2015).[31] J. Beggs and D. Plenz, Journal Neuroscience , 11167 (2003).[32] C. Tetzlaff, S. Okujeni, U. Egert, F. W¨org¨otter, andM. Butz, PLoS Computational Biology , e1001013(2010).[33] G. Hahn, T. Petermann, M. N. Havenith, S. Yu,W. Singer, D. Plenz, and D. Nikoli´c, Journal of Neu-rophysiology , 3312 (2010), pMID: 20631221.[34] T. Ribeiro, M. Copelli, F. Caixeta, H. Belchior,D. Chialvo, M. Nicolelis, and S. Riberio, PLoS ONE ,e14129 (2010).[35] A. Ponce-Alvarez, A. Jouary, M. Privat, G. Deco, andG. Sumbre, Neuron , 1446 (2018), ISSN 0896-6273.[36] M. A. Mu˜noz, Reviews of Modern Physics , 031001(2018).[37] D. Durstewitz, J. K. Seamans, and T. J. Sejnowski, Na-ture Neuroscience , 1184 (2000).[38] A. Kaltenbrunner, V. G´omez, and V. L´opez, NeuralComputation , 3011 (2007).[39] S. Scarpetta, I. Apicella, L. Minati, and A. de Candia,Physical Review E , 062305 (2018).[40] V. Buend´ıa, P. Villegas, R. Burioni, and M. A. Mu˜noz, Hybrid-type synchronization transitions: where marginalcoherence, scale-free avalanches, and bistability live to-gether (2020), arXiv:2011.03263.[41] O. Kinouchi and M. Copelli, Nature Physics , 348(2006).[42] I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. New-man, Chaos: An Interdisciplinary Journal of NonlinearScience , 026103 (2007).[43] B. Sch¨affer, D. Witthaut, M. Timme, and V. Latora,Nature Communications , 1975 (2018).[44] G. ´Odor and B. Hartmann, Physical Review E , 022305(2018).[45] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi, Ad-vances in Physics , 349–476 (2006), ISSN 1460-6976.[46] S. Unicomb, G. Iniguez, and M. Karsai, Scientic Reports , 3094 (2018).[47] M. A. Mu˜noz, R. Juh´asz, C. Castellano, and G. ´Odor,Physical Review Letters , 128701 (2010).[48] G. ´Odor, R. Dickman, and G. ´Odor, Scientific Reports , 14451 (2015).[49] W. Cota, G. ´Odor, and S. C. Ferreira, Scientific Reports , 9144 (2018), ISSN 2045-2322.[50] M. Gastner and G. ´Odor, Scientific Reports , 27249(2016).[51] S. Di Santo, P. Villegas, R. Burioni, and M. Mu˜noz,Proceedings of the National Academy of Sciences of theUnited States of America , E1356 (2018).[52] G. ´Odor, Physical Review E , 056114 (2003).[53] W. Choi, D. Lee, J. Kert´esz, and B. Kahng, PhysicalReview E , 012311 (2018).[54] T. E. Harris, The Annals of Probability , 969 (1974).[55] M. Kaiser and C. Hilgetag, Frontiers in Neuroinformatics (2010).[56] P. E. Latham and N. S., Neural Compuping , 1385(2004).[57] M. Barthelemy, Comptes Rendus Physique , 205–232(2018), ISSN 1631-0705. [58] D. J. Watts and S. H. Strogatz, Nature , 440 (1998),ISSN 1476-4687.[59] A. Fronczak, P. Fronczak, and J. A. Ho lyst, PhysicalReview E , 056110 (2004).[60] M. D. Humphries and K. Gurney, PLOS ONE , 1(2008).[61] J. Palva, A. Zhigalov, J. Hirvonen, O. Korhonen,K. Linkenkaer-Hansen, and S. Palva, Proceedings of theNational Academy of Sciences of the United States ofAmerica , 3585 (2013).[62] A. J. Fontenele, N. A. P. de Vasconcelos, T. Feliciano,L. A. A. Aguiar, C. Soares-Cunha, B. Coimbra, L. DallaPorta, S. Ribeiro, A. J. a. Rodrigues, N. Sousa, et al.,Physical Review Letters , 208101 (2019).[63] P. Moretti and M. A. Mu˜noz, Nature Communications , 2521 (2013).[64] G. ´Odor, Physical Review E , 062411 (2016).[65] G. ´Odor, Physical Review E , 012113 (2019). [66] S. Moosavi, A. Montakhab, and A. Valizadeh, ScientificReports , 7107 (2017).[67] M. Zarepour, J. I. Perotti, O. V. Billoni, D. R. Chialvo,and S. A. Cannas, Physical Review E , 052138 (2019).[68] M. A. Mu˜noz, R. Dickman, A. Vespignani, and S. Zap-peri, Physical Review E , 6175 (1999).[69] R. Burioni and D. Cassi, Physical Review Letter , 1091(1996).[70] R. Burioni, D. Cassi, and A. Vezzani, The EuropeanPhysical Journal B , 665–672 (2000), ISSN 1434-6028.[71] R. Burioni and D. Cassi, Journal of Physics A: Mathe-matical and General , R45 (2005).[72] A. P. Mill´an, J. J. Torres, and G. Bianconi, PhysicalReview E , 022307 (2019).[73] S. Esfandiary, A. Safari, J. Renner, P. Moretti, and M. A.Mu˜noz, Physical Review Research2