PHYSICAL REVIEW E 106, 024304 (2022)

Scaling of avalanche shape and activity power spectrum in neuronal networks

Manoj Kumar Nandi,1 Alessandro Sarracino ,1 Hans J. Herrmann,2,3 and Lucilla de Arcangelis1,*
1
Department of Engineering, University of Campania “Luigi Vanvitelli”, 81031 Aversa, Caserta, Italy
2
PMMH, ESPCI, 7 Quai Saint Bernard, Paris 75005, France
3
Departamento de Fisica, Universidade Federal do Ceará, 60451-970 Fortaleza, Ceará, Brazil

(Received 2 March 2022; revised 17 June 2022; accepted 9 August 2022; published 26 August 2022)

Many systems in nature exhibit avalanche dynamics with scale-free features. A general scaling theory has
been proposed for critical avalanche profiles in crackling noise, predicting the collapse onto a universal avalanche
shape, as well as the scaling behavior of the activity power spectrum as Brown noise. Recently, much attention
has been given to the profile of neuronal avalanches, measured in neuronal systems in vitro and in vivo. Although
a universal profile was evidenced, confirming the validity of the general scaling theory, the parallel study of
the power spectrum scaling under the same conditions was not performed. The puzzling observation is that
in the majority of healthy neuronal systems the power spectrum exhibits a behavior close to 1/ f , rather than
Brown, noise. Here we perform a numerical study of the scaling behavior of the avalanche shape and the power
spectrum for a model of integrate and fire neurons with a short-term plasticity parameter able to tune the system to
criticality. We confirm that, at criticality, the average avalanche size and the avalanche profile fulfill the general
avalanche scaling theory. However, the power spectrum consistently exhibits Brown noise behavior, for both
fully excitatory networks and systems with 30% inhibitory networks. Conversely, a behavior closer to 1/ f noise
is observed in systems slightly off criticality. Results suggest that the power spectrum is a good indicator to
determine how close neuronal activity is to criticality.

DOI: 10.1103/PhysRevE.106.024304

I. INTRODUCTION Gutenberg-Richter [21] for earthquake occurrence in the 50s

and became a general paradigm after the seminal work of Bak
Besides the well-known rich phenomenology, i.e., insta-
[22,23]. Within this context, a real breakthrough was the for-
bilities and metastability transitions [1–3], synchronization
mulation of a general scaling theory for avalanche phenomena
[4,5], the presence of multiple spatiotemporal scales [6,7],
at the critical point, encompassing the scaling behavior of
and long-range temporal correlations [8], spontaneous brain
most relevant properties in the process. This scaling theory,
activity exhibits bursts of activity, named neuronal avalanches.
initially formulated for the Barkhausen noise [24,25], has
These were first detected in vitro in organotypic cultures from
turned out to be extremely general and found in a variety of
coronal slices of rat cortex [9] and in dissociated neurons
different phenomena, from plastic deformations [26] to earth-
from rat hippocampus and cortex [10,11] or leech ganglia
quakes [27]. Among a number of scaling relations for different
[10]. Next, neuronal avalanches were identified in vivo in rat
quantities [25], some scaling laws have been considered in
cortical layers during early postnatal development [12], in
the literature as an indicator to determine if a system acts
the cortex of awake adult rhesus monkeys [13], as well as
at the critical point. In particular, the scaling of the average
in the resting state of the human brain by means of nonin-
avalanche size S versus its duration T [24,25,28], S ∼ T γ ,
vasive techniques such as magnetoencephalography (MEG)
with the exponent γ related to the exponents of the avalanche
[14]. The statistical properties of neuronal avalanches have
distributions:
been the object of intensive investigation both experimentally
and numerically [15–17], focusing mainly on the scaling be- τ −1
γ = . (1)
havior of the avalanche size and duration distributions, P(S) α−1
and P(T ). Consistent evidence indicates a scaling behavior
This exponent γ is also predicted to control the collapse of the
in the universality class of the mean-field branching model
profile of avalanches with different durations, as well as the
[18,19], namely, P(S) ∝ S −α with α  1.5 and P(T ) ∝ T −τ
scaling of the signal power spectral density (PSD), S( f ) ∼
with τ  2.0, even if the existence of a different universality
f −γ [25]. Moreover, other features of the avalanche profile
class has been also proposed in the literature [20].
have received a wide interest. In particular, the avalanche
Avalanching is a widespread phenomenon, occurring in
shape, not necessarily symmetric in the scaling theory [24],
systems where many degrees of freedom interact under
has been found to depend on the interaction kernel; namely,
slow drive, which had, as the first prototype, the model of
asymmetry appears when the interaction is not fully nonlocal,
reflecting broken time-reversal symmetry in the avalanche
dynamics [29]. Symmetric profiles are found in Barkhausen
*
lucilla.dearcangelis@unicampania.it noise [24,25] and in mean-field systems [28]; however,

2470-0045/2022/106(2)/024304(8) 024304-1 ©2022 American Physical Society

