# Theoretical foundations of studying criticality in the brain

Yang Tian  1,2  Zeren Tan  3  Hedong Hou  4  Guoqi Li  5,6  Aohua Cheng  7
Yike Qiu  7  Kangyu Weng  7  Chun Chen  1  Pei Sun  1  1  2  3  4  5  6  7

###### Abstract

Criticality is hypothesized as a physical mechanism underlying efficient
transitions between cortical states and remarkable information processing
capacities in the brain. While considerable evidence generally supports this
hypothesis, non-negligible controversies persist regarding the ubiquity of
criticality in neural dynamics and its role in information processing.
Validity issues frequently arise during identifying potential brain
criticality from empirical data. Moreover, the functional benefits implied by
brain criticality are frequently misconceived or unduly generalized. These
problems stem from the non-triviality and immaturity of the physical theories
that analytically derive brain criticality and the statistic techniques that
estimate brain criticality from empirical data. To help solve these problems,
we present a systematic review and reformulate the foundations of studying
brain criticality, i.e., ordinary criticality (OC), quasi-criticality (qC),
self-organized criticality (SOC), and self-organized quasi-criticality (SOqC),
using the terminology of neuroscience. We offer accessible explanations of the
physical theories and statistic techniques of brain criticality, providing
step-by-step derivations to characterize neural dynamics as a physical system
with avalanches. We summarize error-prone details and existing limitations in
brain criticality analysis and suggest possible solutions. Moreover, we
present a forward-looking perspective on how optimizing the foundations of
studying brain criticality can deepen our understanding of various
neuroscience questions.

###### keywords:

Non-equilibrium criticality, Neural avalanches, Neural dynamics, Directed
percolation

††articletype: Review

Department of Psychology & Tsinghua Laboratory of Brain and Intelligence,
Tsinghua University, Beijing, 100084, China. Laboratory of Advanced Computing
and Storage, Central Research Institute, 2012 Laboratories, Huawei
Technologies Co. Ltd., Beijing, 100084, China. Institute for Interdisciplinary
Information Science, Tsinghua University, Beijing, 100084, China. UFR de
Mathématiques, Université de Paris, Paris, 75013, France. Institute of
Automation, Chinese Academy of Science, Beijing, 100190, China. University of
Chinese Academy of Science, Beijing, 100049, China. Tsien Excellence in
Engineering Program, School of Aerospace Engineering, Tsinghua University,
Beijing, 100084, China.

Pei<EMAIL_ADDRESS>
Brain criticality hypothesis is one of the most focused and controversial
topics in neuroscience and biophysics. This research develops a unified
framework to reformulate the physics theories of four basic types of brain
criticality, i.e., ordinary criticality (OC), quasi-criticality (qC), self-
organized criticality (SOC), and self-organized quasi-criticality (SOqC), into
more accessible and neuroscience-related forms. For the statistic techniques
used to validate brain criticality hypothesis, we also present comprehensive
explanations of them, summarize their error-prone details, and suggest
possible solutions. This framework may help resolve potential controversies in
studying brain criticality hypothesis, especially those arising from the
misconceptions about the theoretical foundations of brain criticality.

##  1 Introduction

Neuroscience is dawning upon revealing physics foundations of the brain Abbott
(2008). Ever since the 1970s, the term neurophysics has been suggested as a
term to indicate the essential role of physics in understanding the brain
Scott Alwyn (1977). More recently, substantial progress has been accomplished
in studying brain connectivity and brain functions with statistical physics
theories Lynn Bassett (2019).

For brain connectivity, physics provides insights for its emergence,
organization, and evolution. Random graphs Betzel Bassett (20171); Betzel .
(2016), percolation Guo . (2021); Breskin . (2006), and other physics theories
of correlated systems Haimovici . (2013); Wolf (2005) are applied to reveal
the underlying mechanisms accounting for the origins of brain network
properties. Complex network theories act as the foundation of characterizing
brain connectivity organizational features (e.g., community Betzel Bassett
(20172); Betzel . (2018); Khambhati . (2018), hub Gong . (2009); Deco .
(2015), and small-world Deco . (2015); Bullmore Sporns (2012) structures) and
embedding attributes into physical space Bassett . (2010); Kaiser Hilgetag
(2006). Network evolution driven by neural plasticity helps to explain the
dynamics of brain connectivity structures during information processing Del
Pozo . (2021); Montague . (1996); Song . (2000); Robert Vignoud (2021); Galván
(2010). For brain functions, physics presents possible explanations for the
origin of information processing capacities from collective neural activities.
From single neuron dynamics models Gerstner . (2014), stochastic network
models of neural populations and circuits Tian Sun (2021); Tian . (2021),
mean-field neural mass models of brain regions Touboul . (2011); David Friston
(2003), eventually to models of entire brain networks Schneidman . (2006);
Hopfield (1982), important efforts have been devoted to characterize
information-processing-related neural dynamics across different scales.
Networks with memory capacities (e.g., Hopfield networks Tyulmankov . (2021)),
which are equivalent to Ising models under specific conditions Lynn Bassett
(2019), have been applied to study neural information storage and recall
Krotov Hopfield (2020); Haldeman Beggs (2005), adaptation to environment
changes Shew . (2015), information transmission optimization Beggs Plenz
(2003), dynamic range maximization Kinouchi Copelli (2006); Shew . (2009), and
neural computation power Bertschinger Natschläger (2004). These models are
further related to maximum entropy models (e.g., specific fine-tuned Ising
models) that predict long-range correlations observed among neurons Schneidman
. (2006); Ganmor . (2011). Moreover, general theories of free-energy principle
Friston (2010, 2009); Guevara (2021) and information thermodynamics Tian Sun
(2022); Capolupo . (2013); Collell Fauquet (2015); Sartori . (2014) are
suggested as the unified foundations of perception, action, and learning in
the brain.

If one needs to specify one of the most focused and controversial topics among
all the works mentioned above, brain criticality may be a potential candidate
Beggs Timme (2012). The hypothesis of the critical brain has received
increasing attention in recent decades, serving as a possible mechanism
underlying various intriguing but elusive phenomena in the brain. In light of
our limited understanding of the complex nature of collective neural dynamics,
these phenomena include, to name a few, efficient transitions between cortical
states Fontenele . (2019), maximal dynamic ranges of neural responses Kinouchi
Copelli (2006); Shew . (2009); Antonopoulos (2016); Gautam . (2015), optimized
information transmission and representation Shew . (2011); X. Li Small (2012);
Antonopoulos (2016), and numerous other issues concerning brain functions that
we have mentioned above. One can see Shew Plenz (2013); Chialvo (2010); Beggs
(2007); Hesse Gross (2014); Cocchi . (2017) for systematic reviews of the
diverse function advantages implied by brain criticality and their
experimental demonstrations. From a Darwinian perspective, one potential
reason for the brain to feature criticality lay in that the most informative
parts of external world principally occur at a borderline between purely
ordered and purely disordered states (information would be trivial in a purely
ordered world while it would be incomprehensible in a purely disordered
world). Becoming critical may be a potential way for the brain to adapt to the
complex world, where non-trivial information has a finite opportunity to occur
Bak (2013); Chialvo (2010). To date, generic features of a critical brain with
the characteristics discussed above, such as divergent correlation length,
neuronal avalanches with power-law behaviours, and long-range correlations on
the microscopic scale (e.g., neural populations), have been extensively
observed in mathematical models in conjunction with experimental data (e.g.,
Fosque . (2021); Shew . (2009); Beggs Plenz (2003); Dalla Porta Copelli
(2019); Poil . (2012); Hardstone . (2014); Petermann . (2009); Shriki .
(2013); Tagliazucchi . (2012); Poil . (2008); Gireesh Plenz (2008); Tkačik .
(2015); G. Scott . (2014); Ponce-Alvarez . (2018)).

Our work does not aim at repeatedly reviewing experimental advances concerning
brain criticality and its biological significance, given that they have been
comprehensively summarized by existing reviews Shew Plenz (2013); Chialvo
(2010); Beggs (2007); Hesse Gross (2014); Cocchi . (2017); Munoz (2018). On
the contrary, our motivation is to present a systematic and accessible review
of the theoretical methods applied to achieve these advances, which have not
received necessary attention yet. These theoretical foundations are initially
thought to be incomprehensible and irrelevant to neuroscience. However,
practice suggests that omitting these physical and mathematical backgrounds
does not significantly improve the accessibility of studies on brain
criticality. Instead, the lack of detailed explanations of theoretical
foundations has frequently misled neuroscientists, leading to diverse
confusions about the precise meaning, identification criteria, and biological
corollaries of brain criticality. As a result, criticality, an analytic
statistical physics theory with solid foundations, unnecessarily becomes an
elusive black box for neuroscientists. To address this issue, we use the
terminology of neuroscience to present a self-contained framework of brain
criticality, reviewing and reformulating (1) physical theories that
analytically derive brain criticality and (2) statistic techniques that
computationally estimate brain criticality from empirical data. Given the
frequent misunderstanding of neural avalanches, our discussions primarily
focus on brain criticality analysis on the microscopic scale of the brain. The
objectives guiding through our review are tripartite: (1) explaining why brain
criticality matters in the brain, (2) understanding what is brain criticality
and what it conveys about the brain, and (3) confirming how to identify
potential brain criticality and ensure the validity of analyses.

##  2 Brain criticality: general concepts

###  2.1 Overview of brain criticality

Brain criticality frequently confuses neuroscientists since too many distinct
phenomena are studied under this name without being properly classified. In
this review, brain criticality refers to a family of critical processes in
neural dynamics where erratic fluctuations appear to reduce dynamic stability.
To present a systematic classification framework, we discuss three fundamental
perspectives concerning brain criticality. Table 1 provides all the necessary
glossaries in comprehensible forms.

####  2.1.1 Being non-equilibrium

First, the brain, similar to other biological systems, generally exhibits
temporal evolution from initial states that are far away from equilibrium Lynn
. (2021); Gnesotto . (2018). These departures from equilibrium arise due to
diverse endogenous causes Perl . (2021); Gnesotto . (2018) to break the
detailed balance to support consciousness, sensing, and adaptation Perl .
(2021); Lynn . (2021). Therefore, potential critical phenomena underlying
neural dynamics, at least in most neural dynamics models and empirical data
sets, are basically non-equilibrium and can not be characterized by
equilibrium statistic mechanics. In Fig. 1A, we illustrate the difference
between equilibrium and non-equilibrium dynamics.

Table 1: Key concepts in describing brain criticality Concept | Meaning
---|---
Equilibrium |  | A case where the system maximizes entropy and conserves energy simultaneously. The stationary
---
probability distribution $\mathcal{P}_{eq}\left(\cdot\right)$ of system states
of a system at equilibrium is the Boltzmann distribution.
At equilibrium, the transition dynamics between system states $c$ and
$c^{\prime}$ satisfies the detailed balance
condition $\mathcal{P}_{eq}\left(c\right)\mathcal{W}\left(c\rightarrow
c^{\prime}\right)=\mathcal{P}_{eq}\left(c^{\prime}\right)\mathcal{W}\left(c^{\prime}\rightarrow
c\right)$, where $\mathcal{W}\left(\cdot\rightarrow\cdot\right)$ denotes the
transition probability.
Non-equilibrium |  | A case where the system is out of equilibrium because the transition dynamics between system states
---
breaks the detailed balance condition. In other words, the transition dynamics
between states becomes
directional rather than symmetric.
Self-organization |  | A process where the internal complexity of a system increases without being tuned by any external
---
mechanism. All potentially emergent properties are created by endogenous
feedback processes or other
internal factors inside the system.
Criticality |  | A kind of phenomena where the systems is generally close to specific critical points separating between
---
multiple system states. Small disturbances are sufficient to make the system
experience dramatic and
sharp transitions between system states.
Quasi-criticality |  | A kind of phenomena where all statistical physics relations required by criticality are principally adhered
---
by the system but slight and inconstant deviations from perfect criticality
can be seen on the actual values
of characteristic variables. These deviations robustly exist and are generally
independent of data noises.
Sub-criticality |  | A kind of system states below criticality. They occur when the order parameter (i.e., the macroscopic
---
observable used to describe system states) remains at zero even with the
addition of derives, corresponding
to disordered system dynamics.
Super-criticality |  | A kind of system states above criticality. They occur when the order parameter is positive, corresponding
---
to ordered system dynamics.

####  2.1.2 Fine tuning versus self-organization

Second, there exist two types of general mechanisms underlying the existence
of brain criticality. One type of mechanisms either arise from the external
manipulations outside the brain (e.g., researchers manipulate the tonic
dopamine D1-receptor stimulation Stewart Plenz (2006, 2008) or adjust network
topology Kaiser Hilgetag (2010); S. Wang Zhou (2012); Rubinov . (2011)) or
belong to the top-down biological processes that globally function on neural
dynamics inside the brain (e.g., anesthesia effects Fontenele . (2019);
Ribeiro . (2010); Hahn . (2017) as well as sleep restoration effects Meisel .
(2013)). Neural dynamics is passively fine tuned towards or away from ordinary
criticality (OC) by these exogenous mechanisms, similar to ordinary critical
phenomena that require the fine tuning of order parameters.

Another type of mechanisms includes all endogenous factors of neural dynamics
(e.g., neural plasticity mechanisms such as spike-timing dependent synaptic
plasticity Shin Kim (2006); Meisel Gross (2009); Effenberger . (2015), short-
term synaptic plasticity Levina . (2007, 2009), retro-synaptic signals
Hernandez-Urbina Herrmann (2017) and Hebbian rules De Arcangelis . (2006); de
Arcangelis Herrmann (2010)), which locally function on neural dynamics as
drive and dissipation components. The interactions between these components
naturally form feedback control loops to support the self-organization of
neural dynamics towards the critical point Beggs (2007); Chialvo (2010). This
spontaneously emerged brain criticality, distinct from ordinary critical
phenomena, is conjectured as a kind of self-organized criticality (SOC)
Chialvo (2010). In Fig. 1B, we present conceptual illustrations of ordinary
criticality and self-organized criticality in the brain.

####  2.1.3 Standard versus non-standard

Third, brain criticality frequently occurs in non-standard forms due to
stimulus derives or endogenous factors. On the one hand, slight and inconstant
deviations from perfect brain criticality can be seen on the actual values of
characteristic variables, differentiating the characterized phenomena from the
standard criticality Williams-García . (2014); Fosque . (2021). On the other
hand, all statistical physics relations required by perfect brain criticality
are still adhered by these actual characteristic variables, distinguishing the
brain from being non-critical Williams-García . (2014); Fosque . (2021).

Figure 1: Please see the caption on the next page.

Figure 1: Conceptual illustrations of brain criticality. A, Difference between
equilibrium and non-equilibrium dynamics in a three-state brain (upper
parallel). Brain states are characterized by three system components. We
illustrate an instance of non-equilibrium dynamics between these states
(bottom parallel). B, Fine tuning with exogenous mechanisms (represented by
animated hands) makes the brain evolve from a non-critical state (upper left)
to the critical state (upper right). Endogenous mechanisms enable the brain to
self-organize from a non-critical state (bottom left) to the critical state
(bottom right). C Increasing stimulus intensity enlarges the quasi-critical
region around the perfect critical point in a quasi-critical system. D The
approaching process to a critical point in a self-organized quasi-critical
system consists of two stages. In the first stage, the brain self-organizes
from a non-critical state to a quasi-critical region based on certain
endogenous mechanisms. In the second stage, additional exogenous mechanisms
are necessary to fine tune the brain to the critical point. Otherwise, the
brain just hovers within the quasi-critical region. E The difference between
four types of brain criticality from the perspective of susceptibility. For
standard brain criticality (e.g., ordinary criticality and self-organized
criticality), susceptibility becomes divergent (i.e., infinite) at the
critical point. For non-standard brain criticality (e.g., quasi-criticality
and self-organized quasi-criticality), susceptibility is always non-divergent
(i.e., finite). The quasi-critical region is defined as a set of all control
parameters where susceptibility values are no less than a specific threshold
(e.g., half-maximum value). F The commonness and difference between four types
of brain criticality.

