Less is More: Wiring-Economical Modular Networks Support Self-Sustained Firing-Economical Neural Avalanches for Efficient Processing

Brain network is remarkably cost-efficient while the fundamental physical and dynamical mechanisms underlying its economical optimization in network structure and activity are not clear. Here we study intricate cost-efficient interplay between structure and dynamics in biologically plausible spatial modular neuronal network models. We find that critical avalanche states from excitation-inhibition balance, under modular network topology with less wiring cost, can also achieve less costs in firing, but with strongly enhanced response sensitivity to stimuli. We derived mean-field equations that govern the macroscopic network dynamics through a novel approximate theory. The mechanism of low firing cost and stronger response in the form of critical avalanche is explained as a proximity to a Hopf bifurcation of the modules when increasing their connection density. Our work reveals the generic mechanism underlying the cost-efficient modular organization and critical dynamics widely observed in neural systems, providing insights to brain-inspired efficient computational designs.


Introduction
The interplay between structure and dynamics of complex networked systems is always a longstanding topic, covering applications in complex systems from diverse scientific fields. Nowadays, its research and applications in brain and neuroscience are experiencing a rapid growth. arXiv: 2007.02511v3 [physics.bio-ph] 28 Jan 2021 Neurons in human brain form a huge complex dynamical network for efficient functional processing with remarkable cost-efficiency. The principles underlying its efficiency have been actively studied during the past years, either from structure or dynamic aspect.
Brain network is globally very sparse: ~100 billion neurons with ~10 14 synaptic connections each so that the overall density is ~10 -8 in human brain [1]. However, the overall low-density connectivity is organized in a hierarchical manner from local circuits, cortical sheets to whole brain connectome [2,3]. Thus, a prominent feature of brain organization is globally sparse with hierarchical relatively dense modular architectures [4][5][6], which is economical in network wiring since most of the connections are in short-range. There is ample evidence for brain networks to achieve local wiring cost minimization from brain structure [7,8] and a trade-off between global wiring cost and processing efficiency [9,10].
Brain activities consume low energy power of only about 20W, remarkably energy-efficient when compared to digital computers [11]. Dynamically, the irregular and sparse firing of neurons [12] can be collectively organized as oscillations and critical avalanches across different scales [13][14][15][16]. Such 'scale-free' dynamic activities are originally explained by critical branching theory [13], where critical avalanches emerge near the transition point between a silent and an overactive phase. Later experimental evidence [14,17] supports that the transition point between an asynchronous and a synchronous phase better explains the observed critical avalanches, especially in terms of the existence of different critical exponents [14,18,19] satisfying scaling relations [20]. Functionally meaningful avalanche dynamics in critical synchronous transition states enable neurons to fire with low rate [21]. Since cortical metabolic energy usage is dominated by action potentials and synaptic transmission [22][23][24][25][26], avalanche dynamic is also energy economical to maintain the sustained spontaneous (resting) state which consumes the majority of brain metabolic cost [27]. Finally, critical states are functionally beneficial to provide a broad dynamic range to stimulations [28,29] and thus a sensitive standing-by state to respond to constantly changing environments for the brain [30].
Though it is recognized that metabolic cost is a unifying principle governing neuronal biophysics [31], the fundamental mechanism underlying the economical interaction between structures and dynamic modes at the neural circuit level is not well understood. Specifically, how the modular network structure and critical dynamics jointly achieve structural and dynamical optimization for energy-efficient processing? Deciphering these mechanisms are also important for developing braininspired efficient computing. Here we address these questions with biologically realistic neural dynamic model of excitation-inhibition (E-I) balanced [32,33] spiking neuronal networks, clustered on two-dimensional (2D) space to represent the resting-state dynamics on a cortical sheet composed of micro-columns. Interestingly, when rewiring the initial globally sparse random network into locally dense modular networks, the firing rates decreases, and the self-sustained dynamics changes from asynchronous state to critical avalanche state and the response sensitivity to weak transient external stimuli is greatly enhanced. Theoretically, we reveal the enhanced response of neurons by clustered firing and elucidate the dynamic transition via a Hopf bifurcation induced by denser connections within modules during rewiring, through a novel mean-field theory. Overall, our integrative study of cost-structure-dynamics-function relationships in neural networks elucidates that locally dense connectivity under E-I balanced dynamics appears to be the key "less-is-more" solution to achieve cost-efficient organization.