MANOJ KUMAR NANDI et al. PHYSICAL REVIEW E 106, 024304 (2022)

considering inertial effects leads to leftward asymmetry (pos- Here we address the issue of scaling for neuronal
itive skewness) [30], whereas rightward asymmetry (negative avalanches within the context of an integrate and fire neuronal
skewness) is observed for increasing interaction range [29]. network model including short- and long-term plasticity [17].
Interestingly, the symmetry of avalanche shape in experimen- The model presents the advantage of being able to tune the
tal neuronal avalanches shows a variety of features: from an system at and far from criticality, allowing the investigation of
almost symmetric shape in cortex slices [31] and in nonhuman the scaling behavior in a wider region of phase space. In par-
primates [32], provided that the modulation of γ oscillations ticular, we study the scaling properties of the avalanche shape
is carefully taken into account, to a leftward asymmetric shape in parallel with the activity PSD, to identify the conditions
in zebrafish larvae [33]. This confirms that the symmetry of for the validity of the avalanche scaling theory. The analysis
the shape is not a necessary requirement for criticality, as is performed for different percentages of inhibitory neurons,
already observed in Ref. [24]. unveiling their crucial role in the PSD scaling.
Interestingly, experimental and numerical studies aimed
at verifying the validity of the avalanche scaling theory for
neuronal avalanches have mainly focused on the scaling of
II. NETWORK MODEL
the average size with duration and the collapse of avalanche
profiles onto a universal curve [31–34]. Little attention has We consider a neuronal network consisting of N neu-
been given to the parallel investigation of the scaling be- rons, randomly placed in a cubic lattice and connected by
havior of the activity power spectrum, which is predicted in directed synapses, with a fraction, Pin , of inhibitory neurons.
Ref. [25] to exhibit Brown noise behavior, i.e., S( f ) ∼ f −β , In order to construct the network of synaptic connections,
with β = γ = 2 [Eq. (1)]. A variety of experimental studies since the out-degree distribution is not fully characterized for
on different signals, such as electroencephalography (EEG), the morphological connections in real brains, the outgoing
MEG, resting-state functional magnetic resonance imaging, connections are assigned according to the scale-free degree
−2
and the local field potentials (LFPs) of spontaneous cortical distribution P(kout ) ∝ kout found experimentally for func-
activity, have evidenced the presence of effective power-law tional networks [45], with kout ∈ [2, 100]. The connections
regimes in the PSD [35–39], which are the background for between two neurons are established using a distance-
peaks at characteristic frequencies corresponding to different dependent probability, P(r) ∝ e−r/r0 , where r is the Euclidean
brain modes. In the majority of studies, a behavior closer distance between two neurons and r0 = 5 is a characteristic
to 1/ f , rather than Brown, noise, has been detected. For connectivity range [46]. We assign to each synaptic connec-
instance, the PSD scaling behavior in the human eyes-closed tion between neuron i and j an initial random strength, gi j ,
and eyes-open resting EEG [38] provides at low frequencies uniformly chosen in the interval [0.4,0.6]. A recent study
(in the range 0.5–8 Hz) β ∼ 1.32 for the eyes-closed condition [47] has shown that in scale-free functional networks the
and β ∼ 1.27 for the eyes-open condition, with the exponent inhibitory neurons are typically hubs; we therefore choose the
varying across brain regions with a standard error of 20%. A inhibitory neurons among the highly connected ones, kout > 5.
similar analysis of both EEG and MEG signals has recently We implement in the model two plasticity mechanisms, the
proposed [37] an average slope of β = 1.06 ± 0.29, varying short- and long-term plasticity. The first one models the dy-
in different brain areas. Conversely, the exponent value β ∼ 2 namics of synaptic resources that are used in all firing events
is typically found in epileptic patients, in the range [2.2, 2.44] and need to be restored in the presynaptic terminal. By tun-
in the awake state and [1.6, 2.87] in the slow wave sleep [40]. ing the efficiency in restoring neurotransmitter resources, the
In addition, a numerical study of a self-organized model for model tunes the network at the critical point. This plasticity
neuronal activity has evidenced that the exponent β depends mechanism is always active during the avalanche dynamics.
on the fraction of inhibitory neurons, crossing over from Conversely, the long-term plasticity is a Hebbian mechanism
Brown noise for fully excitatory networks to 1/ f behavior for that models the experience of the neuronal systems by sculpt-
30% inhibitory neurons [41], typical of mammal brains. ing the synaptic strengths according to their activity. This
Finally, it is important to remark that there is evidence plasticity mechanism mimics the age of the neuronal system
in the literature [42,43] that the observed scaling behavior and modifies the synaptic connections starting from an initial
for avalanches can be also obtained without invoking the random configuration. For this reason, the long-term plasticity
criticality hypothesis. Recent investigation has been triggered is active only during an initial training period, during which
by the results in Ref. [20] showing that data, ranging from avalanche measurements are not monitored, in order to adapt
freely moving to anesthetized mammals, surprisingly show the synaptic strengths according to previous activity. The im-
scaling behavior in a specific activity regime, extending orig- plementation of both plasticity measurements, acting on very
inal results which focused on spontaneous brain activity. It different timescales (milliseconds for the first measurement
has, therefore, been suggested that power laws can emerge and up to years for the second measurement) is motivated by
even in models which are not critical [43] and do not sat- the need to include in the model fundamental experimental ev-
isfy Eq. (1). Therefore, the criteria implemented in Ref. [20] idence for neuronal dynamics. More precisely, each neuron i is
to confirm criticality would be not sufficient to discrimi- characterized by a membrane potential vi , whose initial values
nate critical from noncritical systems. Indeed, the emergence are distributed uniformly in the interval [0.5,1.0]. At each
of power-law distributions has been recently related to time step all neurons with a potential exceeding a threshold
the presence of fluctuating variables in the process, such value of vi  vc = 1.0 fire, leading to the release of a fraction
as input stimuli leading to neuronal firings at different rates of neurotransmitter δu at all synapses and a change in the
[44]. potential of the connected postsynaptic neurons j according