For ordinary criticality, its non-standard form is referred to as quasi-
criticality (qC) Williams-García . (2014); Fosque . (2021). Diverse mechanisms
can force the brain to depart from perfect ordinary criticality and exhibit
quasi-critical neural dynamics, among which, stimulus derive may be the most
common one Williams-García . (2014); Fosque . (2021). In general, sufficiently
strong stimulus drives can capture or even govern neural dynamics. Similar to
the situation where external inputs suppress irregular neural dynamics
Molgedey . (1992), the stimuli that are too strong may evoke intense but less
changeable neural dynamics to make the brain depart from the perfect critical
point Williams-García . (2014); Fosque . (2021). Let us take the qC phenomenon
introduced by Williams-García . (2014); Fosque . (2021) as an instance. Under
specific conditions, the actual brain state may be close to a Widom line in
the three-dimensional space defined by the stimulus intensity $\upsilon$,
refractory period length $\tau$, and branching ratio $\kappa$ (i.e., the time-
dependent average number of subsequent neural activities caused by a single
neuron activation event Haldeman Beggs (2005)). The Widom line is a line of
all the combinations of $\left(\upsilon,\tau,\kappa\right)$ where the
susceptibility of neural dynamics is maximized Williams-García . (2014);
Fosque . (2021). The susceptibility is defined by $\lim_{x\rightarrow
0}\frac{\partial y}{\partial x}$, where $y$ is the neural dynamics state and
$x$ denotes a factor that affects $y$. In general, one can understand
susceptibility as the degree to which fluctuations in the state of each neuron
can propagate to neighbored neurons Williams-García . (2014). Being close to
the Widom line suggests the existence of quasi-criticality in the brain.
Moving along the Widom line as the stimulus intensity increases, the
susceptibility of neural dynamics decreases, and the branching ratio at
maximal susceptibility will decrease as well Williams-García . (2014); Fosque
. (2021). Significant deviations from the Widom line suggest non-criticality
(i.e., the sub-criticality where neural dynamics is disordered and the super-
criticality where neural dynamics is ordered Williams-García . (2014)). In
Fig. 1C, we conceptually illustrate how stimuli imply qC in the brain. In Fig.
2D, the qC phenomenon in Williams-García . (2014); Fosque . (2021) is shown in
details.

As for self-organized criticality (SOC), its non-standard form is defined
according to statistical physics criteria. Perfect self-organized criticality
only exists in conserved neural dynamics (e.g., see integrate-and-fire neurons
analyzed by Levina . (2007)), where system energy (i.e., neural activities)
either conserves within the system and only dissipates at the system boundary,
or dissipates inside the system (i.e., bulk dissipation) with a dissipation
rate vanishing in the system size limit Malcai . (2006). Under more general
conditions where neural dynamics is not conserved (e.g., see leaky integrate-
and-fire neurons analyzed by Millman . (2010); Levina . (2007); Rubinov .
(2011); Stepp . (2015), where neural dynamics dissipates within the system due
to voltage leak), perfect self-organized criticality can be broken by any rate
of bulk dissipation Bonachela Munoz (2009); Bonachela . (2010); Buendía, di
Santo . (2020); de Andrade Costa . (2015). Stronger bulk dissipation implies
larger deviations from perfect self-organized criticality De Arcangelis .
(2006). Consequently, the self-organization process of non-conserved neural
dynamics only make the brain hover around the critical point. Any further
closeness towards the critical point requires the fine tuning of order
parameter by additional exogenous mechanisms, which is different from pure
self-organized criticality Bonachela Munoz (2009); Bonachela . (2010);
Buendía, di Santo . (2020); de Andrade Costa . (2015). This non-conserved
self-organization process is termed as self-organized quasi-criticality (SOqC)
Bonachela Munoz (2009). Similar to SOC in conserved dynamics, neural
plasticity mechanisms, such as spike-timing-dependent synaptic plasticity
Rubinov . (2011), Hebbian rules De Arcangelis . (2006), short-term synaptic
depression in conjunction with spike-dependent threshold increase Girardi-
Schappo . (2021), and inhibitory plasticity in conjunction with network
topology Ma . (2019), can serve as underlying self-organization mechanisms of
SOqC. Because purely conserved neural dynamics is relatively rare in empirical
data (e.g., neural dynamics is conserved for integrate-and-fire neurons Levina
. (2007) and leaky integrate-and-fire neurons whose pre-synaptic inputs are
exactly equal to the sum of voltage leak and potential costs during neural
spiking Bonachela . (2010)), we suggest that SOqC may be more common in the
brain than SOC. In Fig. 1D, we present conceptual instances of the two-stage
approaching process towards the critical point in the brain with SOqC.

####  2.1.4 Classification of brain criticality

The above discussion has presented a classification framework of brain
criticality, i.e., ordinary criticality (OC), quasi-criticality (qC), self-
organized criticality (SOC), and self-organized quasi-criticality (SOqC). In
Fig. 1E, we compare between these four types of brain criticality in term of
susceptibility. In general, susceptibility diverges at the critical point in a
brain with standard criticality (e.g., OC and SOC) while it does not diverge
in the quasi-critical region of a brain with non-standard criticality (e.g.,
qC and SOqC). In Fig. 1F, we summarize the commonness and difference between
these four types of brain criticality discussed in our review. From a
neuroscience perspective, a brain with critical neural dynamics is expected to
be near the critical point and prepared for tremendous changes in cortical
states during a short duration. This intriguing property coincides with the
experimentally observed efficient transitions between cortical states (e.g.,
Jercog . (2017); Holcman Tsodyks (2006); Reimer . (2014); H. Lee . (2020);
Cardin (2019)) and, therefore, interests researchers for the potential
existence of brain criticality. The importance of identifying brain
criticality in neural dynamics is beyond brain criticality itself because it
implies an opportunity to explain and predict brain function characteristics
by various statistical physics theories built on non-equilibrium criticality.

###  2.2 Neural avalanches and their phases

To identify potential non-equilibrium criticality in the brain, researchers
actually characterize neural dynamics as a physical system with absorbing
states and avalanche behaviors Hinrichsen (2000); Lübeck (2004); Larremore .
(2012). In general, one need to consider the propagation of neural dynamics
where neurons are either activated (“on” state) or silent (”off” state) Dalla
Porta Copelli (2019). An silent neuron may be activated with a probability
defined by the number of activated pre-synaptic neurons and the coupling
strength $\theta$ among neurons (e.g., neural correlation Franke . (2016)). An
activated neuron spontaneously becomes silent at a constant rate (e.g., after
the refractory period Squire . (2012); Kinouchi Copelli (2006)). These
definitions naturally support to distinguish between different phases of
neural dynamics. Here we review two kinds of phase partition that are active
in neuroscience.

####  2.2.1 Absorbing versus active

The first group of phases are absorbing and active phases Larremore . (2012).
The absorbing phase refers to cases where couplings between neurons are weak
and all neurons eventually become silent (neural dynamics vanishes). Once a
neural dynamics process vanishes, it can not reappear by itself. The brain
requires new drives (e.g., neurons activated spontaneously or by stimuli) to
trigger new neural dynamics. The active phase, on the other hand, correspond
to cases where “on” state propagates among neurons with strong couplings,
leading to stable self-sustained neural dynamics (e.g., non-zero time- and
ensemble-averaged density of active neurons in the brain). In Fig. 2A, we show
conceptual instances of neural avalanches, self-sustained neural dynamics, and
vanished neural dynamics. Denoting $\rho\left(t\right)$ as the density of
active neurons at moment $t$, we can simply represent the absorbing (Eq. (1))
and active (Eq. (2)) phases of a neural dynamics process triggered by an
active neuron at moment $0$ as

$\displaystyle\rho\left(t\right)=0,\;\exists t>0,$ (1)
$\displaystyle\rho\left(t\right)>0,\;\forall t>0.$ (2)

####  2.2.2 Synchronous versus asynchronous

The second group of phases are synchronous and asynchronous phases Di Santo .
(2018); Fontenele . (2019); Girardi-Schappo . (2021). As their names suggest,
these two phases correspond to the situations where synchronization emerges or
disappears in neural activities, respectively. Synchronization refers to the
cases where “on” states appear in an oscillatory, although not strictly
periodic, manner. To quantify its potential existence, we can measure the
variability of neural dynamics using the coefficient of variation (CV) Di
Santo . (2018); Fontenele . (2019); Girardi-Schappo . (2021) or the Kuramoto
order parameter Arenas . (2008); Acebrón . (2005). CV can be defined from
diverse perspectives, yet the most common definition is the ratio between the
standard deviation and the mean of the inter-spike interval length Di Santo .
(2018); Fontenele . (2019); Girardi-Schappo . (2021). A higher value of CV
implies the reduction of synchronization. For most neural dynamics data, an
empirical choice of the CV threshold that separates between synchronous and
asynchronous phases may be $\simeq 1$ Fontenele . (2019) or
$\simeq\frac{3}{2}$ Fontenele . (2019). The Kuramoto order parameter
$\omega\in\left[0,1\right]$ measures the coherent degree of neural dynamics
based on the Kuramoto model of oscillators (see Arenas . (2008); Acebrón .
(2005) for detailed definitions). Perfect synchronization emerges when
$\omega=1$ and vanishes when $\omega=0$ Arenas . (2008); Acebrón . (2005).

####  2.2.3 Critical point or quasi-critical region

The boundary between these two phases is the critical point, at which the
brain is on the edge of exhibiting self-sustained (for absorbing and active
phases) or synchronous (for synchronous and asynchronous phases) neural
dynamics. Perturbations (e.g., the propagation of “on” state among neurons) to
the absorbing or asynchronous phase do not have characteristic lifetime and
size. These perturbations, referred to as neural avalanches, are expected to
exhibit power-law properties in their lifetime (time difference between the
first and last activation of neurons in between complete quiescent epochs) and
size (number of active neurons along with the excursion) distributions
Hinrichsen (2000); Lübeck (2004); Larremore . (2012); Hesse Gross (2014). In
general, the emergence of neural avalanches implies the slowing down of neural
dynamics, i.e., the brain state recovery process towards the baseline state
after fluctuations changes from fast (exponential) to slow (power-law) Hesse
Gross (2014); Cocchi . (2017). The dynamic stability of neural dynamics is
limited by the slow recovery and, therefore, can not robustly counteract
perturbations. Consequently, small perturbations initiated on the microscopic
scale may still make the brain change sharply on the macroscopic scale Hesse
Gross (2014); Cocchi . (2017). In Fig. 2B, we conceptually illustrate how the
recovery process slows down when the brain is close to the critical point or
the quasi-critical region.

###  2.3 General relations between neural avalanches and brain criticality

The relation between neural avalanches and brain criticality is frequently
neglected or misunderstood. Neural avalanche data alone is not sufficient to
determine the concrete type of brain criticality (i.e., OC, qC, SOC, and SOqC)
unless additional information about the mechanisms underlying neural avalanche
emergence is provided (e.g., if neural dynamics is conserved or self-
organizing). To explore a concrete type of brain criticality, researchers need
to explicitly present its definition depending on different control parameters
(e.g., the balance between excitatory and inhibitory neurons in CROS models
Poil . (2012); Hardstone . (2014)) and order parameters (e.g., active neuron
density and synchronous degree Dalla Porta Copelli (2019)). A brain
criticality hypothesis without strict definitions of control and order
parameters is not informative Cocchi . (2017); Girardi-Schappo (2021). To
present conceptual instances, we illustrate four possible critical phenomena
in Fig. 2, each of which corresponds to a concrete brain criticality type.

####  2.3.1 Instance of ordinary criticality

To produce ordinary criticality (OC), we can control neural dynamics and
manipulate $\langle\theta\rangle$, the expectation of coupling strength
$\theta$ among all neurons (e.g., averaged neural correlation), by some top-
down and global biological effects. These effects, for instance, may be
anesthesia effects (e.g., by ketamine-xylazine Ribeiro . (2010) and isoflurane
Hahn . (2017)) or sleep restoration effects Meisel . (2013). We use the
Kuramoto order parameter $\omega$ Arenas . (2008); Acebrón . (2005) as the
order parameter to define synchronous and asynchronous phases Di Santo .
(2018); Fontenele . (2019). As $\langle\theta\rangle$ increases, we may see
transitions from asynchronous to synchronous phase in some situations (see a
similar instance in Villegas . (2014)). One can see Fig. 2C for conceptual
illustrations.

Figure 2: Please see the caption on the next page.

Figure 2: Conceptual illustrations of the relations between neural avalanches
and brain criticality. A, Instances of neural avalanche, self-sustained neural
dynamics, and vanished neural dynamics. B, The recovery processes of brain
states after the same perturbation in the space of absorbing and active phases
(upper parallel) and the space of synchronous and asynchronous phases (bottom
parallel). The recovery processes after perturbations are relatively fast when
the brain is far from the critical point or the quasi-critical region, These
recovery processes slow down when the brain is close to the critical point or
the quasi-critical region. C The conceptual illustrations of neural dynamics
when the brain state is asynchronous, synchronous, or at ordinary criticality.
D Without stimuli, there initially exist disordered (gray), ordered but
asynchronous (light blue), synchronous (green) phases in the phase space of
the brain. Stimulus inputs imply quasi-criticality in the brain. An increasing
stimulus intensity enlarges the quasi-critical region (purple) around the
Widom line (purple dashed line). E The conceptual illustrations of how
endogenous mechanisms in conserved neural dynamics can function as drive or
dissipation terms to create self-organized criticality between absorbing and
active phases in the brain. F In the self-organized quasi-critical brain,
endogenous mechanisms in non-conserved neural dynamics only support the self-
organization towards a quasi-critical region between asynchronous and
synchronous phases. Extra exogenous mechanisms are required to fine tune the
brain towards the critical point.

####  2.3.2 Instance of quasi-criticality

To produce quasi-criticality (qC), we can manipulate refractory period length
$\tau$, branching ratio $\kappa$, and stimulus intensity $\upsilon$ as control
parameters (e.g., control $\tau$ and $\kappa$ by pharmacological perfusion or
ionic concentration adjustment Shew . (2011); Chiappalone . (2003)). There
exist a disordered phase (sub-critical), an ordered but asynchronous phase
(super-critical), and a synchronous (quasi-periodic) phase in the space of
$\left(\upsilon,\tau,\kappa\right)$ Williams-García . (2014); Fosque . (2021).
These phases can be characterized by specific order parameters related to
synchronization. As $\upsilon$ increases, a qC phenomenon emerges in the
space, where the quasi-critical region is defined by all combinations of
$\left(\upsilon,\tau,\kappa\right)$ whose susceptibility values are at least
half-maximum. Cross-over behaviours (i.e., a generalization of phase
transition with finite susceptibility) emerge when the quasi-critical region
has overlaps with at least two phases Williams-García . (2014); Fosque .
(2021). In Fig. 2D, we show this qC phenomenon in details.

