# On the critical signatures of neural activity

Benedetta Mariani  Department of Physics and Astronomy “Galileo Galilei”,
University of Padova, Padova, Italy  Padova Neuroscience Center, University of
Padova, Padova, Italy  Giorgio Nicoletti  Department of Physics and Astronomy
“Galileo Galilei”, University of Padova, Padova, Italy  Marta Bisio  Padova
Neuroscience Center, University of Padova, Padova, Italy  Department of
Biomedical Science, University of Padova, Padova, Italy  Marta Maschietto
Department of Biomedical Science, University of Padova, Padova, Italy  Stefano
Vassanelli  Padova Neuroscience Center, University of Padova, Padova, Italy
Department of Biomedical Science, University of Padova, Padova, Italy  Samir
Suweis  Department of Physics and Astronomy “Galileo Galilei”, University of
Padova, Padova, Italy  Padova Neuroscience Center, University of Padova,
Padova, Italy

###### Abstract

The critical brain hypothesis has emerged as an attractive framework to
understand neuronal activity, but it is still widely debated. We show that
power-law neuronal avalanches satisfying the crackling-noise relation emerge
in a paradigmatic model as a sole consequence of external modulation, which
induces an extrinsic mutual information in the system. Using data from multi-
electrodes array probes in the rat’s cortex, we find that the spatial
correlations display typical features of critical systems and we show that our
model, with the addition of intrinsic interactions, reproduces these critical
scale-free correlations, coexisting with and disentangled from externally-
driven power-law avalanches.

††preprint: APS/123-QED

The critical brain hypothesis has been much investigated since scale-free
neuronal avalanches were found in 2003 by Beggs and Plenz [1, 2, 3, 4, 5, 6,
7, 3, 8, 9]. In their seminal work they showed that, by analyzing Local Field
Potentials (LFPs) from cortical slices and cultures on chips, cascades of
neuronal activity are power-law distributed in their sizes and lifetimes, with
exponents remarkably close to the ones of a critical branching process. Since
then, such power-laws have been repeatedly observed in experiments [10, 11,
12, 13, 14, 15] suggesting that the brain might be poised near the critical
point of a phase transition. However, this hypothesis is still widely debated
[16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26], and reconciling the different
views and experimental results remains pressing.

On the one hand, many subsequent works showed that the presence of power-law
avalanches is not a sufficient condition for criticality, as they might emerge
from different mechanisms [18, 19, 20, 21]. What is now considered one of the
stringest test for criticality is whether avalanche exponents satisfy the
crackling-noise relation [27, 28], which is believed to hold exactly only at
criticality [29]. Moreover, recent findings [22, 30] suggest that the
avalanche exponents found in different experimental settings lay along the
scaling line defined by the crackling-noise relation with a seemingly
universal exponent $\delta\approx 1.28$. However, this relation can be
fulfilled even in models of independent spiking units, for a range of choices
of the power-law fitting method [24]. Arguably, a more fundamental signature
of criticality is the presence of long-range correlations in space [31], which
manifest themselves in a correlation length that scales with the system size.
Nevertheless, these methods have usually been applied at coarser scales, such
as in whole brain data [32, 33], and not in specific cortical areas [34]. When
such spatial information is not available, phenomenological renormalization
group procedures [35, 36] have been recently proposed.

On the other hand, the debate about the “nature” of the transition is very
much open, and thus its hypothetical universality class is poorly understood.
In particular, recent works have proposed that the observed transition might
be related to a synchronous-asynchronous one [29, 37, 38, 39] or a disorder-
induced transition [33], contrarily to the original, simpler hypothesis of a
critical branching process. All in all, these results still lack a unifying
framework. In this work, we show that signatures of criticality beyond
avalanches are present in spatial LFPs, in the form of a spatial correlation
length that grows with the system size. Then, we introduce an archetypal
framework in which the origin of these different signatures can be
disentangled, shedding a light on the underlying mechanisms from which such
properties of neural activity emerge. In particular, our work suggests that,
whereas avalanches may emerge from extrinsic stochasticity, interactions
between neural populations are the fundamental biological mechanism that gives
rise to seemingly scale-free correlations.

Figure 1: (a) Left: scheme of the array used to obtain the LFPs data from all
the cortical layers of the barrel cortex (modified from [40]); right an
example of the LFPs signals for different layers and the corresponding
discretization. An array of 256 channels organized in a $64\times 4$ matrix is
inserted in a barrel column and the signals from the cortical layers are
collected by $55\times 4$ electrodes. (b-d) Avalanche statistics obtained from
the analysis of LFPs data in a rat. Both the distribution of the avalanches
(b) sizes and (c) duration are power-laws, and (d) the crackling-noise
relation is satisfied. (e) Scaling of the correlation length with the system
size in LFPs data, averaging over four different rats. The error bars are
shown as $5$ standard deviations from the mean for visual ease. The
correlation length scales linearly with the system size with no plateau in
sight, a hallmark of criticality.

We measured the population activity in the primary somatosensory cortex of
four rats using state of the art spatially-extended multi-electrodes arrays.
The cortical activity is recorded through a 256-channels array organized in a
$64$ rows $\times$ $4$ columns matrix with an inter-electrode distance of 32
$\mu m$ [41], see Figure 1a. The avalanche statistics is analyzed in LFPs data
following standardized pipelines for event detection [13, 10, 42, 41], and we
report a full statistical analysis of our dataset in [43]. Briefly, the
distribution $p(s)$ of the avalanches sizes and of the avalanche duration
$p(T)$ are computed and fitted using a corrected maximum likelihood method
[41, 44, 43]. We find that both are statistically compatible with the expected
power-laws $P(s)\sim s^{-\tau}$ and $P(T)\sim T^{-\tau_{t}}$ as we show in
Figure 1b-c. Averaging over four rats, we find an inter-rat variability with
average exponents $\expectationvalue{\tau}=1.75\pm 0.1$ and
$\expectationvalue{\tau_{t}}=2.1\pm 0.3$ [41]. In Figure 1d we show that the
crackling-noise relation [27] holds, by comparing
$\delta_{\mathrm{pred}}=\frac{\tau_{t}-1}{\tau-1}$ with the exponent obtained
from $\expectationvalue{s(T)}\sim T^{\delta_{\mathrm{fit}}}$. Averaging over
each of our rats, we find $\expectationvalue{\delta_{\mathrm{pred}}}=1.47\pm
0.18$ and $\expectationvalue{\delta_{\mathrm{fit}}}=1.46\pm 0.04$. Further
details on the fitting procedure can be found in [43].

We then focus on the spatial correlations of the fluctuations of the measured
LFP activity and whether they display any signature of criticality [45, 32].
By studying the scaling of the correlation length $\xi$ as a function of the
system sizes $L$ [46, 41], we find that $\xi$ scales linearly with $L$,
signaling the presence of long-range correlations. This behavior matches
exactly what would happen at a critical point, where the correlation length
diverges in the thermodynamic limit and thus grows with the size of a finite
system.

Hence, we find that the measured neural activity in the barrel cortex at rest
displays two different signatures of criticality - power-law avalanches and
scale-free spatial correlations. In order to try and unfold the underlying
processes from which these collective properties emerge, we assume that neural
activity may be decomposed in two parts [47, 48]: (i) the _intrinsic activity_
, which is the activity driven by interactions between neurons or populations
of neurons - in our case, the propagation dynamics across the multi-layer
network of the interconnected neurons along the barrel; (ii) the _extrinsic
activity_ , which instead corresponds to the activity modulated by an external
or global unit - in our case, the external inputs triggering or modulating the
propagation (e.g. synaptic current injection from the thalamic inputs).

We introduce a paradigmatic model model of $N$ variables
$(v_{1},\dots,v_{N})$, denoting the activity of $N$ units (e.g. neurons or, as
in our case, distinct populations of neurons as measured by our LFPs). This
model is simple enough to be treated analytically and to provide a clear
physical interpretation, while being complex enough to display non-trivial
behaviors [49]. Focusing first on the extrinsic activity only, we model these
variables as conditionally independent through an Ornstein-Uhlenbeck (OU)
process [50, 51, 52, 53] with time-dependent diffusion

$\frac{dv_{i}(t)}{dt}=-\frac{1}{\gamma_{i}}v_{i}(t)+\sqrt{\mathcal{D}(t)}\eta_{i}(t),$
(1)

where $\gamma_{i}$ is the characteristic time of the $i$-th neuron,
$\eta_{i}(t)$ is a white noise, and $\mathcal{D}(t)$ corresponds to a noise
strength modulation from an external input that is shared among all the units.
In particular, we define the noise modulation $D(t)$ as

$\mathcal{D}(t)=\begin{cases}\mathcal{D}^{*}&\text{if}\quad
D(t)\leq\mathcal{D}^{*}\\\ D(t)&\text{if}\quad
D(t)>\mathcal{D}^{*}\end{cases}$ (2)

where, following [19], $D(t)$ is itself an OU process
$\dot{D}(t)=-D(t)/\gamma_{D}+\sqrt{\theta}\eta_{D}(t)$ [54]. Therefore, the
noise modulation $\mathcal{D}$ can be either constant in time and equal to
$\mathcal{D}^{*}$ or change in time according to an OU process with value
$\mathcal{D}(t)>\mathcal{D^{*}}$.

Figure 2: Avalanche statistics generated by the model at $\mathcal{D}^{*}=0.3$
(a-b, e-g) and at $\mathcal{D}^{*}=5$ (c-d,h-j), with $\gamma_{D}=15$ and
$\theta=1$ and $\gamma_{i}=\gamma=0.05$ for the extrinsic model. (a-b)
Comparison between the trajectories of $\mathcal{D}(t)$, $v_{i}$ and the
corresponding discretization in the low-$\mathcal{D}^{*}$ regime for (a) the
extrinsic model and (b) the interacting one. (c-d) Same, but in the
high-$\mathcal{D}^{*}$ regime. (e-g) If $\mathcal{D}^{*}$ is low, avalanches
are power-law distributed with almost identical exponents in the extrinsic and
interacting model, $\tau^{\mathrm{ext}}=1.60\pm 0.01$,
$\tau^{\mathrm{int}}=1.55\pm 0.01$ and $\tau^{\mathrm{ext}}_{t}=1.77\pm 0.01$,
$\tau^{\mathrm{int}}_{t}=1.74\pm 0.01$. The crackling-noise relation is
verified in both cases. (h-j) Same plots, now in the high-$\mathcal{D}^{*}$
regime. Avalanches are now fitted with an exponential distribution. Notice
that larger events, corresponding to periods in which
$\mathcal{D}(t)>\mathcal{D}^{*}$, show up in the distributions’ tails,
suggesting that the shift between exponentials and power-laws is smooth. (j)
The average avalanche size as a function of the duration scales with an
exponent that, as $D^{*}$ increases, becomes closer to the trivial one
$\delta_{\mathrm{fit}}^{\mathrm{ext}}\approx\delta_{\mathrm{fit}}^{\mathrm{int}}\approx
1$.