024304-2
SCALING OF AVALANCHE SHAPE AND ACTIVITY … PHYSICAL REVIEW E 106, 024304 (2022)

to the following equations [17,48]: III. DISTRIBUTION OF AVALANCHE SIZES AND

DURATIONS
v j (t + 1) = v j (t ) ± vi (t )ui (t )δugi j,
As mentioned in the description of the model, it is possi-
ui (t + 1) = ui (t )(1 − δu), ble to tune the system in different activity states, subcritical,
vi (t + 1) = 0, critical, and supercritical, by adjusting the parameter δurec
[17,48] controlling the efficiency in neurotransmitter recov-
where the ± sign stands for an excitatory or inhibitory ery. Since inhibitory neurons hamper the activity propagation
presynaptic neuron and ui is the amount of releasable neuro- and limit large avalanches, by increasing their fraction the
transmitter at neuron i, which is initially set equal to 1 for all system will move away from the critical state into the sub-
neurons. After a neuron fires, it enters into a refractory state critical region. Here we consider networks in the critical
for the tr = 1 time step, during which it is unable to elicit any state, for which δurec is appropriately increased for increasing
signal. The activity propagates and stops when the potential fractions of inhibitory neurons. More precisely, we slowly
of all neurons is below vc . The present model is the discrete increase the value of δurec , moving from a subcritical regime
time version of the integrate and fire model of Ref. [16]. Here (exponential size distribution) to a regime where a power-
the unit time step is of the order of 10 ms and represents the law behavior is detected in the size distribution. The critical
time elapsed between the elicitation of the action potential and value of δurec is identified as the largest value providing a
the change in membrane potential at the postsynaptic neuron. power law with the cutoff as close as possible to N. We also
The amount of available neurotransmitter at each neuron ui is verify that no bump appears in the tail of the distribution,
called the readily releasable pool [49,50] and only a 5% frac- a sign of an excess of large avalanches. In previous studies
tion is released at each firing [51,52] (δu = 0.05), according to it has been verified that activity can be tuned to be gen-
short-term synaptic plasticity. When all neurons are below the uinely critical; i.e., the cutoff in the scale-free avalanche size
firing threshold, activity stops and can be triggered again by distribution correctly scales with the system size N [17,48].
setting any random neuron to its threshold value, generating The distributions of avalanche sizes P(S) and durations P(T )
another burst of firing neurons, called the neuronal avalanche. are shown in Fig. 1 for different percentages of inhibitory
We measure for each avalanche the size S as the number neurons, each time tuning δurec to the critical state. The
of firing neurons, independently of their firing rate, and the size distribution exhibits a power-law behavior P(S) ∝ S −α ,
duration T as the number of time steps in the activity propaga- with α  1.5, quite independently of the percentage of in-
tion. During an avalanche, because of short-term plasticity, the hibitory neurons Pin , followed by an exponential cutoff which
available neurotransmitter decreases constantly. Therefore to moves towards smaller sizes of S as we increase the per-
sustain further activity a recovery mechanism is needed. Since centage of inhibitory neurons, hindering the occurrence of
the synaptic recovery takes a time of the order of seconds large avalanches. The distribution of durations also shows
[53–55], whereas synaptic transmission acts on the scale of a power-law behavior P(T ) ∝ T −τ , with τ  2.0, followed
milliseconds [9,56], the available neurotransmitter recovery by an exponential cutoff at long avalanche durations de-
is implemented at the end of each avalanche propagation as pending on Pin . Both scaling exponents are in agreement
ui = ui + δurec for all neurons. Here δurec is a tunable param- with experimental values [9,13,14,58] and consistent with
eter which determines whether the network is in a subcritical, the mean-field branching model universality class [18] (see
a critical, or a supercritical state; i.e., in the critical state the Table I).
cutoff in the scale-free avalanche size distribution correctly The critical exponents for the avalanche distributions have
scales with the system size N [17,48], showing that the power been also found in a self-organized neuronal network model
law is due to criticality. without short-term plasticity, namely, in the absence of any
The synaptic network structure is sculpted by the activity- tuning parameter [15,41]. For this model with different Pin
dependent long-term plasticity, following the principles of the size distribution has been found to follow the scaling
Hebbian plasticity [57]. Whenever a neuron i fires, all synap- form P(S) = S −α f (S/Pin−θ ), where α = 1.5 and the cutoff
tic strengths gi j to postsynaptic neurons j are strengthened size scales with Pin with an exponent θ  2.2 for scale-free
proportionally to the potential variation induced in the post- networks. In self-organized models the fraction of inhibitory
synaptic neuron, as gi j (t + 1) = gi j (t ) + δgi j , where δgi j = neurons strongly affects the extension of the scaling regime;
ε|(v j (t + 1) − v j (t ))| and ε = 0.04 is a parameter controlling however, data collapse confirms the value of α by appropri-
the strength of plastic adaptation. Conversely, at the end of ately considering the cutoff scaling. To understand the role
each avalanche, all gi j are reduced by the average increase in of inhibitory neurons in the present model, which is not self-
synaptic strength per bond, gi j = gi j − δgi j /Ns , where Ns organized but tuned at criticality, we first determine the value
is the total number of synapses in the network. Synaptic con- of δurec which sets a purely excitatory system (Pin = 0%) in
nections whose strength becomes smaller than gmin = 10−5 the critical state. Then we fix the value of δurec and increase
are pruned, i.e., permanently removed from the network. In the percentage of inhibitory neurons in the system. By doing
order to modulate the initially random strengths, the long-term so, the system moves away from the critical point, into the
plastic adaptation is either implemented for a fixed number subcritical regime, which leads to a decrease of the scaling
of stimulations Np = 10 000 or stopped at the first synaptic regime in the distributions with the exponential cutoff moving
pruning. We have studied systems with N = 16 000 neurons towards smaller avalanche sizes [Fig. 2(a)]. Similar behavior
in a cube of side L = 100 and data are averaged over 5000 is also observed for the distribution of avalanche durations
network configurations. [Fig. 2(b), see Table I].