####  2.3.3 Instance of self-organized criticality

To study self-organized criticality (SOC), we consider the conserved neural
dynamics generated by integrate-and-fire neurons Levina . (2007). The order
parameter is active neuron density $\rho$, whose dynamics is controlled by
parameter $\langle\theta\rangle_{A}$, the averaged coupling strength $\theta$
between activated neurons and their post-synaptic neurons (here $A$ denotes
the set of activated neurons). In specific cases, the considered neural
dynamics may self-organize to the critical point under the joint effects of
excitatory and inhibitory neurons, neural spiking processes (activation and
silence), as well as neural plasticity. In Fig. 2E, we conceptually illustrate
a case where these endogenous mechanisms enable the brain to self-organize to
the criticality between absorbing and active phases.

####  2.3.4 Instance of self-organized quasi-criticality

To analyze self-organized quasi-criticality (SOqC), we consider the non-
conserved neural dynamics affected by two homeostatic adaptation processes,
i.e., the short-term depression of inhibition and the spike-dependent
threshold increase. These processes are controlled by $\widehat{y}$, the
maximum inhibitory coupling strength, as well as $\tau_{x}$ and $\tau_{y}$,
the decay time scales of neural activation threshold increase and synaptic
depression. These control parameters affect neural activation threshold $x$
and inhibition strength $y$ to shape neural dynamics states (e.g., the active
neuron density $\rho$). With appropriate $x$, $y$, and $\rho$, neural
avalanches with power-law behaviours will occur to indicate the criticality
between an asynchronous phase (stochastic oscillations) and a synchronous
phase (periodic oscillations). According to Girardi-Schappo . (2021), $x$ and
$\rho$ self-organize to their appropriate values through quasi-critical
fluctuations under biologically reasonable conditions (i.e., $\tau_{x}\gg 1$)
while $y$ hovers around the expected value. Additional fine tuning of $y$
based on exogenous mechanisms are necessary to place neural dynamics at the
perfect criticality. Meanwhile, synaptic homeostasis is discovered as
constantly canceled by the variation of the activation threshold, impeding
neural dynamics from self-organizing to perfect criticality. In Fig. 2F, we
conceptually illustrate the defined SOqC phenomenon in a similar manner of
Fig. 2D and Fig. 2E. As for the precise description of quasi-critical
fluctuations, one can see Girardi-Schappo . (2021) for details.

To this point, we have conceptually introduced the phenomenological properties
of brain criticality. To verify the hypothetical brain criticality, one need
to learn about analytic brain criticality theories and the properties of
neural avalanche predicted by them. Below, we present accessible expositions
of these theoretical foundations.

##  3 Brain criticality: physical theories

###  3.1 Mean field and stochastic field theories of brain criticality

One of the main challenges faced by neuroscientists in studying ordinary
criticality (OC), quasi-criticality (qC), self-organized criticality (SOC),
and self-organized quasi-criticality (SOqC) is how to understand their
theoretical relations Girardi-Schappo (2021). Overcoming this challenge is
crucial for understanding why we can verify the existence of different types
of brain criticality with certain theoretical tools. To present a concise and
thorough review, we first focus on brain criticality between absorbing and
active phases, where we generalize the idea in Bonachela Munoz (2009);
Buendía, Di Santo . (2020) to present a possible framework for unification.

####  3.1.1 Langevin formulation of ordinary criticality

In general, brain criticality in the space of absorbing and active phases are
related to directed percolation Dalla Porta Copelli (2019), a universality
class of continuous phase transitions into absorbing states Hinrichsen (2000);
Lübeck (2004). Here a universality class can be understood as the set of all
systems with the same scaling properties Hinrichsen (2000); Lübeck (2004);
Sethna . (2001). Directed percolation theory initially covers OC phenomena
Hinrichsen (2000); Lübeck (2004). Let us begin with a variant of the classic
Reggeon field theory, the simplest description of absorbing phase transitions
Henkel . (2008). The Langevin equation of the activity neuron field
$\rho\left(\vec{x},t\right)$ is defined as

$\displaystyle\frac{\partial}{\partial t}\rho\left(\vec{x},t\right)$
$\displaystyle=\left(a+b\nu\left(\vec{x},t\right)\right)\rho\left(\vec{x},t\right)-c\rho^{2}\left(\vec{x},t\right)+d\nabla^{2}\rho\left(\vec{x},t\right)+e\sqrt{\rho\left(\vec{x},t\right)}\sigma\left(\vec{x},t\right),$
(3) $\displaystyle\frac{\partial}{\partial t}\nu\left(\vec{x},t\right)$
$\displaystyle=\nabla^{2}\nu\left(\vec{x},t\right)+f\left(\vec{x},t\right)-g\left(\vec{x},t\right)\rho\left(\vec{x},t\right),$
(4)

where $\vec{x}$ represents spatial coordinates, $a\in\mathbb{R}$,
$b\in\left(0,\infty\right)$, $c\in\left(0,\infty\right)$, $d\in\mathbb{R}$ is
the diffusion factor, and $e\in\mathbb{R}$ is the noise factor. Function
$\sigma\left(\cdot,\cdot\right)$ defines a zero-mean Gaussian noise with a
spatio-temporal correlation
$\langle\rho\left(\vec{x},t\right)\rho\left(\vec{x}^{\prime},t^{\prime}\right)\rangle=\delta\left(\vec{x}-\vec{x}^{\prime}\right)\delta\left(t-t^{\prime}\right)$,
where $\delta\left(\cdot\right)$ is the delta function. In general,
$\sigma\left(\cdot,\cdot\right)$ reflects the collective fluctuations in
neural activities that vanish in the absorbing phase
$\rho\left(\vec{x},t\right)=0$ under the effects of factor
$\sqrt{\rho\left(\vec{x},t\right)}$. The term
$\nabla^{2}\rho\left(\vec{x},t\right)$ reflects the propagation of neural
dynamics. The function $\nu\left(\vec{x},t\right)$ defines the energy (i.e.,
membrane potential) that propagates according to
$\nabla^{2}\nu\left(\vec{x},t\right)$, increases with external drives
$f\left(\vec{x},t\right)$, and decreases with bulk dissipation
$g\left(\vec{x},t\right)$. Please note that $\rho\left(\vec{x},t\right)\geq 0$
and $\nu\left(\vec{x},t\right)\geq 0$ always hold. The initial active neuron
density and energy are assumed as non-zero. It is clear that
$a+b\nu\left(\vec{x},t\right)<0$ makes the neural dynamics eventually vanish
(i.e., absorbing phase) while $a+b\nu\left(\vec{x},t\right)>0$ does not (i.e.,
active phase). Therefore, we can fine tune the control parameter
$\nu\left(\vec{x},t\right)$ to make the brain exhibit OC dynamics at
$a+b\nu_{c}\left(\vec{x},t\right)=0$, a critical point defined by $\nu_{c}$.
The fine tuning relies on manipulating $f\left(\vec{x},t\right)$ and
$g\left(\vec{x},t\right)$ by exogenous mechanisms.

####  3.1.2 Langevin formulation of quasi-criticality

Then we turn to analyzing qC, whose mean field approximation is initially
derived based on the cortical branching model Williams-García . (2014); Fosque
. (2021). A cortical branching model with no stimulus input belongs to
directed percolation universality class according to the Janssen-Grassberger
conjecture Williams-García . (2014). Non-zero stimulus inputs make the
cortical branching model depart from directed percolation universality class
to create qC Williams-García . (2014). Nevertheless, the above mean field
theory is defined in the space of synchronous and asynchronous phases. To
derive a qC phenomenon between absorbing and active phases, we can
provisionally analyze a mean field approximation of Eqs. (3-4)

$\displaystyle\frac{\partial}{\partial t}\rho\left(\vec{x},t\right)$
$\displaystyle=\left(a+b\nu\left(\vec{x},t\right)\right)\rho\left(\vec{x},t\right)-c\rho^{2}\left(\vec{x},t\right),$
(5) $\displaystyle\frac{\partial}{\partial t}\nu\left(\vec{x},t\right)$
$\displaystyle=f\left(\vec{x},t\right)-g\left(\vec{x},t\right)\rho\left(\vec{x},t\right),$
(6)

where $\nabla^{2}\rho\left(\vec{x},t\right)$,
$\nabla^{2}\nu\left(\vec{x},t\right)$, and $\sigma\left(\vec{x},t\right)$ in
Eqs. (3-4) are neglected under the mean field assumption. We consider the
cases where stimulus inputs vanish, i.e., $f\left(\vec{x},t\right)\equiv 0$.
The critical point between active and absorbing phase becomes
$\nu_{c}=-\frac{a}{b}$. The steady state solutions of Eqs. (5-6) are

$\displaystyle\rho\left(\vec{x},t\right)$ $\displaystyle=0,$ (7)
$\displaystyle\nu\left(\vec{x},t\right)$
$\displaystyle=r\in\left(0,\infty\right),$ (8)

respectively. Therefore, OC is one of the steady states of neural dynamics
when there is no stimulus. In the cases where stimulus inputs become
increasingly strong, there exists no steady state solution of Eqs. (5-6)
unless $\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow
r\in\left(0,\infty\right)$. If
$\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow
r\in\left(0,\infty\right)$ holds, we can derive

$\displaystyle\rho\left(\vec{x},t\right)$
$\displaystyle=\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow
r,$ (9) $\displaystyle\nu\left(\vec{x},t\right)$
$\displaystyle=\frac{1}{b}\left(c\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}-a\right)\rightarrow\frac{1}{b}\left(cr-a\right).$
(10)

Because the critical point $\nu_{c}=-\frac{a}{b}$ is not necessarily a steady
state, it can be disturbed by diverse factors (e.g., by stimuli). Unless there
exist certain ideal exogenous mechanisms that persistently enlarge
$g\left(\vec{x},t\right)$ whenever $f\left(\vec{x},t\right)$ increases, the
fine tuning of neural dynamics can not cancel the effects of
$f\left(\vec{x},t\right)$. Consequently, the fine tuning process may only
enable the brain to reach a quasi-critical region where the susceptibility of
neural dynamics is relatively large. The initial OC vanishes and is replaced
by qC.

####  3.1.3 Langevin formulation of self-organized criticality

Although SOC is treated as a rather isolated concept after its first discovery
in statistical physics Bak . (1987), subsequent analyses demonstrate SOC as
relevant with ordinary continuous phase transitions into infinitely many
absorbing states Narayan Middleton (1994); Sornette . (1995); Dickman . (1998,
2000). Specifically, SOC models can be subdivided into two families, which we
refer to as external dynamics family (e.g., Bak-Sneppen model Bak Sneppen
(1993)) and conserved field family (e.g., sandpile models such as Manna model
Manna (1991) and Bak-Tang- Wiesenfeld model Bak . (1987)). The second family,
being the main theoretical source of studying SOC in neural dynamics,
corresponds to absorbing-state transitions since it can represent any system
with conserved local dynamics and continuous transitions to absorbing states
Dickman . (2000); Lübeck (2004). Although the universality class of the second
family should be precisely referred to as conserved directed percolation, the
explicit behaviours (e.g., avalanche exponents and scaling relations) of
conserved directed percolation are similar to those of directed percolation in
high-dimensional systems (e.g., neural dynamics) Buendía, Di Santo . (2020);
Bonachela Muñoz (2008); Bonachela Munoz (2009). Therefore, SOC and OC share
some identification criteria in practice. To understand the connections
between SOC and OC more precisely, we can consider the cases where
$f\left(\vec{x},t\right)\rightarrow 0$ such that
$\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow 0$ (i.e.,
infinite separation of timescales). The steady state solutions of Eqs. (5-6)
become

$\displaystyle\rho\left(\vec{x},t\right)$
$\displaystyle=\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow
0,$ (11) $\displaystyle\nu\left(\vec{x},t\right)$
$\displaystyle=\frac{1}{b}\left(c\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}-a\right)\rightarrow\nu_{c},$
(12)

respectively. Self-organization properties are reflected by the following
processes: if the brain is in the absorbing phase because neural dynamics
vanishes, i.e, $\rho\left(\vec{x},t\right)\rightarrow 0$, Eq. (6) becomes
$\frac{\partial}{\partial t}\nu\left(\vec{x},t\right)=f\left(\vec{x},t\right)$
to shift the brain towards the active phase; if the brain is in the active
phase, Eq. (6) becomes $\frac{\partial}{\partial
t}\nu\left(\vec{x},t\right)\simeq-g\left(\vec{x},t\right)\rho\left(\vec{x},t\right)$
to reduce neural dynamics since $f\left(\vec{x},t\right)\ll
g\left(\vec{x},t\right)$. These feedback-control loops drive the brain to the
critical point. One may be curious about why energy conservation, i.e.,
$g\left(\vec{x},t\right)\rightarrow 0$, is necessary for SOC since the above
derivations seem to be independent of $g\left(\vec{x},t\right)\rightarrow 0$.
Later we show that the absence of $g\left(\vec{x},t\right)\rightarrow 0$ in
Eq. (14) makes the active phase no longer exist. In other words, the non-
conserved energy implies a kind of continuous phase transition that does not
belong to conserved directed percolation or directed percolation when the
infinite separation of timescales is satisfied. Therefore, energy conservation
is necessary for SOC.

####  3.1.4 Langevin formulation of self-organized quasi-criticality

As for SOqC, non-zero bulk dissipation breaks the conservation law to generate
non-Markovian components in neural dynamics Bonachela Munoz (2009); Buendía,
Di Santo . (2020). In the ideal cases where the drive terms (e.g., stimulus
inputs) of a sufficiently large neural dynamics system occur at an arbitrarily
slow timescale (i.e., only occur in the interval between neural avalanches),
the brain exhibits pure dynamical percolation behaviours Buendía, Di Santo .
(2020). To understand this property, let us consider a variant of Eqs. (5-6)
where the dissipation term $g\left(\vec{x},t\right)$ is non-negligible

$\displaystyle\frac{\partial}{\partial t}\rho\left(\vec{x},t\right)$
$\displaystyle=\left(a+b\nu\left(\vec{x},t\right)\right)\rho\left(\vec{x},t\right)-c\rho^{2}\left(\vec{x},t\right),$
(13) $\displaystyle\frac{\partial}{\partial t}\nu\left(\vec{x},t\right)$
$\displaystyle=-g\left(\vec{x},t\right)\rho\left(\vec{x},t\right).$ (14)

By integrating Eq. (14) and plugging the integral into Eq. (13), we can derive

$\displaystyle\frac{\partial}{\partial t}\rho\left(\vec{x},t\right)$
$\displaystyle=\left(a+b\nu\left(\vec{x},0\right)\right)\rho\left(\vec{x},t\right)-c\rho^{2}\left(\vec{x},t\right)-b\rho\left(\vec{x},t\right)\int_{0}^{t}g\left(\vec{x},t\right)\rho\left(\vec{x},\tau\right)\mathsf{d}\tau,$
(15)

