Symmetry breaking and avalanche shapes in modular neural networks

Modularity is as a key characteristic of structural and functional brain networks across species and spatial scales. We investigate the stochastic Wilson–Cowan model on a modular network in which synaptic strengths differ between intra-module and inter-module connections. The system exhibits a rich phase diagram comprising symmetric (with low and high activity) and “broken symmetry” phases. Symmetric phases are characterized by the same low or high activity in all the modules, while the broken symmetry phases are characterized by a high activity in a subset of the modules and low activity in the remaining ones. There are two lines of critical points, the first between the low activity symmetric phase and the high activity symmetric phase, and the second between the low activity symmetric phase and a broken symmetry phase with one active module. At those lines the system shows a critical behavior, with power law distributions in the avalanches. Avalanche shapes differ systematically along the two lines: they are symmetric or right-skewed at the transition with the symmetric phase, but become left-skewed over intermediate durations along critical line with the broken symmetry phase. These results provide a theoretical framework that accounts for both symmetric and left-skewed neural avalanche shapes observed experimentally, linking modular organization to critical brain dynamics.

1 Introduction

Spontaneous brain activity unfolds over multiple spatial and temporal scales, giving rise to rich dynamical patterns that include intermittent bursts known as neuronal avalanches (Beggs and Plenz, 2003; Mazzoni et al., 2007; Pasquale et al., 2008; Shew et al., 2009; Chialvo, 2010; Massobrio et al., 2015). First identified in rat cortical slices, and later confirmed in vivo across rodents (Gireesh and Plenz, 2008; Xu et al., 2024), non-human primates (Petermann et al., 2009; Yu, 2011, 2017; Miller et al., 2021) and humans (Palva et al., 2013; Shriki et al., 2013; Meisel et al., 2013; Solovey et al., 2012; Arviv et al., 2016; Tagliazucchi et al., 2012; Scarpetta et al., 2023), these cascades display size and lifetime distributions devoid of a characteristic scale, a signature long interpreted as evidence that cortical networks reside near a nonequilibrium critical point (Hengen and Shew, 2025). Scale-free avalanches were reproduced also in many computational models of neural activity (de Arcangelis et al., 2006; Levina et al., 2007; Millman et al., 2010; Poil et al., 2012; Scarpetta and de Candia, 2013; Scarpetta et al., 2018; Dalla Porta and Copelli, 2019), among which in a stochastic version of the Wilson–Cowan model (Benayoun et al., 2010; de Candia et al., 2021; Apicella et al., 2022; Mariani et al., 2022; Piuvezam et al., 2023; Alvankar Golpayegan and de Candia, 2023).

These scaling laws place the brain within the broader realm of crackling noise, in which slowly driven systems respond through impulsive events spanning many orders of magnitude (Sethna et al., 2001; Kuntz and Sethna, 2000), ranging from earthquakes (Bak and Tang, 1989; Olami et al., 1992) to Barkhausen noise in ferromagnets (Perković et al., 1995; Mehta et al., 2002; Durin and Zapperi, 2000). Historically, the power-law distributions of avalanche sizes and durations have been considered hallmarks of criticality, however this could not be a sufficient indicator of criticality, as they can arise from various non-critical mechanisms (Newman, 2005; Stumpf and Porter, 2012).

To overcome the ambiguity of power-law exponents alone, recent research has increasingly focused on the mean temporal profiles (shapes) of neuronal avalanches as a more stringent and reliable test for criticality (Papanikolaou et al., 2011; Laurson et al., 2013; Gleeson and Durrett, 2017). Critical systems predict that when appropriately rescaled, the mean temporal profiles of avalanches of widely varying durations should collapse onto a single universal scaling function, often approximated by an inverted parabola (Sethna et al., 2001).

Experimental observations, however, have yielded mixed results regarding these shapes. A parabolic profile for avalanches in the scaling regime has been observed with two-photon imaging of neurons in the superficial cortex of awake mice (Capek et al., 2023; Srinivasan et al., 2024). Nearly symmetric avalanche shapes have been observed in intracranial depth recordings in humans (Priesemann et al., 2013), as well as in the spiking activity of the mammalian primary visual cortex, both in anesthetized animals (rats and monkeys) and freely moving mice, where leftward asymmetric avalanches tend to occur at lower values of the spiking variability (Fontenele et al., 2019). In nonhuman primates, chronically implanted high-density microelectrode arrays showed a nearly parabolic avalanche shape, modulated by γ–scillations at intermediate raster plot binning (Miller et al., 2019, 2021).

Many other experimental measurements have revealed clear departures from perfect symmetry, displaying leftward skewing and extended tails. Whole-brain light-sheet imaging of transgenic zebrafish larvae has shown that average avalanche profiles collapse into a distinctly leftward-asymmetric form, characterized by rapid initiation and slow decay (Ponce-Alvarez et al., 2018). Similar leftward asymmetry has been observed in high-resolution in vitro recordings of dissociated cortical cultures (Friedman et al., 2012), whereas organotypic cortical cultures exhibit collapse onto approximately inverted parabolic profiles with only minimal leftward skewing. This characteristic leftward asymmetry has also been reported in human EEG recordings, both during recovery from hypoxia (Roberts et al., 2014) and in preterm infants (Iyer et al., 2015). Remarkably, EEG data from preterm infants revealed systematic changes in avalanche shape asymmetry as a function of avalanche duration, with avalanches transitioning from symmetric to leftward-asymmetric forms. Moreover, the degree of asymmetry at intermediate durations correlated significantly with long-term neurodevelopmental outcomes (Iyer et al., 2015). Similar asymmetric profiles are also observed in physical crackling noise systems, such as Barkhausen noise (Dobrinevski et al., 2013) and earthquakes (Houston et al., 1998; Houston, 2001), where they have been linked to complex underlying mechanisms like a negative “effective mass” or transient forces counteracting propagation (Zapperi et al., 2005). The origins of this leftward asymmetry in biological data remain puzzling. A notable contrast arises from the simple (non-modular) stochastic Wilson–Cowan model, where symmetric or even right-skewed avalanche profiles have been found (de Candia et al., 2021; Apicella et al., 2022).

Beyond the intrinsic dynamics of neuronal populations, the structural organization of brain networks plays a fundamental role in shaping neuronal activity. The brain is widely recognized as a complex network characterized by a modular and often hierarchically modular architecture (Meunier et al., 2009; Sporns et al., 2016; Gallos et al., 2012; Bassett and Bullmore, 2009). Modules are typically defined as subsets of nodes, often corresponding to anatomically or functionally related neuronal regions, that are densely and strongly interconnected internally, while maintaining comparatively sparser and weaker connections with nodes in other modules (Bullmore and Sporns, 2009). The role of modular network topology in spontaneous collective dynamics and the emergence of critical behavior has been investigated in previous studies (Moretti and Muñoz, 2013; Russo et al., 2014; Ódor et al., 2015; Cota et al., 2018; Rusch et al., 2025; Zhigalov et al., 2017; Yamamoto et al., 2023; Angiolelli et al., 2025; Myrov et al., 2024; Cirunay et al., 2025), but the impact of modularity on avalanche shape and on the manner in which avalanches propagate across modular regions remains an important open question. Motivated by these unresolved questions and the persistent discrepancies between current models and experimental findings regarding avalanche shapes, our study aims to bridge this gap by investigating neuronal avalanche distributions and, critically, their mean temporal shapes within a stochastic Wilson–Cowan model incorporating a modular network topology.

We first show that, on a modular network, the Wilson–Cowan model shows different phases, depending on the strength of excitatory and inhibitory connections between modules and within each module. Phases can be symmetric (SL and SH) characterized by the same (high and low) mean activity on all the modules, or they can be characterized by a breaking of the symmetry between modules (Bm, with m = 1, 2,…) corresponding to m modules having a high activity, and the others a low activity. Between SL and SH, and between SL and B1, there are two lines of critical points, where the system displays a scale-free distribution of the avalanches.

2 Wilson–Cowan model on a modular network

We consider the stochastic Wilson–Cowan model (Benayoun et al., 2010). In this model, neurons have two possible states, namely active and quiescent, and make transitions from one state to the other with specified rates. The transition from active to quiescent takes place with a constant rate α, while that from quiescent to active with a rate f(s_i), where s_i is the input of the i-th neuron that we are considering, and f(s) is an activation function taken equal to

$\begin{array}{ll}f(s)=\{\begin{array}{ll}\beta \tanh (s)\hfill & ifs>0,\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}& (1)\end{array}$

The input of the i-th neuron is given by

$\begin{array}{ll}{s}_i={\displaystyle \sum _{j=1}^N}{w}_{ij}{a}_j+h_i,& (2)\end{array}$