We now assume that the timescales of the two processes are separated [47],
that is $\gamma_{D}\gg\gamma_{i}$, so that the stationary probability
distribution of the activity is $p(v)=\int
d\mathcal{D}\prod_{i=1}^{N}p(v_{i}|\mathcal{D})p(\mathcal{D})$ and in general
$p(v_{i},v_{j})\neq p(v_{i})p(v_{j})$. Thus, the presence of the unobserved
modulation results in the emergence of an effective dependence between the
units. Remarkably, one can show that the units, although not independent, are
always uncorrelated
$\expectationvalue{v_{i}v_{j}}-\expectationvalue{v_{i}}\expectationvalue{v_{j}}=0\quad\forall
i\neq j$ [41]. At the same time, the dependence between the units is shaped by
the parameter $\mathcal{D}^{*}$. If $\mathcal{D}^{*}$ is high enough, the
modulation is rare and the units are always dominated by noise as we can see
in Figure 2c. On the other hand, if $\mathcal{D}^{*}$ is small, whenever
$\mathcal{D}(t)=\mathcal{D}^{*}$ the noise contribution to the units will be
vanishing and the activity will follow an exponential decay. Therefore, in
this regime, each $v_{i}$ will typically alternate periods of quasi-silence to
periods of activity. In other words, depending on the value of
$\mathcal{D}^{*}$, this model can either reproduce a noise-driven behavior or
a bursty, coordinated one, as shown in Figure 2a.

Most importantly, the low $\mathcal{D}^{*}$ regime is also the onset of power-
law distributed avalanches, as we see by simulating the model at different
$\mathcal{D}^{*}$ and performing the same analysis as in the LFP. We find
that, as $\mathcal{D}^{*}$ decreases, a transition between exponential
decaying avalanches and power-law distributed ones appears. Figures 2e-f-g
show that for $\mathcal{D}^{*}$ small enough, the stochastic modulation
produces scale-free avalanches in both size and time with exponents
$\tau^{\mathrm{ext}}=1.60\pm 0.01$ and $\tau^{\mathrm{ext}}_{t}=1.77\pm 0.01$
satisfying the crackling-noise relation, i.e.
$\delta_{\mathrm{fit}}^{\mathrm{ext}}=1.21\pm
0.01\approx\delta_{\mathrm{pred}}^{\mathrm{ext}}=1.28\pm 0.02$. Moreover, the
rescaled temporal profiles of these avalanches collapse to a single curve, as
we show in the Supplemental Material [41]. On the other hand, for higher
$\mathcal{D}^{*}$ only exponential-decaying avalanches are present as we see
in Figures 2h-i-j.

Although this is an archetypal example, our results are qualitatively
unchanged if we use a more biological plausible model, at the price of losing
analytical tractability. In particular, as shown in the Supplemental Material
[41], we model the firing rates of both the external modulation and of the
neural dynamics through non-interacting Wilson-Cowan units [55]. We set the
external dynamics in a balanced state [56] and we still find that the
avalanches of this extrinsic Wilson-Cowan model do display power law
avalanches and the crackling noise relation is satisfied as well.

We have thus shown that avalanches statistics can be generated by non-critical
mechanisms such as the one presented above. However the correlations as
observed in Figure 1c cannot be explained be the extrinsic model. We thus
introduce an intrinsic component in the OU model, i.e., we add a reconstructed
effective connectivity [53] to Eqs. (S3-2). Namely, we now consider

${\frac{dv_{i}(t)}{dt}=-\sum_{j}A_{ij}v_{j}(t)+\sqrt{\mathcal{D}(t)}\eta_{i}(t)}$
(3)

where $A_{ij}$ is the effective connectivity matrix describing the interaction
between the $i$-th and the $j$-th unit. Since in the extrinsic model the units
are uncorrelated, we show in the Supplemental Material [41] that we can infer
the values of $A_{ij}$ directly from the data, by solving the inverse problem
in such a way that the correlations of Eq. (3) match the experimentally-
measured correlations $\sigma_{ij}$ of our LFPs. The different regimes for the
interacting model are plotted in Figures 2b and 2d.

As shown in Figures 2e-f-g, all avalanches exponents $\tau\approx 1.6$,
$\tau_{t}\approx 1.75$ and $\delta\approx 1.28$ are not changed significantly
by the inclusion of direct interactions among the units, nor the high
$\mathcal{D}^{*}$ regime is changed either as shown in Figures 2h-i-j. Let us
note that these exponents are different from the ones obtained in LFP data,
but this is perhaps not surprising. In fact, beside the simplicity of this
paradigmatic model, such exponents have been found to depend on the
experimental settings [30] and on individual variability [22]. Nevertheless,
our framework reproduces the ubiquitous scaling exponent $\delta\approx 1.28$
[22, 39], which we also find in our experimental settings in MUAs data [43].

The fact that the exponents do not change when we add the interaction network
to our model has a profound implication: if we measure neural avalanches
alone, we might not be able to infer anything on the intrinsic neural dynamics
of our model. However, we might have access to such dynamics if we investigate
the correlations of the fluctuations. This disentangling is deeply related to
the dependency structure of our model, described by the mutual information

$I=\int_{-\infty}^{+\infty}dv_{i}\int_{-\infty}^{+\infty}dv_{j}p(v_{i},v_{j})\log\frac{p(v_{i},v_{j})}{p(v_{i})p(v_{j})},$
(4)

which captures pairwise dependencies in the system that go well beyond simple
correlations. In Figure 3a we show that a non-zero mutual information emerges
in the extrinsic model in the low $\mathcal{D}^{*}$ limit. The onset of the
coordinated behavior between the units, i.e., of power-law distributed
avalanches, is also the onset of a non-zero mutual information. When
interactions are added back, the mutual information is simply shifted
independently of $\mathcal{D}^{*}$ [41], signaling that the effects of
external modulation and the effective interactions are completely disentangled
[57].

Furthermore, with the effective interactions reconstructed from the data we
can also study the scaling of the correlation length $\xi$ as a function of
the system sizes $L$ and we plot the result in Figure 3b. As in our data,
$\xi$ scales linearly with $L$ in the interacting model, displaying thus a
signature of criticality that not only is strictly dependent on the
interaction network, but that is also completely unaffected by avalanches. In
fact, as with the mutual information, the correlations in the interacting
model scale trivially with $\mathcal{D}^{*}$ and thus the correlation length
is not affected by $\mathcal{D}^{*}$ (see the Supplemental Material [41]).

Figure 3: (a) Comparison between the mutual information in the extrinsic model
($\theta=1$, $\gamma_{D}=10$, $\gamma_{1}=0.1$, $\gamma_{2}=0.5$) and, as an
example, in the interacting model with two units (see the Supplemental
Material [41]). Notice that the onset of a non-vanishing mutual information
induced by $\mathcal{D}(t)$ is also the onset of power-law distributed
avalanches, whereas the mutual information arising from interactions is
independent on $\mathcal{D}^{*}$. (b) The correlation length of the
interacting model scales linearly with the system size, as in the data. In the
extrinsic model, as expected, the correlation length of the fluctuations is
constant and equal to 1, i.e., the correlation function drops to zero for
adjacent electrodes.

Overall, our work shows that while power-law avalanches and the crackling-
noise relation both appear in LFPs data from the rat’s barrel cortex, their
relation with criticality is far from stringent. We have developed a
paradigmatic, but of analytical ease, framework where the intrinsic
contributions to the neuronal activity, due to the direct interaction between
the units themselves, and the extrinsic ones, arising from externally-driven
modulated activity, are exactly disentangled. Our model displays a regime in
which power-law avalanches that satisfy the crackling-noise relation emerge
and that are compatible with the universal exponent $\delta\approx 1.28$ that
appears in [22], a general result that holds also when considering a different
extrinsic model, such as Wilson-Cowan units. The fact that this exponent turns
out to be so ubiquitous in a variety of neural systems, from freely moving or
anesthetized mammals to cultured slices, suggests that it could be explained
by some basic phenomena shared by all these systems. We propose that such
phenomena can be represented by a simple time varying extrinsic dynamics that
does not require any fine tuning to criticality [24].

At the same time, we show that in spatially-extended LFPs data seemingly
scale-free spatial correlations are present. While they can coexist with
phenomena like power-law avalanches, these kind of signatures are not
explained by the extrinsic activity alone and thus may be deeply related to
the origin of criticality in the brain, playing a fundamental role in the
advantages it might achieve by being critical [9, 58, 3, 59]. Crucially, our
archetypal model allows to combine this extrinsic dynamics to an intrinsic
interaction matrix, inferred directly from the experimental data to match the
seemingly scale-free spatial correlations we find in our experiments. When we
do so, we show that these two signatures of criticality are completely
disentangled - avalanches appear as a sole consequence of the external
modulation and are not affected by the interactions, and, vice-versa, the
interactions determine the spatial correlations independently on the external
modulation.

Remarkably, it was recently shown that in these kind of models the mutual
information always receives distinct and disentangled contributions from the
internal interactions and from an environmental, external dynamics [57].
Although, clearly, the presence of a non-zero mutual information cannot be a
sufficient condition for power-law avalanches to appear, in our extrinsic
model their emergence does correspond to the onset of a non-vanishing mutual
information. This fact suggests a promising future perspective. By explicitly
considering both the intrinsic activity and the extrinsic contributions, one
might be able to combine all these considerations into a unified information-
theoretic view - and perhaps helping to unfold the underlying biological
mechanisms at the origin of the observed signatures of criticality in neural
activity.