The non-Markovian term
$-b\rho\left(\vec{x},t\right)\int_{0}^{t}g\left(\vec{x},t\right)\rho\left(\vec{x},\tau\right)\mathsf{d}\tau$
in Eq. (15) makes the regions already visited by neural dynamics become more
unlikely to be activated Bonachela Munoz (2009); Buendía, Di Santo . (2020).
Therefore, the pure self-sustained active phase vanishes and is replaced by a
spreading phase, where local perturbations can transiently propagate across
the whole system without reaching a self-sustained state, and a non-spreading
phase, where local perturbations can never span the entire system Bonachela
Munoz (2009); Buendía, Di Santo . (2020). The phase transition and
corresponding critical point $\nu_{d}>\nu_{c}$ between spreading and non-
spreading phases belong to the universality class of dynamical percolation
rather than conserved directed percolation Bonachela Munoz (2009); Buendía, Di
Santo . (2020). The initial neural dynamics can be created by random shifts at
moment $0$ Bonachela Munoz (2009); Buendía, Di Santo . (2020)

$\displaystyle\rho\left(\vec{x}^{*},0\right)$
$\displaystyle\rightarrow\epsilon,$ (16)
$\displaystyle\nu\left(\vec{x}^{*},0\right)$
$\displaystyle\rightarrow\nu\left(\vec{x}^{*},0\right)+h\left(\vec{x}^{*},0\right),$
(17)

where $\vec{x}^{*}$ is a randomly selected coordinate, and function
$h\left(\cdot,\cdot\right)$ is a driving function of energy at moment $0$.
Every time a neural avalanche occurs after random shifts, the strong
dissipation term $g\left(\vec{x},t\right)$ pushes the brain towards the sub-
critical phase. Consequently, the brain can not exactly self-organize to the
perfect criticality. Instead, the brain just hovers around the critical point
$\nu_{d}$ to form a quasi-critical region, exhibiting finite fluctuations to
the both sides of $\nu_{d}$. In the more realistic cases where the drive terms
do not necessarily occur at an arbitrarily slow timescale (i.e., can occur at
an arbitrary moment), however, neural dynamics may be phenomenology controlled
by conserved directed percolation transitions and hover around the critical
point. Let us add a drive term in Eq. (14)

$\displaystyle\frac{\partial}{\partial t}\nu\left(\vec{x},t\right)$
$\displaystyle=f\left(\vec{x},t\right)-g\left(\vec{x},t\right)\rho\left(\vec{x},t\right).$
(18)

Then Eq. (15) becomes

$\displaystyle\frac{\partial}{\partial t}\rho\left(\vec{x},t\right)$
$\displaystyle=\left(a+b\nu\left(\vec{x},0\right)\right)\rho\left(\vec{x},t\right)-c\rho^{2}\left(\vec{x},t\right)-b\rho\left(\vec{x},t\right)\int_{0}^{t}\left(f\left(\vec{x},\tau\right)-g\left(\vec{x},\tau\right)\rho\left(\vec{x},\tau\right)\right)\mathsf{d}\tau,$
(19)

If we can ideally fine tune the drive term $f\left(\vec{x},t\right)$ to ensure
that $\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow
r\in\left(0,\infty\right)$, the steady state solutions of Eqs. (18-19) are

$\displaystyle\rho\left(\vec{x},t\right)$
$\displaystyle=\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}\rightarrow
r,$ (20) $\displaystyle\nu\left(\vec{x},0\right)$
$\displaystyle=\frac{1}{b}\left(c\frac{f\left(\vec{x},t\right)}{g\left(\vec{x},t\right)}-a\right)\rightarrow\frac{1}{b}\left(cr-a\right),$
(21) $\displaystyle\nu\left(\vec{x},t\right)$
$\displaystyle=\nu\left(\vec{x},0\right)+\int_{0}^{t}\left(f\left(\vec{x},\tau\right)-g\left(\vec{x},\tau\right)\rho\left(\vec{x},\tau\right)\right)\mathsf{d}\tau\rightarrow\frac{1}{b}\left(cr-a\right).$
(22)

Eqs. (20-22) correspond to a steady state of the brain with
$\rho\left(\vec{x},t\right)\rightarrow r$ and conserved energy, which is
similar to SOC. Therefore, the brain may self-organize to a quasi-critical
region around $\nu_{c}$, the critical point of SOC. Reaching the critical
point requires ideal fine tuning. These emerged conserved-directed-percolation
behaviours enable scientists to recognize SOqC in a similar manner of SOC in
practice (i.e., when stimulus inputs can occur at any moment) Bonachela Munoz
(2009); Buendía, Di Santo . (2020).

####  3.1.5 Summary of theoretical relations

Taken together, neuroscientists can approximately verify the existence of
brain criticality in the space of absorbing and active phases with specific
tools coming from directed percolation theory. This is because OC, qC, SOC,
and SOqC exhibit or approximately exhibit directed percolation behaviours
under certain conditions. The verification may be inaccurate since the
approximation holds conditionally. As for the brain criticality between
asynchronous and synchronous phases, however, the universality class
properties become rather elusive because an analytic and complete theory of
synchronous phase transitions in the brain remains absent yet (see Di Santo .
(2018); Buendía . (2021) for early attempts). Although some behaviours of
absorbing phase transitions can be observed in synchronous phase transitions
(e.g., see Di Santo . (2018); Girardi-Schappo . (2021); Buendía . (2021);
Fontenele . (2019)), there also exist numerous differences between them (e.g.,
see Fontenele . (2019); Buendía . (2021); Girardi-Schappo . (2021)). As
suggested by Dalla Porta Copelli (2019), it remains elusive if directed
percolation properties are applicable, at least conditionally applicable, to
analyzing synchronous phase transitions. More explorations are necessary in
the future.

There are numerous properties of brain criticality predicted by directed
percolation theory, among which, neural avalanche exponents (the power-law
exponents of lifetime and size distributions), scaling relation, universal
collapse shape, and slow decay of auto-correlation are applicable in both
analytic derivations and statistical estimations from empirical data. These
properties are our main focuses. For convenience, we summarize important
glossaries and symbol conventions before we discuss theoretical details (Table
2).

Table 2: Glossaries and symbol conventions. Please note that Table 2 mainly contains important glossaries with fixed symbol definitions. There are many symbols uncovered by Table 2 since they are only used for mathematical derivations. Variable | Meaning
---|---
$T$ | The lifetime of the neural avalanche
$S$ | The size of the neural avalanche
$A$ | The area of the neural avalanche
$\langle S\left(T\right)\rangle$ | The averaged size of neural avalanches with lifetime $T$
$\langle S\left(t\mid T\right)\rangle$ | The averaged time-dependent avalanche size at moment $t$ during neural avalanches with the lifetime $T$
$\mathcal{P}_{T}\left(t\right)$ | The probability distribution of neural avalanche lifetime
$\mathcal{P}_{S}\left(s\right)$ | The probability distribution of neural avalanche size
$\alpha$ | Power-law exponent of the neural avalanche lifetime distribution $\mathcal{P}_{T}\left(t\right)\propto t^{-\alpha}$
$\beta$ | Power-law exponent of the neural avalanche size distribution $\mathcal{P}_{S}\left(s\right)\propto s^{-\beta}$
$\gamma$ | Power exponent of the neural avalanche area $A\propto T^{\gamma}$
$\mathcal{H}\left(\cdot\right)$ | Universal scaling function
$\operatorname{Cov}\left(\cdot,\cdot\right)$ | Auto-correlation function
$\chi$ | Power-law decay rate of auto-correlation
$\xi$ | Exponential decay rate of auto-correlation

###  3.2 Neural avalanche exponents

As we have mentioned above, neural avalanches are expected to exhibit power-
law properties in their lifetime and size distributions when the brain is at
the critical point Hinrichsen (2000); Lübeck (2004); Larremore . (2012).
Therefore, it is pivotal to confirm the detailed values of neural avalanche
exponents. To analytically derive these exponents, one can consider critical
branching process García-Pelayo . (1993); Harris Edward (1963); Otter (1949);
di Santo . (2017); Gros (2010), neural field theory Robinson (2021), mean-
field Abelian sandpile models Janowsky Laberge (1993); DS. Lee . (2004), and
avalanches in networks Larremore . (2012). The key idea to derive neural
avalanche exponents shared by these existing theories is to confirm the
explicit forms of $\mathcal{P}_{T}\left(t\right)$ and
$\mathcal{P}_{S}\left(s\right)$, the probability distributions of the lifetime
and size of neural avalanches, under ideal conditions (e.g., when the maximum
lifetime and size are unlimited and can be infinitely large). In real cases
where lifetime and size are restricted because the brain is a finite system,
slight deviations from idea values may be observed but theoretical derivations
of neural avalanche exponents principally hold.

To present accessible expositions, we consider a critical branching process in
Eqs. (23-34) to describe related backgrounds. More importantly, we present a
novel and simple idea to calculate target exponents in the context of
neuroscience in Box 1. Abstractly, one can define
$\mathcal{P}\left(n,t\right)$ as the probability for an active neuron at
moment $t$ to activate $n$ post-synaptic neurons subsequently and define
$\mathcal{Z}\left(n,t\right)$ as the probability of finding $n$ active neurons
at moment $t$. Meanwhile, one denote

$\displaystyle\mathcal{F}\left(x,t\right)=\sum_{n=0}^{\infty}\mathcal{P}\left(n,t\right)x^{n},$
(23)
$\displaystyle\mathcal{G}\left(x,t\right)=\sum_{n=0}^{\infty}\mathcal{Z}\left(n,t\right)x^{n}$
(24)

as the corresponding generating functions Fristedt Gray (2013); Rao Swift
(2006). Then, one can readily see the recursion relation

$\displaystyle\mathcal{G}\left(x,t\right)$
$\displaystyle=\sum_{n=0}^{\infty}\mathcal{Z}\left(n,t-\delta
t\right)\mathcal{F}\left(x,t-\delta t\right)^{n},$ (25)
$\displaystyle=\mathcal{G}\left(\mathcal{F}\left(x,t-\delta t\right),t-\delta
t\right),$ (26)

where $\delta t$ denotes the minimum time step. Eq. (26) implies that
branching processes are Markovian. Similarly, one can measure the expectations

$\displaystyle\mu\left(t\right)$ $\displaystyle=\frac{\partial}{\partial
x}\mathcal{F}\left(x,t\right)\Big{|}_{x=1},$ (27)
$\displaystyle\phi\left(t\right)$ $\displaystyle=\frac{\partial}{\partial
x}\mathcal{G}\left(x,t\right)\Big{|}_{x=1}$ (28)

to derive another recursion relation

$\displaystyle\phi\left(t\right)$ $\displaystyle=\frac{\partial}{\partial
x}\mathcal{F}\left(x,t-1\right)\Big{|}_{x=1}\frac{\partial}{\partial
x}\mathcal{G}\left(x,t-1\right)\Big{|}_{x=1},$ (29)
$\displaystyle=\phi\left(t-\delta\right)\mu\left(t-\delta\right),$ (30)
$\displaystyle=\prod_{\tau=0}^{t-\delta}\mu\left(\tau\right).$ (31)

Note that Eq. (31) is derived from the fact that $\phi\left(0\right)=1$ (one
neuron is activated at moment $0$ to trigger neural avalanches). Please see
Marković Gros (2014) for more explanations of Eqs. (25-31). Assuming that
$\phi\left(t\right)$ scales as $\exp\left(\lambda t\right)$ for large $t$, we
know that $\phi\left(t\right)$ converges to $0$ given a negative Lyapunov
exponent $\lambda$ (the branching process is sub-critical García-Pelayo .
(1993); Harris Edward (1963); Otter (1949); di Santo . (2017); Gros (2010))
and diverges with a positive Lyapunov exponent $\lambda$ (the branching
process is super-critical García-Pelayo . (1993); Harris Edward (1963); Otter
(1949); di Santo . (2017); Gros (2010)). Here $\lambda$ can be defined
according to Eq. (32)

$\displaystyle\lambda$
$\displaystyle=\lim_{t\rightarrow\infty}\ln\left(\frac{1}{t}\phi\left(t\right)\right)=\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{\tau=0}^{t-\delta}\ln\left(\mu\left(\tau\right)\right).$
(32)

If the branching process is homogeneous, namely
$\mathcal{P}\left(n,t\right)=\mathcal{P}\left(n\right)$,
$\mathcal{Z}\left(n,t\right)=\mathcal{Z}\left(n\right)$,
$\mu\left(\tau\right)=\mu$, and $\phi\left(\tau\right)=\phi$ for every moment
$\tau$, then $\mu=1$ is the condition for the branching process to be
critical. To relate these results with neural avalanches, one only need to
consider the avalanche size $S=\sum_{t}z\left(t\right)$, where
$z\left(t\right)\sim\mathcal{Z}$ denotes the number of active neurons at
moment $t$, and the avalanche life time $T=\min\\{t\mid
z\left(t\right)>0\;\text{and}\;z\left(t+\delta t\right)=0\\}$. It has been
analytically proved that in terms of fixed environments and a Poisson
generating function $\mathcal{F}$ one can derive Otter (1949)

$\displaystyle\mathcal{P}_{S}\left(s\right)$ $\displaystyle\sim
s^{-3/2}\mu^{s-1}\exp\left[s\left(1-\mu\right)\right],$ (33)
$\displaystyle\mathcal{P}_{T}\left(t\right)$ $\displaystyle\sim
t^{-2}\mu^{t-1}\exp\left[t\left(1-\mu\right)\right].$ (34)

In the case with $\mu=1$, one can obtain $\mathcal{P}_{S}\left(s\right)\sim
s^{-3/2}$ and $\mathcal{P}_{T}\left(t\right)\sim t^{-2}$, the power-law
distributions of neural avalanche size and neural avalanche life time García-
Pelayo . (1993); Harris Edward (1963); Otter (1949); di Santo . (2017); Gros
(2010); Robinson (2021); Janowsky Laberge (1993); DS. Lee . (2004); Larremore
. (2012); Lombardi . (2017); Jung . (2020), from Eqs. (33-34).

The derivations of avalanche exponents $\alpha=2$ and $\beta=\frac{3}{2}$ are
non-trivial. However, few neuroscience studies elaborate on these details,
impeding researchers from understanding the theoretical foundations of brain
criticality in the brain. The importance of these derivations is beyond the
detailed values of avalanche exponents since they reveal the fundamental
properties of neural dynamics di Santo . (2017); Cocchi . (2017); Girardi-
Schappo (2021). In Box 1, we sketch an original idea to derive these avalanche
exponents in the terminology of neuroscience. In Fig. 3A, we present graphical
illustrations of our idea in Box 1.

Derivations of neural avalanche exponents Consider a time-continuous neural
dynamics process, where an active neuron implies three possibilities: becoming
absorbed with probability $\varsigma$, activating another neuron with
probability $\eta$, or remaining effect-free with probability
$1-\left(\varsigma+\eta\right)$. In critical states, we have $\varsigma=\eta$
García-Pelayo . (1993). We define $\mathcal{A}_{n}\left(t\right)$ as the
probability for $n$ active neurons to exist at $t^{*}+t$ given that $1$ active
neuron exists at $t^{*}$. Assuming the independence of neuron activation, we
have

$\displaystyle\mathcal{A}_{n}\left(t\right)=\sum_{\begin{subarray}{c}n_{1}+\ldots+n_{k}=n\end{subarray}}\mathcal{A}_{n_{1}}\left(t\right)\ldots\mathcal{A}_{n_{k}}\left(t\right).$
(35)