where a_j = 0, 1 if the j-th neuron is quiescent or active, w_ij is the strength of the connection going from j to i, and h_i is the external input to neuron i. In the simplest case, one takes N_E excitatory and N_I inhibitory neurons, with N_E+N_I = N, and the connections w_ij having only two values, namely w_E/N_E if the pre-synaptic neuron j is excitatory, or −w_I/N_I if j is inhibitory.

It was shown in de Candia et al. (2021) that, when the external input h_i goes to zero, the dynamics of the model undergoes a transition (bifurcation) at a critical value of the imbalance w₀ = w_E−w_I, namely $w_{0c}={\beta }^{-1}\alpha$. For w₀<w_0c, the fraction of active neurons goes to zero when h_i → 0, while for w₀>w_0c it is finite also for h_i → 0, so that there is a self-sustained activity. Exactly at w₀ = w_0c, fluctuations are maximized, and the dynamics is characterized by bursts of activity (avalanches), whose distribution is scale-free, with the critical exponents of the branching model. The transition can be understood as a “transcritical bifurcation,” where the quiescent fixed point of the dynamics becomes unstable, and the self-sustained fixed point becomes stable. As such, it can be found looking when one of the eigenvalues of the stability matrix becomes equal to zero.

In this paper, we consider a version of the model in which the network is composed by M modules, each containing N_E excitatory and N_I inhibitory neurons as schematically illustrated in Figure 1. Each neuron in module k receives input from excitatory neurons of module l with synaptic strengths $w_{kl}^E/N_E$, and from inhibitory neurons of module l with synaptic strengths $-w_{kl}^I/N_I$, so that the input s_k of neurons in module k is

$\begin{array}{ll}{s}_k={\displaystyle \sum _{l=1}^M}(w_{kl}^E{x}_l-w_{kl}^I{y}_l)+h_k,& (3)\end{array}$

where x_l and y_l are the fractions of excitatory and inhibitory active neurons in module l (activities), h_k is the external input to module k. In the Gaussian noise approximation (van Kampen, 2007; Benayoun et al., 2010), the dynamics of the network is determined by the equations

$\begin{array}{ll}\frac{d{x}_k}{dt}=-\alpha x_k+(1-x_k)f(s_k)+N_E^{-1/2}\sqrt{\alpha x_k+(1-x_k)f(s_k)}{\eta }_{E,k}(t),& (\mathrm{4a})\end{array}$

$\begin{array}{ll}\frac{d{y}_k}{dt}=-\alpha y_k+(1-y_k)f(s_k)+N_I^{-1/2}\sqrt{\alpha y_k+(1-y_k)f(s_k)}{\eta }_{I,k}(t),& (\mathrm{4b})\end{array}$

where η_{E, k}(t), η_{I, k}(t) are normal, Gaussian independent white noises. The variables x_k and y_k can be written as $x_k=X_k+N_E^{-1/2}{\xi }_k^E$, $y_k=Y_k+N_I^{-1/2}{\xi }_k^I$, where X_k and Y_k are the “deterministic” components of the activities, that obey the deterministic equations

$\begin{array}{ll}\frac{d{X}_k}{dt}=-\alpha X_k+(1-X_k)f(S_k),& (\mathrm{5a})\end{array}$

$\begin{array}{ll}\frac{d{Y}_k}{dt}=-\alpha Y_k+(1-Y_k)f(S_k),& (\mathrm{5b})\end{array}$

with $S_k=\sum _{l=1}^M(w_{kl}^E{X}_l-w_{kl}^I{Y}_l)+h_k$. The deterministic components evolve until they reach a stable (attractive) fixed point, that is a configuration of X_k and Y_k that set to zero the right hand sides of Equation 5. We can then define ${\Sigma }_k=\frac{X_k+Y_k}{2}$, ${\Delta }_k=\frac{X_k-Y_k}{2}$, where Σ_k represents the total activity of the k-th module, and Δ_k the “imbalance” between the activity of excitatory and inhibitory neurons. In terms of these variables, fixed point equations become

$\begin{array}{ll}-\alpha {\Sigma }_k+(1-{\Sigma }_k)f(S_k)=0,& (\mathrm{6a})\end{array}$

$\begin{array}{ll}-\left[\alpha +f(S_k)\right]{\Delta }_k=0.& (\mathrm{6b})\end{array}$

At a fixed point it is therefore Δ_k = 0, which means that the activity of excitatory and inhibitory neurons inside each module is equal. This is a consequence of the structure of the connections, namely the fact that excitatory and inhibitory neurons of a given module receive the same input. The “stochastic” components ${\xi }_k^E$ and ${\xi }_k^I$ obey Langevin equations that can be directly obtained subtracting Equation 5 from Equation 4. Expanding in powers of N^−1/2, at leading order one finds 2M coupled linear Langevin equations, namely

$\begin{array}{l}\frac{d}{dt}\left(\begin{array}{c}{\xi }_k^E\\ {\xi }_k^I\end{array}\right)=\\ {\displaystyle \sum _{l=1}^M}\left(\begin{array}{cc}-(\alpha +f(S_k^{*})){\delta }_{kl}+f^{\prime }(S_k^{*})w_{kl}^E& -f^{\prime }(S_k^{*})w_{kl}^I\\ f^{\prime }(S_k^{*})w_{kl}^E& -(\alpha +f(S_k^{*})){\delta }_{kl}-f^{\prime }(S_k^{*})w_{kl}^I\end{array}\right)\end{array}$

$\begin{array}{ll}\left(\begin{array}{c}{\xi }_l^E\\ {\xi }_l^I\end{array}\right)+\sqrt{2\alpha {\Sigma }_k^{*}}\left(\begin{array}{c}{\eta }_{E,k}(t)\\ {\eta }_{I,k}(t)\end{array}\right),& (7)\end{array}$

where δ_kl is the Kroneker delta, δ_kl = 1 if k = l, δ_kl = 0 otherwise, and ${\Sigma }_k^{*}$, $S_k^{*}$ are evaluated at the fixed point. We then make the assumption that

$\begin{array}{l}{w}_{kl}^E=w_0^E{\delta }_{kl}+w_1^E(1-{\delta }_{kl}),\\ w_{kl}^I=w_0^I{\delta }_{kl}+w_1^I(1-{\delta }_{kl}).\end{array}$

In this case, Equation 6a will be given by F_k({Σ_l}) = 0, where

$\begin{array}{ll}{F}_k(\left\{{\Sigma }_l\right\})=-\alpha {\Sigma }_k+(1-{\Sigma }_k)f(w_0{\Sigma }_k+w_1{\displaystyle \sum _{l\ne k}}{\Sigma }_l+h_k),& (8)\end{array}$

where $w_0=w_0^E-w_0^I$, $w_1=w_1^E-w_1^I$. In the following, we consider the external input h_k = 0 for all the modules, and α = β. Equation 8 always has the solution with all the activities Σ_k = 0. This is an absorbing quiescent fixed point, in the sense that the dynamics ceases completely in that state, and has to be restarted by activating one or more neurons by an external input. Depending on the values of w₀ and w₁, this fixed point can be attractive or repulsive. When it is attractive, we call the phase “SL,” which stands for “symmetric low” activity. In this phase, for a very low external input, the activity of the network is correspondingly low and incoherent, with sparse and uncorrelated spikes. For values of w₀ and w₁ for which the quiescent fixed point is repulsive, different attractive fixed points appear, in which the activity is different from zero also in absence of external input, that is the dynamics is self-sustained. We look for fixed points where modules k = 1, …, m have an activity Σ_k = Σ > 0, while module k = m+1, …, M have Σ_k = 0, so that the fixed point equations become, for h_k = 0

$\begin{array}{ll}-\alpha \Sigma +(1-\Sigma )f\left[(w_0+(m-1)w_1)\Sigma \right]=0,& (\mathrm{9a})\end{array}$

$\begin{array}{ll}f\left(m{w}_1\Sigma \right)=0,& (\mathrm{9b})\end{array}$

respectively for k ≤ m and k > m. By Equation 9a, as the function f(s) is increasing and f(s) ≤ βs for s > 0, there is a solution with Σ > 0 only for $w_0+(m-1)w_1>{\beta }^{-1}\alpha$. On the other hand, if m<M, by Equation 9b it has to be w₁ ≤ 0. Therefore, for w₁ > 0 and $w_0>{\beta }^{-1}\alpha -(M-1)w_1$, there is only a solution with all the Σ_k greater than zero. We call this phase “SH,” which stands for “symmetric high” activity. For w₁ ≤ 0 and $w_0\ge {\beta }^{-1}\alpha -(m-1)w_1$, there will be also the solution with m modules having Σ_k > 0 and M−m modules having Σ_k = 0. We call this solution “Bm,” where the B stands for breaking of symmetry. Of course there are $\genfrac{}{}{0.0pt}{}{M}{m}$ such fixed points, corresponding to which modules are chosen to be active.

Figure 1