###### Acknowledgements.

S.S. acknowledges DFA and UNIPD for SUWE_BIRD2020_01 grant, and INFN for
LINCOLN grant. S.V acknowledges support from the European Commission, FET
Proactive, SYNCH project (GA number 824162).

## References

  * [1] J. M. Beggs and D. Plenz., “Neuronal avalanches in neocortical circuits,” Journal of Neuroscience, vol. 12, p. 23(35):11167–11177, 2003.
  * [2] L. de Arcangelis, C. Perrone-Capano, and H. J. Herrmann, “Self-organized criticality model for brain plasticity,” Physical Review Letters, vol. 96, no. 2, p. 028107, 2006.
  * [3] O. Kinouchi and M. Copelli, “Optimal dynamical range of excitable networks at criticality,” Nature physics, vol. 2, no. 5, pp. 348–351, 2006.
  * [4] W. L. Shew and D. Plenz, “The functional benefits of criticality in the cortex.,” The Neuroscientist, vol. 19,1, pp. 88–100, 2013.
  * [5] J. Hesse and T. Gross, “Self-organized criticality as a fundamental property of neural systems,” Frontiers in systems neuroscience, vol. 8, p. 166, 2014\.
  * [6] G. Tkačik, T. Mora, O. Marre, D. Amodei, S. E. Palmer, M. J. Berry, and W. Bialek, “Thermodynamics and signatures of criticality in a network of neurons,” Proceedings of the National Academy of Sciences, vol. 112, no. 37, pp. 11508–11513, 2015.
  * [7] R. P. Rocha, L. Koçillari, S. Suweis, M. Corbetta, and A. Maritan, “Homeostatic plasticity and emergence of functional networks in a whole-brain model at criticality,” Scientific reports, vol. 8, no. 1, pp. 1–15, 2018.
  * [8] J. Hidalgo, J. Grilli, S. Suweis, M. A. Munoz, J. R. Banavar, and A. Maritan, “Information-based fitness and the emergence of criticality in living systems,” Proceedings of the National Academy of Sciences, vol. 111, no. 28, pp. 10095–10100, 2014.
  * [9] M. A. Muñoz, “Colloquium: Criticality and dynamical scaling in living systems,” Review of Modern Physics, vol. 90, p. 031001, 2018.
  * [10] T. Petermann, T. C. Thiagarajan, M. A. Lebedev, M. A. L. Nicolelis, D. R. Chialvo, and D. Plenz, “Spontaneous cortical activity in awake monkeys composed of neuronal avalanches,” Proceedings of the National Academy of Sciences, vol. 106, no. 37, pp. 15921–15926, 2009.
  * [11] S. Yu, H. Yang, H. Nakahara, G. S. Santos, D. Nikolić, and D. Plenz, “Higher-order interactions characterized in cortical activity,” Journal of Neuroscience, vol. 31, no. 48, pp. 17514–17526, 2011.
  * [12] G. Hahn, T. Petermann, M. Havenith, and al., “Neuronal avalanches in spontaneous activity in vivo,” Journal of Neurophysiology., vol. 104(6), p. 3312‐3322, 2010.
  * [13] E. D. Gireesh and D. Plenz, “Neuronal avalanches organize as nested theta- and beta/gamma-oscillations during development of cortical layer 2/3,” Proceedings of the National Academy of Sciences, p. 7576‐7581, 2008.
  * [14] A. Mazzoni, F. D. Broccard, E. Garcia-Perez, P. Bonifazi, M. Ruaro, and V. Torre, “On the dynamics of the spontaneous activity in neuronal networks,” PLoS One, vol. 2(5):e439, 2007.
  * [15] V. Pasquale, P. Massobrio, L. Bologna, M. Chiappalone, and S. Martinoia, “Self-organization and neuronal avalanches in networks of dissociated cortical neurons,” Neuroscience, vol. 153(4), p. 1354‐1369, 2008.
  * [16] C. Bedard, H. Kroeger, and A. Destexhe, “Does the 1/f frequency scaling of brain signals reflect self-organized critical states?,” Physical Review Letters, vol. 97, no. 11, p. 118102, 2006.
  * [17] J. M. Beggs and N. Timme, “Being critical of criticality in the brain,” Frontiers in physiology, vol. 3, p. 163, 2012.
  * [18] J. Touboul and A. Destexhe, “Can power-law scaling and neuronal avalanches arise from stochastic dynamics?,” PLOS ONE, vol. 5, pp. 1–14, 02 2010\.
  * [19] J. Touboul and A. Destexhe, “Power-law statistics and universal scaling in the absence of criticality,” Physical Review E, vol. 95, p. 012413, 2017.
  * [20] M. Martinello, J. Hidalgo, A. Maritan, S. di Santo, D. Plenz, and M. A. Muñoz, “Neutral theory and scale-free neural dynamics,” Physical Review X, vol. 7, p. 041071, 2017.
  * [21] A. Faqeeh, S. Osat, F. Radicchi, and J. P. Gleeson, “Emergence of power laws in noncritical neuronal systems,” Physical Review E, vol. 100, no. 1, p. 010401, 2019.
  * [22] A. J. Fontenele, N. A. P. de Vasconcelos, T. Feliciano, L. A. A. Aguiar, C. Soares-Cunha, B. Coimbra, L. Dalla Porta, S. Ribeiro, A. J. a. Rodrigues, N. Sousa, P. V. Carelli, and M. Copelli, “Criticality between cortical states,” Physical Review Letters, vol. 122, p. 208101, 2019.
  * [23] J. Wilting and V. Priesemann, “25 years of criticality in neuroscience — established results, open controversies, novel concepts,” Current Opinion in Neurobiology, vol. 58, pp. 105 – 111, 2019.
  * [24] A. Destexhe and J. D. Touboul, “Is there sufficient evidence for criticality in cortical systems?,” eNeuro, vol. 8, no. 2, 2021.
  * [25] J. Ruseckas and B. Kaulakys Chaos, vol. 23, p. 023102, 2013.
  * [26] B. Kaulakys and J. Ruseckas, “Stochastic nonlinear differential equation generating 1/f noise,” Phys. Rev. E, vol. 70, no. 2, p. 020101, 2004.
  * [27] J. Sethna, K. Dahmen., and C. Myers, “Crackling noise,” Nature, vol. 410, p. 242–250, 2001.
  * [28] N. Friedman, S. Ito, B. A. W. Brinkman, M. Shimono, R. E. L. DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler, “Universal critical dynamics in high resolution neuronal avalanche data,” Physical Review Letters, vol. 108, p. 208102, 2012.
  * [29] S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz, “Simple unified view of branching process statistics: Random walks in balanced logarithmic potentials,” Physical Review E, vol. 95, p. 032115, 2017.
  * [30] L. J. Fosque, R. V. Williams-García, J. M. Beggs, and G. Ortiz, “Evidence for Quasicritical Brain Dynamics,” Physical Review Letters, vol. 126, no. 9, p. 098101, 2021.
  * [31] J. J. Binney, N. J. Dowrick, A. J. Fisher, and M. Newman, The Theory of Critical Phenomena: An Introduction to the Renormalization Group.  USA: Oxford University Press, Inc., 1992.
  * [32] A. Haimovici, E. Tagliazucchi, P. Balenzuela, and D. R. Chialvo, “Brain organization into resting state networks emerges at criticality on a model of the human connectome,” Physical Review Letters, vol. 110, p. 178101, 2013\.
  * [33] A. Ponce-Alvarez, A. Jouary, M. Privat, G. Deco, and G. Sumbre, “Whole-brain neuronal activity displays crackling noise dynamics,” Neuron, vol. 100, no. 6, pp. 1446 – 1459.e6, 2018.
  * [34] T. L. Ribeiro, S. Yu, D. A. Martin, D. Winkowski, P. Kanold, D. R. Chialvo, and D. Plenz, “Trial-by-trial variability in cortical responses exhibits scaling in spatial correlations predicted from critical dynamics,” bioRxiv, 2020\.
  * [35] L. Meshulam, J. L. Gauthier, C. D. Brody, D. W. Tank, and W. Bialek, “Coarse graining, fixed points, and scaling in a large population of neurons,” Physical Review Letters, vol. 123, no. 17, p. 178103, 2019.
  * [36] G. Nicoletti, S. Suweis, and A. Maritan, “Scaling and criticality in a phenomenological renormalization group,” Physical Review Research, vol. 2, p. 023144, 2020.
  * [37] L. Dalla Porta and M. Copelli, “Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic M/EEG results,” PLOS Computational Biology, vol. 15, no. 4, pp. 1–26, 2019.
  * [38] S. S. Poil, R. Hardstone, H. D. Mansvelder, and K. Linkenkaer-Hansen, “Critical-state dynamics of avalanches and oscillations jointly emerge from balanced excitation/inhibition in neuronal networks,” Journal of Neuroscience, vol. 32, no. 29, pp. 9817–9823, 2012.
  * [39] V. Buendia, P. Villegas, R. Burioni, and M. A. M. noz, “Hybrid-type synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together,” Phys. Rev. Research, vol. 3, no. 2, p. 023224, 2021.
  * [40] Z.-W. Zhang and M. Deschênes, “Intracortical axonal projections of lamina vi cells of the primary somatosensory cortex in the rat: A single-cell labeling study,” Journal of Neuroscience, vol. 17, no. 16, pp. 6365–6379, 1997.
  * [41] See Supplemental Material for experimental procedures, details on the statistical analyses and on the numerical simulations, analytical derivations and the study of more realistic models.
  * [42] O. Shriki, J. Alstott, F. Carver, T. Holroyd, R. N. Henson, M. L. Smith, R. Coppola, E. Bullmore, and D. Plenz, “Neuronal avalanches in the resting meg of the human brain,” Journal of Neuroscience, vol. 33, no. 16, pp. 7079–7090, 2013.
  * [43] B. Mariani, G. Nicoletti, M. Bisio, M. Maschietto, R. Oboe, A. Leparulo, S. Suweis, and S. Vassanelli, “Neuronal avalanches across the rat somatosensory barrel cortex and the effect of single whisker stimulation,” Frontiers in Systems Neuroscience, vol. 15, p. 89, 2021.
  * [44] M. Gerlach and E. G. Altmann, “Testing statistical laws in complex systems,” Physical Review Letters, vol. 122, p. 168301, 2019.
  * [45] D. Fraiman and D. R. Chialvo, “What kind of noise is brain noise: anomalous scaling behavior of the resting brain activity fluctuations,” Frontiers in Physiology, vol. 3, p. 307, 2012.
  * [46] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale, “Scale-free correlations in starling flocks,” Proceedings of the National Academy of Sciences, vol. 107, no. 26, pp. 11865–11870, 2010.
  * [47] V. Priesemann and O. Shriki, “Can a time varying external drive give rise to apparent criticality in neural systems?,” PLoS computational biology, vol. 14, no. 5, p. e1006081, 2018.
  * [48] U. Ferrari, S. Deny, M. Chalk, G. Tkačik, O. Marre, and T. Mora, “Separating intrinsic interactions from extrinsic correlations in a network of sensory neurons,” Physical Review E, vol. 98, no. 4, p. 042410, 2018\.
  * [49] E. Nozari, J. Stiso, L. Caciagli, E. J. Cornblath, X. He, M. A. Bertolero, A. S. Mahadevan, G. J. Pappas, and D. S. Bassett, “Is the brain macroscopically linear? a system identification of resting state dynamics,” arXiv preprint arXiv:2012.12351, 2020.
  * [50] M. L. Saggio, P. Ritter, and V. K. Jirsa, “Analytical operations relate structural and functional connectivity in the brain,” PLOS ONE, vol. 11, pp. 1–25, 08 2016.
  * [51] L. Arbabyazd, K. Shen, Z. Wang, M. Hofmann-Apitius, P. Ritter, A. R. McIntosh, D. Battaglia, and V. Jirsa, “Virtual connectomic datasets in alzheimer’s disease and aging using whole-brain network dynamics modelling,” eNeuro, vol. 8, no. 4, 2021.
  * [52] M. Gilson, N. E. Kouvaris, G. Deco, J.-F. Mangin, C. Poupon, S. Lefranc, D. Rivière, and G. Zamora-López, “Network analysis of whole-brain fmri dynamics: A new framework based on dynamic communicability,” NeuroImage, vol. 201, p. 116007, 2019.
  * [53] M. Gilson, R. Moreno-Bote, A. Ponce-Alvarez, P. Ritter, and G. Deco, “Estimation of directed effective connectivity from fmri functional connectivity hints at asymmetries of cortical connectome,” PLoS Comput Biol, vol. 12(3), 2016.
  * [54] One can set $\theta=1$ without loss of generality, hence we often do so to reduce the number of free parameters of the model.
  * [55] M. Benayoun, J. D. Cowan, W. van Drongelen, and E. Wallace, “Avalanches in a stochastic model of spiking neurons,” PLOS Computational Biology, vol. 6, no. 7, pp. 1–13, 2010.
  * [56] S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz, “Non-normality, reactivity, and intrinsic stochasticity in neural dynamics: a non-equilibrium potential approach,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2018, p. 073402, jul 2018.
  * [57] G. Nicoletti and D. M. Busiello, “Mutual information disentangles interactions from changing environments,” Physical Review Letters, 2021 (accepted).
  * [58] W. Bialek and T. Mora, “Are biological systems poised at criticality?,” Journal of Statistical Physics, vol. 144, no. 2, 2011.
  * [59] D. R. Chialvo, “Emergent complex neural dynamics,” Nature Physics, vol. 6, no. 744, 2010.
  * [60] P. Kloeden and E. Platen, Numerical Solution of Differential Stochastic Equations.  Springer Berlin, 1999.
  * [61] L. Swanson, Brain maps: structure of the rat brain.  London, Academic press, third ed., 2003.
  * [62] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences, vol. 13 of Springer Series in Synergetics.  Springer-Verlag, third ed., 2004.
  * [63] Let us note that the process is simulated through its discrete version by employing the Euler–Maruyama method [60], and for this purpose a finite $dt$ has to be chosen. The deterministic part of the discrete version of the process is a discrete map of the form $\bm{v}(t+dt)=K\bm{v}(t)$, where $K=\mathbb{I}+(W-\frac{\mathbb{I}}{\gamma})dt$. When choosing $dt$ one has to ensure that $K$ has spectral radius less than unity, in order to guarantee stability.
  * [64] R. H. Bartels and G. W. Stewart, “Solution of the matrix equation AX + XB = C,” Commun. ACM, vol. 15, no. 9, p. 820–826, 1972.
  * [65] D. A. Martin, T. L. Ribeiro, S. A. Cannas, T. S. Grigera, D. Plenz, and D. R. Chialvo, “Box scaling as a proxy of finite size correlations,” Scientific Reports, vol. 11, no. 1, pp. 1–9, 2021.
  * [66] W. Shew, W. Clawson, J. Pobst, and al, “Adaptation to sensory input tunes visual cortex to criticality.,” Nature Phys., vol. 11, p. 659–663, 2015\.
  * [67] G. Buzsáki, C. Anastassiou, and C. Koch, “The origin of extracellular fields and currents — EEG, ECoG, LFP and spikes,” Nature Review Neuroscience, vol. 13, pp. 407–420, 2012.
  * [68] A. Clauset, C. R. Shalizi, and M. E. J. Newman, “Power-law distributions in empirical data,” SIAM Rev., p. 661–703, 2009.
  * [69] K. Bansal, J. O. Garcia, N. Lauharatanahirun, S. F. Muldoon, P. Sajda, and J. M. Vettel, “Scale-specific dynamics of large-amplitude bursts in eeg capture behaviorally meaningful variability,” 2020.
  * [70] N. Marshall, N. M. Timme, N. Bennett, M. Ripp, E. Lautzenhiser, and J. M. Beggs, “Analysis of power laws, shape collapses, and neural complexity: New techniques and matlab support via the ncc toolbox,” Frontiers in Physiology, vol. 7, p. 250, 2016.
  * [71] S. Scarpetta, I. Apicella, L. Minati, and A. de Candia, “Hysteresis, neural avalanches, and critical behavior near a first-order transition of a spiking neural network,” Phys. Rev. E, vol. 97, p. 062305, 2018.
  * [72] L. M. Ricciardi and L. Sacerdote, “The ornstein-uhlenbeck process as a model for neuronal activity,” Biological cybernetics, vol. 35, no. 1, pp. 1–9, 1979.
  * [73] N. Brunel, “Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons,” Journal of Computational Neuroscience, vol. 8, p. 183–208, 2000.