If $\mathcal{A}_{n}(t),n\in\mathbb{N}^{+}$ admits a Maclaurin expansion
$\mathcal{A}_{n}\left(t\right)=a_{n}t+o\left(t^{2}\right)$ (when $n\neq 1$) or
$\mathcal{A}_{n}\left(t\right)=a_{n}t+1+o\left(t^{2}\right)$ (when $n=1$)
where $a_{n}=\mathrm{d}\mathcal{A}_{n}\left(0\right)/\mathrm{d}t$, we can
readily derive $a_{0}=a_{2}=\varsigma$ and $a_{1}=-2\varsigma$ García-Pelayo .
(1993). Meanwhile, we can know

$\displaystyle\mathcal{A}_{n}\left(t+\mathrm{d}t\right)-\mathcal{A}_{n}\left(t\right)=\sum_{k=0}^{\infty}a_{k}\mathcal{A}_{n-k}\left(t\right)\mathrm{d}t.$
(36)

Eqs. (15-16) readily lead to

$\displaystyle\frac{\partial}{\partial
t}\mathcal{W}\left(x,t\right)=\sum_{k=0}^{\infty}a_{k}\sum_{n=0}^{\infty}\left(\sum_{\begin{subarray}{c}n_{1}+\ldots+n_{k}=n-k\end{subarray}}\prod_{i=1}^{k}\mathcal{A}_{n_{i}}\left(t\right)\right)x^{n}=\sum_{k=0}^{\infty}a_{k}\mathcal{W}\left(x,t\right)^{k},$
(37)

where
$\mathcal{W}\left(x,t\right)=\sum_{n=0}^{\infty}\mathcal{A}_{n}\left(t\right)x^{n},\;x\in\left[0,1\right]$
denotes the generating function. Applying a trick introduced in García-Pelayo
. (1993), we define $\mathcal{H}\left(x\right)=\frac{\partial}{\partial
t}\mathcal{W}\left(x,0\right)$, which naturally leads to
$\frac{\partial}{\partial
t}\mathcal{W}\left(x,t\right)=\mathcal{H}\left(\mathcal{W}\left(x,t\right)\right)$.
Meanwhile, $\mathcal{H}\left(x\right)=\varsigma\left(1-x\right)^{2}$ can be
derived based on $a_{0}$, $a_{1}$, and $a_{2}$ García-Pelayo . (1993). Taken
together, we have

$\displaystyle\frac{\partial}{\partial
t}\mathcal{W}\left(x,t\right)=\varsigma\left(1-\mathcal{W}\left(x,t\right)\right)^{2}.$
(38)

Note that the initial condition is $\mathcal{W}\left(x,0\right)=x$ since one
neuron is activated at $t^{*}$. Solving Eq. (38), we derive that

$\displaystyle\mathcal{W}\left(x,t\right)=\frac{\varsigma\left(1-x\right)t}{\varsigma\left(1-x\right)t+1}.$
(39)

Therefore, we have
$\mathcal{A}_{0}\left(t\right)=\mathcal{W}\left(0,t\right)=\frac{\varsigma
t}{\varsigma t+1}$, supporting a calculation of lifetime distribution
$\mathcal{P}_{T}\left(t\right)$

$\displaystyle\lim_{t\rightarrow\infty}\mathcal{P}_{T}\left(t\right)=\lim_{t\rightarrow\infty}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{W}\left(0,t\right)\sim
t^{-2}.$ (40)

Following García-Pelayo . (1993); Harris Edward (1963); Otter (1949), one can
similarily calculate

$\displaystyle\lim_{s\rightarrow\infty}\mathcal{P}_{S}\left(s\right)\sim
s^{-\frac{3}{2}}.$ (41)

There are three important things to remind. First, the lifetime exponent
$\alpha=2$ and size exponent $\beta=\frac{3}{2}$ can only be treated as ideal
exponents under mean field assumptions of directed percolation. There are
numerous factors, such as granularity, network topology, and neural dynamics
variability, can be considered in derivations to affect the detailed values of
avalanche exponents Girardi-Schappo (2021); Bonachela Munoz (2009). In Table
3, we summarize the possible intervals of $\alpha$ and $\beta$ in empirical
neural data. Second, $\alpha$ and $\beta$ alone are not sufficient to verify
the existence of brain criticality. Even when the actual values of $\alpha$
and $\beta$ in empirical data are exactly equal to theoretical predictions,
they may still not satisfy the scaling relation and universal collapse.
Meanwhile, as we shall discuss later, estimating $\alpha$ and $\beta$ in
practice is statistically error-pone. Third, one can not confirm or disprove a
detail type of brain criticality based on $\alpha$ and $\beta$ unless
additional information is provided. Although four types of brain criticality
exhibit (e.g., OC) or approximatively exhibit (e.g., qC, SOC, and SOqC)
directed percolation behaviours under certain conditions, these preconditions
are difficult to verify in practice.

###  3.3 Scaling relation

In the previous section, we discuss how the neural avalanche lifetime and size
distributions exhibit power-law properties when the brain is at the critical
point Hinrichsen (2000); Lübeck (2004); Larremore . (2012). Apart from
lifetime $T$ and size $S$, there are several other quantities that
characterize neural avalanches, such as area $A$ (number of distinct active
neurons, measured as $A\simeq\langle S\left(T\right)\rangle$ where the
expectation $\langle\cdot\rangle$ is averaged across all neural avalanches
with the same lifetime $T$) and radius exponent $R$ (radius of gyration)
Lübeck Heger (2003); Lübeck (2004). In general, the corresponding probability
distributions of these four quantities decay algebraically

$\displaystyle\mathcal{P}_{X}\left(x\right)\propto x^{-\lambda_{X}},$ (42)

where random variable $X\in\\{S,T,A,R\\}$ can be an arbitrary quantity to
characterize neural avalanches. The avalanche exponent $\lambda_{X}$ is
defined according to the selected meaning of $X$ (e.g., $\lambda_{T}=2$ and
$\lambda_{T}=\frac{3}{2}$ under mean filed assumptions). Assuming that
variables $\\{S,T,A,R\\}$ scale as a power of each other

$\displaystyle X^{\prime}\propto X^{\psi_{X^{\prime}X}},\;\forall
X,\;X^{\prime}\in\\{S,T,A,R\\},$ (43)

we can derive the scaling relation from Eqs. (42-43)

$\displaystyle\psi_{X^{\prime}X}=\frac{\lambda_{X}-1}{\lambda_{X^{\prime}}-1}.$
(44)

If we let $X^{\prime}=A$ and $X=T$, we can specify Eq. (44) as

$\displaystyle\gamma$ $\displaystyle=\frac{\alpha-1}{\beta-1},$ (45)

where $\mathcal{P}_{T}\left(t\right)\propto t^{-\alpha}$,
$\mathcal{P}_{S}\left(s\right)\propto s^{-\beta}$, and $A\propto T^{\gamma}$.
Eq. (45) leads to $\gamma=2$ in the mean field theory of directed percolation.
In Table 3, one can see the possible interval of $\gamma$ in empirical neural
data. Eq. (45) is widely used as a criterion to verify if the brain is at the
critical point in neuroscience studies (e.g., Dalla Porta Copelli (2019); N.
Friedman . (2012); Fontenele . (2019); Ponce-Alvarez . (2018)). Once the
scaling relation is confirmed among observed neural avalanche exponents, it
indicates key features of the universality class (please note that $\alpha$,
$\beta$, and $\gamma$ should be derived independently). For neuroscientists,
the importance of Eq. (45) lies in that it provides extra verification of the
validity of estimated neural avalanche exponents. This verification is
necessary given that neural avalanche exponent estimation is frequently
inaccurate Fontenele . (2019). In Fig. 3B, we illustrate the scaling relation
in Eq. (45) under mean field assumptions.

In Lübeck (2004), one can further learn about how brain criticality is mapped
to an directed percolation transition characterized by ordinary critical
exponents. Meanwhile, one can see how to connect these neural avalanche
exponents with second order phase transition exponents Lübeck Heger (2003).

Figure 3: Conceptual illustrations of the neural avalanche properties predicted by analytic theories of brain criticality. A, Illustrations of the framework to derive neural avalanche exponents in Box 3. B, Illustrations of the scaling relation satisfied by neural avalanches under mean field assumptions. C, Illustrations of the universal collapse shape of neural avalanches. The un-scaled plot of $t$ vs. $\langle S\left(t\mid T\right)\rangle$ (upper left) and the scaled plot $\langle S\left(t\mid T\right)\rangle T^{1-\gamma}$ vs. $\frac{t}{T}$ (upper right) are shown for comparison. Here terms $\frac{1}{T}$ and $T^{1-\gamma}$ respectively serve as scale factors on x-axis and y-axis to create a universal collapse shape. Meanwhile, the symmetric collapse shape in Eq. (51) (bottom left) and the asymmetric collapse shape controlled by skewness $c$ in Eq. (52) (bottom right) are also presented. D, Auto-correlations and their decays in critical (left) and non-critical (right) cases are shown. Auto-correlations are calculated after a $t_{i}\in\left[0,T\right)$ is randomly selected as a reference (upper left and right). Meanwhile, the auto-correlation decays measured on $t_{j}\in\left[t_{i},T\right]$ is fitted to derive the coefficients $\chi$ (x-axis corresponds to $\ln\left(\frac{t_{j}-t_{i}}{T}\right)$) and $\xi$ (x-axis corresponds to $\frac{t_{j}-t_{i}}{T}$) in Eq. (53) and Eq. (55) (bottom left and right). One can see that auto-correlations in the critical case have slower decays (smaller $\chi$ and $\xi$) than those in the non-critical case. Table 3: Neural avalanche exponents with scaling relation in empirical data. The data is acquired from Girardi-Schappo (2021), where 45 experimental observations of neuronal avalanches reported by 30 studies are summarized. These observations can be classified according to the recording techniques of neural avalanches. Detailed data classification criteria (e.g., details of spike sorting and thresholding) can be seen in Girardi-Schappo (2021). We only include the data where $\alpha$, $\beta$, and $\gamma$ are all recorded and satisfy the scaling relation in Eq. (45). For LFP recordings filtered by spike sorting, included observations are reported by Carvalho . (2021); Fontenele . (2019); Senzai . (2019); Fosque . (2021); Ma . (2019); Mariani . (2021). For LFP recordings with thresholding, included observations come from Mariani . (2021); Shew . (2015). For Ca and voltage imaging, observations are provided by Ponce-Alvarez . (2018); Yaghoubi . (2018). Although numerous studies report neural avalanches in whole-brain imaging (e.g., MEG, M/EEG, and invasive ECoG), these studies either do not report three exponents together Shriki . (2013); Zhigalov . (2015); Palva . (2013) or have not observed the scaling relation in Eq. (45) Varley . (2020). One can see Girardi-Schappo (2021) for a summary of these results. Data type | Observed interval of $\alpha$ | Observed interval of $\beta$ | Observed interval of $\gamma$
---|---|---|---
LFP recordings filtered by spike sorting | $\alpha\in\left[1.35,2.67\right]$ | $\beta\in\left[1.3,2.5\right]$ | $\gamma\in\left[1.16,1.48\right]$
LFP recordings with thresholding | $\alpha\in\left[1.82,2.84\right]$ | $\beta\in\left[1.57,2.59\right]$ | $\gamma\in\left[1.12,1.39\right]$
Ca and voltage imaging | $\alpha\in\left[2.15,3.5\right]$ | $\beta\in\left[1.5,2.3\right]$ | $\gamma\in\left[1.75,2.5\right]$

###  3.4 Universal collapse shape

####  3.4.1 Universal collapse with an implicit scaling function

Apart from the scaling relation discussed above, the average temporal shape of
bursts, a fundamental signature of avalanches Baldassarri . (2003); Laurson .
(2013); Papanikolaou . (2011), can also be used to verify the existence of
brain criticality in a more precise manner. This approach has been previously
applied on diverse physical systems, such as plastically deforming crystals
Laurson Alava (2006) and Barkhausen noise Mehta . (2002); Papanikolaou .
(2011), and is recently introduced into neuroscience Ponce-Alvarez . (2018);
Fontenele . (2019); Pausch . (2020); Dalla Porta Copelli (2019); N. Friedman .
(2012). To understand this approach, let us step back to the power relation in
Eq. (43) and specify that $X^{\prime}=A$ and $X=T$. These settings naturally
lead to

$\displaystyle A\simeq\langle S\left(T\right)\rangle\equiv\int_{0}^{T}\langle
S\left(t\mid T\right)\rangle dt\propto T^{\gamma},$ (46)

where $\langle S\left(t\mid T\right)\rangle$ measures the averaged time-
dependent avalanche size during an avalanche and the expectation
$\langle\cdot\rangle$ is averaged across all neural avalanches with the same
lifetime $T$. Eq. (46) can be readily reformulated as

$\displaystyle\langle S\left(t\mid T\right)\rangle\propto T^{\gamma-1}.$ (47)

The general form of Eq. (47) is usually given by Baldassarri . (2003); Laurson
. (2013); Papanikolaou . (2011)

$\displaystyle\langle S\left(t\mid
T\right)\rangle=T^{\gamma-1}\mathcal{H}\left(\frac{t}{T}\right),$ (48)

where $\mathcal{H}\left(\cdot\right)$ denotes a universal scaling function.
When the brain is at the critical point, all data of $\langle S\left(t\mid
T\right)\rangle T^{1-\gamma}$ is expected to collapse onto
$\mathcal{H}\left(\cdot\right)$ with reasonable errors Baldassarri . (2003);
Laurson . (2013); Papanikolaou . (2011). Here the terminology “collapse onto”
means that all data generally exhibits a similar pattern in a plot of $\langle
S\left(t\mid T\right)\rangle T^{1-\gamma}$ vs. $\frac{t}{T}$ (e.g., all data
follows function $\mathcal{H}\left(\cdot\right)$). Meanwhile, scaling function
$\mathcal{H}\left(\cdot\right)$ is expected to be a parabolic function
Baldassarri . (2003); Laurson . (2013); Papanikolaou . (2011). By testing
these properties, neuroscientists can verify whether the brain is at
criticality (e.g., Dalla Porta Copelli (2019); N. Friedman . (2012); Fontenele
. (2019); Ponce-Alvarez . (2018)). In Fig. 3C, we graphically illustrate these
properties.

####  3.4.2 Universal collapse with an explicit scaling function

Under specific conditions, researchers can further consider an explicit form
of scaling function $\mathcal{H}\left(\cdot\right)$ Laurson . (2013). Assuming
that the early-time growth of neural avalanches averagely follows a power-law
of time, one can derive that $\langle S\left(t\mid T\right)\rangle\propto
t^{\kappa}$ for certain $\frac{t}{T}\leq\varepsilon\ll 1$. Meanwhile, one
knows that $\langle S\left(\varepsilon T\mid T\right)\rangle\propto
T^{\gamma-1}$ should hold according to Eq. (28). To ensure these two
properties, one needs to have $\langle S\left(\varepsilon T\mid
T\right)\rangle\propto\left(\varepsilon T\right)^{\kappa}\propto
T^{\gamma-1}$, which readily leads to $\kappa=\gamma-1$. Based on these
derivations, one can know

$\displaystyle\langle S\left(t\mid T\right)\rangle\propto t^{\gamma-1},\;t\ll
T.$ (49)

To find an explicit form of $\mathcal{H}\left(\cdot\right)$ that satisfies
Eqs. (48-49), one can consider a possible answer Laurson . (2013)