Figure 1. Connectivity of the network. Only the first two modules, and the connections outgoing from neurons of the first module, are shown. Synaptic strengths are $w_0^E$ and $w_0^I$ for excitatory and inhibitory connections within the same module (intra-module), and $w_1^E$ and $w_1^I$ for excitatory and inhibitory connections between different modules (intermodule), defining the effective imbalances $w_0=w_0^E-w_0^I$ and $w_1=w_1^E-w_1^I$. The connectivity matrix is symmetric under permutation of the modules, so that the connections outgoing from any module to any other module are equal.

Note that the order parameter of the transition, namely the mean activity Σ of the neurons, which is zero in the phase SL, is continuous along the two critical lines separating SL from SH and B1, and therefore goes to zero while approaching the lines from above, with an exponent β = 1 (not to be confused with the rate appearing in Equation 1). The exponent can be derived by making a Taylor expansion of f(s) around s = 0 in Equation 9a, and taking only the first linear term. The fluctuations of the order parameter can be also computed (de Candia et al., 2021), and it is found that they diverge with the exponent γ = 1. Interestingly these exponents are the same as the ones of the branching model, or mean-field directed percolation.

Note also that the present model does not support Hopf-bifurcations and oscillations, due to the structure of the connections. Indeed, the connections depend only on the pre-synaptic neuron, so that all the neurons of a given module receive the same input, and the connections are invariant for a permutation of the modules. As a consequence, the eigenvalues of the Jacobian matrix appearing in Equation 7 are always real, and no oscillations are possible.

To look for the stability of the fixed points, one has to look at the sign of the eigenvalues of the Jacobian matrix. If all the eigenvalues are negative, then the fixed point is attractive. If on the other hand one or more eigenvalues become positive, then the fixed point becomes repulsive. In Figure 2 we plot the phase diagram in the case M = 3, showing which fixed point is stable depending on the values of w₀ and w₁. The two solid lines, respectively $w_0={\beta }^{-1}\alpha$ for w₁ ≤ 0, and $w_0={\beta }^{-1}\alpha -(M-1)w_1$ for w₁≥0, represent the limit of stability of the low activity fixed point, that is of the phase SL. Exactly on those lines, one of the eigenvalues of the Jacobian matrix becomes equal to zero, so that the stability of the low activity phase is only marginal. If one or more neurons are activated by an external input, then an “avalanche” is generated, that is a very long sequence of activations and deactivations, before the system comes back to the quiescent state. These avalanches have a power law distributions of the durations and sizes, as it will be shown in the next section.

Figure 2

Figure 2. Phase diagram of the Wilson–Cowan model on a modular network, in the case M = 3 and for α = β. In each region of parameter space, the phases that are stable (that is whose fixed point is attractive) are listed. In the phase SL the activities of all the modules are low (proportional to h) for h → 0. In the phase SH the activities of all the modules are high (finite in the limit h → 0). In phase B1 (B2) one (two) modules have high activity and the others low activity. The values of the activity are continuous at the boundary of phase SL, as the activity in phases SH and B1 tends to zero when approaching the separation lines, that is transitions are “second order” in the language of statistical physics. The fixed point of the dynamics is “marginally stable” there, so that avalanches have a scale-free distribution. Blue dots (A–C) are the points where avalanche distributions and their average shape are investigated in the following section. On the other hand, at dashed lines one or more fixed points lose stability in a discontinuous manner.

3 Critical dynamics at the edge of the low activity phase

In this section, we study the distribution of bursts of activity (avalanches), and their average shape, varying the strength of the excitatory and inhibitory connections, both between different modules (inter-module connections, labeled by index 1) and within each module (intra-module connections, labeled by index 0). While the fixed points of the dynamics, defined by Equation 6, only depend on the “imbalances” $w_0=w_0^E-w_0^I$, $w_1=w_1^E-w_1^I$, that is the differences between excitatory and inhibitory connections, we will see that the probability distributions of the avalanches and their average shape depend separately on all the parameters $w_0^E$, $w_0^I$, $w_1^E$ and $w_1^I$. We therefore, for given values of w₀ and w₁, also change the values of the sums $w_0^S=w_0^E+w_0^I$ and $w_1^S=w_1^E+w_1^I$. As we have defined $w_0^E$, $w_0^I$, $w_1^E$ and $w_1^I$ to be non-negative, it has to be $w_0^S\ge \left|w_0\right|$ and $w_1^S\ge \left|w_1\right|$.

In the following we consider M = 3, α = β = 0.1 ms⁻¹, $N_E=N_I=10^5$. We moreover use h = 0. This means that the state where all the neurons are quiescent is an absorbing one, that is the dynamics of the network cease completely. Therefore, each time the network reaches such an absorbing state, we start the dynamics again by setting one neuron of one of the modules (chosen randomly) active “by hand.” The dynamics therefore corresponds to the limit h → 0, where after the avalanche is started, neurons can become active only due to the input of the network itself, since there is no external input.

After setting the first neuron active, we simulate the system by means of the event-driven Gillespie algorithm (Gillespie, 1977). The steps of the algorithm are the following: (1) for each neuron i compute the transition rate r_i: r_i = α if neuron i is active, or r_i = f(s_i) if it is quiescent; (2) compute the sum over all neurons $r=\sum _i{r}_i$; (3) draw a time interval dt from an exponential distribution with rate r; (4) choose the i−th neuron with probability r_i/r and change its state; and (5) update the time to t+dt. The avalanche ends when all the neurons become quiescent. The size of the avalanche is defined as the total number of spikes, that is transitions of a neuron from the quiescent to the active state. The duration is the difference in time between the first spike (the first neuron that is activated “by hand”), and the last deactivation of a neuron.

3.1 Disconnected modules

In this subsection, we start by considering the case w₁ = 0, $w_1^S=0$, that corresponds to $w_1^E=w_1^I=0$. This is therefore the case of disconnected modules, that coincides with the case of a single module M = 1, that is with a homogeneous (non-modular) network, already studied in (de Candia et al. 2021). In this case we have just two phases (two attractive fixed points), corresponding to a phase of low incoherent activity for w₀<w_0c, and a phase of high self-sustained activity for w₀ > w_0c, where $w_0=w_0^E-w_0^I$ is the “imbalance” between excitatory and inhibitory connections. These phases correspond to the phases SL and SH of the modular network. The critical point is given by $w_{0c}={\beta }^{-1}\alpha$, that is w_0c = 1 as we consider here α = β. For w₀<w_0c, the attractive fixed point corresponds to a quiescent network, while for w₀ > w_0c the attractive fixed point corresponds to a finite value of the activity, so that the dynamics of the network is self-sustained. At w₀ = w_0c = 1 the quiescent state is marginally stable, which means that an avalanche started by setting just one neuron active will have a power law distribution of sizes and durations, with a cut-off only determined by the finite number of neurons. The point w₀ = 1, w₁ = 0 corresponds to the point “A” in Figure 2. Note that in the case of disconnected modules (w₁ = 0, $w_1^S=0$), and in general in the case w₁≥0 (balanced or “excitatory imbalanced” inter-module connections), there can be no “broken symmetry” phases Bm, that appear only for “inhibitory imbalanced” (w₁ < 0) inter-module connections.

In Figure 3 we show the avalanche statistics for different values of $w_0^S=w_0^E+w_0^I=1,2,5,10$. The case $w_0^S=1$ corresponds to a value $w_0^I=0$ of the inhibitory connections, that is to a purely excitatory network. In Figures 3A, B we show the probability distribution of avalanche sizes and durations, while in Figure 3C we show the average size of the avalanche as a function of its duration. Dashed lines correspond to the predictions of the branching model, or mean-field directed percolation (MFDP). It is interesting to note that the larger the value of $w_0^S$, the smaller the deviations with respect to the branching model predictions. The purely excitatory network, while converging to the predicted exponents for very large avalanches, has the most significant deviations for not very large ones. In Figure 3D we show the mean skewness of the avalanches as a function of their duration. The skewness is strictly zero for the purely excitatory network, and increases for durations around T = 100 ms, when increasing the strength of both inhibitory and excitatory connections (note that the difference $w_0=w_0^E-w_0^I$ has to be kept constant, in order to remain in the critical point). Finally in Figures 3E, F we show the average shape of the avalanches, that is the average number of spikes per millisecond as a function of the time since the start of the avalanche, for avalanches of duration T = 100 and T = 1, 000 ms. Notably, the single module only produces symmetric avalanches (in the purely excitatory case), or avalanches with a rightward asymmetry for an intermediate range of durations, with the skewness increasing with the increase of the strength of the connections (and thus of the relative strength of inhibitory connections).

Figure 3

