Controlling extended criticality via modular connectivity

Criticality has been conjectured as an integral part of neuronal network dynamics. Operating at a critical threshold requires precise parameter tuning and a corresponding mechanism remains an open question. Recent studies have suggested that topological features observed in brain networks give rise to a Griffiths phase, leading to power-laws in brain activity dynamics and the operational benefits of criticality in an extended parameter region. Motivated by growing evidence of neural correlates of different states of consciousness, we investigate how topological changes affect the expression of a Griffiths phase. We analyze the activity decay in modular networks using a Susceptible-Infected-Susceptible propagation model and find that we can control the extension of the Griffiths phase by altering intra- and intermodular connectivity. We find that by adjusting system parameters, we can counteract changes in critical behavior and maintain a stable critical region despite changes in network topology. Our results give insight into how structural network properties affect the emergence of a Griffiths phase and how its features are linked to established topological network metrics. We discuss how those findings can contribute to understand the observed changes in functional brain networks. Finally, we indicate how our results could be useful in the study of disease spreading.


Introduction
The criticality hypothesis states that biological neuronal networks are poised to operate at the critical threshold of a phase transition [1,2,3,4]. It offers an explanation to characteristic scaling of power-law activity dynamics observed in such networks [5,6,7]. This sheds light on the brain's information processing capabilities, as critical operation has been conjectured to optimize computational capability [8,9], information transmission and storage [10,11,12], and signal sensitivity and range [13,14]. While evidence for critical neuronal dynamics has been increasing [7,15,16], an explanation of how the brain could self-regulate at a precise critical point remains an open question [17,18].
It has been shown that certain topological structures present in neuronal networks can cause the emergence of a critical region [19,20], substituting a single critical point, which would relax the necessity for fine-tuning parameters. Quenched disorder in networks can induce rare-region effects, resulting in critical behavior in an entire parameter region below the critical point, i.e., a Griffiths phase [21,22]. Griffiths phases have been observed in synthetic hierarchical modular networks as well as in empirical and biologically inspired networks [20,23,24]. Thereafter, it has been shown that sufficiently heterogeneous modular networks can support a Griffiths phase without hierarchy [25].
This study was driven by the following question: Given a self-regulating system poised at criticality, how would changes in its network topology affect its dynamics? Topological changes have been observed in functional brain networks of individuals in diverse states of consciousness, such as induced by psychedelics or anesthetics, in sleep or in coma [26,27,28,29,30,31,32,33]. The critical properties of a system are strongly determined by its topology [34,35] and changes in topology can lead to an altered critical point. If the topology of a self-regulating system that only operates at criticality is changed, it would be forced to adapt. Maintaining critical operation could be achieved by either adjusting its parameters to the altered critical region, e.g. the rate of activity spread in brain networks, or by modifying its structure to revert the critical region to its previous parameter range.
In this paper, we investigate how topological properties influence dynamical processes in modular networks featuring a Griffiths phase. We study which network features are responsible for the emergence of a Griffiths phase and how one can manipulate its properties. We find a connection between the Griffiths phase width, i.e., the range of system parameter values that lead to power-law decay, and the network's topological properties. In short, we find that the Griffiths phase can persist in a changing topology and its width can be controlled via both intra-and intermodular connectivity. Alterations in the critical region that stem from a change in either connectivity can be counteracted by tuning the opposing structure. We argue that this could provide a mechanism of self-regulation in modular systems that operate at criticality and add to the functional benefits of modularity. Our results give insight into how the structural properties of modular networks lead to the emergence of a Griffiths phase and connect it to established topological metrics. We highlight the importance of low global efficiency and propose that it is a central feature in this context. We suggest further inquiry into other artificial modular networks or real-world networks, such as empirical brain networks.
Finally, we discuss how an altered Griffiths phase could be connected to the topological changes observed in functional brain networks during altered states of consciousness. If consciousness relates to critical operation, could an increase in Griffiths phase width be connected to the reported changes in conscious quality under the influence of mind altering substances? Alongside, we briefly discuss how our results can help in understanding the persistence of a pandemic disease. This paper is structured as follows: In the Methods section, we introduce the modular networks, the epidemiological model and our approach to analyze the Griffiths phase. In the Results section, we visualize the topological disorder in the modular networks. We show how a change in inter-and intramodular connectivity affects the Griffiths phase and topological network metrics. We conclude with a discussion of our results.