$\displaystyle\mathcal{H}\left(\frac{t}{T}\right)=\left[\frac{t}{T}\left(1-\frac{t}{T}\right)\right]^{\gamma-1},$
(50)

which can be analytically derived by multiplying Eq. (49) by
$\left(1-\frac{t}{T}\right)^{\gamma-1}$. Here
$\left(1-\frac{t}{T}\right)^{\gamma-1}$ is a term to characterize the
deceleration at the ends of neural avalanches Laurson . (2013). Because
$\gamma=2$ is expected for critical neural avalanches under mean field
assumptions, Eq. (48) and Eq. (50) imply that

$\displaystyle\langle S\left(t\mid T\right)\rangle\propto
t\left(1-\frac{t}{T}\right).$ (51)

This result is consistent with the prediction by the ABBM model in the limit
of vanishing drive rate and demagnetizing factor N. Friedman . (2012); Laurson
. (2013).

A potential limitation of Eq. (51) in applications lies in its internal
symmetry property Laurson . (2013). Although avalanches under mean-field
frameworks have a symmetric average shape N. Friedman . (2012), it does not
mean that symmetry generally holds in real complex systems Laurson . (2013).
Applying Eq. (48) on neural data, researchers may observe a nonstandard
parabolic function $\mathcal{H}\left(\cdot\right)$ with specific skewness.
This does not necessarily mean that neural dynamics is not at criticality.
When neural avalanches are time-irreversible (this is generally true in the
brain since the detailed balance of neural dynamics is frequently broken Lynn
. (2021)), one can consider small temporal asymmetry in the collapse shape
Laurson . (2013). To characterize potential asymmetry, one can add a
correction term controlled by skewness degree $c$ into Eq. (51)

$\displaystyle\langle S\left(t\mid T\right)\rangle\propto
t\left(1-\frac{t}{T}\right)\left[1-c\left(\frac{t}{T}-\frac{1}{2}\right)\right].$
(52)

If $c=0$, then Eq. (52) reduces to Eq. (51). Otherwise, neural avalanches can
have a temporally asymmetric collapse shape with a positive ($c>0$) or
negative ($c<0$) skewness Laurson . (2013). We suggest that Eq. (52) may be
more applicable to real data of neural dynamics. In Fig. 3C, we show examples
of Eqs. (51-52).

###  3.5 Slow decay of auto-correlation

In applications, researchers can also consider a more practical verification
of the potential brain criticality. When the brain is at the critical point, a
slow decay of auto-correlation is expected to occur in neural avalanches,
corresponding to long-range correlations Schaworonkow . (2015); Smit . (2011);
Dalla Porta Copelli (2019); Erdos . (2018). This slow decay property is
initially derived from the power-law decay of auto-correlation, which can be
analytically derived as a part of the scaling relation if ordinary critical
exponents of directed percolation transition are considered (see Girardi-
Schappo (2021) for details). The power-law decay is expressed as

$\displaystyle\ln\left[\frac{\operatorname{Cov}\left(S\left(t_{i}\mid
T\right),S\left(t_{j}\mid
T\right)\right)}{\operatorname{Cov}\left(S\left(t_{i}\mid
T\right),S\left(t_{i}\mid
T\right)\right)}\right]=-\chi\ln\left(\frac{t_{j}-t_{i}}{T}\right)+r,$ (53)

where $t_{i}\in\left[0,T\right)$ is used as a reference and
$t_{j}\in\left[t_{i},T\right]$ traverses the entire interval Schaworonkow .
(2015); Smit . (2011). According to the Wiener–Khinchin theorem, coefficient
$\chi$ is related to $\mathcal{S}\left(f\right)$, the power spectrum of neural
avalanches (notion $f$ denotes frequency) Bak . (1987); Linkenkaer-Hansen .
(2001); Girardi-Schappo (2021). One may expect $\mathcal{S}\left(f\right)\sim
f^{-\upsilon}$ at the critical point, where $\chi=1-\upsilon$ Bak . (1987);
Linkenkaer-Hansen . (2001); Girardi-Schappo (2021). The power-law decay of
auto-correlation in Eq. (53) breaks down when $\upsilon=1$, leading to
infinitely long temporal correlations. Therefore,
$\chi\in\left[0,\infty\right)$ in Eq. (53) is expected to be sufficiently
small. Certainly, the actual value of $\chi$ may not be perfectly zero in
empirical data. For instance, $\chi\in\left[0.58\pm 0.23,0.73\pm 0.31\right]$
and $\chi\in\left[0.52\pm 0.35,0.81\pm 0.32\right]$ are observed in
spontaneous alpha oscillations in MEG and EEG data, respectively Linkenkaer-
Hansen . (2001).

Apart from verifying power-law decay directly, one can also consider the
exponential decay, which is active in neuroscience as well Pausch . (2020);
Wilting Priesemann (2019); Miller Wang (2006). The exponential decay can
described by

$\displaystyle\frac{\partial}{\partial
t}\operatorname{Cov}\left(S\left(t_{i}\mid T\right),S\left(t_{j}\mid
T\right)\right)=-\xi\operatorname{Cov}\left(S\left(t_{i}\mid
T\right),S\left(t_{j}\mid T\right)\right).$ (54)

Eq. (54) directly leads to

$\displaystyle\ln\left[\frac{\operatorname{Cov}\left(S\left(t_{i}\mid
T\right),S\left(t_{j}\mid
T\right)\right)}{\operatorname{Cov}\left(S\left(t_{i}\mid
T\right),S\left(t_{i}\mid
T\right)\right)}\right]=-\xi\left(\frac{t_{j}-t_{i}}{T}\right)+r.$ (55)

The exponential decay can be seen in the dynamics with short-term correlations
(i.e., correlations have a characteristic time scale). Mathematically, the
exponential decay can be related to power-law decay in a form of
$x^{-y}=\Gamma\left(y\right)\int_{0}^{\infty}z^{y-1}\exp\left(-xz\right)\mathsf{d}z$,
where $\Gamma\left(\cdot\right)$ denotes the Gamma function. When
$\xi\in\left[0,\infty\right)$ is sufficiently small, Eq. (55) can be treated
as a looser criterion that approximately verifies the slow decay of auto-
correlation and may be more applicable to non-standard brain criticality
(e.g., qC and SOqC) Wilting Priesemann (2019). Despite of its practicality,
this looser criterion should be used with cautions since it is not
analytically derived from criticality theories.

In Fig. 3D, we illustrate examples of auto-correlation slow decay in critical
cases and compare them with non-critical cases. Compared with other properties
previously discussed, a slow auto-correlation decay can be readily verified by
conventional data fitting. However, we need to note that one should not
confirm or reject the possibility of brain criticality only based on the decay
characteristic of auto-correlation in Eqs. (53-55). This is because Eqs.
(53-55) only serve as the approximate descriptions of long-range correlations
at criticality. The strict criterion $\chi,\xi\rightarrow 0$ is rarely seen in
empirical data while the determination of whether $\chi$ and $\xi$ are
sufficiently small in the looser criterion is relatively subjective.

In summary, we have reviewed the physical foundations of identifying and
characterizing criticality in the brain. Based on these analytic derivations,
we attempt to present systematic explanations of what is brain criticality and
how to identify potential criticality in neural dynamics. Nevertheless,
physical theories alone are not sufficient to support neuroscience studies
because the implementation of these theories on empirical data is even more
challenging than the theories themselves. To overcome these challenges, one
needs to learn about statistic techniques to computationally estimate brain
criticality from empirical data.

##  4 Brain criticality: statistic techniques

While most properties of neural avalanches analytically predicted by the
physical theories of brain criticality can be estimated by conventional
statistic techniques, there exist several properties that frequently imply
serious validity issues and deserve special attention. Below, we discuss them
in detail.

###  4.1 Estimating neural avalanche exponents

Perhaps the estimation of neural avalanche exponents from empirical data is
the most error-prone step in brain criticality analysis. The least-square
approach is abused in fitting power-law data and frequently derives highly
inaccurate results Clauset . (2009); Virkar Clauset (2014). To derive neural
avalanche exponents $\alpha$ and $\beta$ in Eq. (25) with reasonable errors,
one need to consider the maximum likelihood estimation (MLE) approach and
corresponding statistic tests (see MLE on un-binned data Clauset . (2009) and
binned data Virkar Clauset (2014)). Taking the avalanche size distribution as
an instance, the estimator $\widehat{\beta}$ of distribution exponent $\beta$
is expected to maximize the log-likelihood function

$\displaystyle\mathcal{L}\left(\beta\right)$
$\displaystyle=-n\ln\left[\zeta\left(\beta,s^{\prime}\right)\right]-\beta\sum_{i=1}^{n}\ln\left(s_{i}\right),$
(56) $\displaystyle\mathcal{L}\left(\beta\right)$
$\displaystyle=n\left(\beta-1\right)\ln
b^{\prime}+\sum_{i=1}^{k}h_{i}\ln\left(b_{i}^{1-\beta}-b_{i+1}^{1-\beta}\right).$
(57)

Here Eq. (56) and Eq. (57) denote the log-likelihood functions on un-binned
and binned data, respectively. Function $\zeta\left(\cdot,\cdot\right)$ is the
generalized zeta function Clauset . (2009); Bauke (2007). Notion $s$ denotes
avalanche size samples in Eq. (33) and Eq. (41) Clauset . (2009). Notion $b$
denotes bin boundaries defined on these samples and $h$ counts the number of
samples within each bin Virkar Clauset (2014). Notions $s^{\prime}$ and
$b^{\prime}$ are the lower cutoffs of un-binned and binned power-law
distributions Clauset . (2009); Virkar Clauset (2014). They are necessary
because few empirical data exhibits power-law properties on the entire
distribution Clauset . (2009). Notions $n$ and $k$ measure the numbers of
samples and bins above cutoffs, respectively Clauset . (2009); Virkar Clauset
(2014). To estimate $\widehat{\beta}$ precisely, researchers are suggested to
follow several indispensable steps Clauset . (2009); Virkar Clauset (2014):
(1) for each potential choice of $s^{\prime}$ or $b^{\prime}$, estimate the
power-law model on the distribution tail above the cutoff. Compute the
Kolmogorov–Smirnov (KS) goodness-of-fit statistic between the cumulative
probability distributions of power-law model and empirical data. Find the
ideal choice of $s^{\prime}$ or $b^{\prime}$ that minimizes KS statistic; (2)
derive the corresponding estimator $\widehat{\beta}$ and KS statistic based on
the chosen cutoff; (3) use the semi-parametric bootstrap to generate numerous
synthetic data distributions that follow the estimated power-law model above
the cutoff but follow the empirical data distribution below the cutoff.
Estimate new power-law models on these synthetic data distributions and
measure the goodness-of-fit by KS statistic. Define a $p$-value, the fraction
of these KS statistics in step (3) that are no less than the KS statistic in
step (2). Rule out the estimated power-law model in steps (1-2) if $p<0.1$
(conservative criterion). Apart from these necessary steps, one can further
consider Vuong’s likelihood ratio test for alternative distribution checking
Vuong (1989); Clauset . (2009); Virkar Clauset (2014) and information loss
measurement of binning approach Virkar Clauset (2014). During the above
process, we measure the goodness-of-fit by KS statistic instead of the well-
known $\chi^{2}$ statistic because the latter has less statistic power Bauke
(2007); Clauset . (2009); Virkar Clauset (2014). Meanwhile, KS statistic is
measured on cumulative probability distributions rather than probability
distributions to control the effects of extreme values in empirical data
Clauset . (2009); Virkar Clauset (2014). Except for the above approach Clauset
. (2009); Virkar Clauset (2014), one can also consider the BIC method (for un-
binned data) Schwarz (1978) and the RT method (for binned data) Reiss Thomas
(2007) for comparisons. In practice, the approaches proposed by Clauset _et
al._ are more robust Clauset . (2009); Virkar Clauset (2014) and have
attracted numerous follow-up studies for improvements (e.g., Yu . (2014);
Marshall . (2016); Deluca Corral (2013)).

###  4.2 Estimating universal collapse shape

Another error-prone step is the calculation and evaluation of the universal
collapse shape, which is closely related to scaling relation analysis.
Deriving the collapse shape from empirical data may be problematic because the
goodness evaluation of collapse shape is rather subjective (e.g., depends on
personal opinions about whether all data follows function
$\mathcal{H}\left(\cdot\right)$ in Eq. (48)) in most cases Marshall . (2016).
Although important efforts have been devoted to quantify if a given data set
exhibits shape collapse Shaukat Thivierge (2016); Bhattacharjee Seno (2001),
common approaches in practice still depend on specific shape collapse
algorithms that search potential scaling parameters (e.g., $\gamma$ in Eq.
(48)) in a data-driven manner Marshall . (2016). In these algorithms,
thresholding on neural avalanches before analyzing the shape collapse is a
standard pre-processing scheme to control noises (e.g., set an avalanche size
threshold and remove all data below the threshold) Marshall . (2016);
Papanikolaou . (2011). While experimental noises are partly limited,
unexpected excursions of scaling parameters away from theoretical predictions
may occur after thresholding as well Laurson . (2009). To our best knowledge,
the effects of thresholding on brain criticality analysis are non-negligible.
Although being highly practical, thresholding may lead to significant
transient effects to cloud the true scaling property Villegas . (2019).
Therefore, any qualitative evaluation of collapse shape after thresholding is
questionable regardless of its practicability. Although an ideal approach
requires further explorations, we suggest researchers to consider following
methods: (1) estimate $\gamma$ by area fitting (e.g., follow Eq. (47) in
scaling relation analysis) and collapse shape fitting (e.g., follow Eq. (48)
in collapse shape analysis), respectively; (2) compare between $\gamma$
derived by these two kinds of fitting and measure the difference. Search for a
threshold that minimizes the difference (e.g., makes variation amplitude
$<1\%$) and maintains a reasonable sample size (e.g., maintains $>80\%$
samples); (3) Given the chosen threshold and corresponding $\gamma$, measure
the difference (e.g., the dynamic time warping Keogh Pazzani (2001)) between
$\langle S\left(t\mid T\right)\rangle T^{1-\gamma}$ derived on neural
avalanches with different lifetime $T$ in the plot of $\langle S\left(t\mid
T\right)\rangle T^{1-\gamma}$ vs. $\frac{t}{T}$. Denote the shape collapse
error as the averaged difference. Combining these three steps, researchers may
partly avoid the errors implied by subjective judgment. Similar ideas can be
seen in Marshall . (2016).

###  4.3 Estimating the slow decay of auto-correlation

Finally, the analysis of slow decay of auto-correlation is also error-prone in
practice. Although this approach is practical and has been extensively applied
(e.g., Pausch . (2020); Wilting Priesemann (2019)), the criterion to determine
if the decay is truly slow (e.g., $\chi>0$ in Eq. (54) and $\xi>0$ in Eq. (54)
are sufficiently small) remains ambiguous. A fixed criterion (e.g.,
$\chi,\xi<0.5$) may serve as an explicit condition of a slow decay. However,
this presupposed criterion may deviate from real situations. For instance, the
baseline of decay rate in a non-critical brain may be essentially high (e.g.,
$\chi,\xi>10$). Even though the decay rate drops significantly when the brain
becomes critical (e.g., $\chi,\xi\simeq 1$), the presupposed criterion is
still unsatisfied and leads to unnecessary controversies on criticality
hypothesis. Given that $\xi$ is principally independent from spatial sub-
sampling on neurons or brain regions at criticality Pausch . (2020), we
suggest researchers to consider the following approaches: (1) do spatial sub-
sampling in both critical and non-critical brains to derive two groups of
$\chi$ or $\xi$ (one group for criticality and another group for non-
criticality); (2) use appropriate statistic tests (e.g., choose $t$-test Kanji
(2006), Kolmogorov-Smirnov test Berger Zhou (2014), or Wilcoxon-Mann-Whitney
test Fay Proschan (2010) according to sample distribution properties) to
verify if two groups of $\chi$ or $\xi$ belong to different distributions.
Test if the expectation and variance of $\chi$ or $\xi$ drops significantly
from the non-critical group to the critical group according to certain effect
sizes.