Figure 3. Avalanche statistics in the case of disconnected modules, or of a single module, for $w_0=w_0^E-w_0^I=1$ (critical value), and different values of the parameter $w_0^S=w_0^E+w_0^I=1,2,5,10$. The case $w_0^S=1$ corresponds to a value $w_0^I=0$ of the inhibitory connections (purely excitatory network), and is characterized by a perfect symmetry of the avalanches, with skewness equal to zero. (A) Probability distribution of the sizes; (B) Probability distribution of the durations; (C) Mean size as a function of the duration; (D) Mean skewness of the avalanches as a function of the duration; (E) Average shape of the avalanches as a function of time, for avalanches of duration T = 100 ms; (F) Same as (E) but for durations T = 1, 000 ms.

In conclusion, in a range of intermediate durations 10 < T < 1, 000, for small values of $w_0^S$ we find nearly symmetric avalanches but strong deviations from MFDP exponents in the distribution of avalanches. On the other hand, for high values of $w_0^S$, we find very good agreement with MFDP exponents, but right-skewed avalanches. This is quite counterintuitive, as avalanches in MFDP have a symmetric (inverted parabolic) shape.

3.2 Modular network with balanced inter-module connections

In this subsection, we consider the case in which the system is at the point w₀ = 1, w₁ = 0 (as in the previous subsection), corresponding to the crossing of all the lines in Figure 2 (point “A”), on the boundary between phase SL and the different phases of high activity, but this time $w_1^S>0$. The modules that constitute the network are therefore connected by excitatory and inhibitory connections greater than zero but exactly balanced, that is with $w_1^E=w_1^I$. We fix the value $w_0^S=10$, and consider different values of $w_1^S$. As previously remarked, the fixed points of the dynamics, and their stability, depend only on the “imbalances” $w_0=w_0^E-w_0^I$ and $w_1=w_1^E-w_1^I$, so that once w₀ and w₁ are fixed to the critical point, changing the values of $w_0^S=w_0^E+w_0^I$ and $w_1^S=w_1^E+w_1^I$ cannot bring the system away from criticality. Moreover, it can be expected that, for a very large system and for very long avalanches, the exponents of power law distributions, and their average shape, do not depend on the values of $w_0^S$ and $w_1^S$. However, their value can change the probability distributions and the shape of not very large avalanches. This was already observed in Figure 3 for a disconnected network (or a single module) with $w_1^S=0$, and is confirmed in Figure 4 for a connected network formed by M = 3 modules and $w_1^S>0$. Deviations from MFDP behavior are largest for small values of $w_1^S$ (small values of both excitatory and inhibitory inter-module connections), and tend to reduce for larger values of $w_1^S$. Moreover, while very large avalanches tend to be symmetric, intermediate ones with durations 10 < T < 1, 000 are generally right-skewed (negative skewness), apart from the case $w_1^S=5$ that shows avalanches slightly left-skewed for durations around T = 1, 000.

Figure 4

Figure 4. Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w₀ = 1, w₁ = 0, so that inter-module excitatory and inhibitory connections are exactly balanced, $w_1^E=w_1^I$. The value of $w_0^S=w_0^E+w_0^I$ is fixed to $w_0^S=10$, while $w_1^S=w_1^E+w_1^I$ varies between 1 and 10. For reference, we include in the plots also the case $w_1^S=0$ (red lines), that corresponds to the yellow lines in Figure 3. (A) Probability distribution of the sizes; (B) Probability distribution of the durations; (C) Mean size as a function of the duration; (D) Mean skewness of the avalanches as a function of the duration; (E) Average shape of avalanches as a function of time, for avalanches of duration T = 100 ms; (F) Same as (E) but for durations T = 1, 000 ms.

3.3 Modular network with “excitatory imbalanced” inter-module connections

If we shift the balance of inter-module connections making them excitatory imbalanced, that is take w₁ > 0, we have to decrease correspondingly the imbalance of intra-module connections to remain at the critical point. Indeed as discussed after Equation 9 and shown in Figure 2, the critical line for w₁≥0 is given by the equation

$w_0+(M-1)w_1={\beta }^{-1}\alpha =1,$

having chosen α = β. In Figure 5 we show the case M = 3, w₁ = 0.1, w₀ = 0.8, corresponding to the point “B” in Figure 2. As can be seen by comparing Figures 4, 5, the behavior for w₁ > 0 is quite similar to that for w₁ = 0. There are deviations from the MFDP predictions at low values of $w_1^S$, that are reduced increasing $w_1^S$. Moreover, for large values of $w_1^S$, avalanches become right-skewed for intermediate values of the duration.

Figure 5

Figure 5. Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w₀ = 0.8, w₁ = 0.1 (point “B” in Figure 2), that is for “excitatory imbalanced” inter-module connections. The value of $w_0^S=w_0^E+w_0^I$ is fixed to $w_0^S=10$, while $w_1^S=w_1^E+w_1^I$ varies between 0.1 and 10. (A) Probability distribution of the sizes; (B) Probability distribution of the durations; (C) Mean size as a function of the duration; (D) Mean skewness of the avalanches as a function of the duration; (E) Average shape of avalanches as a function of time, for avalanches of duration T = 100 ms; (F) Same as (E) but for durations T = 1, 000 ms.

3.4 Modular network with “inhibitory imbalanced” inter-module connections

As shown in Figure 2, if inter-module connections are “inhibitory imbalanced,” that is $w_1^E<w_1^I$ and w₁ < 0, then “broken symmetry” phases appear, where a subset of the modules have high activity, while others have low activity. In particular, symmetric phase SL becomes unstable at the critical line w₀ = 1 toward phase B1, where only one of the modules is active. In Figure 6 we show avalanche statistics for w₀ = 1, w₁ = −0.15, that is at point “C” in Figure 2, for $w_0^S=10$ and different values of $w_1^S$ between 0.5 and 10.

Figure 6

Figure 6. Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w₀ = 1, w₁ = −0.15 (point “C” in Figure 2), that is for “inhibitory imbalanced” inter-module connections. The value of $w_0^S=w_0^E+w_0^I$ is fixed to $w_0^S=10$, while $w_1^S=w_1^E+w_1^I$ varies between 0.5 and 10. (A) Probability distribution of the sizes; (B) Probability distribution of the durations; (C) Mean size as a function of the duration; (D) Mean skewness of the avalanches as a function of the duration; (E) Average shape of avalanches as a function of time, for avalanches of duration T = 50 ms; (F) Same as (E) but for durations T = 400 ms; (G) Average shape of avalanches of different durations, for $w_1^S=7$; (H) Same as (G) but for $w_1^S=10$; (I) Mean skewness of the avalanches as a function of the duration, for different number of modules, $w_1^S=10$.

As shown in Figures 6A–C, in this case, differently from the case w₁ > 0, the largest deviations from MFDP are observed at large values of $w_1^S$, in particular between T = 100 ms and T = 1, 000 ms, where both the distributions of sizes P(S) and of durations P(T) decay with an exponent larger (steeper decay) than the one of MFDP, while the average size as a function of duration 〈S〉(T) increases with a smaller exponent. Notably, in the same range of avalanche durations, that is between T = 100 ms and T = 1, 000 ms, the average shape of the avalanches tends to become leftward skewed (positive skewness) as the strength of inter-module interactions increases, that is for larger values of $w_1^S$, as shown in Figure 6D. Figures 6E, F show the average shape for avalanches of duration 50 ms and 400 ms. For $w_1^S=10$ and durations around 400 ms, the skewness assumes a quite large value of almost 0.4. In Figures 6G, H we plot the average shape for different avalanche durations, respectively for $w_1^S=7$ and $w_1^S=8$. Finally in Figure 6I we plot the skewness as a function of the number of modules, showing that the effect is enhanced for greater number of modules. In the followinf subsection, we investigate the origin of such a leftward asymmetry, given that similar results have been reported in several experimental systems, as discussed in the Introduction. Specifically, we analyze more in detail the differences between avalanches for “excitatory imbalanced” and “inhibitory imbalanced” inter-module connections, that is at points B and C of Figure 2.

3.5 Comparison between “excitatory imbalanced” and “inhibitory imbalanced” inter-module connections