Methods
We explore how topological changes influence the Griffiths phase by simulating the Susceptible-Infectious-Susceptible (SIS) propagation model [36,37] on synthetic modular networks. We choose the modular topology introduced in [25] because, to our knowledge, it is the simplest modular structure reported in the literature that induces extensive Griffiths phase effects with different propagation models. The modular networks consist of a loosely connected ensemble of modules and offer a direct way to manipulate intra-and intermodular connectivity individually. An illustration of the networks can be seen in Figure 1.

Constructing the modular networks
We construct the networks by generating and randomly interconnecting M modules of size N , each module drawn from the same power-law degree distribution P intra (k) ∼ k −γ . For simplicity, we generate the networks with modules of equal size and intermodular connections, leading to a random regular intermodular structure of degree k inter . This architecture is referred to as a monodisperse modular network [25]. Detailed instructions to generate the networks are given in Appendix A.1.
We consider networks with M = 1000 modules consisting of N = 1000 nodes each. The minimum degree in each module is k min = 3, and the cut-off is set to k max , corresponding to the average maximal degree k max (N ) ∼ N 1 γ−1 of a power-law network generated by the configuration model [38]. This cut-off leads to a distribution of critical points in individual module realizations, creating topological disorder within the modular networks. A discussion of how the cut-off choice impacts SIS dynamics in power-law networks is given in [39].

Dynamical spreading process
For the activity density decay analysis, we utilize the SIS spreading process. It was originally introduced as an epidemiological model for diseases that do not confer any immunity [40]. A population is compartmentalized into susceptible and infectious members. Infectious members spread a disease to susceptible members with rate λ and recover with rate µ. After recovery an infectious member is again susceptible to reinfection. The SIS model features an absorbing state phase transition: A critical spreading rate λ c separates a stationary from an absorbing phase. Above λ c , the system converges to a stable density of infected/active members ρ. Below λ c , the disease/activity eventually dies out. Due to its minimal assumptions, the SIS process is readily applicable in contexts that go beyond epidemiology, such as the spread of computer viruses or, as in our case, the activity in neuronal networks. In a network model each node represents a member of the population and infected nodes spread activity to every susceptible neighbor node, which results in a high susceptibility to degree variations. The SIS process is a continuous-time Markov chain and its dynamics in a network can be simulated with the statistically exact Gillespie algorithm [41]. In the present study, we use an optimized version of the algorithm that reduces the computational load of the simulation [42]. Our implementation follows the description in [25] and is detailed in Appendix A.2.

Susceptibility
An important quantity in the analysis of complex, coupled systems is the susceptibility. It diverges when a system undergoes a phase transition in dynamical spreading processes and can be utilized to calculate the critical threshold. We consider the susceptibility defined as follows [43,44]: The average activity density ρ was computed via the quasistationary (QS) method, where the system is kept in a QS state by returning it to a previous state whenever the activity dies out [45]. During the initial m time steps of a simulation the state of the system is saved. At each subsequent time step, a randomly chosen saved system state is overwritten by the current state with probability p QS . If the process reaches the absorbing state with N inf = 0, the system is returned to a randomly chosen saved state. We used m = 70 saved states and an overwriting probability of p QS = 0.01, for which the system converges to a QS state. Then, the n-th moment of activity density is estimated by taking the respective temporal average of the steady state where T denotes the observation period, and t is set large enough to discard the initial dynamical transient before the QS state is reached.

