Symmetry breaking and avalanche shapes in modular neural networks

TL;DR

This paper explores how modular neural networks exhibit different activity patterns and critical behaviors, explaining experimental observations of symmetric and left-skewed neural avalanches.

Contribution

The study introduces a theoretical framework linking modular network structure to critical dynamics and distinct avalanche shapes.

Findings

The system shows symmetric and broken symmetry phases with distinct activity patterns.

Avalanche shapes vary systematically along critical lines, showing symmetric, right-skewed, or left-skewed distributions.

The model accounts for both symmetric and left-skewed neural avalanches observed in experiments.

Abstract

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.

Figures

Connectivity of the network. Only the first two modules, and the connections outgoing from neurons of the first module, are shown. Synaptic strengths are w0E and w0I for excitatory and inhibitory connections within the same module (intra-module), and w1E and w1I for excitatory and inhibitory connections between different modules (intermodule), defining the effective imbalances w0=w0E-w0I and w1=w1E-w1I. The connectivity matrix is symmetric under permutation of the modules, so that the connections outgoing from any module to any other module are equal.

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.

Avalanche statistics in the case of disconnected modules, or of a single module, for w0=w0E-w0I=1 (critical value), and different values of the parameter w0S=w0E+w0I=1,2,5,10. The case w0S=1 corresponds to a value w0I=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.

Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w0 = 1, w1 = 0, so that inter-module excitatory and inhibitory connections are exactly balanced, w1E=w1I. The value of w0S=w0E+w0I is fixed to w0S=10, while w1S=w1E+w1I varies between 1 and 10. For reference, we include in the plots also the case w1S=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.

Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w0 = 0.8, w1 = 0.1 (point “B” in Figure 2), that is for “excitatory imbalanced” inter-module connections. The value of w0S=w0E+w0I is fixed to w0S=10, while w1S=w1E+w1I 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.

Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w0 = 1, w1 = −0.15 (point “C” in Figure 2), that is for “inhibitory imbalanced” inter-module connections. The value of w0S=w0E+w0I is fixed to w0S=10, while w1S=w1E+w1I 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 w1S=7; (H) Same as (G) but for w1S=10; (I) Mean skewness of the avalanches as a function of the duration, for different number of modules, w1S=10.

(A–C) Examples of avalanches of duration T = 400 ms in the case of “excitatory imbalanced” inter-module connections (parameters w0 = 0.8, w1 = 0.1). (D–F) The same for “inhibitory imbalanced” inter-module connections (parameters w0 = 1, w1 = −0.15). In both cases w0S=w1S=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).