Dynamic transition from random network to modular network
We study a model of = 5 × 10 4 neurons spread on a 2D plane. Considering that other tissue like vessels may separate micro-columns of neurons, neurons are randomly placed on = 100 square regions (modules). Modules are separated by blank space (Fig. 1 top) and each of the modules contains 500 neurons (80% excitatory and 20% inhibitory). Initially, we construct a random network (RN) by randomly connecting each neuron pair with a probability ( = 0.0017 in the main text).
To build a modular network (MN), inter-modular links are rewired, with a probability , into the same module to become intra-modular links. The rewiring method is equivalent to construct smallworld networks, an essential feature of brain networks [34]. Here the essential structural property captured in our model is that the network consists of coupled modules embedded in space [3,5,9]. The voltage (membrane potential) of a neuron in the E-I network is governed by conductance-based (COB) leaky integrate-and-fire (IF) dynamic [35], As the initial RN is rewired into MN, it is interesting to find that the spike-generating and spiketransmission costs of the network are significantly decreased by orders of magnitude ( Fig. 2(a)). Here, the spike-generating cost is defined as the average firing rate, and the spike-transmission cost of a neuron is defined as the product of its firing rate and the total length of its outgoing synapses (the average spike-transmission cost shown in Fig. 2(a) is normalized by the value of RN). These two measures for running costs mimic the energy cost for generating spikes and transmitting spikes [36,37]. Besides, since many links become local short-ranged connections, the normalized wiring cost (defined as the total Euclidian length of all links, normalized by the value of RN) is decreased by orders of magnitude and the connection density within modules increases, approaching to a value ≈ 0.17 ( times of the whole network density , refer to Eq. (4) below). Such smaller wiring length is desirable as the reduced membrane areas of the fibers can reduces metabolic cost and can reduce the transmission delay (although synaptic delay is not considered in our model for simplicity).
The wiring cost reduction is more pronounced for larger networks with more modules (see Fig. S1 and detailed analysis in Supplementary Notes I. Network Setting). Thus, modular network structures reduce both the wiring and running cost.
The initial large and sparse RN can self-sustain asynchronous activity without external input, giving a sustain probability = 1.0, which maintains as the network is rewired into MN ( Fig. 1(a)).
This self-sustained activity [38] resembles the resting states of brain and thus may play functional role.
The results are similar when the overall connection density changes (Fig. S4). However, a denser MN (larger ) with too weak inter-modular connections ( → 1) may not maintain self-sustained activities (Fig. S4). This breakdown of the self-sustainability can be understood later by dynamic analysis of a separate module. Very interestingly, the dynamics modes of networks also co-vary during rewiring. RNs exhibit classical E-I balanced asynchronous state with Poisson-like neuronal spiking [32]. We measure the balance by the net synaptic input current rescaled by the excitatory synaptic current ( / ) averaged over time and neurons. It maintains around 0 when RN is changed into MN ( Fig. 2(a)), suggesting the maintenance of overall balance. The asynchronous state in RN has a noisy fluctuation of the mean voltages around an equilibrium value ( Fig. 3(a) Importantly, MN supports critical neuronal avalanches in modules. Here, the time bin for measuring avalanches is the average inter-spike interval (ISI) of the merged spiking train [13] in a module. When rewiring RN into strongly MN (e.g., =0.995), the avalanche size and duration distribution of a module changes from exponential-decay to power-law ( Fig. 3(b)). Statistical test and estimation of critical exponents are made by accustomed truncation algorithm [39].  [20] approximately holds (error ~ 0.15). The size distribution is fitted into a power-law function [39] and its mean square deviation (MSD) from the fitted curve in Fig. 2(a) bottom shows that modules in MN with ≥ 0.99 has avalanches with power law distributions, exhibiting features of criticality. Other measurements of avalanche in MN associated with the threshold of the average membrane potential are presented in Fig. S5, which also exhibit power-law distribution. This transition from asynchronous spiking to critical avalanches dynamics is the approach to a continuous synchronous transition point, as seen from the increase of CV of activity ( Fig. 2(a)). The selfsustained activity of coupled modules in the critical states provides the ideal scheme that networks can work with low firing rate. The reduced firing rate at criticality is the feature of the critical synchronous transition model [21], whereas traditional branching process model does not exhibit this propertythe firing rate of branching processes at critical states should be larger than that at subcritical states.
Critical states also induce greatly enhanced response sensitivity to transient stimuli. Here, the stimulus is modeled by raising the voltage to −40 of a proportion of the neurons in all modules. We call the stimulus strength. These neurons driven above the firing threshold emit spikes immediately, similar to optogenetic stimulation in experiments. The response sensitivity of a system can be reflected by the returning process of a signal to its baseline value after a transient perturbation.
We measure the response in membrane potential and firing rate of the network modules. The size of response is defined as ∫ | ( ) − | 0 + 0 , which is the area between the signal ( ) (the trialaverage voltage or firing rate of the network) and its resting value , within a window of = 250 , beginning from stimulus onset at 0 (see also details in Supplementary Notes II. Neural Dynamics). As shown by the stimulus-induced trial-average voltage and firing rate of a module ( Fig.   3(c)), the response of MN is much larger and pronounced than that of RN. Interestingly, MN shows a stronger damped oscillation-like response pattern, which is a characteristic of the event-related potential in electroencephalogram signal of brain's response to stimuli [40]. Importantly, both in membrane potential and in firing rate, MN with critical avalanches exhibits response sensitivity much higher than RN with asynchronous spiking activity for small stimulus strength ( Fig. 2(c)).
The above structure-dynamics relationships are robust with respect to the overall connection density (see Fig. S4) and also hold in an extended modeling procedure where the number of intermodular links decay with distance (see Fig. S6). In a word, MN can support cost-efficient critical dynamical modes with greatly enhanced response sensitivity to encode variable input strength, whereas globally sparse RN is both costly in architecture and in running, and cannot properly respond to weak input signals.

Structural correlation and dynamic transition of a single separate module
The key features in the structure-dynamics relationship can be understood from an isolated module separate from the whole network, but subjected to a background excitatory Poisson input train with rate . Here, we use = 50 , to approximate the weak input received by a neuron from other modules in the highly rewired region in the original network. As the rewiring probability increases, local connection within a module become denser ( Fig. 1 and Fig. 2(b)). For a separate module, as its connection density increases, neurons tend to have more common neighbors in the module, the common signal received by a pair of neurons becomes stronger, and their output spikes can be more correlated, as shown in Fig. 4

(a).
This correlation in spiking changes the internal interactions in the network. Fig. 4(b) shows the synaptic current of a randomly selected neuron. For low density, the net input current fluctuates slightly around zero due to strong E-I balance ( Fig. 4(b), upper panels), and the distribution is close to a normal distribution ( Fig. 4(c)). With higher density where spike correlation becomes prominent, correlated excitatory spikes induces quick activation of the network, followed by the activation of inhibitory neurons after an effective delay (due to slower inhibitory synaptic time) and then the activity is depressed. Thus, the net current exhibits oscillations around zero (thus the network maintains E-I balance on average), as shown in Fig. 4(b) lower panels, and its distribution has a large tail at the positive side ( Fig. 4(c)). The dynamic pattern is an alternation between synchronized firing and quiescent state with no spikes ( Fig. 4(b), lower panels). Furthermore, as the module becomes denser, the ability of self-sustaining activity of the module decreases ( Fig. 4(a)). Here, the self-sustainability is tested by tuning off the external inputs, i.e. letting = 0 after initial 200ms. The activity almost cannot sustain when > 0.17, which is the module density in the original MN when → 1 ( Fig. 2(b), also see Eq. (4) below). The weaker sustainable ability of denser network results from clustered firing dynamic mode ( Fig. 4(b), lower panel). The silent period during which no neuron fires increases with ( Fig. 4(a)). If this period is too long, all recurrent inputs drop out and the network activity dies out since there is no external driving.
Under fixed weak external background inputs = 50 Hz, as the density increases, the network dynamics undergoes a transition from asynchronous firing pattern ( Fig. 4(b) upper panel) to critical avalanches ( Fig. 4(b) lower panel) with reduced firing rate ( Fig. 6(b)). The MSD of avalanche size distribution from its best-fitted power-law function in Fig. 4

Effect of correlation: insight from a simplified model
To quantitatively illustrate the impact of input correlation on response sensitivity under E-I balanced dynamic, we can consider a simplified model as follows. A single neuron receives spike inputs from other K=200 Poisson excitatory spike trains. Each of the received spike generates a unit of post-synapse current lasting for = 0.01 . The input signal of the neuron is the summation of these arriving spikes minus a constant equal to the mean current generated by these spikes, in order to mimic the E-I balance. The input correlation is introduced by copying a common Poisson spike train into all input trains [41], see Fig. 5(a) for a paradigm illustration. To construct spike trains with rate and correlation , the common spike train has a rate = and independent spike trains have rates = (1 − ) . Assuming a threshold ( = 20) of the input signal above which the neuron fires a spike, we can numerically obtain the input-output rate response curves (Fig. 5(c)). Compared with independent input trains ( = 0), correlation in inputs induces a positive tail in the distribution of input signal ( Fig. 5(b)), qualitatively capturing the feature of the IF module ( Fig. 4(c)). In the simulation example of Fig. 5(b), the common spike train is copied into a portion of randomly selected input synapses at different time (such that the probability that more synapses receive spikes simultaneously is lower), resulting a decaying positive tail in the distribution of input signal when = 0.05 , which quantitatively resembles the observation from the spiking neural network simulations (Fig. 4(c)). We can see that the correlation increases the output rate when the input rate is the same (Fig. 5(c)). Moreover, this simplified model allows an analytic treatment to explain the effect of the correlation on the response rate (see Methods for details). The theoretical results (black dashed line for = 0 and red solid line for = 0.05 in Fig. 5(c)) fit well to the simulation results of this simplified model.
Hence, the correlation in spikes injected from different recurrent synapses improves the responsiveness of neuron. With input correlation and the response sensitivity, each neuron can maintain generating spikes when the overall firing rate is low.

Mean-field theory of single module dynamics
To further understand the dynamical mechanism underlying the transition of dynamics modes together with the reduction of firing rate, we derive the equations of average neural activity in each module and the interaction among the modules by a novel mean-field technique [19] where , are the average E, I voltages, ( ) = IF dynamics challenges an analytical (self-consistent) estimation of the effective parameters [19].
Taking different fixed , the field equations can qualitatively predict the decay of rate with connection density (Fig. 6(b)). To achieve the best prediction, we numerically estimate the effective parameters through the formula from simulations of the single module IF spiking network to numerically obtain the steady-state mean voltage and mean firing rate of neurons. The results of from modules with different density are shown in Fig. 4(a) bottom panel. Under this setting, the field equations can well quantitatively predict the decrease of firing rate as increases ( Fig. 6(b)). Note that qualitative prediction can already be achieved by fixing the effective parameters value in Eq. (2) (Fig. 6(b)).
Importantly, the field equations reveal that the change of dynamics is associated to a (supercritical) Hope bifurcation. The dominant eigenvalue of the equilibrium in Eq. (2) is complex and its real part approaches zero as increases ( Fig. 6(b)). Thus, the firing rate oscillation emerges through approaching to the Hopf bifurcation under noise-perturbation, which induces the critical avalanches [19]. However, the finite-size effect in a small module (500 neurons) hinders the precision of a meanfield theory. Thus, in the spiking IF model the MSD achieve a minimal at around = 0.17 ( Fig.   4(a)), whereas the field equations do not reach the Hopf bifurcation point and the dynamic is perturbed to bifurcation by noise. Note that a Hopf bifurcation indicates a periodic motion emerges from zero amplitude, corresponding to continuously increasing synchrony as the spiking network ( Fig.   2(a)). The CV of activity, measured by the firing rate series of field equations, grows as the increase of ( Fig. 6(b)). Finally, response size of voltage computed from the field equations ( Fig. 6(b)) also qualitatively predicts the increase of response sensitivity for denser module (examples of = 0.05 and = 0.25 are shown in Fig. 6(a), compared to Fig. 3(c)). This is because when approaching a bifurcation point, the system will respond more sensitively and will take longer time to damp back to the fixed point after perturbation, a phenomenon called critical slowing down [42].

Mean-field theory of the modular network
The above investigation of separated modules with various connection densities under weak external background driving provides understanding of the change of dynamics modes and firing rates with respect to the rewiring probability in the original MN (Fig. 1). First, there is a correspondence between the density in a module , and the rewiring probability , that (refer to Eq. (S1.5)), as shown in Fig. 6 with , Φ , , the corresponding quantities of neurons in the k-th module (see Methods for more details). Thus, the whole MN can be considered as coupled identical neural oscillators.
During the rewiring process, the coupling strength between different modules (~ 1 − ) reduces whereas self-coupling strength (~ 1 + ( − 1) ) increases. In this process, although different modules become less affecting each other, the increase of their internal density significantly shapes the dynamic properties of each module as revealed in separated modules (Fig. 4). Here, the effective parameters , to construct in Eq. (5) depend on the rewiring probability , through their optimal dependence on in separated modules shown in Fig. 4(a) bottom panel and the relationship between and (Eq. (4)). Numerical results in Fig. 6(c) show that the coupled field equations qualitatively predict the decrease of firing rate, increase of CV of activity and response sensitivity to transient stimuli for increasing rewiring probability as observed in the spiking neural model in Fig.   2. Furthermore, Eq. (5) with = 0 also predicts the decrease of nonzero firing rate (Fig. 6(c)), which is a qualitative prediction of the self-sustainability during the rewiring ( Fig. 2(a)). Note, however, that the mean-field analysis here does not capture the effect of changes of input patterns (e.g. increased input correlation) of a module during the rewiring process of the MN. This is a source of prediction errors that lead to the difference between single-module field equations and a module in the coupled-modules field equations when the latter is constructed by the parameters of the former (refer to Fig. S8). An improvement in the future may be made by assuming oscillatory input in the coupled-modules field model where the oscillatory amplitude increases with the rewiring. To conclude, the mean-field theory predicts the dynamical transition (approaching to a Hopf bifurcation) of a module with increasing internal density and this emergent behavior is maintained for the whole MN with mutually coupled modules when rewiring the inter-modular links to intra-modular links.

Conclusion and Discussion
In summary, we have unveiled the principle of neural networks allowing cost-efficient optimization in both structure and dynamics simultaneously. There were many studies on either of the sides, considering the optimization of brain network structure [7][8][9], or the energy-efficient neural dynamics [22][23][24][25][26]. For example, the energy efficient cortical action potentials is facilitated by body temperature [25], and cellular ion channel expression are optimized to achieve function while minimizing metabolic cost of action potential [31]. However, most of previous studies considered the efficiency of network structure (wiring-cost) or the efficiency of dynamics (running-cost) separately. Here, considering both structure and dynamics at the circuit level, we show that wiring-economical modular network can support response-sensitive critical dynamics with much lower running cost while maintaining self-sustainability. This is a notable counterintuitive "less is more" result, because we obtained greatly enhanced functional values with the significant decreases of cost rather than a tradeoff between them.
In our model, the efficiency of activity is achieved by critical avalanche states. Different from traditional critical branching region, the critical dynamics in the synchronous transition region simultaneously achieves greater response sensitivity and lower firing rate. Previous studies showed that critical avalanches can appear under various network topology, for instance, scale-free networks with small-world features [43]. Here we show that locally dense while globally sparse MN is the efficient organization of the network structure that enables both the low global wiring cost and the response-sensitive critical dynamics with low running cost. It would be interesting to further explore its dynamic advantages on specific cognitive tasks such as working memory recalling and decision making.
The origin and mechanism of functionally sensitive critical dynamics in neural systems [13][14][15][16][28][29][30]] is a long-standing, challenging and controversial topic. Considering the physical mechanism that supports such a co-optimization of structure and dynamic, here we reveal that with increasing topological correlation in E-I balanced network, the spike correlation increases, and so does the fluctuation of the inputs received by neurons. In this case, neurons can be activated by lower firing rate, and the network has higher sensitivity. From the perspective of nonlinear dynamics, these features are captured by a novel mean-field analysis, which reduces the whole modular network into coupled oscillators describing the macroscopic dynamics of each module. We elucidate the dynamical mechanism for producing avalanche as the proximity to a Hopf bifurcation in the mean-field. Close to the bifurcation point, the resulting synchronized spikes in each module are temporally organized as critical avalanches. This stronger collective firing rate variability allows greater computation and coding power [44]. In the highly (yet not totally) rewired MN, the sparse inter-modular connections can provide weak external input to a module from other modules. Meanwhile, as modules are dense enough to be around the response-sensitive critical dynamic states, these weak inputs are enough to maintain the whole MN in self-sustained states with low rates.
In principle, the analytical theory for treating biological plausible COB IF neuronal network is still an open question [45]. Our approximation semi-analytical mean-field technique serves as an effective theory to study the macroscopic dynamics of such realistic networks. It is important to stress that our work put several important features of neural systems into an integrated framework. Spatial embedding of neural circuits under the wring cost constraint gives rise to local dense connections and modular organization [7][8][9]. The E-I balance is a fundamental property of the neural circuits [32,33].
Collective activities such as critical avalanches and oscillations are pronounced dynamical features of neural networks [13][14][15][16][28][29][30]. Our modeling and theoretical analysis framework reveals intricate interactions among wiring and running costs, modular network topology, critical avalanche dynamic modes and sensitive response to weak stimulations. Thus, it provides an integrative principle for structural-dynamic cost-efficient optimization in neural systems. Our integrative studies with generic network manipulation and novel mean-field theory with realistic neural dynamics can be extended to coupled cortical areas to offer understanding of critical dynamics across the whole brain [46,47] based on a hierarchical modular connectome [48]. Furthermore, spatial networks with connections to only the nearest neighbors can exhibit propagating waves with critical dynamic properties [49]. It would be interesting to generalize our model to such nearest-neighbors coupling scenario and explore the effect of extra sparse long-range connection in such models. This type of models may share similar principles revealed in this study, as in our extended model where short-distance links are dominant (see Fig. S6). The physical principles revealed in our work can guide further development of braininspired efficient computing [50].