024304-3
MANOJ KUMAR NANDI et al. PHYSICAL REVIEW E 106, 024304 (2022)

(a) 0 (a) 0
10 10
Pin% Pin%
-1 -1
10 0 10 0
5 5
P(S)-2 10 P(S) 10
-2
10 15 10 15
20 20
-3 25 -3 25
10 30 10 30
-1.5 -1.5
S S
-4 -4
10 10
-5 -5
10 10

-6 -6
10 10

-7 -7
10 10
1 2 3 4 1 2 3 4
10 10 S 10 10 10 10 S 10 10
(b) 0 (b) 100
10
Pin% Pin%
-1 -1
10 0 10 0
5 5
-2 10 -2 10
10 15 10 15
P(T) 20 P(T) 20
25
-3 25 -3
10 30 10 30
-2 -2
T T
-4 -4
10 10

-5 -5
10 10

-6 -6
10 10

-7 -7
10 10
1 2 1 2
10 10 10 T 10
T
FIG. 1. Avalanche size (a) and duration (b) distributions for 5000 FIG. 2. Avalanche size (a) and duration (b) distributions for 5000
configurations of a network of N = 16 000 neurons for Pin = 0, 5, 10, configurations of a network of N = 16 000 neurons for Pin = 0, 5, 10,
15, 20, 25, and 30% in the critical state. 15, 20, 25, and 30% for fixed δurec = 0.001 value corresponding to
the critical state for Pin = 0%.

As for the self-organized critical model [15,41], upon

implementation of short-term plasticity but with no tuning where δurec (Pin ) and δurec (Pin = 0) are the parameter values
parameter, the avalanche size distribution exhibits (Fig. 3) the for the system with a fraction Pin of inhibitory neurons and
scaling behavior for a fully excitatory system to be in the critical state, re-
spectively. The argument of the scaling function f (x), x =
P(S) ∝ S −α f {S[Pin δurec (Pin )/δurec (Pin = 0)]θ }, (2) S[Pin δurec (Pin )/δurec (Pin = 0)]θ , now takes into account that

TABLE I. Avalanche exponents for different Pin . The exponents τ and α have been obtained by fitting the distributions with a power law
in the critical case and by a power law truncated by an exponential in the subcritical case. The error on γ is obtained by propagation of
uncertainty. The exponent of the spectrum S( f ) is obtained by fitting with a power law in the high-frequency regime. The error on γcoll is
obtained by exploring the range of exponent values providing good collapse.
τ −1
Pin δurec α τ γ = α−1
γcoll γS γS( f )