(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.

Equations

Full text

Introduction

1

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.

Wilson–Cowan model on a modular network

2

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(si), where si is the input of the i-th neuron that we are considering, and f(s) is an activation function taken equal to

$$[eqn]$$

The input of the i-th neuron is given by

$$[eqn]$$

where aj = 0, 1 if the j-th neuron is quiescent or active, wij is the strength of the connection going from j to i, and hi is the external input to neuron i. In the simplest case, one takes NE excitatory and NI inhibitory neurons, with NE+NI = N, and the connections wij having only two values, namely wE/NE if the pre-synaptic neuron j is excitatory, or −wI/NI if j is inhibitory.

It was shown in de Candia et al. (2021) that, when the external input hi goes to zero, the dynamics of the model undergoes a transition (bifurcation) at a critical value of the imbalance w0 = wE−wI, namely $[eqn]$ . For w0<w0c, the fraction of active neurons goes to zero when hi → 0, while for w0>w0c it is finite also for hi → 0, so that there is a self-sustained activity. Exactly at w0 = w0c, 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 NE excitatory and NI inhibitory neurons as schematically illustrated in Figure 1. Each neuron in module k receives input from excitatory neurons of module l with synaptic strengths $[eqn]$ , and from inhibitory neurons of module l with synaptic strengths $[eqn]$ , so that the input sk of neurons in module k is

$$[eqn]$$

where xl and yl are the fractions of excitatory and inhibitory active neurons in module l (activities), hk 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

$$[eqn]$$ $$[eqn]$$

where η_E, k(t), ηI, k(t) are normal, Gaussian independent white noises. The variables xk_ and yk can be written as $[eqn]$ , $[eqn]$ , where Xk and Yk are the “deterministic” components of the activities, that obey the deterministic equations

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ . The deterministic components evolve until they reach a stable (attractive) fixed point, that is a configuration of Xk and Yk that set to zero the right hand sides of Equation 5. We can then define $[eqn]$ , $[eqn]$ , 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

$$[eqn]$$ $$[eqn]$$

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 $[eqn]$ and $[eqn]$ 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

$$[eqn]$$ $$[eqn]$$

where δ_kl_ is the Kroneker delta, δ_kl_ = 1 if k = l, δ_kl_ = 0 otherwise, and $[eqn]$ , $[eqn]$ are evaluated at the fixed point. We then make the assumption that

$$[eqn]$$

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

$$[eqn]$$

where $[eqn]$ , $[eqn]$ . In the following, we consider the external input hk = 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 w0 and w1, 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 w0 and w1 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 hk = 0

$$[eqn]$$ $$[eqn]$$

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 $[eqn]$ . On the other hand, if m<M, by Equation 9b it has to be w1 ≤ 0. Therefore, for w1>0 and $[eqn]$ , there is only a solution with all the Σ_k_ greater than zero. We call this phase “SH,” which stands for “symmetric high” activity. For w1 ≤ 0 and $[eqn]$ , 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 $[eqn]$ such fixed points, corresponding to which modules are chosen to be active.

Connectivity of the network. Only the first two modules, and the connections outgoing from neurons of the first module, are shown. Synaptic strengths are w0E and w0I for excitatory and inhibitory connections within the same module (intra-module), and w1E and w1I for excitatory and inhibitory connections between different modules (intermodule), defining the effective imbalances w0=w0E-w0I and w1=w1E-w1I. 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 w0 and w1. The two solid lines, respectively $[eqn]$ for w1 ≤ 0, and $[eqn]$ for w1≥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.

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.

Critical dynamics at the edge of the low activity phase

3

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” $[eqn]$ , $[eqn]$ , 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 $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . We therefore, for given values of w0 and w1, also change the values of the sums $[eqn]$ and $[eqn]$ . As we have defined $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ to be non-negative, it has to be $[eqn]$ and $[eqn]$ .

In the following we consider M = 3, α = β = 0.1 ms^−1^, $[eqn]$ . 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 ri: ri = α if neuron i is active, or ri = f(si) if it is quiescent; (2) compute the sum over all neurons $[eqn]$ ; (3) draw a time interval dt from an exponential distribution with rate r; (4) choose the i−th neuron with probability ri/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.

Disconnected modules

3.1

In this subsection, we start by considering the case w1 = 0, $[eqn]$ , that corresponds to $[eqn]$ . 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 w0<w0c, and a phase of high self-sustained activity for w0>w0c, where $[eqn]$ 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 $[eqn]$ , that is w0c = 1 as we consider here α = β. For w0<w0c, the attractive fixed point corresponds to a quiescent network, while for w0>w0c the attractive fixed point corresponds to a finite value of the activity, so that the dynamics of the network is self-sustained. At w0 = w0c = 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 w0 = 1, w1 = 0 corresponds to the point “A” in Figure 2. Note that in the case of disconnected modules (w1 = 0, $[eqn]$ ), and in general in the case w1≥0 (balanced or “excitatory imbalanced” inter-module connections), there can be no “broken symmetry” phases Bm, that appear only for “inhibitory imbalanced” (w1 < 0) inter-module connections.

In Figure 3 we show the avalanche statistics for different values of $[eqn]$ . The case $[eqn]$ corresponds to a value $[eqn]$ 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 $[eqn]$ , 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 $[eqn]$ 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).

Avalanche statistics in the case of disconnected modules, or of a single module, for w0=w0E-w0I=1 (critical value), and different values of the parameter w0S=w0E+w0I=1,2,5,10. The case w0S=1 corresponds to a value w0I=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 $[eqn]$ we find nearly symmetric avalanches but strong deviations from MFDP exponents in the distribution of avalanches. On the other hand, for high values of $[eqn]$ , 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.

Modular network with balanced inter-module connections

3.2