Analysis of the simplified model with correlated inputs
We measure the topological correlation in a RN by the ratio between the number of common neighbors and total number of distinct neighbors of a pair of neurons in the network, that is,  Fig. 4(c)) with the mean = 0 and the variance 2 ( 0 ( )) = 2 ( ( )).
For each synapse, as = 0 or 1 randomly, we have 〈 〉 = 〈 2 〉 = and the variance is 2 ( ) = − ( ) 2 . Thus, 2 ( 0 ) = [ − ( ) 2 ] . The output firing rate is determined by the probability that the input signal is above the threshold . Using the error function the firing rate [41] is obtained as as shown by the red dashed line ( = 0) in the Fig. 5(c).
As the correlation of spikes is present in the input synapses ( > 0) when a common spike train is added into all spike trains, the correlated spikes at the time t give the input signal ( ) = + ∑ [ ( ) − ]and its mean strength is , which can active the neuron at time t with probability 1 because ≫ . As the rate of the correlated spikes is = , the firing rate of the neuron is since the correlated spikes will always induce spiking of the neuron in this simplified model. Here in Eq. (7), the is given by Eq. (6) with replaced by = (1 − ) .

Mean-field theory of IF neural dynamics
In this section, we present the outline of the mean-field theory for deriving the field equations Eq.
Taking 〈 〉 ∈ , or 〈 〉 ∈ , to the second and third equations of Eq. (8), we get the second equation of Eq. (5), which finishes the construction of the coupled field equations.
In the limit of → 1 (all rewired), modules are almost separated. Let = 1 in Eq. (5) and we get the field equations corresponding to one separate module with additional external excitatory inputs, i.e.
Eq. (2), with = being the connection density of the module.