Activity density decay analysis
We study the Griffiths phase by performing an activity density decay analysis, as described in the following. Starting with a fully active network, we monitor the density of active nodes ρ over time, averaged over multiple runs and network realizations. We use the spreading rate λ as a control parameter and ρ serves as the order parameter. By exploring a range of λ values we observe an extended region showing power-law decay of ρ(t) in the transition from subcritical absorbing states -characterized by exponential decay -to supercritical steady states. This extended region of slow decay is the signature of a present Griffiths phase [20,25]. The range ∆GP of the power-law decay region is determined by the margin between the critical point λ c , defined as the highest value of λ that does not lead to a steady state, and the lowest spreading rate that shows power-law decay λ low , determined by the topological properties of the modules: Increasing the intramodular connectivity of a modular network to the limit at which it becomes non-modular, its critical point converges to the lower boundary of the Griffiths phase λ c ≈ λ low , annulling any Griffiths phase effects. We therefore identify λ low as the critical point of this non-modular network that has the same structure as a single module, but the size of all modules combined.

Determining λ c and λ low
Susceptibility diverges when a system undergoes phase transition in dynamical spreading processes and can be utilized to calculate the critical threshold [44]. A Griffiths phase is accompanied by an extended region of high susceptibility, which makes this approach challenging. We therefore measure λ c by increasing λ in the density decay simulations until a value is reached that shows a steady state, and is therefore above the critical point. λ c is then taken as the value midway between the first value above the critical point and the last value below it. The range between these rates is taken as the error.
It should be noted that the error is not a standard deviation, since the likelihood of finding the true critical point within the error region does not have a Gaussian profile. At the lower end of the Griffiths phase the power-law decay transitions into exponential decay continuously. To determine a distinct λ low that separates the regions within and below the Griffiths phase, we additionally take into account how the decay behavior is influenced by changes in intermodular connectivity. By increasing the intermodular connectivity, we lower λ c until we reach a non-modular structure that is equivalent to a single module of size M · N and has a critical point λ c . The decay behavior for any λ below λ c is not significantly affected by the change in intermodular connectivity (see supplementary material). However, any λ above λ c , that lies in the Griffiths phase at low k inter , will lead to a steady state when k inter is increased. Therefore λ c = λ low is the natural lower limit of the Griffiths phase.
In short, to determine λ low , we generate non-modular networks of size M · N and determine their critical point via susceptibility peaks. Note that the modular networks can be seen as diluted power-law networks, similar to the diluted Ising lattice, in which the Griffiths phase was originally proposed [21].

Averaging over multiple networks
The intermodular links are assigned at random in each network realization, which leads to a slightly varying λ c and differing decay behavior, when k inter is increased. Above k inter = 10 it is necessary to average over multiple networks to observe consistent powerlaw decay within the Griffiths phase. However, by increasing the number of modules to M = 10 4 we can observe a consistent Griffiths phase in single network realizations up to k inter = 100. If we increase the intermodular connectivity beyond k inter = 100 the decay transitions into the decay of a non-modular power-law network even for very large modular networks. Further details can be found in Figure A2 of the supplementary material.