In this subsection, we consider the case in which the system is at the point w0 = 1, w1 = 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 $[eqn]$ . The modules that constitute the network are therefore connected by excitatory and inhibitory connections greater than zero but exactly balanced, that is with $[eqn]$ . We fix the value $[eqn]$ , and consider different values of $[eqn]$ . As previously remarked, the fixed points of the dynamics, and their stability, depend only on the “imbalances” $[eqn]$ and $[eqn]$ , so that once w0 and w1 are fixed to the critical point, changing the values of $[eqn]$ and $[eqn]$ 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 $[eqn]$ and $[eqn]$ . 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 $[eqn]$ , and is confirmed in Figure 4 for a connected network formed by M = 3 modules and $[eqn]$ . Deviations from MFDP behavior are largest for small values of $[eqn]$ (small values of both excitatory and inhibitory inter-module connections), and tend to reduce for larger values of $[eqn]$ . 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 $[eqn]$ that shows avalanches slightly left-skewed for durations around T = 1, 000.

Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w0 = 1, w1 = 0, so that inter-module excitatory and inhibitory connections are exactly balanced, w1E=w1I. The value of w0S=w0E+w0I is fixed to w0S=10, while w1S=w1E+w1I varies between 1 and 10. For reference, we include in the plots also the case w1S=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.

Modular network with “excitatory imbalanced” inter-module connections

3.3

If we shift the balance of inter-module connections making them excitatory imbalanced, that is take w1>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 w1≥0 is given by the equation

$$[eqn]$$

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

Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w0 = 0.8, w1 = 0.1 (point “B” in Figure 2), that is for “excitatory imbalanced” inter-module connections. The value of w0S=w0E+w0I is fixed to w0S=10, while w1S=w1E+w1I 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.

Modular network with “inhibitory imbalanced” inter-module connections

3.4

As shown in Figure 2, if inter-module connections are “inhibitory imbalanced,” that is $[eqn]$ and w1 < 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 w0 = 1 toward phase B1, where only one of the modules is active. In Figure 6 we show avalanche statistics for w0 = 1, w1 = −0.15, that is at point “C” in Figure 2, for $[eqn]$ and different values of $[eqn]$ between 0.5 and 10.

Avalanche statistics in the case of a network constituted by M = 3 modules, at the critical point w0 = 1, w1 = −0.15 (point “C” in Figure 2), that is for “inhibitory imbalanced” inter-module connections. The value of w0S=w0E+w0I is fixed to w0S=10, while w1S=w1E+w1I 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 w1S=7; (H) Same as (G) but for w1S=10; (I) Mean skewness of the avalanches as a function of the duration, for different number of modules, w1S=10.

As shown in Figures 6A–C, in this case, differently from the case w1>0, the largest deviations from MFDP are observed at large values of $[eqn]$ , 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 $[eqn]$ , as shown in Figure 6D. Figures 6E, F show the average shape for avalanches of duration 50 ms and 400 ms. For $[eqn]$ 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 $[eqn]$ and $[eqn]$ . 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.

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

3.5

In Figure 7 we show different examples of avalanches of duration T = 400 ms, in the case of “excitatory imbalanced” (upper row, parameters w0 = 0.8, w1 = 0.1) and “inhibitory imbalanced” (lower row, parameters w0 = 1, w1 = −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〉 = 〈(xk−yk)/2〉 between the excitatory activity xk_ (fraction of active excitatory neurons) and the inhibitory activity yk, on the k-th module, and the average total activity 〈Σ_k〉 = 〈(xk+yk)/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 w0 = 0.8, w1 = 0.1 (w0 = 1, w1 = −0.15). From these, one can compute the “intra-module” and “inter-module” input of neurons. The total input sk_ of a neuron in module k, given by Equation 3, can indeed be written as sk = s0, k+s1, k+hk, where

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

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

$$[eqn]$$

The two contributions 〈s0, k〉 and 〈s1, 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.

(A–C) Examples of avalanches of duration T = 400 ms in the case of “excitatory imbalanced” inter-module connections (parameters w0 = 0.8, w1 = 0.1). (D–F) The same for “inhibitory imbalanced” inter-module connections (parameters w0 = 1, w1 = −0.15). In both cases w0S=w1S=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).

(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.

Conclusions

4

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.