In Figure 7 we show different examples of avalanches of duration T = 400 ms, in the case of “excitatory imbalanced” (upper row, parameters w₀ = 0.8, w₁ = 0.1) and “inhibitory imbalanced” (lower row, parameters w₀ = 1, w₁ = −0.15) inter-module connections, highlighting the contributions of the individual modules to the avalanche. The most apparent feature of the plots is the fact that, in the first case where inter-module connections are “excitatory imbalanced,” there is a strong positive correlation between the profiles of the avalanche across all three modules. It is quite predictable then that the profile is not very different from the one observed in the case of a single module, or disconnected modules, that is with a rightward asymmetry. On the other hand, in the case of “inhibitory imbalanced” inter-module connections, the behavior is markedly different. We observe a first interval of time, of about 100 ms, where the avalanche spreads in a nearly uniform way in the three modules, while for subsequent times one of the modules prevails over the others, and remains almost the only active module until the end of the avalanche. To understand the origin of this behavior, we plot in Figure 8A the average imbalance 〈Δ_k〉 = 〈(x_k−y_k)/2〉 between the excitatory activity x_k (fraction of active excitatory neurons) and the inhibitory activity y_k, on the k-th module, and the average total activity 〈Σ_k〉 = 〈(x_k+y_k)/2〉 (Inset of Figure 8A). As 〈⋯ 〉 represents the average over all the avalanche of fixed duration, 〈Δ_k〉 and 〈Σ_k〉 are equal for all the modules k = 1, …M. Red (orange) line corresponds respectively to the excitatory (inhibitory) imbalanced case, that is to parameters w₀ = 0.8, w₁ = 0.1 (w₀ = 1, w₁ = −0.15). From these, one can compute the “intra-module” and “inter-module” input of neurons. The total input s_k of a neuron in module k, given by Equation 3, can indeed be written as s_k = s_{0, k}+s_{1, k}+h_k, where

$\begin{array}{ll}{s}_{0,k}=w_0{\Sigma }_k+w_0^S{\Delta }_k,\langle s_{0,k}\rangle =w_0\langle {\Sigma }_k\rangle +w_0^S\langle {\Delta }_k\rangle ,& (\mathrm{10a})\end{array}$

$\begin{array}{l}{s}_{1,k}={\displaystyle \sum _{j=1}^M}(1-{\delta }_{jk})(w_1{\Sigma }_j+w_1^S{\Delta }_j),\end{array}$

$\begin{array}{ll}\langle s_{1,k}\rangle =(M-1)(w_1\langle {\Sigma }_k\rangle +w_1^S\langle {\Delta }_k\rangle ).& (\mathrm{10b})\end{array}$

The input s_{0, k} comes from neurons belonging to the same module, while s_{1, k} from neurons belonging to different modules. In the case w₁ < 0 (“inhibitory imbalanced” inter-module connections) the input coming from module k to neurons in different modules can still be positive, if

$\begin{array}{ll}{\Delta }_k>\frac{\left|w_1\right|}{w_1^S}{\Sigma }_k.& (11)\end{array}$

The two contributions 〈s_{0, k}〉 and 〈s_{1, k}〉 are plotted in Figures 8B, C. It can be noticed that, in the case of “inhibitory imbalanced” inter-module connections (orange lines), the external input (Figure 8C) is positive in a first interval of around 50 ms, due to the prevalence of the activity of excitatory neurons in the beginning of the avalanche, so that Equation 5 is satisfied. Therefore, although inter-module inhibitory connections are stronger than inter-module excitatory ones, in the beginning of the avalanche one observes a “mutual reinforcement” between the modules, as the one observed in the case of excitatory imbalanced connections (albeit of course weaker). Only for times larger than 100–150 ms mutual inhibition, and therefore the symmetry breaking between modules, takes hold. This explains why leftward avalanches are observed only for durations between 100 and 1,000 ms. For shorter durations, symmetry breaking has not developed yet. On the other hand, for durations larger than 1,000 ms, the initial time window where all the modules are engaged in the avalanche can be neglected. For intermediate durations, the avalanche is formed by a first interval where all the modules are engaged, and a second one where only one module is left, giving rise therefore to a leftward asymmetry.

Figure 7

Figure 7. (A–C) Examples of avalanches of duration T = 400 ms in the case of “excitatory imbalanced” inter-module connections (parameters w₀ = 0.8, w₁ = 0.1). (D–F) The same for “inhibitory imbalanced” inter-module connections (parameters w₀ = 1, w₁ = −0.15). In both cases $w_0^S=w_1^S=10$. The instantaneous firing rates of the M = 3 modules are shown with red, green, and blue lines, and three different realizations of an avalanche are shown in the three columns. Black lines show the average firing rate (over all avalanches with duration 400 ms), on a single module in the first case (upper row), and the total average rate on the three modules in the second case (lower row).

Figure 8

Figure 8. (A) Average half-difference 〈Δ_k〉 between the excitatory and inhibitory activity (fraction of active neurons) on the k-th module during an avalanche of duration T = 400 ms, in the case of “excitatory imbalanced” (red line) and “inhibitory imbalanced” (orange line) inter-module connections. Inset: average half-sum 〈Σ_k〉 of the activities; (B) Average input to the neurons from the same module to which they belong; (C) Average input to the neurons from the modules to which they do not belong.

4 Conclusions

Modularity, likely shaped by evolutionary and functional constraints, represents a fundamental organizing principle of complex networks across social, technological, and biological systems (Girvan and Newman, 2002; Newman, 2006). This principle is observed across levels of biological organization, from genetic regulatory and protein interaction networks to the structure of ecosystems. A modular community architecture constitutes a core feature of both functional and structural brain networks (Chen et al., 2008; Sporns, 2010; Bullmore and Sporns, 2009; Bassett and Bullmore, 2009), observed at both whole-brain and cellular scales. Here, we have studied a stochastic Wilson–Cowan model incorporating a modular architecture, to investigate how such structural organization influences spontaneous dynamics, exploring the phase diagram and the avalanche shape near the analytically derived critical lines.

Experimentally, avalanche shapes can exhibit marked left-skewed asymmetry, that are not reproduced in stochastic Wilson–Cowan models lacking modular network organization, which, both under all-to-all coupling and within two-dimensional network structures, display either symmetric or rightward avalanche shapes in the vicinity of the critical point (de Candia et al., 2021; Apicella et al., 2022).

To clarify the role of network architecture, while capturing the modular organization observed in the brain, we investigated a stochastic modular Wilson–Cowan model composed of M modules, coupled by intra and inter-module connections, under a vanishing external drive. By analyzing the fixed-point equations and linearizing the dynamics around them, we derived a phase diagram comprising symmetric phases with low (SL) or high (SH) activity in all modules, and broken symmetry phases Bm, where only m<M modules are active. Along two critical lines, between phases SL and SH, and between SL and B1, the distribution of the avalanches is scale-free.

We then characterized the average shape of the avalanches along these lines. At the critical line between phases SL and SH, avalanches show a rightward shape, as in the case of a non-modular architecture. On the other hand, at the critical line between SL and B1, that is at the edge of the broken symmetry phase, avalanches show a leftward asymmetry for an intermediate range of durations. We have analyzed in detail the dynamics of the network during an avalanche in this latter condition, finding that it proceeds in two stages: an initially cooperative regime, where excitatory activity is prevalent, followed by inhibitory competition that selects one dominant module and suppresses the others. This is the relevant mechanism leading to a fast rise of the avalanche, followed by a slower decay, and therefore to leftward asymmetry.

It has been shown that the measured values of the exponents and avalanche shapes in a system at criticality may depend on experimental conditions, like for example subsampling (Girardi-Schappo, 2021; Carvalho et al., 2021; Conte and de Candia, 2025). We stress that the reported asymmetry is an intrinsic dynamical outcome of the fully simulated network, and is not due to measurement artifacts like subsampling. Therefore, while subsampling may in principle distort measured skewness in experiments, it is not required to generate the skewness mechanism reported here.

Our results are consistent with the observation that the location of criticality in E–I networks is not uniquely tied to a canonical “balanced” condition (Girardi-Schappo et al., 2020). In our model, criticality is defined by the marginal stability of the low-activity (absorbing) fixed point, and this does not correspond to a strict balance between excitation and inhibition. However, the skewed avalanche shapes we report are not a generic consequence of being “unbalanced.” They arise at a standard absorbing-state critical boundary, specifically the SL–B1 line, and reflect an additional dynamical mechanism present there: symmetry-breaking inter-module competition following an initial cooperative recruitment stage.

A symmetric avalanche profile typically arises in branching-process models, i.e. mean-field directed percolation, and in some neural models lacking modular organization (Poil et al., 2012; Dalla Porta and Copelli, 2019). In the non-modular fiber-bundle model avalanche profiles are right-skewed, and become symmetric in the mean-field limit of high dimensionality (Danku et al., 2018). Symmetric or right-skewed profiles are observed also in non-modular stochastic Wilson–Cowan models (see Figure 3). On the other hand, we have shown here that introducing a modular structure, with inhibitory-biased inter-module connections, can give rise to a left-skewed avalanche shape. Notably, a left-skewed profile has also been reported in a variable-threshold model implemented on a modular architecture derived from the KK-18 human connectome (Ódor, 2016), further suggesting that network topology may contribute to left-skewed avalanche dynamics.

The role of inhibition in producing leftward-skewed shapes has also been recently found (Zaccariello et al., 2025) in an integrate-and-fire model with short- and long-term plasticity and more highly connected inhibitory neurons, which highlighted the role of inhibition, and of differences in synaptic recovery rates between excitatory and inhibitory neurons, in generating leftward asymmetry.