In summary, statistic techniques bridge between brain criticality theory and
empirical data. However, misconception and misuse of statistic analyses of
neural avalanche properties still occasionally appear in practice. Although
existing techniques remain imperfect in brain criticality analysis, we wish
that our discussion may inspire future studies.

##  5 Brain criticality and other neuroscience theories

Ever since brain criticality is introduced into neuroscience, it is frequently
speculated as contradictory with other traditional neuroscience hypotheses,
such as the conjectured hierarchical processing characteristic of neural
information Felleman Van Essen (1991) and the asynchronous-irregular
characteristic of neural dynamics (e.g., neurons spike independently in
Poisson manners Burns Webb (1976); Softky Koch (1993); Stein . (2005)).
Meanwhile, the differences between brain criticality and scale-free neural
dynamics Chialvo (2010); Martinello . (2017); He (2014) are frequently
neglected. Before we put an end to our review, we discuss the relations
between brain criticality and these neuroscience theories.

###  5.1 Brain criticality and hierarchical processing

The hierarchical processing of neural information Felleman Van Essen (1991) is
initially speculated to contradict critical neural dynamics since hierarchical
topology has not been used as a explicit condition to analytically derive
criticality (e.g., see derivations in García-Pelayo . (1993); Harris Edward
(1963); Otter (1949); di Santo . (2017); Gros (2010); Robinson (2021);
Janowsky Laberge (1993); DS. Lee . (2004); Larremore . (2012)). On the
contrary, random graphs without strict hierarchical structures seem to be more
widespread in criticality derivations. Recently, this speculation has been
challenged by numerous discoveries of the facilitation effects of hierarchical
modular structures on critical dynamics EJ. Friedman Landsberg (2013); Kaiser
Hilgetag (2010); S. Wang Zhou (2012); Rubinov . (2011). Meanwhile,
computational analysis suggests that information transmission in standard
feed-forward networks is maximized by critical neural dynamics Beggs Plenz
(2003). Parallel to neuroscience, a recent machine learning study empirically
observes and analytically demonstrates that artificial neural networks, a kind
of hierarchical structure, self-organize to criticality during learning
Katsnelson . (2021). Therefore, brain criticality is not necessarily
contradictory with hierarchical information processing, yet more analyses are
required to understand how brain criticality affects hierarchical processing
schemes.

###  5.2 Brain criticality and asynchronous-irregular characteristic

Brain criticality and the asynchronous-irregular (AI) characteristic may
correspond to distinct encoding schemes in the brain Wilting Priesemann
(2019); Girardi-Schappo . (2021). While AI characteristic can minimize
redundancy Atick (1992); Bell Sejnowski (1997); Van Hateren van der Schaaf
(1998); Barlow . (1961) to improve neural encoding Van Vreeswijk Sompolinsky
(1996), brain criticality may optimize encoding performance utilizing a series
of reverberations of neural activities Bertschinger Natschläger (2004);
Haldeman Beggs (2005); Kinouchi Copelli (2006); XR. Wang . (2011); Boedecker .
(2012); Shew Plenz (2013); Del Papa . (2017). The coexistence of empirical
evidence of AI and brain criticality characteristics initially confuses
researchers since these characteristics are hypothesized as contradictory with
each other Wilting Priesemann (2019); Girardi-Schappo . (2021). In
experiments, AI characteristic is supported by small correlations between the
spike rates of different neurons in cortical microcircuits Ecker . (2010);
Cohen Kohn (2011) and exponential distributions of inter-spike intervals Kara
. (2000); Carandini Stevens (2004) while brain criticality characteristic is
observed in neural dynamics recorded from multiple species (e.g., awake
monkeys Petermann . (2009), anesthetized rats Gireesh Plenz (2008), slices of
rat cortices Shew . (2009); Beggs Plenz (2003), and humans Poil . (2008)). A
recent study demonstrates that cortical spikes may propagate at somewhere
between perfect criticality (e.g., OC or SOC depending on whether underlying
mechanisms are exogenous or endogenous) and full irregularity Wilting
Priesemann (2019), similar to the cases of qC and SOqC. Meanwhile, it is known
that stimulus drives suppress irregularity in neural activities Molgedey .
(1992). These results imply that brain criticality may not necessarily
contradict AI characteristic. On the contrary, they may coexist when stimulus
drives are too weak to disturb brain criticality (e.g., OC or SOC) and
suppress AI characteristic. In our previous discussions, we have analytically
proven that neural avalanche exponents, the fundamental properties of brain
criticality, can still be derived under the condition of independent neuron
activation, a key feature of AI characteristic Wilting Priesemann (2019). This
result suggests that brain criticality and AI characteristic do not contradict
each other. As for the case where stimulus drives are non-negligible, a recent
study presents an elegant theory to prove that two homeostatic adaptation
mechanisms (i.e., the short-term depression of inhibition and the spike-
dependent threshold increase) enable synaptic excitation/inhibition balance,
AI characteristic, and SOqC to appear simultaneously in the same neural
dynamics Girardi-Schappo . (2021). Similarly, it is suggested that neural
dynamics with criticality or with AI characteristic can be generated by the
same neural populations if the synaptic excitation/inhibition balance is fine
tuned appropriately J. Li Shew (2020).

###  5.3 Brain criticality and power-law behaviours in neural dynamics

Neural dynamics with power-law behaviours is a necessary but insufficient
condition of brain criticality. This property is frequently neglected in
practice. Power-law behaviours are widespread in the nature because it can be
generated by diverse mechanisms, such as exponential curve summation and
preferential attachment Reed Hughes (2002); Mitzenmacher (2004). It has been
reported that the aggregate behaviours of non-critical stochastic systems may
also create scale-free dynamics within a limited range Touboul Destexhe (2010,
2017). In the brain, the generic scale-free properties can be implied by
neutral dynamics, a kind of dynamics where the population size of neutral
individuals (or dynamically homogeneous individuals) does not tend to increase
or decrease after adding a new individual that is neutral to existing ones
(see neutral theories for further explanations Blythe McKane (2007); Liggett
(2006)). This generic property can generate power-law neural avalanches
without criticality Martinello . (2017). Meanwhile, bistability phenomena, a
kind of fine tuned or self-organized discontinuous phase transitions with
limit cycles rather than critical points, can also create neural dynamics with
power-law properties Buendía, Di Santo . (2020); di Santo . (2016); Cocchi .
(2017). Consequently, we emphasize that neural avalanches exponents alone are
insufficient to prove or disprove any brain criticality hypothesis. These
power-law exponents are meaningless for brain criticality hypothesis unless
they satisfy the scaling relation.

##  6 Brain criticality: conclusions on current progresses and limitations

Given what have been reviewed above, we arrive at a point to conclude on the
current progresses and limitations in establishing theoretical foundations of
different types of brain criticality, i.e., ordinary criticality (OC), quasi-
criticality (qC), self-organized criticality (SOC), and self-organized quasi-
criticality (SOqC). As we have suggested, an inescapable cause of various
controversies is the non-triviality of physical theories that analytically
derive brain criticality and statistic techniques that estimate brain
criticality from empirical data. Immoderate omitting of these theoretical
foundations, especially their imperfection, in practice may lead to confusions
on the precise meaning, identification criteria, and biological corollaries of
brain criticality. To address these problems, we have introduced mainstream
theoretical foundations of brain criticality, reformulated them in the
terminology of neuroscience, and discussed their mistakable details.

Thanks to the increasing efforts devoted to improving theoretical frameworks
of criticality in the brain, researchers have seen substantial progresses in
explaining various important neuroscience problems, including but not limited
to efficient cortical state transitions Fontenele . (2019), dynamic range
maximization in neural responses Kinouchi Copelli (2006); Shew . (2009), and
optimization of information transmission and representation Shew . (2011).
These advances have been comprehensively reviewed by existing works Shew Plenz
(2013); Chialvo (2010); Beggs (2007); Hesse Gross (2014); Cocchi . (2017);
Munoz (2018) and are no longer discussed in details in our review. The
benefits of studying brain criticality, as we have suggested, lay in the
possibility to analyze brain function characteristics with numerous
statistical physics theories relevant to brain criticality, such as directed
percolation Hinrichsen (2000); Lübeck (2004), conserved directed percolation
Bonachela Muñoz (2008); Bonachela . (2010), and dynamical percolation theories
Steif (2009); Bonachela . (2010). These theories characterize the brain as a
physical system with avalanche behaviors, enabling researchers to analyze
various propagation, synchronization, and correlation properties of neural
dynamics (e.g., continuous phase transitions). These properties intrinsically
shape neural information processing (e.g., encoding Bertschinger Natschläger
(2004); Haldeman Beggs (2005); Kinouchi Copelli (2006); XR. Wang . (2011);
Boedecker . (2012); Shew Plenz (2013); Del Papa . (2017), transmission Shew .
(2011), and memory Krotov Hopfield (2020); Haldeman Beggs (2005)) and can be
readily recorded in neuroscience experiments. Therefore, the non-equilibrium
dynamic processes and potential criticality defined by statistical physics
theories are highly applicable to characterizing brain functions. As we have
discussed in Fig. 2, researchers can consider diverse brain criticality
phenomena in neural dynamics by defining different control (e.g., the balance
between excitatory and inhibitory neurons Poil . (2012); Hardstone . (2014))
and order (e.g., active neuron density Dalla Porta Copelli (2019)) parameters,
corresponding to multifarious biological mechanisms underlying neural dynamics
(e.g., synaptic depression Levina . (2007)). Meanwhile, the definition of
neural avalanches can flexibly change from neural spikes, local field
potentials, to global cortical oscillations. The flexibility of brain
criticality and neural avalanche definitions enables researchers to analyze
different functional properties on distinct organizational levels in the
brain.

The limited theoretical foundations of brain criticality in the brain,
however, have become irreconcilable with their increasingly widespread
applications. Although the analytic theories of brain criticality have solid
physics backgrounds, they needlessly become black boxes for neuroscientists in
practice. On the one hand, the details of brain criticality theory frequently
experience immoderate neglecting in neuroscience studies. On the other hand,
to our best knowledge, there is no accessible and systematic introduction of
the statistical physics foundations of brain criticality in the terminology of
neuroscience yet. These obstacles severely impede neuroscientists from
comprehensively understanding brain criticality, eventually motivating us to
present this review. When we turn to bridging between brain criticality
theories and experiments, one can find non-negligible gaps separating between
theories and experiments. Although numerous biological factors (e.g., neural
plasticity Levina . (2007, 2009); De Arcangelis . (2006), membrane potential
leakage Millman . (2010); Levina . (2007); Rubinov . (2011); Stepp . (2015),
retro-synaptic signals Hernandez-Urbina Herrmann (2017), spatial heterogeneity
Moretti Munoz (2013); Girardi-Schappo . (2016), and refractory period
Williams-García . (2014); Fosque . (2021)) have been considered in brain
criticality characterization, existing theories more or less suffer from
deviations from actual neural system properties. For instance, the
requirements of conserved neural dynamics and an infinite time scale
separation between the dissipation and drive processes required by SOC may not
be biologically realistic Munoz (2018). The implicit requirement of a
sufficiently large system size by the mean field theories of brain criticality
may not always be satisfied during neural avalanche recording, implying non-
negligible finite size effects Girardi-Schappo (2021). Meanwhile, precisely
verify the existence of a detailed type of brain criticality (e.g., confirm
the actual universality class) in empirical neural data is principally
infeasible. As we have explained, the common criteria used for brain
criticality hypothesis verification, such as neural avalanche exponents
Clauset . (2009); Bauke (2007); Yu . (2014); Marshall . (2016); Deluca Corral
(2013), scaling relation Lübeck Heger (2003); Lübeck (2004), universal
collapse shape Marshall . (2016); Papanikolaou . (2011); Laurson . (2009);
Bhattacharjee Seno (2001), and slow decay of auto-correlation Pausch . (2020);
Wilting Priesemann (2019), are derived according to directed percolation
theory under mean field assumptions. Among four types of brain criticality in
absorbing phase transitions, only OC originally belongs to directed
percolation universality class while qC, SOC, and SOqC conditionally exhibit
directed percolation behaviours. In most cases, one can only verify if the
brain is plausibly at criticality (e.g., whether neural avalanches obey
universal collapse and have the power-law exponents that satisfy the scaling
relation). When observed neural avalanche exponents depart from their mean
field approximation results but still satisfy the scaling relation, there may
exist an OC phenomenon affected by non-mean-field factors (e.g., network
topology Girardi-Schappo (2021)) or exist a certain qC, SOC, or SOqC
phenomenon caused by diverse mechanisms. Additional information of neural
dynamics properties is inevitably required to determine the category belonging
of the hypothesized brain criticality, which poses daunting challenges to
neuroscience experiment designs. Moreover, the potential validity issues of
applying the theoretical tools derived from directed percolation theory to
verify brain criticality in synchronous phase transitions deserve special
attention (see Dalla Porta Copelli (2019) for similar opinions). It remains
controversial if absorbing and synchronous phase transitions robustly share
specific features (see reported similarities Di Santo . (2018); Girardi-
Schappo . (2021); Buendía . (2021); Fontenele . (2019) and differences
Fontenele . (2019); Buendía . (2021); Girardi-Schappo . (2021)). Any
speculated relations between these two kinds of critical phenomena should be
tested with cautions. Furthermore, statistic techniques to estimate and verify
brain criticality from empirical data are yet imperfect. The estimation of
some properties of neural avalanches is error-prone in practice and may lead
to serious validity issues. Although we suggest compromised solutions to these
issues, more optimal approaches are required in future studies.

##  7 Brain criticality: suggestions of future direction

We pursue that this review not only summarizes latest developments in the
field of studying criticality in the brain, but also serves as a blueprint for
further explorations. Below, we offer concrete recommendations of future
directions.

First, we suggest researchers to carefully rethink the theoretical foundations
of criticality in the brain. Immoderately omitting these foundations in
neuroscience needlessly muddies an already complex scientific field and leads
to potential validity issues. While we have presented a self-contained
framework of brain criticality to characterize neural dynamics as a physical
system with avalanches, plentiful details are uncovered in this article (e.g.,
the Landau–Ginzburg theory Di Santo . (2018)) because the statistical physics
theories of brain criticality are essentially grand. We recommend researchers
to further improve our work and explore a more accessible and systematic
reformulation of related physics theories, such as directed percolation,
conserved directed percolation, dynamic percolation, non-equilibrium dynamics,
in the context of neuroscience. Moreover, we note that these theories are not
initially proposed for brain analysis. It is normal to see gaps between these
theories and real situations of the brain. We urge researchers to develop new
variants of criticality formalism that is more applicable to the brain or even
explore new universality classes of continuous phase transitions.