Supplemental Material for “On the critical signatures of neural activity”

## A. Experimental setting

### 1\. Surgical procedures

LFPs recordings are performed on Wistar rats, which are maintained under
standard environmental conditions in the animal research facility of the
Department of Biomedical Sciences of the University of Padova. All the
procedures are approved by the local Animal Care Committee (O.P.B.A.) and the
Italian Ministry of Health (authorization number 522/2018-PR). Young adult
rats aged 36 to 43 days and weighting between 150 and 200 g are anesthetized
with an intra-peritoneal induction mixture of tiletamine-xylazine (2 mg and
1.4 g/100 g body weight, respectively), followed by additional doses (0.5 mg
and 0.5 g/100 g body weight) every hour. The anesthesia level is constantly
monitored by testing the absence of eye and hind-limb reflexes and whiskers’
spontaneous movements. Each animal is positioned on a stereotaxic apparatus
where the head is fixed by teeth- and ear-bars. To expose the cortical area of
interest, an anterior-posterior opening in the skin is made in the center of
the head and a window in the skull is drilled over the somatosensory barrel
cortex at stereotaxic coordinates $-1\div-4$ AP, $+4\div+8$ ML referred to
bregma [61]. A slit in the meninges is then carefully made with fine forceps
at coordinates $-2.5$ AP, $+6$ ML for the subsequent insertion of the
recording probe. As a reference, the depth is set at 0 µm when the electrode
proximal to the chip tip touches the cortical surface. The neuronal activity
is recorded from the entire barrel cortex (from 0 to $-1750$ $\mu m$), which
is constantly bathed in Krebs’ solution (in mM: NaCl 120, KCl 1.99, NaHCO3
25.56, KH2PO4 136.09, CaCl2 2, MgSO4 1.2, glucose 11). An Ag/AgCl electrode
bathed in the extracellular solution in proximity of the probe is used as
reference.

### 2\. Recordings

LFPs are recorded through a custom-made needle which integrates a high density
array, whose electrodes are organized in a $64\times 4$ matrix. The operation
principle of the multi-electrode-arrays used to record LFPs is an extended
CMOS based EOSFET (Electrolyte Oxide Semiconductor Field Effect Transistor).
The recording electrodes are 7.4$\mu$m in diameter size and the needle is 300
$\mu$m in width and 10$mm$ long. The x- and y-pitch (i.e. the distance between
adjacent recording sites) are 32$\mu$m. The multiplexed signals are then
digitized by a NI PXIe-6358 (National Instruments) up to 1.25MS/s at 16bit
resolution and saved to disk by a custom LabVIEW acquisition software. The LFP
signal is sampled at 976.56 Hz and band-pass filtered (2-300 Hz). The dataset
analyzed for this work consist in 20 trials of basal activity lasting 7.22
seconds, that are recorded from 4 rats.

### 3\. Barrel cortex

The barrel cortex is the region of the primary somatosensory cortex (S1) that
encodes tactile sensory inputs from the rodents’ whiskers. The barrel presents
a high degree of segmentation into vertical columns and horizontal layers. The
whisker-related barrel column is a cylindrical structure spanning vertically
the six layers of the barrel cortex, although its border is defined
exclusively by spatially aligned sub-cellular structures in layer 4 called
barrels.

Each barrel-column is composed of an archetypal circuit that is repeated in
each column. Importantly, the barrels are laid out in a pattern that is nearly
identical to the whiskers on the rat’s snout that, together with the
underlying neural circuit, suggests a highly specific correspondence between
whiskers and barrel columns.

Figure S1: Comparison between the stationary distributions obtained in the
main text and the results $10^{3}$ simulations. Semi-transparent lines
represent different simulations. Filled areas of the plots represent one
standard deviation from the mean distribution. (a) Probability distribution of
$\mathcal{D}$. (b) Probability distribution of a single $v_{i}$.

## B. The extrinsic model