In our framework, the key ingredient behind symmetry breaking and avalanche shape distortion is the emergence of an effective competition between modules (inhibitory-imbalanced inter-module coupling), preceded by an initial cooperative recruitment stage. An interesting extension of our study would be the introduction of a hierarchically modular topology, with “small scale” modules within “large scale” ones. Introducing such a hierarchy would not remove the main mechanism behind avalanche skewness; rather, it would naturally introduce multiple coupling scales, so that we expect the symmetry-breaking scenario to become multi-stage, with sequential recruitment/competition across hierarchical levels. This would likely (i) shift the location of the transition (since it is controlled by the relevant eigenmodes of the block-structured coupling matrix); and (ii) further enhance and/or structure the avalanche-shape distortion by adding additional time scales (e.g., a longer cooperative phase followed by stronger late-stage selection).

In conclusion, this model demonstrates that the leftward shape of avalanches emerges as a direct consequence of the symmetry-breaking transition driven by competitive inter-module connections, where inhibition dominates over excitation between modules. Notably, this framework enables systematic manipulation of inter-module interactions, biasing them toward either competitive (inhibitory) or cooperative (excitatory) dynamics, while maintaining critical behavior by adjusting the E/I balance of intra-module connections. Together, these findings highlight how the interplay between intra- and inter-module connectivity shapes both the nature of the symmetry-breaking transition and the emergence of asymmetric avalanches, providing a mechanistic link between network architecture and spontaneous neural dynamics.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

AC: Conceptualization, Investigation, Writing – original draft, Writing – review & editing. DC: Conceptualization, Data curation, Investigation, Writing – original draft, Writing – review & editing. HG: Conceptualization, Data curation, Investigation, Writing – original draft, Writing – review & editing. SS: Conceptualization, Investigation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. SS acknowledge financial support under the National Recovery and Resilience Plan, Mission 4, Component 2, Investment 1.1, Call for tender No. 1409 published on 14.9.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU – PRIN 2022 PNRR P2022BS4EH – CUP D53D23016350001.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alvankar Golpayegan, H., and de Candia, A. (2023). Bistability and criticality in the stochastic Wilson–Cowan model. Phys. Rev. E 107:034404. doi: 10.1103/PhysRevE.107.034404

PubMed Abstract | Crossref Full Text | Google Scholar

Angiolelli, M., Stefano, M., Filippi, S., Chiodo, L., Scarpetta, S., Palmieri, V., et al. (2025). Role of criticality in the structure-function relationship in the human brain. Phys. Rev. Res. 7:043153. doi: 10.1103/PhysRevResearch.7.043153

Crossref Full Text | Google Scholar

Apicella, I., Scarpetta, S., de Arcangelis, L., Sarracino, A., and de Candia, A. (2022). Power spectrum and critical exponents in the 2D stochastic wilson cowan model. Sci. Rep. 12:21870. doi: 10.1038/s41598-022-26392-8

PubMed Abstract | Crossref Full Text | Google Scholar

Arviv, O., Medvedovsky, M., Sheintuch, L., Goldstein, A., and Shriki, O. (2016). Deviations from critical dynamics in interictal epileptiform activity. J. Neurosci. 36, 12276–12292. doi: 10.1523/JNEUROSCI.0809-16.2016

PubMed Abstract | Crossref Full Text | Google Scholar

Bak, P., and Tang, C. (1989). Earthquakes as a self-organized critical phenomenon. J. Geophys. Res. Solid Earth 94, 15635–15637. doi: 10.1029/JB094iB11p15635

Crossref Full Text | Google Scholar

Bassett, D. S., and Bullmore, E. T. (2009). Human brain networks in health and disease. Curr. Opin. Neurol. 22, 340–347. doi: 10.1097/WCO.0b013e32832d93dd

PubMed Abstract | Crossref Full Text | Google Scholar

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23:11167. doi: 10.1523/JNEUROSCI.23-35-11167.2003

PubMed Abstract | Crossref Full Text | Google Scholar

Benayoun, M., Cowan, J. D., van Drongelen, W., and Wallace, E. (2010). Avalanches in a stochastic model of spiking neurons. PLoS Comput. Biol. 6:e1000846. doi: 10.1371/journal.pcbi.1000846

PubMed Abstract | Crossref Full Text | Google Scholar

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

PubMed Abstract | Crossref Full Text | Google Scholar

Capek, E., Ribeiro, T. L., Kells, P., Srinivasan, K., Miller, S. R., Geist, E., et al. (2023). Parabolic avalanche scaling in the synchronization of cortical cell assemblies. Nat. Comm. 14:2555. doi: 10.1038/s41467-023-37976-x

PubMed Abstract | Crossref Full Text | Google Scholar

Carvalho, T. T. A., Fontenele, A. J., Girardi-Schappo, M., Feliciano, T., Aguiar, L. A. A., Silva, T. P. L., et al. (2021). Subsampled directed-percolation models explain scaling relations experimentally observed in the brain. Front. Neural Circuits 14:576727. doi: 10.3389/fncir.2020.576727

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, Z. J., He, Y., Rosa-Neto, P., Germann, J., and Evans, A. C. (2008). Revealing modular architecture of human brain structural networks by using cortical thickness from MRI. Cereb. Cortex 18, 2374–2381. doi: 10.1093/cercor/bhn003

PubMed Abstract | Crossref Full Text | Google Scholar

Chialvo, D. (2010). Emergent complex neural dynamics. Nat. Phys. 6:744. doi: 10.1038/nphys1803

Crossref Full Text | Google Scholar

Cirunay, M. T., Batac, R. C., and Ódor, G. (2025). Learning and criticality in a self-organizing model of connectome growth. Sci. Rep. 15:31890. doi: 10.1038/s41598-025-16377-8

PubMed Abstract | Crossref Full Text | Google Scholar

Conte, D., and de Candia, A. (2025). Inferring global exponents in subsampled neural systems. iScience 28:113049. doi: 10.1016/j.isci.2025.113049

PubMed Abstract | Crossref Full Text | Google Scholar

Cota, W. Ódor, G., and Ferreira, S. C. (2018). Griffiths phases in infinite-dimensional, non-hierarchical modular networks. Sci. Rep. 8:9144. doi: 10.1038/s41598-018-27506-x

PubMed Abstract | Crossref Full Text | Google Scholar

Dalla Porta, L., and Copelli, M. (2019). Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: continuously varying exponents mimic M/EEG results. PLoS Comput. Biol. 15:e1006924. doi: 10.1371/journal.pcbi.1006924

PubMed Abstract | Crossref Full Text | Google Scholar

Danku, Z., Ódor, G., and Kun, F. (2018). Avalanche dynamics in higher-dimensional fiber bundle models. Phys. Rev. E 98:042126. doi: 10.1103/PhysRevE.98.042126

Crossref Full Text | Google Scholar

de Arcangelis, L., Perrone-Capano, C., and Herrmann, H. J. (2006). Self-organized criticality model for brain plasticity. Phys. Rev. Lett. 96:028107. doi: 10.1103/PhysRevLett.96.028107

PubMed Abstract | Crossref Full Text | Google Scholar

de Candia, A., Sarracino, A., Apicella, I., and de Arcangelis, L. (2021). Critical behavior of the stochastic wilson-cowan model. PLoS Comput. Biol. 17:e1008884. doi: 10.1371/journal.pcbi.1008884

Crossref Full Text | Google Scholar

Dobrinevski, A., Le Doussal, P., and Wiese, K. J. (2013). Statistics of avalanches with relaxation and barkhausen noise: a solvable model. Phys. Rev. E 88:032106. doi: 10.1103/PhysRevE.88.032106

PubMed Abstract | Crossref Full Text | Google Scholar

Durin, G., Papanikolaou, S., Bohn, F., Sethna, J. P., Sommer, R. L., and Zapperi, S. (2011). Universality beyond power laws and the average avalanche shape. Nat. Phys. 7:316. doi: 10.1038/nphys1884

Crossref Full Text | Google Scholar

Durin, G., and Zapperi, S. (2000). Scaling exponents for barkhausen avalanches in polycrystalline and amorphous ferromagnets. Phys. Rev. Lett. 84, 4705–4708. doi: 10.1103/PhysRevLett.84.4705

PubMed Abstract | Crossref Full Text | Google Scholar

Fontenele, A. J., de Vasconcelos, N. A. P., Feliciano, T., Aguiar, L. A. A., Soares-Cunha, C., Coimbra, B., et al. (2019). Criticality between cortical states. Phys. Rev. Lett. 122:208101. doi: 10.1103/PhysRevLett.122.208101

PubMed Abstract | Crossref Full Text | Google Scholar

Friedman, N., Ito, S., Brinkman, B. A., Shimono, M., DeVille, R. E., Dahmen, K. A., et al. (2012). Universal critical dynamics in high resolution neuronal avalanche data. Phys. Rev. Lett. 108:208102. doi: 10.1103/PhysRevLett.108.208102