Topological metrics
The structural properties of the modular networks change when generated in different configurations. To characterize these changes we utilize various topological metrics.
Global efficiency [46] is defined as with total number of nodes N and geodesic distance d(i, j) from node i to j in graph G. It measures how well information can be exchanged within a network. It scales inverse to the characteristic path length and is high in integrated networks with low diameter and low when a network is segregated. We also calculate the geodesic entropy of the networks. Geodesic entropy [47] is a measure for how distributed the geodesic distances for a given node to all other nodes are. We calculate it by first determining the probabilities p i (r) that a given node i ∈ G is of distance d(i, j) = r to a randomly chosen node j ∈Ḡ withḠ = {j|j ∈ G \ {i}} where r lies in the interval 1 ≤ r ≤ r max with r max = max j∈Ḡ (d(i, j)). The geodesic entropy is then given by and the characteristic geodesic entropy by taking the average over all nodes The geodesic entropy allows us to quantify entropic changes due to structural properties, such as altering the intermodular connectivity in the modular networks. This leaves the degree distribution unchanged and can not be measured in the entropy of the degree distribution.
To measure entropic changes connected to the intramodular connectivity, we use the degree entropy Another quantity we evaluate is the extended, local clustering coefficient [48]. It reveals neighbor relations that go beyond direct connectivity: with the set of neighbors N i of node i. It measures the ratio between the number of pairs in N i whose distance is d in G(V \ {i}) and the total number of pairs of neighbors. At d = 1 it returns the standard clustering coefficient. We calculate the extended clustering coefficient with a function from the graph-tool library in python [49].  Figure 2 shows the susceptibility of individual network modules. Individual module realizations have large variations in the number and connectivity of high degree outlier nodes [39]. This leads to varying critical points with the SIS, which can be seen in shifting peaks of dynamical susceptibility. An extended region of high susceptibility is revealed by averaging the susceptibility of an ensemble of modules. This expresses how locally varying dynamics, caused by topological disorder, build the basis of the Griffiths phase. Figure 3a shows how the Griffiths phase width depends on intermodular connectivity. We randomly distribute a fixed number k inter of new intermodular connections per module, while keeping the degree distribution and intramodular connectivity unchanged. The number of added intermodular links is small enough to not significantly change the networks modularity. The increase in k inter leads to a reduction of the Griffiths phase width ∆GP. Figure 3b-d display detailed activity density decay simulations for modular networks with different numbers of intermodular connections (highlighted in Figure 3a). The reduction of ∆GP stems from the reduction of λ c , since λ low remains constant.

Changing the intramodular connectivity
Introducing or removing intermodular links has a consistent influence on dynamic behavior, because of the regular intermodular structure. The influence of intramodular  links depends on where the links are attached to: High-degree nodes have a stronger individual influence on the SIS dynamics than low-degree nodes. We therefore lower intramodular connectivity via two approaches: The first is to reduce each node degree k > 3 by a constant value, maintaining the minimal degree of k min = 3, which creates a shift in the degree distribution. We named this approach the offset method. This method affects every node in the network and since most nodes are of low degree, it changes the intramodular connectivity via the low-degree nodes.
The reduction of intramodular connectivity via the offset method increases λ c , such that the whole curve λ c versus k inter moves upwards (cf. Figure 4). λ low is also increased, but less than λ c which leads to an increase in ∆GP with lowered intramodular connectivity. One can observe that at low intermodular connectivity, λ c and ∆GP are more susceptible to changes in inter-than intramodular links. This behavior can be better understood by considering how topological metrics are affected by the changes in connectivity, in particular global efficiency.
Our second approach to lower the intramodular connectivity is to increase the power-law exponent of the degree distribution γ. This leads to a lower chance to draw high-degree nodes, which reduces connectivity via the outlier nodes of the modules. This method leads to an increase in λ c and λ low at the same rate, which moves the Griffiths phase to a different parameter region (cf. Figure 5b). Both methods increase    Figure 4 and 5a. The networks have intra-and intermodular connectivity tuned contrary to each other and have the same λ c and similar power-law decay regions. A higher offset leads to an increased λ low and decreased ∆GP, a higher k inter leads to faster decay to steady states above λ c and shorter lifetimes below. From left to right the networks global efficiency increases and entropy decreases, as marked in Figure 7 and 8a. Figure 4 shows that network configurations with different connectivity can have equal values of either λ c or ∆GP, if k inter and offset are tuned accordingly. In Figure 6 we see the density decay of three networks with varying topology tuned to equal λ c . The networks at higher offsets have an increased λ low , and therefore reduced ∆GP, but the three networks still have a significant overlap in power-law decay region, despite highly varying global efficiency and both geodesic and degree entropy (cf. Figure 7, 8).