As explained in the main text, we introduce a model of $N$ variables
$(v_{1},\dots,v_{N})$ that are conditionally independent given the state of
the external input. In particular, let us assume that there are some
parameters $\bm{\xi}=(\xi_{1},\dots,\xi_{M})$ controlling the extrinsic
modulation, that is

$p(v_{i},v_{j},t\,|\,\bm{\xi})=p(v_{i},t\,|\,\bm{\xi})\,p(v_{j},t\,|\,\bm{\xi})\quad\quad\forall
i,j.$ (S1)

Since at this moment we assume no direct interaction is happening, we think
about the probability distribution $p(v_{i}|\bm{\xi})$ as shaped by extrinsic
activity alone. That is, once we specify $\bm{\xi}$, no interaction happens
between the $i$-th variable and the $j$-th variable. However, we now suppose
that we do not have access on the states of the external parameters. Hence, we
can only hope of describing to the joint probability distribution

$p(v_{1},\dots,v_{N},t)=\int
d\bm{\xi}\prod_{i=1}^{N}p(v_{i},t\,|\,\bm{\xi})p(\bm{\xi},t)$ (S2)

which is what we typically observe in an experimental setting. In general, due
to the marginalization over the external parameters this probability
distribution is not factorizable, and thus in this sense it is not trivial.

In general, we want to work in the stationary limit and thus there is a
complication we need to take into account. In particular, since the integrand
of Equation (S2) can be arbitrarily complicated, we would like to be able to
perform the stationary limit before the marginalization, i.e,

$\displaystyle\lim_{t\to\infty}p(v_{1},\dots,v_{N},t)=\int
d\bm{\xi}\prod_{i=1}^{N}$
$\displaystyle\lim_{t\to\infty}p(v_{i},t\,|\,\bm{\xi})\cdot\lim_{t\to\infty}p(\bm{\xi},t).$

This limit, of course, is not always true, but it does hold if we assume that
the time-scales of the two processes - the process for $\bm{v}$ and for
$\bm{\xi}$ \- are separated. In particular, we are interested in the limit in
which the timescale of $\bm{\xi}$ is slower than $\bm{v}$, so that $\bm{v}$
relaxes to its stationary state in a time-frame in which the external
modulation can be considered constant.

Figure S2: Comparison between the analytical expressions of the joint
probability distribution in Eq. (S8) and its factorization $p(v_{i},v_{j})$.
(a) The factorized distribution. (b) The joint distribution. (c) Comparison
between the two in a log-plot. (d) Relative difference between the two with
respect to the factorized distribution. (e-f) Comparison between the joint and
the factorized distributions, with the same parameters except for
$\mathcal{D}^{*}=5$.

Let us now describe the particular choices of the main text, where

$\frac{dv_{i}(t)}{dt}=-\frac{1}{\gamma_{i}}v_{i}(t)+\sqrt{\mathcal{D}(t)}\eta(t)$
(S3)

is the process that generates $\bm{v}$ and the single-parameter external
modulation is

$\mathcal{D}(t)=\begin{cases}\mathcal{D}^{*}&\text{if}\quad
D(t)\leq\mathcal{D}^{*}\\\ D(t)&\text{if}\quad
D(t)>\mathcal{D}^{*}\end{cases}$ (S4)

where

$\frac{dD(t)}{dt}=-\frac{1}{\gamma_{D}}D(t)+\sqrt{\theta}\eta(t).$ (S5)

For the above considerations to hold, we assume that
$\gamma_{D}\gg\gamma_{i}$. This is simply a time-scale separation limit. In
this limit, the process of $v_{i}$ reaches stationarity much faster than the
process of $\mathcal{D}$, thus the overall stationary distribution is the
stationary distribution $p(v_{i}|\mathcal{D})$ averaged over the stationary
distribution $p(\mathcal{D})$. Hence we have

$\displaystyle
p(\mathcal{D})=\frac{1}{2}\left[1+\text{Erf}\left(\frac{\mathcal{D}^{*}}{\sqrt{\theta\gamma_{D}}}\right)\right]\delta(\mathcal{D}-\mathcal{D}^{*})+\frac{H(\mathcal{D}-\mathcal{D}^{*})}{\sqrt{\pi\theta\gamma_{D}}}e^{-\frac{\mathcal{D}^{2}}{\theta\gamma_{D}}},$
(S6)

where $H$ is the Heaviside step function. Then, the computation of the
stationary probability distributions is quite easy. The single-unit
probability is given by

$\displaystyle p(v_{i})$ $\displaystyle=\int
d\mathcal{D}\,p(v_{i}|\mathcal{D})\,p(\mathcal{D})=\int
d\mathcal{D}\,\frac{1}{\sqrt{\pi\gamma_{i}\mathcal{D}}}\,e^{-\frac{v_{i}^{2}}{\gamma_{i}\mathcal{D}}}\,p(\mathcal{D})$
$\displaystyle=\frac{1+\text{Erf}\left(\frac{\mathcal{D}^{*}}{\sqrt{\theta\gamma_{D}}}\right)}{2\sqrt{\pi\mathcal{D}^{*}\gamma_{i}}}e^{-\frac{v_{i}^{2}}{\mathcal{D}^{*}\gamma_{i}}}+\frac{1}{\sqrt{\pi^{2}\theta\gamma_{D}\gamma_{i}}}\int_{\mathcal{D}^{*}}^{\infty}\frac{dD}{\sqrt{D}}e^{-\left[\frac{D^{2}}{\theta\gamma_{D}}+\frac{v_{i}^{2}}{D\gamma_{i}}\right]}$
(S7)

and, in principle, we can compute the probability distributions for any number
of variables $\\{v_{i}\\}$ in the same way.

We note that although before the marginalization the two-point probability
distribution is factorizable

$p(v_{i},v_{j},\mathcal{D})=p(v_{i},v_{j}|\mathcal{D})p(\mathcal{D})=p(v_{i}|\mathcal{D})p(v_{j}|\mathcal{D})p(\mathcal{D}),$

the marginalization itself breaks the factorization, namely

$\int d\mathcal{D}\,p(v_{i},v_{j},\mathcal{D})\neq\left[\int
d\mathcal{D}\,p(v_{i},\mathcal{D})\right]\left[\int
d\mathcal{D}\,p(v_{j},\mathcal{D})\right].$

Hence, we can write down the two-unit joint probability as

$\displaystyle
p(v_{i},v_{j})=\frac{1+\text{Erf}\left(\frac{\mathcal{D}^{*}}{\sqrt{\theta\gamma_{D}}}\right)}{2\pi\mathcal{D}^{*}\sqrt{\gamma_{i}\gamma_{j}}}e^{-\frac{1}{\mathcal{D}^{*}}\left(\frac{v_{i}^{2}}{\gamma_{i}}+\frac{v_{j}^{2}}{\gamma_{i}}\right)}+\frac{1}{\sqrt{\gamma_{i}\gamma_{j}\gamma_{D}\pi^{3}\theta}}\int_{\mathcal{D}^{*}}^{\infty}\frac{dD}{D}e^{-\frac{1}{D}\left(\frac{v_{i}^{2}}{\gamma_{i}}+\frac{v_{j}^{2}}{\gamma_{i}}\right)}e^{-\frac{D^{2}}{\theta\gamma_{D}}}$
(S8)

where the second term is not factorizable.

From these probabilities we can immediately show that

$\expectationvalue{v_{i}v_{j}}-\expectationvalue{v_{i}}\expectationvalue{v_{j}}=0\quad\forall
i\neq j$ (S9)

and, in general, all the expectation values where a variable $v_{i}$ appears
an odd number of times vanish. This follows from the simple fact that, for
instance,

$\displaystyle\expectationvalue{v_{i}v_{j}}=\int
dv_{i}\,dv_{j}\,p(v_{i},v_{j})v_{i}\,v_{j}=\int
d\mathcal{D}\,p(\mathcal{D})\left(\int
dv_{i}\,v_{i}\,p(v_{i}|\mathcal{D})\right)\left(\int
dv_{j}\,v_{j}\,p(v_{j}|\mathcal{D})\right)=0\quad\forall i\neq j$

since $\int dv_{i}\,v_{i}\,p(v_{i}|\mathcal{D})=0$. In fact, let us note that
this is not a general feature of these kind of models, but it is rather a
property inherited from the choice of an OU process for $v_{i}$ given
$\mathcal{D}$, for which all the odds moment vanish. The diagonal entries of
the covariance matrix, on the other hand, are given by

$\langle
v_{i}^{2}\rangle=\frac{\gamma_{i}\mathcal{D}^{*}}{2}\left[1+\text{Erf}\left(\frac{\mathcal{D}^{*}}{\sqrt{\theta\gamma_{D}}}\right)\right]+\sqrt{\frac{\gamma_{i}^{2}\theta\gamma_{D}}{16\pi}}e^{-\frac{(\mathcal{D}^{*})^{2}}{\theta\gamma_{D}}}$

Thus, in this particular setting, the units are uncorrelated and have a
variance that is independent on the other units. At the same time, however, it
is clear that in general $p(v_{i},v_{j})\neq p(v_{i})p(v_{j})$ so while
uncorrelated the variables are not independent.

Figure S3: Comparison between the analytical expressions of the joint
probability distribution $p(v_{i},v_{j})$ given by Eq. (S8) and the one
obtained from with $10^{3}$ simulations. We also show the corresponding
results for the factorized probabilities. The blue line corresponds to the
analytical expression of $p(v_{i},v_{j})$. The corresponding dots represent
the histogram of the distribution obtained from $10^{3}$ simulations, and the
semitransparent filled areas represent one standard deviation from this
estimate. Similarly, the gray dashed line represent the analytical expression
of $p(v_{i})p(v_{j})$. (a) Section along the $v_{j}$ direction for small
$v_{j}$, so that we are looking at the tails of the distribution. Even though
the estimate along the tails is noisy, we clearly see that the estimate from
the simulations lies along the analytical prediction. (b) As before, but for
higher $v_{j}$. (c-d) As before, but with values of $v_{j}$ close to zero so
we look at the bulk of the distribution near its peak. Even though joint
probability and its factorization now are more similar, once more the estimate
from the simulation match the analytical expression $p(v_{i},v_{j})$.

## C. Testing the predictions of the extrinsic model