Supplementary Notes I. Network Setting
The whole system consists of neurons, 80% are excitatory neurons and 20% are inhibitory neurons. To model a local cortical surface, neurons are placed on square regions which are separated by blank space on a two-dimensional (2D) plane, as illustrated in Fig. 1 top panel in the main text. To measure the distance between different neurons, we first assume the length of each square and the width of the blank space are set as 1 ( Fig. 1 top panel). Each neuron , with the coordinate of its position ( , ), is randomly distributed in each square. In all cases, we assume each module (square) consists of = 500 neurons and the ratio between excitatory and inhibitory neurons in each module is kept as 4:1. The example in Fig.  1 top panel contains = 10 × 10 modules and in this case the whole network size is = × = 50,000. Throughout the work, the default setting is = 100, = 500 so that = 50,000 and = 0.0017 unless extra specifying. The random network (RN) is built by connecting each possible neuron pair with a probability . To build a modular network (MN), links between modules are rewired, with a probability , into square to become intra-modular links. For example, for an inter-module link → whose source neuron and target neuron are in different square modules, we replace the target node with a randomly selected neuron in the same module of (initially a link between and is absent). This rewiring makes the connection probability in a module higher than that between different modules, forming modular structure in the network while maintaining connection density of the whole network unchanged, that is, the total number of links on average is = ( − 1) ≈ 2 , which is fixed during the rewiring process.
In the initial random network, the number of intra-module links, denoted by , is = ( − 1) ≈ .
(S1.3) In the MN with the rewiring probability , the number of intra-modular links becomes (S1.4) Thus the intra-module density is The number of inter-modular links in MN becomes (S1.6) Thus the inter-module density is The length of a link from neuron to neuron is the Euclidean distance between the pair of neurons in the 2D plane = √ ( − ) 2 + ( − ) 2 . The normalized wiring cost of the whole network is defined as the summation of the length of all links rescaled by that value in initial RN without rewiring. Thus, the normalized wiring cost of initial RN is 1. Fig. S1 shows the relation between the normalized wiring cost and the rewiring probability . As the rewiring probability increases, the wiring cost decreases. When the value of is larger than 0.99, the normalized wiring cost tends to a constant. The results also hold for different module numbers (Fig. S1). This normalized wiring cost in MN is independent of the overall connection density , as shown in Fig. S4 bottom panels where the costs with connection density = 0.001 , 0.0015 , 0.002 are presented. This property can be understood by an analytic treatment as follows. We assume that the mean length of the intra-modular links is , whereas the mean length of the inter-modular links is . The normalized wiring cost is . (S1.8) Thus, it is independent of the connection density of the network, the same as the results shown in Fig. S4 bottom panels. As the rewiring probability tends to 1.0, the normalized wiring cost is (S1.9) Furthermore, the mean intra-modular link length does not change with the size of the 2D plane, while the mean inter-modular link length increases with the size of plane. Therefore, the normalized wiring cost decreases with the size of the system, as shown in Fig. S1. are independent of each other for different , . Here, we use noise strength = 10. The properties of the self-sustained dynamics is independent of this noise strength.
In the model of a separate module simulation, there is an external input spike train { ( ), ≥ 1} with rate = 50 . We use 100 realizations of the network in computing the probability of sustained activity . If the network still spikes at 1 second after the noise or external input has been removed, it is regarded to exhibit self-sustained activity.