Measuring the topological changes
Increasing intermodular connectivity leads to a decrease in the average shortest path length between nodes and an increase in global efficiency. Figure 7 depicts the relation between global efficiency, intermodular connectivity and the Griffiths phase. The increase of global efficiency due to the increase in intermodular connectivity is particularly strong at low intermodular connectivity, when the modular networks are close to segregation. The reason for that lies in the presence of a percolation phase transition at k inter = 3 [51]. We observe that global efficiency can serve as the order parameter of the percolation phase transition. Close to k inter = 3 global efficiency is therefore highly sensitive to changes in k inter . The inset in Figure 7 shows that this sensitivity extends to λ c and ∆GP, as they scale linearly with global efficiency.
Geodesic entropy is also connected to the average path length and decreases with increasing intermodular connectivity, as the variability of shortest paths decreases (cf. Figure 8a). In the inset plot we can see that in modular networks a higher geodesic entropy correlates with a larger Griffiths phase width. Degree entropy on the other hand is higher with larger intramodular connectivity (cf. Figure 8b).
Despite a reduction in average degree, when intramodular connectivity is reduced,  Figure 7. Global efficiency in modular powerlaw networks with exponent γ = 2.7 with M = N = 10 3 modules and nodes, with increasing intermodular connectivity. At k inter = 3 the networks shift from disconnected rings of modules to a connected small-world, which leads to a drastic reduction in average path length. Global efficiency scales inverse to average path length and increases non-analytically from zero to a finite value. This marks the presence of a percolation phase transition [51] with global efficiency as the order parameter. The coloring distinguishes the values at which extended power-law decay was observed (blue) and where it transitions into exponential decay (green). The density decay at the orange triangles in is shown in Figure 6. local clustering and efficiency remain unchanged. This indicates that local clustering is not responsible for the changes observed in the Griffiths phase. We therefore evaluate the extended clustering coefficient [48], which is a generalization of the traditional local clustering coefficient. While the local clustering coefficient measures connectivity in the direct neighborhood of a node, the extended clustering coefficient can additionally detect clusters of greater distance. Figure 9 shows that the reduction of intramodular connectivity increases the distance of extended clustering and therefore leads to a less clustered structure. A change in intermodular connectivity does not affect extended clustering. On the other hand, intramodular connectivity does not affect overall Extended clustering in modular power-law networks with exponent γ = 2.7 with M = N = 10 3 modules and nodes at k inter = 3 and decreasing intramodular connectivity. The standard averaged local clustering coefficient (d = 1) implies no clustering and remains unchanged with rising degree offset. Yet deeper measures (d = 3, 4) reveal a clustered structure that shifts to a more distant clustering (d = 5, 6) with rising degree offset. efficiency.