We now show that simulations of the extrinsic model described by Equation S3
agree with the analytical results. Whenever not specified, we assume that the
parameters of the model are given by $\mathcal{D^{*}}=0.3$, $\theta=1$,
$\gamma_{D}=10$ together with $\gamma_{i}=0.1$, $\gamma_{j}=0.5$. Thus, we are
in the limit of timescale separation considered in the main text.

Let us begin with Figure S1. Albeit trivial, we first check in Figure S1a that
the stationary distribution of $\mathcal{D}$ is indeed the one of Eq. (S6) and
in Figure S1b that the stationary distribution of a single $v_{i}$ does
correspond to the analytical expression of Eq. (B. The extrinsic model). If we
compare this distribution to a standard distribution of an Ornstein-Uhlenbeck
process with a diffusion coefficient equal to the mean
$\expectationvalue{\mathcal{D}(t)}$ [54] we immediately see that the
distribution of our model is considerably more peaked around zero and displays
longer tails. Indeed, one expects that due to the fact that
$\mathcal{D}^{*}<\expectationvalue{\mathcal{D}(t)}$ the system tends to wander
more easily close to zero, especially in the time windows where the diffusion
coefficient is constant and equal to $\mathcal{D}^{*}$. At the same time, the
fact that $\mathcal{D}(t)$ can change in time favors the presence of values of
$v$ that are larger in absolute value, which is the mechanism at the origin of
the bursty behavior seen in the main text.

In Figure S2 we look instead at the properties of the joint probability
distribution $p(v_{i},v_{j})$. The most natural quantity to compare this
distribution with is its factorization $p(v_{i})p(v_{j})$, which is equivalent
to ignoring the feedback effects between $v_{i}$ and $v_{j}$ due to the shared
extrinsic modulation of $\mathcal{D}(t)$. Since we are setting
$\mathcal{D}^{*}=0.5$, we expect that these effects are going to be
particularly relevant for the dynamics of the model. In particular, in Figure
S2c-d we see that the most important differences between the two occur in the
tails of the two-dimensional distribution, with the joint distribution
typically showing dramatically longer tails. This translates to the fact that
far-from-zero values of the two variables can occur more easily at the same
time.

The situation is completely reversed when we increase $\mathcal{D}^{*}$. In
Figure S2e-f we see that for $\mathcal{D}^{*}=8$ the joint probability
distribution and the factorized distribution are almost indistinguishable.
Hence, this example shows explicitly that if $\mathcal{D}^{*}$ is high enough
the dependence induced by the extrinsic modulation vanishes.

Let us keep focusing on the case $\mathcal{D}^{*}=0.5$ for the time being. In
Figure S3a-d we compare the one-dimensional sections of the analytical
expression of the joint probability distribution with the results of $10^{3}$
simulations of the model, together with the sections of the factorized
distribution. The joint distribution estimated from the simulation matches
particularly well the analytical prediction. Once more, and perhaps more
clearly, in Figure S3a-b we see the stark difference that emerges along the
tails between the joint probability distribution and its factorization.
Interestingly, panel (c) and panel (d) show that the situation in the bulk of
the distribution is reversed with respect to the tails and now the joint
probability distribution is more peaked with respect to is factorization,
albeit only slightly. That is, the modulation in the low $\mathcal{D}^{*}$
regime favors both large values of $v_{i}$ and $v_{j}$ and values very close
to zero.

Overall, this brief analysis showed us how the bursty behavior that we see for
small values of $\mathcal{D}^{*}$ emerges from the underlying probability
distributions, which in turn emerge from the simple marginalization that
occurs in Eq. (S2). Similar arguments, albeit impractical, could be carried
out for the probability distributions beyond the two-point ones. In a sense,
one could argue that the fundamental properties of the model are inherited
from the fact that there are some unobserved physical quantities, and these
are the quantities that drive the global response of the single variables.

For the specific case of a double Ornstein-Uhlenbeck process we chose, the net
effect of the marginalization is the widening of the tails of both the one-
point $p(v)$ and the two-point $p(v_{i},v_{j})$ probability distributions when
$\mathcal{D}^{*}$ is small enough. As we increase $\mathcal{D}^{*}$, this
effect becomes less and less important until it is completely negligible. In
this sense, we can effectively think of $\mathcal{D}^{*}$ as a control
parameter that changes the qualitative behavior of the system. Most
importantly, the fact that the tails of the joint probability distribution are
wider when $\mathcal{D}^{*}$ is small reflects dynamically in the emergence of
a non trivial coordination between the variables, from which in turn power-law
avalanches emerge.

Figure S4: Comparison between (a) the correlations of the data and (c) the
correlations of the interacting model after solving the inverse problem. Panel
(b) shows the inferred matrix $\tilde{A}$, which is independent on the
parameters of the model.

## D. The multivariate Ornstein-Uhlebeck process

The process studied in the main text is a multivariate Ornstein-Uhlebeck
process [62] of the form

$d\bm{v}(t)=-A\bm{v}(t)dt+B(t)d\bm{W}(t),$ (S10)

where in our case $B(t)$ is a diagonal matrix whose diagonal elements are
given by $\sqrt{\mathcal{D}(t)}$. In the case of decoupled units, which we use
to model the extrinsic activity, the matrix $A$ is again diagonal with entries
$A_{ij}=\delta_{ij}/\gamma_{i}$. The formal solution of this stochastic
process is given by

$\displaystyle\bm{v}(t)=$
$\displaystyle\exp\left(-At\right)\bm{v}(0)+\int_{0}^{t}\exp\left[-A\left(t-t^{\prime}\right)\right]B\left(t^{\prime}\right)d\bm{W}\left(t^{\prime}\right).$

In particular, since in our case the matrix $B$ is on itself a stochastic
variable, from now on we will focus on the original variables conditioned on
$B$, which in turn we assume follows its own stationary distribution $p(B)$.
We call these conditioned variables $\bm{v}(t,B)$. If the eigenvalues of $A$
have all positive real parts, a stationary solution of Eq. S10 exists and it
is of the form [63]

$\displaystyle\bm{v}_{s}(t,B)=\int_{-\infty}^{t}\exp\left[-A\left(t-t^{\prime}\right)\right]Bd\bm{W}(t^{\prime}).$

Hence, the stationary covariance matrix conditioned on $B$ reads as

$\displaystyle\sigma(B)$
$\displaystyle=\left\langle\bm{v}_{s}(t,B),\bm{v}_{\mathbf{s}}^{\mathrm{T}}(t,B)\right\rangle$
$\displaystyle=\int_{-\infty}^{t}dt^{\prime}\exp\left[-A(t-t^{\prime})\right]BB^{\mathrm{T}}\exp\left[A^{\mathrm{T}}(t-t^{\prime})\right]$

which solves the algebraic equation

$\displaystyle A\sigma(B)+\sigma(B)A^{\mathrm{T}}=M,$

where the matrix $M$ is given by

$\displaystyle M=$
$\displaystyle\int_{-\infty}^{t}A\exp\left[-A\left(t-t^{\prime}\right)\right]BB^{\mathrm{T}}\exp\left[-A^{\mathrm{T}}\left(t-t^{\prime}\right)\right]dt^{\prime}+\int_{-\infty}^{t}\exp\left[-A\left(t-t^{\prime}\right)\right]BB^{\mathrm{T}}\exp\left[-A^{\mathrm{T}}\left(t-t^{\prime}\right)\right]A^{\mathrm{T}}dt^{\prime}.$

We can rewrite it as

$\displaystyle
M=\int_{-\infty}^{t}\frac{d}{dt^{\prime}}\left\\{\exp\left[-A\left(t-t^{\prime}\right)\right]BB^{\mathrm{T}}\exp\left[-A\left(t-t^{\prime}\right)\right]\right\\}dt^{\prime}$

Carrying out the integral, we find that the lower limit vanishes due to the
assumed positivity of the eigenvalues of $A$ and hence only the upper limit
remains, giving

$\displaystyle A\sigma(B)+\sigma(B)A^{\mathrm{T}}=BB^{\mathrm{T}},$ (S11)

which is a continuous Lyapunov equation [64] for the covariance matrix
$\sigma(B)$ [53]. Then, we only need to marginalize over $B$ and we obtain the
equation for the covariance matrix of our original variables $\bm{v}_{s}(t)$

$A\sigma+\sigma A^{\mathrm{T}}=Q,$ (S12)

where for the sake of brevity we call $\sigma$ the covariance matrix of the
original variables and $Q$ is a diagonal matrix whose elements are given by

$Q_{ij}=\delta_{ij}\int_{\mathcal{D}^{*}}^{\infty}\mathcal{D}\,p(\mathcal{D})\,d\mathcal{D}:=\delta_{ij}f(\mathcal{D}^{*},\gamma_{D},\theta).$
(S13)

Then, taking the transpose of Eq. (S12), and assuming that the matrix A is
symmetric, we end up with a Lyapunov equation for the matrix A. Thus, the
inverse problem of reconstructing the connectivity matrix from the data
reduces to solving numerically $\sigma A+A\sigma=Q$, starting from the
correlation matrix of the data and given a set of parameters for the model.

Hence we end up with a model

$\dot{v}_{i}(t)=-\sum_{j}A_{ij}v_{j}(t)+\sqrt{\mathcal{D}(t)}\xi_{i}(t)$ (S14)

where $A_{ij}$ depends on the the parameters of the stochastic modulation
$(\mathcal{D}^{*},\gamma_{D},\theta)$. Notice that if we write
$\tilde{A}_{ij}=A_{ij}/f(\mathcal{D}^{*},\gamma_{D},\theta)$ we need to solve
the Lyapunov equation $\sigma\tilde{A}+\tilde{A}\sigma=\mathbb{1}$ that only
depends on $\sigma$, the correlation matrix of the data. If we introduce
$\tilde{\mathcal{D}}=\mathcal{D}/f$ and $\tilde{v}_{i}=v_{i}/\sqrt{f}$ we can
write

$\dot{\tilde{v}}_{i}(t)=-\sum_{j}\tilde{A}_{ij}\tilde{v}_{j}(t)+\sqrt{\tilde{\mathcal{D}}(t)}\xi_{i}(t)$

and clearly