Measure the response sensitivity in spiking networks
In experimental studies on the response of the brain to the cognitive events, it is necessary to average the measured brain activity over trails of experiments. Here, the stimulus method is that proportion of the randomly selected neurons in all modules are activated by increasing their membrane potential to −40 (above the threshold) in one simulation step. is termed the stimulus strength in our study.
We measure the response sensitivity in terms of the membrane potential properties and the spiking properties, using the average membrane potential and firing rate of the network modules. To study the sensitivity, we can consider the returning process of a signal to its baseline value after a transient perturbation.
The collective behavior of a module responding to additional stimulus can be reflected by its mean membrane potential and its mean firing rate: An example of the averaged signal is shown in Fig. S2. In general, the ongoing spontaneous fluctuation is almost eliminated by averaging over realizations, but the response behavior manifests itself in the averaged signal. If the response of the network is weak, the network returns to its baseline ongoing activity quickly, whereas strong response causes a waveform of damped oscillation as illustrated in Fig. S2 (see also Fig. 3(c) in the maintext for the cases of a module in RN and MN).

Measure the response sensitivity in the field models
We can apply similar analysis in the field equations simulation. In the single-module field equation Eq. (S2.10), the stimulus is to raise the E, I membrane potential ( and ) above the threshold to −40 . Then, we measure the corresponding area index calculated from the average membrane potential of the module = 0.8 + 0.2 . In the coupled version of the field equations Eq. (S2.6), we exert the above stimulus in a selected module and detect the response by the signal of this module as described above. Fig. S2. The illustration of computing the response size. The colored area between the signal curve (the trial-average voltage or firing rate) and its pre-stimulus baseline value (dashed line) in 0~250ms after stimuli is used to measure the response size. Stimuli are applied at the time marked by the arrows. Refer to Fig. 3(c) for the cases of RN and MN in the network model. Eq. (S2.3) essentially captures the sub and supra threshold microscopic dynamics of a spiking network, that is, ( ) represents the proportion of type neurons that spike between and + ∆ (∆ is an infinite small quantity) as well as the mean firing rate of type neurons at time with unit per ms [2]. Here, are effective parameters to construct the voltage-dependent mean population firing rate. Note that this approximation scheme based only on the first-order statistics neglects several factors that affect the accurate firing rate, including higher order statistics, noise correlation and refractory time. Thus, does not have an analytical form and should be estimated numerically. A complete analytical approach for conductance-based integrate-and-fire neural network is still an open issue [3,4]. Thus, the steady-state firing rate