Discussion
The detection of a Griffiths phase in complex networks appears to be an important step towards understanding critical behavior in biological systems [20]. An extended critical region in brain-like networks supports the hypothesis that the brain operates at criticality [1] and relaxes the necessity for fine-tuning parameters around a precise critical point. Previous work demonstrates that structural heterogeneity and modularity, both features of human brain networks [52,53,54,55], are sufficient conditions to enable a Griffiths phase [25]. Functional brain networks have been shown to change significantly in different mind conditions, such as sleep, coma, anesthesia or under the influence of psychedelics [47,27,31,32,33]. Studies indicate a neural correlate between brain network integration and consciousness states [29,28,56,57,58] which some quantify using topological network metrics such as the clustering coefficient, local and global efficiency and entropy [26,30,47]. We argue whether the observed changes in topological metrics, especially during psychedelic states, can influence critical behavior by altering the expression of an existing Griffiths phase.
We construct modular networks with varying intra-and intermodular connectivity and explore numerically how the connectivity structure affects the critical behavior of a spreading process and evaluate the link between the observed changes and topological network metrics. We show that an increase in either connectivity leads to a decline of the networks critical point, hence a reduction of the Griffiths phase width. In addition, a decrease in intramodular connectivity leads to a slight increase in the lower limit of the Griffiths phase. Intra-and intermodular connectivity therefore offer two independent ways of controlling critical dynamics and can be used as a tuning mechanism for criticality. If one connectivity structure is changed, we can adapt the opposing structure, leading to networks with differing topology, but identical critical points and similar Griffiths phase regions.
We observe that global efficiency is a central measure in the emergence of a Griffiths phase. The Griffiths phase width scales linearly with global efficiency in the modular networks. Global efficiency captures the reduced information exchange at low intermodular connectivity which enables rare regions effects.
We find that a decreased intermodular connectivity leads to an increase in entropy, even more so when intramodular connectivity is adapted to counter changes in criticality. Recent researches suggests that psychedelics disrupt the hierarchy of brain network topology [59,60] and increase network segregation, indicated by increased shortest path length [26,58] and decreased global efficiency [26]. We argue whether these networks would exhibit a larger Griffiths phase width, similar to the synthetic networks we present here, and whether this increase could explain the enhanced entropy observed in functional brain networks of individuals under psychedelic influence [26,47,56,61].
We hypothesize that the topological changes observed in functional human brain networks during altered states of consciousness may be connected to a mechanism that ensures critical operation. Either by tuning topological parameters to maintain a stable critical region or by extending the critical region. An increased Griffiths phase width via a heightened critical point moves activity rates that were previously supercritical to the critical region. This increases the dynamical parameter range in the networks and could explain the increase of entropy observed during psychedelic states. Further studies should clarify our hypothesis by reproducing these results in networks with more empirical parameters.
We highlight some limitations of our work. The networks are dynamically heterogeneous with the SIS due to the power-law distribution within the modules. Different propagation models, such as the contact process, would require a different intramodular structure to achieve a distribution of critical points among the modules [25] and the structure manipulation would have to be adapted. We expect, however, that the modulation of the critical point is reproducible in any network with distinct intraand intermodular structures and a spreading process that can be separately affected by intramodular properties and changes in average path length created through alterations in intermodular connectivity. Intramodular connectivity can be altered via either low or high degree nodes. Decreasing the degree of each by node by a chosen value affects predominantly low degree nodes and leads to an increase in the Griffiths phase width. This method is only appropriate for small off-sets because large values can distorts the power-law distribution. By increasing the power-law exponent of the degree distribution we focus the connectivity reduction to high degree nodes. This keeps the Griffiths phase width constant and moves the Griffiths phase to a different parameter region. Further research could be focused on how both width and parameter region of the Griffiths phase could be continuously controlled via the intramodular distribution.
An interesting point of further study is if our results are applicable to the spread of infectious diseases, such as the current outbreak of SARS-CoV-2. If we can model the social networks of international communities as sufficiently heterogeneous modular networks, there is a possibility for the emergence of Griffiths effects. This could explain unexpectedly long lifetimes of spreading processes, despite a sub critical spreading rate.

Conclusion
To conclude, we show that the extension of a critical region in modular networks can be controlled via both intra-and intermodular connectivity individually. Changes in critical dynamics that stem from a change in either connectivity can be counteracted by tuning the opposing structure. We find that low global efficiency is central in the emergence of a Griffiths phase and connect our results to the observed changes in functional brain networks in altered states of consciousness.  Figure A2. (a) The peak in susceptibility corresponds to the critical point λ c ≈ 0.10 of a single powerlaw module with γ = 2.7 of size N = 10 6 . (b) We compare how the density decay of modular networks with M = 10 3 modules develops below and above λ low with increasing k inter . Below the lower bound of the Griffiths phase λ c = λ low ≈ 0.10 the lifetime of activity is not significantly affected by the increase in intermodular connectivity. Above λ low the decay gradually shifts from power-law to slow-decay with an exponential cut-off and finally to a super critical steady state when the network resembles a single module. (c) At k inter = 400 we do not see extended power-law decay anymore.