Second, neuroscience is in urgent need of new physical theories and statistic
techniques to bridge between brain criticality hypotheses and experiments.
Although existing theories and techniques have become increasingly booming and
covered most of the pivotal details of brain criticality, there remain various
limitations as we have suggested. Specifically, we suggest five potential
directions to resolve these problems: (1) combine brain criticality theories
with large-scale neural dynamics recording or computation to include more
realistic biological details into brain criticality theories and establish a
closer connection with experimental observations; (2) try to summarize,
standardize, and subdivide these theories according to the concrete biological
meanings of brain criticality phenomena, prerequisites of model definitions,
and scopes of application. Try to avoid abusing or misusing of different brain
criticality theories; (3) develop open-source toolboxes of theoretical models
and statistic techniques to routinize brain criticality analysis in
neuroscience studies (one can see existing efforts to achieve this objective
Marshall . (2016)); (4) establish open-source, multi-species, and large-scale
data sets of neural dynamics recorded from both critical and non-critical
brains. Validate different statistic techniques of brain criticality
estimation and testing on these data sets and, more importantly, confirm
appropriate baselines to define the criteria of brain criticality
identification (see notable contributions in Girardi-Schappo (2021)); (5)
explore new non-equilibrium statistical physics theories for synchronous phase
transitions or analytically verify the theoretical validity of directed
percolation formulation of synchronous phase transitions.

Third, parallel to neuroscience, the discoveries of critical phenomena in
other learning and computation systems also merit attention. Learning or
computing at the edge of chaos has been proven as a mechanism to optimize the
performance of learners (e.g., recurrent neural networks Bertschinger
Natschläger (2004)). The well-known residual connections can control the
performance degradation of artificial neural networks because they enable
networks to self-organize to criticality between stability and chaos to
preserve gradient information flows Yang Schoenholz (2017). It is recently
demonstrated that any artificial neural network generally self-organizes to
criticality during the learning process Katsnelson . (2021). In the future, it
would be interesting to explore whether information processing processes in
brains and artificial neural networks can be universally characterized by a
unified criticality theory.

Overall, we anticipate the potential of well-validated studies of criticality
in the brain to greatly deepen our understanding of neural dynamics
characteristics and their roles in neural information processing. Laying solid
theoretical foundations of studies is the most effective and indispensable
path to contributing to this booming research area.

###### Acknowledgements.

Correspondence of this paper should be addressed to P.S. Author Y.T.
conceptualizes the idea, develops theoretical frameworks, and writes the
manuscript. Authors Z.R.T. and H.D.H. contribute equally to developing
mathematics of neural avalanche exponents. Author G.Q.L. contributes to
reviewing and revising the manuscript. Authors A.H.C., Y.K.Q., and K.Y.W.
contribute equally to the proof reading of Langevin formulation. Author C.C.
contributes to literature collection and summarizing. Author P.S. contributes
to idea conceptualization, manuscript writing, and project supervision. This
project is supported by the Artificial and General Intelligence Research
Program of Guo Qiang Research Institute at Tsinghua University (2020GQG1017)
as well as the Tsinghua University Initiative Scientific Research Program.
Authors are grateful for discussions and assistance of Drs. Yaoyuan Wang and
Ziyang Zhang from the Laboratory of Advanced Computing and Storage, Central
Research Institute, 2012 Laboratories, Huawei Technologies Co. Ltd., Beijing,
100084, China.

## References

  * Abbott (2008) abbott2008theoreticalAbbott, LF.  2008\.  Theoretical neuroscience rising Theoretical neuroscience rising.  Neuron603489–495.
  * Acebrón . (2005) acebron2005KuramotoAcebrón, JA., Bonilla, LL., Pérez Vicente, CJ., Ritort, F. Spigler, R.  2005Apr.  The Kuramoto model: A simple paradigm for synchronization phenomena The kuramoto model: A simple paradigm for synchronization phenomena.  Rev. Mod. Phys.77137–185.  https://link.aps.org/doi/10.1103/RevModPhys.77.137 10.1103/RevModPhys.77.137
  * Antonopoulos (2016) antonopoulos2016dynamicAntonopoulos, CG.  2016\.  Dynamic range in the C. elegans brain network Dynamic range in the c. elegans brain network.  Chaos: An Interdisciplinary Journal of Nonlinear Science261013102.
  * Arenas . (2008) arenas2008synchronizationArenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y. Zhou, C.  2008\.  Synchronization in complex networks Synchronization in complex networks.  Physics reports469393–153.
  * Atick (1992) atick1992couldAtick, JJ.  1992\.  Could information theory provide an ecological theory of sensory processing? Could information theory provide an ecological theory of sensory processing?  Network: Computation in neural systems32213–251.
  * Bak (2013) bak2013natureBak, P.  2013\.  How nature works: the science of self-organized criticality How nature works: the science of self-organized criticality.  Springer Science & Business Media.
  * Bak Sneppen (1993) bak1993punctuatedBak, P. Sneppen, K.  1993\.  Punctuated equilibrium and criticality in a simple model of evolution Punctuated equilibrium and criticality in a simple model of evolution.  Physical review letters71244083.
  * Bak . (1987) bak1987selfBak, P., Tang, C. Wiesenfeld, K.  1987\.  Self-organized criticality: An explanation of the 1/f noise Self-organized criticality: An explanation of the 1/f noise.  Physical review letters594381.
  * Baldassarri . (2003) baldassarri2003averageBaldassarri, A., Colaiori, F. Castellano, C.  2003\.  Average shape of a fluctuation: Universality in excursions of stochastic processes Average shape of a fluctuation: Universality in excursions of stochastic processes.  Physical review letters906060601.
  * Barlow . (1961) barlow1961possibleBarlow, HB. .  1961\.  Possible principles underlying the transformation of sensory messages Possible principles underlying the transformation of sensory messages.  Sensory communication101.
  * Bassett . (2010) bassett2010efficientBassett, DS., Greenfield, DL., Meyer-Lindenberg, A., Weinberger, DR., Moore, SW. Bullmore, ET.  2010\.  Efficient physical embedding of topologically complex information processing networks in brains and computer circuits Efficient physical embedding of topologically complex information processing networks in brains and computer circuits.  PLoS computational biology64e1000748.
  * Bauke (2007) bauke2007parameterBauke, H.  2007\.  Parameter estimation for power-law distributions by maximum likelihood methods Parameter estimation for power-law distributions by maximum likelihood methods.  The European Physical Journal B582167–173.
  * Beggs (2007) beggs2007buildBeggs, JM.  2007\.  How to build a critical mind How to build a critical mind.  Nature Physics312835–835.
  * Beggs Plenz (2003) beggs2003neuronalBeggs, JM. Plenz, D.  2003\.  Neuronal avalanches in neocortical circuits Neuronal avalanches in neocortical circuits.  Journal of neuroscience233511167–11177.
  * Beggs Timme (2012) beggs2012beingBeggs, JM. Timme, N.  2012\.  Being critical of criticality in the brain Being critical of criticality in the brain.  Frontiers in physiology3163.
  * Bell Sejnowski (1997) bell1997independentBell, AJ. Sejnowski, TJ.  1997\.  The “independent components” of natural scenes are edge filters The “independent components” of natural scenes are edge filters.  Vision research37233327–3338.
  * Berger Zhou (2014) berger2014kolmogorovBerger, VW. Zhou, Y.  2014\.  Kolmogorov–smirnov test: Overview Kolmogorov–smirnov test: Overview.  Wiley statsref: Statistics reference online.
  * Bertschinger Natschläger (2004) bertschinger2004realBertschinger, N. Natschläger, T.  2004\.  Real-time computation at the edge of chaos in recurrent neural networks Real-time computation at the edge of chaos in recurrent neural networks.  Neural computation1671413–1436.
  * Betzel . (2016) betzel2016generativeBetzel, RF., Avena-Koenigsberger, A., Goñi, J., He, Y., De Reus, MA., Griffa, A.others  2016\.  Generative models of the human connectome Generative models of the human connectome.  Neuroimage1241054–1064.
  * Betzel Bassett (20171) betzel2017generativeBetzel, RF. Bassett, DS.  20171\.  Generative models for network neuroscience: prospects and promise Generative models for network neuroscience: prospects and promise.  Journal of The Royal Society Interface1413620170623.
  * Betzel Bassett (20172) betzel2017multiBetzel, RF. Bassett, DS.  20172\.  Multi-scale brain networks Multi-scale brain networks.  Neuroimage16073–83.
  * Betzel . (2018) betzel2018diversityBetzel, RF., Medaglia, JD. Bassett, DS.  2018\.  Diversity of meso-scale architecture in human and non-human connectomes Diversity of meso-scale architecture in human and non-human connectomes.  Nature communications911–14.
  * Bhattacharjee Seno (2001) bhattacharjee2001measureBhattacharjee, SM. Seno, F.  2001\.  A measure of data collapse for scaling A measure of data collapse for scaling.  Journal of Physics A: Mathematical and General34336375.
  * Blythe McKane (2007) blythe2007stochasticBlythe, RA. McKane, AJ.  2007\.  Stochastic models of evolution in genetics, ecology and linguistics Stochastic models of evolution in genetics, ecology and linguistics.  Journal of Statistical Mechanics: Theory and Experiment200707P07018.
  * Boedecker . (2012) boedecker2012informationBoedecker, J., Obst, O., Lizier, JT., Mayer, NM. Asada, M.  2012\.  Information processing in echo state networks at the edge of chaos Information processing in echo state networks at the edge of chaos.  Theory in Biosciences1313205–213.
  * Bonachela . (2010) bonachela2010selfBonachela, JA., De Franciscis, S., Torres, JJ. Munoz, MA.  2010\.  Self-organization without conservation: are neuronal avalanches generically critical? Self-organization without conservation: are neuronal avalanches generically critical?  Journal of Statistical Mechanics: Theory and Experiment201002P02015.
  * Bonachela Muñoz (2008) bonachela2008confirmingBonachela, JA. Muñoz, MA.  2008\.  Confirming and extending the hypothesis of universality in sandpiles Confirming and extending the hypothesis of universality in sandpiles.  Physical Review E784041102.
  * Bonachela Munoz (2009) bonachela2009selfBonachela, JA. Munoz, MA.  2009\.  Self-organization without conservation: true or just apparent scale-invariance? Self-organization without conservation: true or just apparent scale-invariance?  Journal of Statistical Mechanics: Theory and Experiment200909P09009.
  * Breskin . (2006) breskin2006percolationBreskin, I., Soriano, J., Moses, E. Tlusty, T.  2006\.  Percolation in living neural networks Percolation in living neural networks.  Physical review letters9718188102.
  * Buendía, Di Santo . (2020) buendia2020feedbackBuendía, V., Di Santo, S., Bonachela, JA. Muñoz, MA.  2020\.  Feedback mechanisms for self-organization to the edge of a phase transition Feedback mechanisms for self-organization to the edge of a phase transition.  Frontiers in physics8333.
  * Buendía, di Santo . (2020) buendia2020selfBuendía, V., di Santo, S., Villegas, P., Burioni, R. Muñoz, MA.  2020\.  Self-organized bistability and its possible relevance for brain dynamics Self-organized bistability and its possible relevance for brain dynamics.  Physical Review Research21013318.
  * Buendía . (2021) buendia2021hybridBuendía, V., Villegas, P., Burioni, R. Munoz, MA.  2021\.  Hybrid-type synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together Hybrid-type synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together.  Physical Review Research32023224.
  * Bullmore Sporns (2012) bullmore2012economyBullmore, E. Sporns, O.  2012\.  The economy of brain network organization The economy of brain network organization.  Nature reviews neuroscience135336–349.
  * Burns Webb (1976) burns1976spontaneousBurns, BD. Webb, A.  1976\.  The spontaneous activity of neurones in the cat’s cerebral cortex The spontaneous activity of neurones in the cat’s cerebral cortex.  Proceedings of the Royal Society of London. Series B. Biological Sciences1941115211–223.
  * Capolupo . (2013) capolupo2013dissipationCapolupo, A., Freeman, WJ. Vitiello, G.  2013\.  Dissipation of ‘dark energy’by cortex in knowledge retrieval Dissipation of ‘dark energy’by cortex in knowledge retrieval.  Physics of life reviews10185–94.
  * Carandini Stevens (2004) carandini2004amplificationCarandini, M. Stevens, C.  2004\.  Amplification of trial-to-trial response variability by neurons in visual cortex Amplification of trial-to-trial response variability by neurons in visual cortex.  PLoS biology29e264.
  * Cardin (2019) cardin2019functionalCardin, JA.  2019\.  Functional flexibility in cortical circuits Functional flexibility in cortical circuits.  Current opinion in neurobiology58175–180.
  * Carvalho . (2021) carvalho2021subsampledCarvalho, TT., Fontenele, AJ., Girardi-Schappo, M., Feliciano, T., Aguiar, LA., Silva, TP.Copelli, M.  2021\.  Subsampled directed-percolation models explain scaling relations experimentally observed in the brain Subsampled directed-percolation models explain scaling relations experimentally observed in the brain.  Frontiers in neural circuits83.
  * Chialvo (2010) chialvo2010emergentChialvo, DR.  2010\.  Emergent complex neural dynamics Emergent complex neural dynamics.  Nature physics610744–750.
  * Chiappalone . (2003) chiappalone2003networksChiappalone, M., Vato, A., Marcoli, M., Davide, F., Martinoia, S. .  2003\.  Networks of neurons coupled to microelectrode arrays: a neuronal sensory system for pharmacological applications Networks of neurons coupled to microelectrode arrays: a neuronal sensory system for pharmacological applications.  Biosensors and Bioelectronics185-6627–634.
  * Clauset . (2009) clauset2009powerClauset, A., Shalizi, CR. Newman, ME.  2009\.  Power-law distributions in empirical data Power-law distributions in empirical data.  SIAM review514661–703.
  * Cocchi . (2017) cocchi2017criticalityCocchi, L., Gollo, LL., Zalesky, A. Breakspear, M.  2017\.  Criticality in the brain: A synthesis of neurobiology, models and cognition Criticality in the brain: A synthesis of neurobiology, models and cognition.  Progress in neurobiology158132–152.
  * Cohen Kohn (2011) cohen2011measuringCohen, MR. Kohn, A.  2011\.  Measuring and interpreting neuronal correlations Measuring and interpreting neuronal correlations.  Nature neuroscience147811–819.
  * Collell Fauquet (2015) collell2015brainCollell, G. Fauquet, J.  2015\.  Brain activity and cognition: a connection from thermodynamics and information theory Brain activity and cognition: a connection from thermodynamics and information theory.  Frontiers in psychology6818.
  * Dalla Porta Copelli (2019) dalla2019modelingDalla Porta, L. Copelli, M.  2019\.  Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic M/EEG results Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic m/eeg results.  PLoS computational biology154e1006924.
  * David Friston (2003) david2003neuralDavid, O. Friston, KJ.  2003\.  A neural mass model for MEG/EEG:: coupling and neuronal dynamics A neural mass model for meg/eeg:: coupling and neuronal dynamics.  NeuroImage2031743–1755.
  * de Andrade Costa . (2015) de2015cande Andrade Costa, A., Copelli, M. Kinouchi, O.  2015\.  Can dynamical synapses produce true self-organized criticality? Can dynamical synapses produce true self-organized criticality?  Journal of Statistical Mechanics: Theory and Experiment20156P06004.