Mean-field reduction of neuronal network to coupled neural oscillators
are independent of , which cannot capture the effect of the firing rate reduction during rewiring (Fig. 2 in the main text).

Field equations of a single module
To understand the overall dynamic principle of the modular network, we can focus on analyzing the dynamics of single separate module. In the limit of → 1 (all rewired), modules are almost separated. Let = 1 in Eq. (S2.6) and we get the field equations corresponding to one separate module with additional external excitatory inputs:  Results in the main text are obtained by = 0.0017. Conclusions are the same when is around this value. However, MN (for → 1) may not be able to self-sustain when is larger, as shown in (c). This can be understood from the self-sustain property of a single isolated module ( Fig. 4(a)). However, since we do not employ synaptic scaling (i.e. let the synaptic strength decreases with neighbor numbers) in our model, the network cannot maintain E-I balance for too large (data not shown). Fig. S5. Additional measurements of avalanches. In the main text, we measure avalanches from neuronal spikes ( Fig. 3 and 4). Here, we present other measurements of avalanches using the average membrane potential associated with a threshold 0 . (a) Avalanches defined by the above-threshold-event of the average membrane potential 〈 〉 of a module, as used in [5]. The size of an avalanche is defined as the area below the curve 〈 〉 and above the threshold in the event.  In (a, b), a threshold 0 = −60 is used. The threshold is approximately the mean (−66.7 ) plus 2 times of the standard deviation (2.9 ) of 〈 〉, as in previous studies. The green lines on the top indicate the ranges of estimated power-law distributions. The results suggest that our model can also exhibit features of critical avalanches in mesoscopic scale. Their properties deserve further study, for example, in terms of critical exponents and in association to the underlying spiking avalanches.