$\expectationvalue{\tilde{v}_{i}\tilde{v}_{j}}-\expectationvalue{\tilde{v}_{i}}\expectationvalue{\tilde{v}_{j}}=\frac{\expectationvalue{v_{i}v_{j}}-\expectationvalue{v_{i}}\expectationvalue{v_{j}}}{f}$

so the correlation between $\tilde{v}_{i}$ and $\tilde{v}_{j}$ is proportional
to the correlation between $v_{i}$ and $v_{j}$. This means that at different
$(D^{*},\theta,\gamma_{D})$ we simply find a rescaled interaction matrix
$A_{ij}$, but the scaling of correlation length does not change.

The stationary probability distribution solution of the interacting model can
be written, in general, as

$p(v_{1},\dots,v_{N})=\frac{1+\mathrm{Erf}\left[\frac{\mathcal{D}^{*}}{\sqrt{\theta\gamma_{D}}}\right]}{2\sqrt{(\pi\mathcal{D}^{*})^{N}\det\Sigma}}e^{-\frac{1}{\mathcal{D}^{*}}\mathbf{v}^{T}{\Sigma}^{-1}\mathbf{v}}+\frac{1}{\sqrt{(\gamma_{D}\theta)^{N}\pi^{N+1}\det\Sigma}}G_{N}\left(\frac{\mathcal{D}^{*}}{\sqrt{\theta\gamma_{D}}},\frac{\mathbf{v}^{T}{\Sigma}^{-1}\mathbf{v}}{\sqrt{\theta\gamma_{D}}}\right)$
(S15)

where $\left({A}{\Sigma}+{\Sigma}{A}^{T}\right)/2=\mathbb{1}$ and

$G_{N}(\alpha,\beta)=\int_{\alpha}^{\infty}\frac{dx}{x^{N/2}}e^{-\frac{\beta}{x}-x^{2}}.$
(S16)

In general we could define a multivariate information between these $N$
variables and the results of the main text would not change. In practice,
however, it is very hard to perform the related numerical integration if $N$
is large. Therefore, in the main text we show an example of two interacting
units that interact through the matrix element $\tilde{A}_{12}=2$. The mutual
information of this model receives a distinct (and constant) contribution from
the interaction matrix and from the modulation induced by $\mathcal{D}^{*}$.
In fact, if we take $A=\mathrm{diag}(\gamma_{1},\dots,\gamma_{N})$, i.e., we
consider the extrinsic model, the only the contribution to the mutual
information comes from the external modulation, whereas the addition of
interactions simply shift the mutual information by a constant value at all
$\mathcal{D}^{*}$.

## E. Scaling of the correlation length

To study if the correlations of the system under study exhibit critical-like
properties we determine the correlation length $\xi$ of the system at various
system’s sizes. The correlation length can be defined as the average distance
at which the correlations of the fluctuations around the mean crosses zero
[46], and it is known to diverge at criticality in the thermodynamic limit
[31].

For finite systems, however, this behavior can be demonstrated by showing that
the correlation length grows with system size. Thus, we first compute for each
time series their fluctuations around the mean, namely

$\tilde{v}_{i}(t)=v_{i}(t)-\frac{\sum_{i=1}^{N}v_{i}(t)}{N}$ (S17)

where $1/N\sum_{i=1}^{N}v_{i}(t)$ is the mean activity, i.e. the mean at time
$t$ computed over the $N$ channels, and $v_{i}$ is the activity at channel
$i$. By definition the mean of the fluctuations vanishes, i.e.
$\frac{\sum_{i=1}^{N}\tilde{v}_{i}(t)}{N}=0,\forall t$. As said, different
sizes (portions of the array) of the system are selected and, importantly, the
mean activity is computed for each system size, considering the channels
inside the observation window.

In particular, since the maximum system size (corresponding to $N$ channels)
both in our experimental data and in our model is fixed, we investigate how
$\xi$ changes with system sizes corresponding to different subsamples from the
multi-array probe [34]. In fact, through simulations on control models that
display a critical point [65], subsampling has been shown to be practically
equivalent to consider systems of different sizes. We assume that the units of
our model have the same topology of our data, i.e. we assume that the units
are placed as the channels in the $55\times 4$ array of our experimental
setup.

Figure S5: Avalanche statistics in four different rats. (1a-4a) The
distribution of the avalanches’ sizes is consistently a power-law, with an
exponent that slightly depends on the single rat. (1b-4b) The avalanche
durations are once more power-law distributed in all rats with some
variability in the exponents, even though the range accessible with the
experimental setup only covers two decades. (1c-4c) The crackling-noise
relation, however, is consistently satisfied in each rat.

Since in our case the array shape is rectangular, we consider the number of
rows as the relevant dimension and build subsampled systems of size $L\times
4$, with $L$ that decreases from the maximum of $55$ channels up to $5$
channels. Next, for each system’s subset, we compute the average correlation
function of the fluctuations between all pairs of channels separated by a
distance $r$,

$\expectationvalue{C(r)}=\expectationvalue{\frac{\expectationvalue{\left(\tilde{v}_{i}-\overline{\tilde{v}}_{i}\right)\left(\tilde{v}_{j}-\overline{\tilde{v}}_{j}\right)}_{t}}{\sigma_{\tilde{v}_{i}}\sigma_{\tilde{v}_{j}}}}_{i,j}$
(S18)

where $\expectationvalue{\cdot}_{t}$ stands for the average over time,
$\expectationvalue{\cdot}_{i,j}$ the average over all pairs of channels
separated by a distance $r$ and

$\displaystyle\overline{\tilde{v}}_{i}=\frac{1}{T}\sum_{t=1}^{T}\tilde{v_{i}}\left(t\right)$
$\displaystyle\sigma_{\tilde{v}_{i}}^{2}=\frac{1}{T}\sum_{t=1}^{T}\left(\tilde{v_{i}}\left(t\right)-\overline{\tilde{v}}_{i}\right)^{2}$

with $T$ is the length of the time series. Then $\xi$ is computed as the zero
of the correlation function $C(r=\xi)=0$. To reduce the noise effects, results
were averaged across all possible sub-regions for any given size. Then the
$\xi$ are plotted against the relative system size $L$ and the slope of the
fit is obtained through a linear regression.

## F. Avalanches statistics in LFPs data

### 1\. LFPS peaks detection and avalanche analysis

For the detection of LFP events, the standard deviation (SD) and the mean of
the signal was computed for each channel. In order to distinguish real events
from noise, a three SD threshold was chosen basing on the distribution of the
signal amplitudes which significantly deviated from a Gaussian best fit above
that threshold. Both negative and positive LFPs (i.e., nLPFs adn pLFPs,
respectively) were considered as events in accordance with previous works
[66]. Within our specific experimental settings, one reason is that across the
depth of the cortex there are polarity changes of the LFP signal because of
compensatory capacitive ionic currents particularly along dendrites of
pyramidal cells [67]. Since in our experiments electrodes span multiple
cortical layers, both nLFPs and pLFPs were found and detected. Moreover,
alternatively, pLFPs can be related to activation of populations of inhibitory
neurons leading to inhibitory outward postsynaptic currents, which also
justifies their inclusion in the events count. For detection, each deflection
was considered terminated only after it crossed the mean of the signal. After
defining events of activity, the avalanche statistics are analyzed through the
usual methods of the literature. In particular, the data are temporally
binned, avalanches are defined as sequences of bins that present activity, and
an avalanche ends once an empty bin is found - the temporal bin chosen is the
average inter-event interval, as it is typically done in avalanche analysis
[1]. Then, the distribution $p(s)$ of the avalanches sizes - the number of
events in each avalanche - and of the avalanche duration $p(T)$ are computed
and fitted using a corrected maximum likelihood method [44, 43]. In
particular, the power law distribution is tested following the method proposed
in [68, 44]. Avalanche sizes and lifetimes are fitted with discrete power laws
$p(y;\alpha)=\frac{y^{-\alpha}}{\sum_{x=x_{min}}^{x=x_{max}}x^{-\alpha}}$. The
parameter $x_{max}$ is set to the maximum observed size or duration. $x_{min}$
is selected as the one that minimizes the Kolmogorov-Smirnov distance (KS)
between the cumulative distribution function (CDF) of the data and the CDF of
the theoretical distribution fitted with the parameter that best fits the data
for $y\geq x_{min}$ [68]. To assess goodness-of-fit we compared the
experimental data against 1000 surrogate datasets drawn from the best-fit
power law distribution with the same number of samples as the experimental
dataset. The p-value of the power-law fit was defined as the fraction of the
KS statistics of the surrogates which were greater than the KS statistic for
the experimental data. The data were considered power law distributed if the
the p-value was greater than 0.1. We also take into account the fact that
while maximum likelihood methods rely on the independence assumption, actual
data often display correlations, and this may lead to false rejection of the
statistical laws [44]. As the author of [44] suggest, before performing the
fit and assessing p-values, we undersample the data in order to decorrelate
them, by estimating the time $\tau^{*}$ after which two observations (e.g.,
the avalanche sizes) are independent from each other, as done in [44, 43].

### 2\. Avalanches statistics

Here we report the avalanches statistics from the other rats that we analyzed.
The avalanches statistics is computed by considering all available the $20$
trials of basal activity for each rat, that are $7.22\mathrm{s}$ long. Results
are reported in Figure S5. Inter-rat variability is present with respect to
avalanche exponents, and is expected as found in previous experiments [22, 66,
69]. Moreover, a theoretical explanation for the difference in avalanche
exponents has been also recently proposed as signature of quasi-criticality
[30]. However, the fundamental point is that the crackling-noise relation is
always verified compatibly with the experimental errors, a feature that is
usually considered a hallmark of criticality [22].

Figure S6: Collapse of the average profile of avalanches of varying duration
in the extrinsic model, for the low $\mathcal{D}^{*}$ regime. (a) Profile of
the avalanches before the rescaling. (b) If we rescale with an exponent
$\delta\approx 1.33$, which is remarkably close to the one found in the main
text through the crackling-noise relation, we obtain an optimal collapse onto
the same scaling function.

## G. Additional avalanches statistics in the model