PubMed Abstract | Crossref Full Text | Google Scholar

Gallos, L., Makse, H., and Sigmanb, M. (2012). A small world of weak ties provides optimal global integration of self-similar modules in functional brain networks. Proc. Nat. Acad. Sci. U.S.A. 109:2825. doi: 10.1073/pnas.1106612109

PubMed Abstract | Crossref Full Text | Google Scholar

Gillespie, D. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81:2340. doi: 10.1021/j100540a008

Crossref Full Text | Google Scholar

Girardi-Schappo, M. (2021). Brain criticality beyond avalanches: open problems and how to approach them. J. Phys. Complex. 2:031003. doi: 10.1088/2632-072X/ac2071

Crossref Full Text | Google Scholar

Girardi-Schappo, M., Brochini, L., Costa, A. A., Carvalho, T. T. A., and Kinouchi, O. (2020). Synaptic balance due to homeostatically self-organized quasicritical dynamics. Phys. Rev. Res. 2:012042(R). doi: 10.1103/PhysRevResearch.2.012042

Crossref Full Text | Google Scholar

Gireesh, E. D., and Plenz, D. (2008). Neuronal avalanches organize as nested theta- and beta/gamma-oscillations during development of cortical layer 2/3. Proc. Nat. Acad. Sci. U.S.A. 105, 7576–7581. doi: 10.1073/pnas.0800537105

PubMed Abstract | Crossref Full Text | Google Scholar

Girvan, M., and Newman, M. E. J. (2002). Community structure in social and biological networks. Proc. Nat. Acad. Sci. U.S.A. 99, 7821–7826. doi: 10.1073/pnas.122653799

PubMed Abstract | Crossref Full Text | Google Scholar

Gleeson, J. P., and Durrett, R. (2017). Temporal profiles of avalanches on networks. Nat. Comm. 8:1227. doi: 10.1038/s41467-017-01212-0

PubMed Abstract | Crossref Full Text | Google Scholar

Hengen, K. B., and Shew, W. L. (2025). Is criticality a unified setpoint of brain function? Neuron 113, 2582–2598.e2. doi: 10.1016/j.neuron.2025.05.020

PubMed Abstract | Crossref Full Text | Google Scholar

Houston, H. (2001). Influence of depth, focal mechanism, and tectonic setting on the shape and duration of earthquake source time functions. J. Geophys. Res. Solid Earth 106, 11137–11150. doi: 10.1029/2000JB900468

Crossref Full Text | Google Scholar

Houston, H., Benz, H. M., and Vidale, J. E. (1998). Time functions of deep earthquakes from broadband and short-period stacks. J. Geophys. Res. Solid Earth 103, 29895–29913. doi: 10.1029/98JB02135

Crossref Full Text | Google Scholar

Iyer, K. K., Roberts, J. A., Hellström-Westas, L., Wikström, S., Hansen Pupp, I., Ley, D., et al. (2015). Cortical burst dynamics predict clinical outcome early in extremely preterm infants. Brain 138, 2206–2218. doi: 10.1093/brain/awv129

PubMed Abstract | Crossref Full Text | Google Scholar

Kuntz, M. C., and Sethna, J. P. (2000). Noise in disordered systems: the power spectrum and dynamic exponents in avalanche models. Phys. Rev. B 62, 11699–11708. doi: 10.1103/PhysRevB.62.11699

Crossref Full Text | Google Scholar

Laurson, L., Illa, X., Santucci, S., Tore Tallakstad, K., Måløy, K. J., and Alava, M. J. (2013). Evolution of the average avalanche shape with the universality class. Nat. Comm. 4:2927. doi: 10.1038/ncomms3927

PubMed Abstract | Crossref Full Text | Google Scholar

Levina, A., Herrmann, J. M., and Geisel, T. (2007). Dynamical synapses causing selforganized criticality in neural networks. Nat. Phys. 3:857. doi: 10.1038/nphys758

Crossref Full Text | Google Scholar

Mariani, B., Nicoletti, G., Bisio, M., Maschietto, M., Vassanelli, S., and Suweis, S. (2022). Disentangling the critical signatures of neural activity. Sci. Rep. 12:10770. doi: 10.1038/s41598-022-13686-0

PubMed Abstract | Crossref Full Text | Google Scholar

Massobrio, P., de Arcangelis, L., Pasquale, V., Jensen, H. J., and Plenz, D. (2015). Criticality as a signature of healthy neural systems. Front. Syst. Neurosci. 9:22. doi: 10.3389/fnsys.2015.00022

PubMed Abstract | Crossref Full Text | Google Scholar

Mazzoni, A., Broccard, F. D., Garcia-Perez, E., Bonifazi, P., Ruaro, M. E., and Torre, V. (2007). On the dynamics of the spontaneous activity in neuronal networks. PLoS ONE 2:e439. doi: 10.1371/journal.pone.0000439

PubMed Abstract | Crossref Full Text | Google Scholar

Mehta, A. P., Mills, A. C., Dahmen, K. A., and Sethna, J. P. (2002). Universal pulse shape scaling function and exponents: critical test for avalanche models applied to barkhausen noise. Phys. Rev. E 65:046139. doi: 10.1103/PhysRevE.65.046139

PubMed Abstract | Crossref Full Text | Google Scholar

Meisel, C., Olbrich, E., Shriki, O., and Achermann, P. (2013). Fading signatures of critical brain dynamics during sustained wakefulness in humans. J. Neurosci. 33, 17363–17372. doi: 10.1523/JNEUROSCI.1516-13.2013

PubMed Abstract | Crossref Full Text | Google Scholar

Meunier, D., Lambiotte, R., Fornito, A., Ersche, K., and Bullmore, E. (2009). Hierarchical modularity in human brain functional networks. Front. Neuroinfo. 3:37. doi: 10.3389/neuro.11.037.2009

PubMed Abstract | Crossref Full Text | Google Scholar

Miller, S., Yu, S., and Plenz, D. (2019). The scale-invariant, temporal profile of neuronal avalanches in relation to cortical γ-oscillations. Sci. Rep. 9:16403. doi: 10.1038/s41598-019-52326-y

PubMed Abstract | Crossref Full Text | Google Scholar

Miller, S. R., Yu, S., Pajevic, S., and Plenz, D. (2021). Long-term stability of avalanche scaling and integrative network organization in prefrontal and premotor cortex. Netw. Neurosci. 5, 505–526. doi: 10.1162/netn_a_00188

PubMed Abstract | Crossref Full Text | Google Scholar

Millman, D., Mihalas, S., Kirkwood, A., and Niebur, E. (2010). Self-organized criticality occurs in non-conservative neuronal networks during ‘up' states. Nat. Phys. 6:801. doi: 10.1038/nphys1757

Crossref Full Text | Google Scholar

Moretti, P., and Muñoz, M. A. (2013). Griffiths phases and the stretching of criticality in brain networks. Nat. Commun. 4:2521. doi: 10.1038/ncomms3521

PubMed Abstract | Crossref Full Text | Google Scholar

Myrov, V., Suleimanova, A., Knapič, S., Liu, W., Palva, J. M., Partanen, P., et al. (2024). Hierarchical whole-brain modeling of critical synchronization dynamics in human brains. bioRxiv [Preprint]. doi: 10.1101/2024.05.08.593146

Crossref Full Text | Google Scholar

Newman, M. E. J. (2005). Power laws, pareto distributions and zipf's law. Contemp. Phys. 46, 323–351. doi: 10.1080/00107510500052444

Crossref Full Text | Google Scholar

Newman, M. E. J. (2006). Modularity and community structure in networks. Proc. Natl. Acad. Sci. U.S.A. 103, 8577–8582. doi: 10.1073/pnas.0601602103

PubMed Abstract | Crossref Full Text | Google Scholar

Ódor, G. (2016). Critical dynamics on a large human open connectome network. Phys. Rev. E 94:062411. doi: 10.1103/PhysRevE.94.062411

PubMed Abstract | Crossref Full Text | Google Scholar

Ódor, G., Dickman, R., and Ódor, G. (2015). Griffiths phases and localization in hierarchical modular networks. Sci. Rep. 5:14451. doi: 10.1038/srep14451

PubMed Abstract | Crossref Full Text | Google Scholar

Olami, Z., Feder, H. J. S., and Christensen, K. (1992). Self-organized criticality in a continuous, non-conservative cellular automaton modeling earthquakes. Phys. Rev. Lett. 68:1244–1247. doi: 10.1103/PhysRevLett.68.1244

Crossref Full Text | Google Scholar

Palva, J. M., Zhigalov, A., Hirvonen, J., Korhonen, O., Linkenkaer-Hansen, K., and Palva, S. (2013). Neuronal long-range temporal correlations and avalanche dynamics are correlated with behavioral scaling laws. Proc. Natl. Acad. Sci. U.S.A. 110, 3585–3590. doi: 10.1073/pnas.1216855110