0% 0.0010 1.50 ± 0.05 2.05 ± 0.05 2.1 ± 0.2 2.08 ± 0.02 2.10 ± 0.05 1.98 ± 0.05
10% 0.0014 1.50 ± 0.05 2.07 ± 0.05 2.1 ± 0.2 2.14 ± 0.01 2.14 ± 0.05 1.95 ± 0.05
20% 0.0023 1.50 ± 0.05 2.09 ± 0.05 2.2 ± 0.2 2.16 ± 0.02 2.3 ± 0.2 1.95 ± 0.05
30% 0.0037 1.50 ± 0.05 2.15 ± 0.05 2.3 ± 0.3 2.32 ± 0.03 2.4 ± 0.2 2.00 ± 0.05
10% 0.0010 1.50 ± 0.05 2.1 ± 0.1 2.2 ± 0.3 2.11 ± 0.02 2.15 ± 0.1 1.95 ± 0.05
20% 0.0010 1.5 ± 0.1 1.9 ± 0.2 1.8 ± 0.5 2.18 ± 0.04 2.15 ± 0.1 1.80 ± 0.05
30% 0.0010 1.7 ± 0.2 1.7 ± 0.2 1.0 ± 0.5 2.08 ± 0.04 2.1 ± 0.1 1.3 ± 0.1

024304-4
SCALING OF AVALANCHE SHAPE AND ACTIVITY … PHYSICAL REVIEW E 106, 024304 (2022)

0
10

-1
10
S P(S)

Pin
1.5

-2
10 0.15
0.18
0.20
0.23
0.25
0.28
-3 0.30
10

0 1 2 3
10 10 10 10
2.4
S(PinGurec(Pin)/Gurec(Pin=0))
FIG. 4. Average avalanche size S vs duration T . The slope for
FIG. 3. Collapse of avalanche size distributions off
all values of Pin is close to 2. Inset: The same quantity for systems
criticality onto a universal curve by plotting S α P(S) vs
not tuned at criticality, i.e., for fixed δurec = 0.001 and Pin = 10, 20,
S[Pin δurec (Pin )/δurec (Pin = 0)]θ , with α = 1.5 and θ = 2.4.
and 30%.