Another signature of criticality is the collapse of the average profile of
avalanches of widely varying duration onto a single scaling function. For
avalanches of duration $T$ we can write down the average number of firing at
time $t$ as $s(t,T)=T^{\delta-1}F(t/T)$ where F is a universal scaling
function that determines the shape of the average temporal profile.
$\expectationvalue{S(T)}$ and $s(t,T)$ are related by
$\expectationvalue{S(T)}=\int_{0}^{T}s(t,T)dt$. At the critical point we
expect that plots of $t/T$ versus $s(t,T)T^{1-\delta}$ for different $T$ will
collapse onto the same universal scaling function [27].

Thus, finding the exponent for which the goodness of the collapse is higher
provides another way to estimate $\delta$. For testing the avalanche shape
collapse, we used the methodology introduced in [70]. To determine the quality
of the collapse, the averaged and rescaled avalanche profiles of different
lifetimes $F(t/T)=T^{1-\delta}s(t/T,T)$ are first linearly interpolated at
$1000$ points along the scaled duration. The variance across the different
$F(t/T)$ is calculated at each interpolated point, and the shape collapse
error $\epsilon(\delta)$ is then defined as the mean variance divided by the
squared span of the avalanche shapes, where the span equals the maximum minus
the minimum value of all rescaled avalanche profiles. In the presented
analysis, avalanche shapes of $T>10$ bins with at least 10 samples were used.

The collapse has been tested on the extrinsic model in the case of low
$\mathcal{D}^{*}$, and the results are plotted in Figure S6. We find that the
exponent that minimizes $\epsilon(\gamma)$ turns out to be $\approx 1.33$,
close to the estimates of $\delta$ found through the linear fit of average
size given duration and through the prediction of the crackling-noise relation
in the main text. Again, it is also close to the apparently super-universal
exponent found in [22] and in [28].

## H. The crackling noise relation and its link with criticality

The so called “crackling noise” relation among avalanches’ distributions’
exponents [27] is usually considered a strong sign of criticality [22],
however it is possible to show that it can be valid under general conditions
that do not necessarily require criticality. Indeed, following arguments
similar to the one proposed in [71], if we express the size $S$ of an
avalanche as a function of its duration $T$ by assuming that $S\sim
T^{\delta}$, i. e. we assume that fluctuations in the size fixed the duration
can be disregarded, the following equivalence holds:

$p(S(T))\left|\frac{dS}{dT}\right|dT=p(T)dT$

Then, if we assume that both $p(S(T))$ and $p(T)$ are described by power laws,
with exponents respectively $\tau$ and $\tau_{t}$, we see that the above
equation reduces to

$(T^{\delta})^{-\tau}\delta T^{\delta-1}=T^{-\tau_{t}}$

that leads ideed to the relation

$\delta=\frac{\tau_{t}-1}{\tau-1}.$

Thus, to derive the above relation we had to assume that $p(S)$ and $p(T)$ are
well described by power laws, and that fluctuations in the size given the
duration are negligible ($S\sim T^{\delta}$). These features are certainly
proper of critical points (where moreover the exponent $\delta$ can be
expressed as a combination of other critical exponents [27]), but can be found
to be valid also away from criticality, as our modeling framework shows.

## I. Alternative models

The multivariate Ornstein-Uhlenbeck (MOU) process has been previously used in
the literature [50, 51, 52, 53] to model neural activity at large scale (e.g.
fMRI signals). However, it is a paradigmatic model which, albeit simple enough
to allow for an analytical treatment, does not account for many biological
aspects. Thus, in order to show that our results generalize to more complex
situations, here we introduce a more biologically sound model that can
describe the population-level neural activity measured by our LFPs data, the
Wilson-Cowan model. It includes both excitatory and inhibitory synapses and
non-linearities in the transfer function. We extend the framework of the
extrinsic model presented in the main text, and we also provide a possible
biological generative mechanism for the external time varying modulation.

Figure S7: (a) External input $h=\frac{E^{(h)}+I^{(h)}}{2}$ in a balanced
state ($\omega^{(h)}_{E}=50.05$ and $\omega^{(h)}_{I}=49.05$) with high level
of noise amplitude $\sigma^{(h)}$ (in orange), and low level of noise
amplitude (in blue). When the noise level is high enough bursty behavior is
present, that is reflected also on the firing rate $\frac{E_{i}+I_{i}}{2}$ of
the units that receive $h$ as input (Figure (b), in orange). Instead, if the
noise level is low, $h$ converges to the up state predicted by the model in
the absence of noise (Fig. (a), blue), hence also the firing rate of the units
(Figure (b), blue) doesn’t display transient behaviors.

The model considered is a stochastic Wilson-Cowan model [56, 55], which
describes the dynamics of the activity, i.e., the density of active neurons,
for the two subpopulations of excitatory (E) and inhibitory (I) neurons. Each
unit $i$, i.e., each neural population, is now represented by a subpopulation
of excitatory neurons and a subpopulation of inhibitory neurons, whose
activity evolves according to the following set of stochastic equations:

$\begin{cases}\frac{dE_{i}}{dt}=-\alpha
E_{i}+(1-E_{i})f(\omega_{E}E_{i}-\omega_{I}I_{i}+h)+\sqrt{(\alpha
E_{i}+(1-E_{i})f(\omega_{E}E_{i}-\omega_{I}I_{i}+h))}\eta_{E_{i}}\\\
\frac{dI_{i}}{dt}=-\alpha
I_{i}+(1-I_{i})f(\omega_{E}E_{i}-\omega_{I}I_{i}+h)+\sqrt{(\alpha
I_{i}+(1-I_{i})f(\omega_{E}E_{i}-\omega_{I}I_{i}+h))}\eta_{I_{i}}\end{cases}$
(S19)

where $\alpha$ is the rate of spontaneous activity decay, $\omega_{E,I}$ are
the synaptic efficacies, $s=\omega_{E}E_{i}-\omega_{I}I_{i}+h$ is the averaged
incoming current, $\eta_{E,I}$ are uncorrelated Gaussian white noises, with
network-size dependent amplitude $\sigma\sim\frac{1}{\sqrt{\mathcal{N}}}$
[55], interpreted in the Ito sense. The incoming current is modeled as the sum
of all synaptic inputs plus an external input current $h$, filtered by a
sigmoid response function

$\begin{cases}f(s)=\tanh(s)\hskip 6.25963pts\geq 0\\\ f(s)=0\hskip
34.14322pts<0.\end{cases}$ (S20)

We consider a case in which the units are inhibition dominated, i.e. when
$\omega_{I}>\omega_{E}$ ($\omega_{I}=7$, $\omega_{E}=6.8$). In this situation
the noise level is not going to produce significant changes in the dynamics,
and the noise amplitude of $\eta_{E_{i}}$, $\eta_{I_{i}}$ is set to a low
value, $\sigma=10^{-4}\,\forall i$. Similarly, the value of $\alpha$ is not
going to affect the inhibition dominated dynamics, and it is set to
$\alpha=1$. For each unit, we consider the quantity
$\Sigma_{i}=(E_{i}+I_{i})/2$, that describes the firing rate of the neural
population $i$.

As in the main text, we consider an external, time-varying input that enters
the model through $h$. As a potential candidate for a biological realization
of this external, stochastic driving, we consider the effective input that
comes from other, yet unobserved, neural populations. Hence, we model $h$ as
the firing rate $h=(E^{(h)}+I^{(h)})/2$ of another Wilson-Cowan model. In
particular, we take it to be in the balanced state
$0\ll\omega^{(h)}_{0}\ll\omega^{(h)}_{\rm
sum}=\omega^{(h)}_{E}+\omega^{(h)}_{I}$ and with a sufficiently high level of
noise, so that the resulting firing rate spends most of the time close to the
down state, while showing frequent bursts of activity [55]. Note that, since
the balanced state requires that $\omega_{0}\ll\omega_{S}$, it does not
coincide with the “critical” state $\omega_{0}=0$ \- being balanced is not a
sufficient condition and neither necessary for being critical.

We end up with a time-varying modulation

$\frac{dh}{dt}=\frac{d}{dt}\left[\frac{E^{(h)}+I^{(h)}}{2}\right]$ (S21)

that is qualitatively similar to the one considered in the main text,
alternating periods of silence and bursts, but which is more realistic from a
biological point of view. In particular, we choose $\omega^{(h)}_{E}=50.05$,
$\omega^{(h)}_{I}=49.95$, $h^{(h)}=10^{-3}$, $\alpha^{(h)}=0.1$, and,
importantly, $\sigma^{(h)}$, the amplitude of the noise for both the
excitatory and the inhibitory populations, is increased to a relatively high
value, $1.2\times 10^{-3}$ (i e. corresponding to a finite but still large
network size), so that the up state can be destabilized by the noise.
Intuitively, we are considering a case in which $i=1,\dots,N$ populations of
neurons evolve according to the WC model in the inhibition dominated phase and
are all receiving the same input by another population of neurons in a
balanced state.

From each firing rate $\Sigma_{i}$ we generate spike trains and analyze
avalanches in the usual way, by temporal binning the events with the average
inter-event interval (see the main text). Again, we find that avalanche sizes
and durations are power law distributed, satisfying the crackling noise
relation (Fig. S8). We note that we do not expect to find the exponents of
avalanches’ sizes and durations found in [56], even though $h$ itself is a WC
model. In fact, in [56] avalanches were found by considering the excursions of
$h$ over a small threshold, in contrast to the study of spike trains carried
out here. Moreover, and most importantly, the external input $h$ that we are
considering feeds the units $i$ passing through the non linearity of $f(s)$,
hence $h$ and $\Sigma_{i}$ do not coincide.

Thus, we have shown that the results presented in the main text, and in
particular the ones of the extrinsic model, extend well beyond the
paradigmatic example of a multivariate Ornstein-Uhlenbeck process with a
shared external input. With respect to the approach presented in this section,
the MOU has the significant advantage of being amenable to analytical
treatments and to the inclusion of interactions reconstructed from the data.
However, the WC model has a deeper biological relevance and proves that the
mechanism proposed here for avalanche generation appears to be quite general.
Finally, we note that $h=(E^{(h)}+I^{(h)})/2$ is only one of the possible
choices that could be considered for the external modulation. Other candidates
include a simple Ornstein Uhlenbeck process, that is often used to model noisy
neural activity [72] or the Brunel’s model in the Synchronous Irregular phase
[73], which is known to be able to reproduce avalanching behavior [19].

Figure S8: Avalanches with Wilson Cowan units.