PubMed Abstract | Crossref Full Text | Google Scholar

Pasquale, V., Massobrio, P., Bologna, L. L., Chiappalone, M., and Martinoia, S. (2008). Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience 153:1354. doi: 10.1016/j.neuroscience.2008.03.050

PubMed Abstract | Crossref Full Text | Google Scholar

Perković, O., Dahmen, K. A., and Sethna, J. P. (1995). Avalanches, barkhausen noise, and plain old criticality. Phys. Rev. Lett. 75, 4528–4531. doi: 10.1103/PhysRevLett.75.4528

PubMed Abstract | Crossref Full Text | Google Scholar

Petermann, T., Thiagarajan, T. C., Lebedev, M. A., Nicolelis, M. A., Chialvo, D. R., and Plenz, D. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Nat. Acad. Sci. U.S.A. 106:15921. doi: 10.1073/pnas.0904089106

PubMed Abstract | Crossref Full Text | Google Scholar

Piuvezam, H. C., Marin, B., Copelli, M., and Muñoz, M. A. (2023). Unconventional criticality, scaling breakdown, and diverse universality classes in the wilson-cowan model of neural dynamics. Phys. Rev. E 108:034110. doi: 10.1103/PhysRevE.108.034110

PubMed Abstract | Crossref Full Text | Google Scholar

Poil, S.-S., Hardstone, R., Mansvelder, H. D., and Linkenkaer-Hansen, K. (2012). Critical-state dynamics of avalanches and oscillations jointly emerge from balanced excitation/inhibition in neuronal networks. J. Neurosci. 32, 9817–9823. doi: 10.1523/JNEUROSCI.5990-11.2012

PubMed Abstract | Crossref Full Text | Google Scholar

Ponce-Alvarez, A., Jouary, A., Privat, M., Deco, G., and Sumbre, G. (2018). Whole-brain neuronal activity displays crackling noise dynamics. Neuron 100:1446. doi: 10.1016/j.neuron.2018.10.045

PubMed Abstract | Crossref Full Text | Google Scholar

Priesemann, V., Valderrama, M., Wibral, M., and Le Van Quyën, M. (2013). Neuronal avalanches differ from wakefulness to deep sleep - evidence from intracranial depth recordings in humans. PLoS Comput. Biol. 9:e1002985. doi: 10.1371/journal.pcbi.1002985

Crossref Full Text | Google Scholar

Roberts, J., Iyer, K., Finnigan, S., Vanhatalo, S., and Breakspear, M. (2014). Scale-free bursting in human cortex following hypoxia at birth. J. Neurosci. 34:6557. doi: 10.1523/JNEUROSCI.4701-13.2014

PubMed Abstract | Crossref Full Text | Google Scholar

Rusch, F. R., Kinouchi, O., and Roque, A. C. (2025). Influence of topology on the critical behavior of hierarchical modular neuronal networks. Commun. Phys. 8:168. doi: 10.1038/s42005-025-02074-5

Crossref Full Text | Google Scholar

Russo, R., Herrmann, H., and de Arcangelis, L. (2014). Brain modularity controls the critical behavior of spontaneous activity. Sci. Rep. 4:4312. doi: 10.1038/srep04312

PubMed Abstract | Crossref Full Text | Google Scholar

Scarpetta, S., Apicella, I., Minati, L., and de Candia, A. (2018). Hysteresis, neural avalanches, and critical behavior near a first-order transition of a spiking neural network. Phys. Rev. E 97:062305. doi: 10.1103/PhysRevE.97.062305

PubMed Abstract | Crossref Full Text | Google Scholar

Scarpetta, S., and de Candia, A. (2013). Neural avalanches at the critical point between replay and non-replay of spatiotemporal patterns. PLoS ONE 8:e64162. doi: 10.1371/journal.pone.0064162

PubMed Abstract | Crossref Full Text | Google Scholar

Scarpetta, S., Morisi, N., Mutti, C., Azzi, N., Trippi, I., Ciliento, R., et al. (2023). Criticality of neuronal avalanches in human sleep and their relationship with sleep macro- and micro-architecture. iScience 26:107840. doi: 10.1016/j.isci.2023.107840

PubMed Abstract | Crossref Full Text | Google Scholar

Sethna, J. P., Dahmen, K. A., and Myers, C. R. (2001). Crackling noise. Nature 410, 242–250. doi: 10.1038/35065675

Crossref Full Text | Google Scholar

Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009). Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci. 29, 15595–15600. doi: 10.1523/JNEUROSCI.3864-09.2009

PubMed Abstract | Crossref Full Text | Google Scholar

Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N., Smith, M. L., et al. (2013). Neuronal avalanches in the resting meg of the human brain. J. Neurosci. 33:7079. doi: 10.1523/JNEUROSCI.4286-12.2013

PubMed Abstract | Crossref Full Text | Google Scholar

Solovey, G., Miller, K. J., Ojemann, J. G., Magnasco, M. O., and Cecchi, G. A. (2012). Self-regulated dynamical criticality in human ecog. Front. Integr. Neurosci. 6:44. doi: 10.3389/fnint.2012.00044

PubMed Abstract | Crossref Full Text | Google Scholar

Sporns, O. (2010). Networks of the Brain. Cambridge, MA: MIT Press. doi: 10.7551/mitpress/8476.001.0001

Crossref Full Text | Google Scholar

Sporns, O., and Betzel, R. F. (2016). Modular brain networks. Annu. Rev. Psychol. 67, 613–640. doi: 10.1146/annurev-psych-122414-033634

Crossref Full Text | Google Scholar

Srinivasan, K., Ribeiro, T. L., Kells, P., and Plenz, D. (2024). The recovery of parabolic avalanches in spatially subsampled neuronal networks at criticality. Sci. Rep. 14:19329. doi: 10.1038/s41598-024-70014-4

PubMed Abstract | Crossref Full Text | Google Scholar

Stumpf, M. P. H., and Porter, M. A. (2012). Critical truths about power laws. Science 335, 665–666. doi: 10.1126/science.1216142

Crossref Full Text | Google Scholar

Tagliazucchi, E., Balenzuela, P., Fraiman, D., and Chialvo, D. R. (2012). Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis. Front. Physiol. 3:15. doi: 10.3389/fphys.2012.00015

PubMed Abstract | Crossref Full Text | Google Scholar

van Kampen, N. G. (2007). Stochastic Processes in Physics and Chemistry. Amsterdam: North Holland. doi: 10.1016/B978-044452965-7/50006-4

Crossref Full Text | Google Scholar

Xu, Y., Schneider, A., Wessel, R., and Hengen, K. B. (2024). Sleep restores an optimal computational regime in cortical networks. Nat. Neurosci. 27, 328–338. doi: 10.1038/s41593-023-01536-9

PubMed Abstract | Crossref Full Text | Google Scholar

Yamamoto, H., Spitzner, F. P., Takemuro, T., Buendía, V., Murota, H., Morante, C., et al. (2023). Modular architecture facilitates noise-driven control of synchrony in neuronal networks. Sci. Adv. 9:eade1755. doi: 10.1126/sciadv.ade1755

PubMed Abstract | Crossref Full Text | Google Scholar

Yu, S., Ribeiro, T. L., Meisel, C., Chou, S., Mitz, A., Saunders, R., et al. (2017). Maintained avalanche dynamics during task-induced changes of neuronal activity in nonhuman primates. eLife 6:e27119. doi: 10.7554/eLife.27119.012

PubMed Abstract | Crossref Full Text | Google Scholar

Yu, S., Yang, H., Nakahara, H., Santos, G. S., Nikolić, D., and Plenz, D. (2011). Higher-order interactions characterized in cortical activity. J. Neurosci. 31, 17514–17526. doi: 10.1523/JNEUROSCI.3127-11.2011

PubMed Abstract | Crossref Full Text | Google Scholar

Zaccariello, R., Herrmann, H. J., Sarracino, A., Zapperi, S., and de Arcangelis, L. (2025). Inhibitory neurons and the asymmetric shape of neuronal avalanches. Phys. Rev. E 111:024133. doi: 10.1103/PhysRevE.111.024133

PubMed Abstract | Crossref Full Text | Google Scholar

Zapperi, S., Castellano, C., Colaiori, F., and Durin, G. (2005). Signature of effective mass in crackling-noise asymmetry. Nat. Phys. 1:46. doi: 10.1038/nphys101

Crossref Full Text | Google Scholar

Zhigalov, A., Arnulfo, G., Nobili, L., Palva, S., and Palva, J. M. (2017). Modular co-organization of functional connectivity and scale-free dynamics in the human brain. Netw. Neurosci. 1, 143–165. doi: 10.1162/NETN_a_00008

PubMed Abstract | Crossref Full Text | Google Scholar