the activity is driven in the subcritical regime not only by in Fig. 4). Interestingly, the scaling S ∼ T 2 is still observed
increasing the fraction of inhibitory neurons but also by not but there is no agreement with the Sethna relation Eq. (1).
tuning δurec : The ratio δuδurecrec(P(P in )
in =0)
is a measure of how far the
system is from the critical state. Interestingly, the values of
IV. SHAPE OF AVALANCHE PROFILE
the scaling exponents are in good agreement with those found
for the self-organized model, namely, α  1.50 ± 0.01 and Next, we analyze the shape of the avalanche profiles
θ  2.4 ± 0.1. This result suggests that the extension of the for fixed avalanche duration in the critical state for dif-
scaling regime depends solely on the distance from the critical ferent percentages of inhibitory neurons. To this end, we
state and not on the different mechanisms implemented in the consider avalanches, whose sizes are in the scaling regime
neuronal dynamics. and have a duration in the range T ± 3, and we plot the
Within the context of the crackling noise, one has the average number of firing neurons during the avalanche as
scaling relation S ∼ T γ , with γ given in Eq. (1). A simple function of the rescaled time t/T . At criticality, the avalanche
argument to justify this relation is the following [59]: Con- profiles are nearly parabolic for all durations [Fig. 5(a)]. Ac-
sider an avalanche with size S  and duration T  in the scaling cording to the scaling of the average avalanche size with
regime of the distributions. The probability to observe an duration, it is possible to rescale the axes to collapse the
avalanche smaller than S  is P(S < S  ) ∼ S 1−α , analogously curves for different T onto a universal profile, according
P(T < T  ) ∼ T 1−τ . If S  is the average size of an avalanche to S(t, T ) = T γ −1 f (t/T ). The best collapse would then
of duration T  and for a narrow distribution of avalanche provide an independent evaluation of the exponent γ and
sizes with fixed duration, then P(S < S  ) ≈ P(T < T  ) and therefore a validation of the avalanche exponents at criticality.
T γ (1−α) ∼ T 1−τ , leading to Eq. (1). This relation holds for In Figs. 5(b)–5(f) we show the collapse of avalanche profiles
avalanche sizes and durations in the critical scaling regime with different durations for systems with various percentages
and is very robust. It has been extended to a variety of other of inhibitory neurons. The collapse is obtained in all cases
avalanche processes [29,60] and verified experimentally in with an exponent γcoll  2 (see Table I) and the shape of the
neuronal systems in vitro and in vivo [31–33]. Numerically, profile is well fitted by a parabolic function for all values of
Eq. (1) has been recently validated in networks of integrate Pin . Interestingly, the same analysis performed off criticality
and fire neurons [59] and for the Wilson-Cowan model with by fixing δurec = 0.001, although confirming the value γcoll 
different populations of excitatory and inhibitory neurons 2, exhibits a profile not compatible with a parabolic function.
[34]. We have then fitted the profiles with the asymmetric function
In Fig. 4 we show the scaling of the average size with [29]
fixed duration vs T for avalanches in the power-law regime
S/T γ −1 ≈ [t/T (1 − t/T )]γ −1 [1 − a(t/T − 0.5)], (3)
for systems with different percentage of inhibitory neurons.
The separate fit of the different curves shows that the slope where a is the measure of asymmetry. The best collapse
for all systems is very close to 2, the value expected on the [Fig. 5(f)] is obtained for a  −0.25 for 20% inhibitory
basis of Eq. (1) for α = 1.5 and τ = 2.0 (see Table I). Inter- neurons, with rightward asymmetry increasing for larger Pin
estingly, this scaling is also verified for systems off criticality, (a  −0.5 for Pin = 30%). Interestingly, a similar rightward
except for fractions of inhibitory neurons larger than 20%. We asymmetric avalanche shape is obtained in the Wilson-Cowan
observe in this case that avalanches with a long duration are model for different percentages of inhibitory neurons if the
not observed, which strongly reduces the scaling regime (inset system is set in the critical state by tuning an appropriate

024304-5
MANOJ KUMAR NANDI et al. PHYSICAL REVIEW E 106, 024304 (2022)

FIG. 5. (a) Average avalanche size S vs time for different durations in the critical state for Pin = 0. The behavior of the plots for other
inhibitory percentages are similar (not shown here). (b)–(e) The quantity S/T γcoll −1 versus t/T for different Pin . The average avalanche shape
for different durations collapses onto a universal profile. The green dashed line is a parabolic fit to the data. (f) The quantity S/T γcoll −1 versus
t/T for fixed δurec = 0.001 and Pin = 20%. The dashed green line shows the fit with Eq. (3).

parameter [34]. We have tabulated the exponents obtained inset, the effective exponent β continuously varies with Pin ,
from different methods in Table I. from β  2.0 for purely excitatory systems to β  1.0 for
30% inhibitory neurons, the fraction typically found in mam-
V. POWER SPECTRAL DENSITY mal brains (see Table I). Moreover, the crossover to white
Next we analyze the behavior of the activity PSD in critical noise moves to larger frequencies for increasing Pin , corre-
and subcritical regimes, searching for the best agreement with sponding to the inverse of the longest avalanche duration in
experimental data. We consider the temporal sequence of the the system. The behavior observed in Fig. 6(b) has been also
number of firing neurons a(t ) and define the power spectrum evidenced in the self-organized model for neuronal activity
as S( f ) = â( f )â∗ ( f ), where â( f ) is the discrete Fourier trans- [41]. To fully characterize the properties of the PSD in the
form of a(t ): entire phase diagram, we have systematically tuned the pa-
T −1
rameter δurec over a range of values at and off criticality for
 different Pin and measured the exponent β. The contour plot
â( f ) = a(t )e−2iπ f t/T . (4)
in Fig. 7 suggests an interesting behavior: At criticality (con-
t=0
tinuous black line) and in the supercritical phase, Brown noise
We first analyze the networks at criticality, where δurec is (β ∼ 2) is always found, which is in fact the value measured in
adjusted to the values in Table I for different Pin . Figure 6(a) epileptic patients. Conversely, moving away from the critical
shows the scaling behavior S( f )  f −β which exhibits a sta- line, into the subcritical region, provides β values closer to
ble critical exponent β  2.0 for all percentages of inhibitory experimental data for healthy patients. The model, therefore,
neurons, with a crossover towards white noise at small fre- suggests that systems in healthy conditions act slightly off
quencies. This result is fully in agreement with the scaling criticality, in the subcritical regime.
expected for crackling noise [24,25], where the theoretical
prediction is S( f )  f −γ for avalanches at criticality. Un-
VI. CONCLUSIONS
fortunately, this result does not comply with experimental
evidence, which rather exhibits an exponent β close to unity In neuronal systems inhibition has a crucial role in keeping
for a variety of healthy neuronal systems. activity balanced under healthy conditions. Inhibitory neurons
Next we evaluate the PSD for the activity signal in net- act as sinks, hampering activity propagation; therefore, their
works not tuned at criticality but fixing the parameter δurec presence affects avalanche dynamics. In particular, the open
at the value for the critical state of Pin = 0%. Results shown question is how their dissipative role can affect the predictions
in Fig. 6(b) confirm an effective power-law behavior over a by the general scaling theory for avalanches, initially formu-
scaling regime decreasing with Pin and a crossover towards lated for crackling noise [24]. Among a variety of scaling
white noise at small frequencies. Moreover, as shown in the relations we investigate the robustness of the predictions on

024304-6
SCALING OF AVALANCHE SHAPE AND ACTIVITY … PHYSICAL REVIEW E 106, 024304 (2022)

(a) 0.003

0.0025
-0.8

-1
0.002

Gurec
-1.2

-1.4
0.0015
-1.6

-1.8
0.001

-2

0.0005
0 5 10 15 20 25 30 35

(b) Pin

FIG. 7. Contour plot of the β exponent values for different δurec

and Pin . The value of the exponent is given by the color code. The
continuous black line is the critical line, with the supercritical (sub-
critical) region for larger (smaller) δurec , respectively. Irregularities in
the contour shape are an artifact due to the numerical discretization.

corrections due to dissipative effects introduced by inhibitory

neurons, but these do not affect the scaling of the average size
and the avalanche shape, except for the emergence of a slight
asymmetry. Interestingly, since the system is off criticality, the
scaling relation, Eq. (1), in this case no longer provides the
value of the exponent describing how the average size scales
with duration and avalanche profiles collapse.
Moreover, the avalanche scaling theory predicts the same
exponent for the power spectrum, S( f ) ∼ f −γ [25], which is
in stark contrast with experimental data for healthy brains,
FIG. 6. (a) Power spectral density for systems at criticality pro- which rather exhibit effectively 1/ f noise. We evidence that
vides an exponent β  2.0 for all fractions of inhibitory neurons. inhibition is responsible for the crossover in the scaling be-
(b) Power spectral density for different percentages of inhibitory neu- havior of the power spectrum: Without tuning the parameter
rons and fixed δurec corresponding to the critical state for Pin = 0%. setting the system at criticality, a fraction of 30% inhibitory
The effective exponent β decreases toward unity and the scaling neurons leads to a behavior closer to 1/ f noise, the same
regime shrinks for increasing Pin → 30%. Inset: The effective ex- evidence found in self-organized models in the absence of
ponent β as function of Pin . a tuning parameter [41]. Not compensating inhibition by in-
creasing the parameter δurec amounts to driving the system
the scaling of the average size with duration, the universal away from criticality, in the subcritical regime, as the per-
profile, and the activity spectra. All these properties, except centage of inhibitory neurons increases. In this case, effective
the last one, have been intensively investigated experimentally exponents close to 1/ f noise can be measured. Therefore,
and numerically in different contexts, confirming the general numerical results suggest that deviations of β from the Brown
validity of the scaling theory. Here we performed a numerical noise value can provide an indicator of how far from criticality
study of a neuronal network model, able to simulate systems the system operates.
at and off criticality and therefore to monitor the role of inhi-
bition. We evidence that the expected scaling behavior for the
ACKNOWLEDGMENTS
average size with duration and the collapse onto a universal
profile are robust properties verified in systems with different L.d.A. would like to thank MUR (Italian Ministry of
fractions of inhibitory neurons at and near criticality, provided University and Research) Project No. PRIN2017WZFTZP
that we consider avalanches in the scaling regime. Therefore, for support. A.S. acknowledges support from MUR (Ital-
the arguments of the Sethna scaling appear to be very robust ian Ministry of University and Research) Project No.
also for neuronal avalanches, confirming the existence of a PRIN201798CZLJ. L.d.A. and A.S. thank the Program VAn-
universal exponent γ ∼ 2 even when the system is not tuned at viteLli pEr la RicErca: VALERE 2019 for financial support.
criticality. In this case the exponents for the size and duration H.J.H. thanks FUNCAP and the University of Campania for
distributions show effective values because of the exponential support through a visiting professorship.

024304-7
MANOJ KUMAR NANDI et al. PHYSICAL REVIEW E 106, 024304 (2022)

[1] W. J. Freeman, M. D. Holmes, G. A. West, and S. Vanhatalo, [31] N. Friedman, S. Ito, B. A. W. Brinkman, M. Shimono, R. E.
Int. J. Intell. Syst. 21, 881 (2006). DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler, Phys.
[2] W. J. Freeman and M. D. Holmes, Neural Networks 18, 497 Rev. Lett. 108, 208102 (2012).
(2005). [32] S. R. Miller, S. Yu, and D. Plenz, Sci. Rep. 9, 16403 (2019).
[3] G. Deco and V. K. Jirsa, J. Neurosci. 32, 3366 (2012). [33] A. Ponce-Alvarez, A. Jouary, M. Privat, G. Deco, and G.
[4] P. L. Nunez, B. M. Wingeier, and R. B. Silberstein, Hum. Brain Sumbre, Neuron 100, 1446 (2018).
Mapp. 13, 125 (2001). [34] A. de Candia, A. Sarracino, I. Apicella, and L. de Arcangelis,
[5] C. G. Bénar, C. Grova, V. K. Jirsa, and J.-M. Lina, J. Comput. PLoS Comput. Biol. 17, e1008884 (2021).
Neurosci. 47, 31 (2019). [35] E. Novikov, A. Novikov, D. Shannahoff-Khalsa, B. Schwartz,
[6] G. Deco, V. K. Jirsa, and A. R. McIntosh, Nat. Rev. Neurosci. and J. Wright, Phys. Rev. E 56, R2387 (1997).
12, 43 (2011). [36] C. Bedard, H. Kröger, and A. Destexhe, Phys. Rev. Lett. 97,
[7] P. Robinson et al., Neuropsychopharmacology 28, S74 (2003). 118102 (2006).
[8] K. Linkenkaer-Hansen, V. V. Nikouline, J. M. Palva, and R. J. [37] N. Dehghani, C. Bedard, S. Cash, E. Halgren, and A. Destexhe,
Ilmoniemi, J. Neurosci. 21, 1370 (2001). J. Comput. Neurosci. 29, 405 (2010).
[9] J. M. Beggs and D. Plenz, J. Neurosci. 23, 11167 (2003). [38] W. Pritchard, Int. J. Neurosci. 66, 119 (1992).
[10] A. Mazzoni et al., PloS One 2, e439 (2007). [39] E. Zarahn, G. Aguirre, and M. D. Esposito, NeuroImage 5, 179
[11] V. Pasquale, P. Massobrio, L. Bologna, M. Chiappalone, and S. (1997).
Martinoia, Neuroscience 153, 1354 (2008). [40] B. J. He, J. M. Zempel, A. Z. Snyder, and M. E. Raichle, Neuron
[12] E. D. Gireesh and D. Plenz, Proc. Natl. Acad. Sci. USA 105, 66, 353 (2010).
7576 (2008). [41] F. Lombardi, H. J. Herrmann, and L. de Arcangelis, Chaos 27,
[13] T. Petermann et al., Proc. Natl. Acad. Sci. USA 106, 15921 047402 (2017).
(2009). [42] J. Touboul and A. Destexhe, Phys. Rev. E 95, 012413 (2017).
[14] O. Shriki et al., J. Neurosci. 33, 7079 (2013). [43] A. Destexhe and J. Touboul, eNeuro 8, 0551 (2021).
[15] L. de Arcangelis, C. Perrone-Capano, and H. J. Herrmann, [44] D. J. Schwab, I. Nemenman, and P. Mehta, Phys. Rev. Lett. 113,
Phys. Rev. Lett. 96, 028107 (2006). 068102 (2014).
[16] A. Levina, J. M. Herrmann, and T. Geisel, Nat. Phys. 3, 857 [45] V. M. Eguíluz, D. R. Chialvo, G. A. Cecchi, M. Baliki, and A. V.
(2007). Apkarian, Phys. Rev. Lett. 94, 018102 (2005).
[17] L. Michiels van Kessenich, M. Lukovic, L. de Arcangelis, and [46] B. Roerig and B. Chen, Cereb. Cortex 12, 187 (2002).
H. J. Herrmann, Phys. Rev. E 97, 032312 (2018). [47] P. Bonifazi et al., Science 326, 1419 (2009).
[18] S. Zapperi, K. B. Lauritsen, and H. E. Stanley, Phys. Rev. Lett. [48] D. Raimo, A. Sarracino, and L. de Arcangelis, Phys. A
75, 4071 (1995). (Amsterdam, Neth.) 565, 125555 (2021).
[19] L. Michiels van Kessenich, L. de Arcangelis, and H. Herrmann, [49] C. Rosenmund and C. F. Stevens, Neuron 16, 1197 (1996).
Sci. Rep. 6, 32071 (2016). [50] P. S. Kaeser and W. G. Regehr, Curr. Opin. Neurobiol. 43, 63
[20] A. Fontenele et al., Phys. Rev. Lett. 122, 208101 (2019). (2017).
[21] B. Gutenberg and C. Richter, Seismicity of the Earth and Asso- [51] S. O. Rizzoli and W. J. Betz, Nat. Rev. Neurosci. 6, 57 (2005).
ciated Phenomena (Princeton University, Princeton, NJ, 1954). [52] K. Ikeda and J. M. Bekkers, Proc. Natl. Acad. Sci. USA 106,
[22] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 2945 (2009).
(1987). [53] H. Markram and M. Tsodyks, Nature (London) 382, 807
[23] H. Jensen, Self-Organized Criticality (Cambridge University, (1996).
Cambridge, England, 1998). [54] T. Abrahamsson, B. Gustafsson, and E. Hanse, J. Physiol. 569,
[24] J. J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature 737 (2005).
(London) 410, 242 (2001). [55] M. Armbruster and T. A. Ryan, Nat. Neurosci. 14, 824 (2011).
[25] M. C. Kuntz and J. P. Sethna, Phys. Rev. B 62, 11699 (2000). [56] F. Lombardi, H. J. Herrmann, C. Perrone-Capano, D. Plenz, and
[26] L. Laurson and M. J. Alava, Phys. Rev. E 74, 066106 (2006). L. de Arcangelis, Phys. Rev. Lett. 108, 228703 (2012).
[27] A. P. Mehta, K. A. Dahmen, and Y. Ben-Zion, Phys. Rev. E 73, [57] D. O. Hebb, The Organization of Behavior: A Neuropsycholog-
056104 (2006). ical Approach (Wiley & Sons, New York, 1949).
[28] S. Papanikolaou et al., Nat. Phys. 7, 316 (2011). [58] P. Fransson et al., Cereb. Cortex 23, 638 (2013).
[29] L. Laurson et al., Nat. Commun. 4, 2927 (2013). [59] S. Scarpetta, I. Apicella, L. Minati, and A. de Candia, Phys.
[30] G. Durin, F. Colaiori, C. Castellano, and S. Zapperi, J. Magn. Rev. E 97, 062305 (2018).
Magn. Mater. 316, 436 (2007). [60] J. B. Mallinson et al., Sci. Adv. 5, eaaw8438 (2019).

024304-8