slightly subcritical state [151, 152] (but see [39]) with occasional
incursions into a supercritical state [153, 154, 36]. This suggests that the
purpose itself that defines perception (reducing the mismatch between
environmental and internal states) is what induces criticality in perceptive
systems [155], similarly to SOC or quasi-SOC systems. In PC, the exogenous
input rate is infinitesimal relative to discharge rates (i.e. the next input
waits for all chain elements to “decide” whether to discharge ot not) and
energy is bulk conserved (within the chain). This entails that, if we neglect
threshold fluctuations (FTA, Section IV.4), PC is analogous to a SOC model
operating near criticality [2]. The finite number of levels $n_{l}$ together
with the boundary (top level) dissipation induce finite-size scaling. True
scale invariance would be achieved in the limit of infinite number of levels
[2]. Conversely, allowing exogenous input injection before the cascade ends
would render the system supercritical (cf. sandpile model with finite driving
rate [2]). When threshold fluctuations are allowed, the system hovers around
the critical point, in a stretched critical phase grounded on threshold
adaptability, similarly to how neural models with dynamic synapses manifest a
Griffiths phase or stretched criticality range [156, 157, 158, 159]. Finally,
introducing surprisal (energy) bulk dissipation would lead to a pseudo-
critical phenomenon called self-organized quasi-criticality [3, 154].

Empirically, the hypothesis that the brain operates near criticality [13, 160]
has also gained much (if not conclusive) support [161, 17, 162]. Power-law
temporal correlations have been found in EEG [23, 24] and MEG signals [163]
across a broad range of frequencies [120, 164], in functional neuroimaging
data [165], and real neuronal avalanches (stereotyped population-size firing
events) have been repeatedly observed in animals in multielectrode arrays, in
vitro in acute cortical slices and organotypic cultures [15], and in vivo
[129, 17], displaying spatio-temporal and neuronal group size scale invariance
[166].

The continuous-input Bak-Tang-Wiesenfeld lattice sandpile model [4, 46] shares
many traits with PC: discrete pulse density peaks (which may be related to
DSI, see section VI.4), energy conservation, and power law spectra for
avalanches and stored energy, and Lorentzian temporal spectra [46]. When
energy is injected from one end and dissipation is added, long range
correlations are quenched, yet temporal spectra become $\sim\nu^{-1}$. This
results from a linear superposition of the (now) uncorrelated independent
lattice elements, which are composed of exponentials with decay times that
depend on their location in the lattice. This is like the mechanism that
underpins the PC power spectrum, where uncorrelated energy (surprisal)
accumulates sequentially at increasingly slow levels, setting up independently
excitable elements of different intrinsic timescales (also to McWhorterś model
of shot noise [167]). Their common key features are separation of timescales
(achieved by decorrelating dissipation in [46] and by hierarchical design in
PC) and unidirectional propagation of energy. In another similar model, input
is injected into the center of a square lattice, constituted by elements
characterized by threshold firing, refractory period, and synaptic plasticity
[40]. This model displays power law distributions of avalanches (with exponent
$\approx-1.2$) and spectrum ($\approx-.8$) presumably because its built-in
plasticity in conjunction with random initial conditions carves a complex
landscape of synapses with a near-critical branching parameter; however, it
has many parameters. Adjusting the ratio of excitation and inhibition can
yield similar results [39], which furnishes an explanation for the wide range
of temporal scaling exponents found in the literature.

###  VI.8 Neuronal avalanches

Cortical neurons are likely to be segregated into layers by their duties as
predictive coding hardware. Specifically, prediction error (pulse) neurons
would be allocated to supragranular and state neurons to infragranular
cortical layers [19, 168]. Thus, we can relate PC pulses to prediction error
spikes in supragranular cortical layers. Neuronal avalanches presumably
originate from activity in superficial (supragranular) layers 2 and 3 that
propagates across multiple cortical columns in all regions of the cortex [169,
129, 170]. This can be taken one more step further, by relating pulses to
local field potentials (LFPs) [171, 32] (the electric potentials recorded in
extracellular space in brain tissue), which are the signals picked up by
micro-electrodes and by EEG —the latter after filtering and diffusion across
head tissues.

Somewhat atypically, here avalanche denotes a volley of pulses triggered
within a singe time step, and its size is proportional to the sum of its
constituting pulse sizes. The mixture density $f_{\mathcal{E}_{mix}}$ (Section
III.3) then corresponds to recordings with high enough temporal (or spatial,
if levels were spatially segregated) resolution to tell apart each individual
pulse or avalanche within a single time step —akin to high resolution
multielectrode recordings in electrophysiology. In contrast, the sum density
$f_{\Sigma\mathcal{E}}$ results from registering all the pulses or avalanches
generated in one time step as their linear superposition, i.e. lumping them
together [172, 116, 173]. The sum density is analogous to the linear
superposition of temporally overlapping post-synaptic potentials on micro-
electrodes or LFPs on surface electrodes [117], where temporal subsampling
[151, 173] or limited spatial resolution may blend separate avalanches into a
single time-varying medley. However, note that the sequential nature of
avalanches is lost in our system due to time discretization, which merges all
pulses within a time step, and the lack of spatial extent within levels. At
any rate, recall that in our system pulses (and avalanches) cannot be naively
interpreted as actual electrical potentials or spike counts (Section VI.5);
for example, the analogous to the branching parameter [15] would be the gain
$g$, but the lack of a within-level spatial dimension precludes sensible
comparison.

Empirical estimates from microelectrode arrays suggest that neuronal avalanche
size and duration bear scaling exponents 3/2 and 2 respectively [15, 169, 163]
(cf. 3/2 with 1.77 in Fig. 7; or with 1.64 in Fig. 16). The number 3/2 matches
the avalanche size scaling exponent in a model of globally coupled (mean-
field) threshold elements with random stimulation location [37] similar to the
Olami-Feder-Christensen earthquake model [7], mean-field self-organized
branching processes [174], many flavors of mean-field zero-temperature Ising
models [175], mean-field directed percolation and mean-field dynamical
percolation [176]. This is notable because the critical behavior of directed
percolation belongs to a robust universality class that is widespread in
systems exhibiting cascading transitions from active to absorbing phases [174,
176, 177, 178]. The closely related (but different for low dimensions) Manna
class or conserved directed percolation class is particularly relevant to SOC
and is associated to conserving dynamics [178, 3]. However, the significance
of 3/2 is moot because there are similar systems with absorbing-state phase
transitions that have different scaling exponents [176], mean-field
approximations are typically incorrect for low dimensional systems (e.g.
2-dimensional sandpile model [4, 2], but see [179]) and the exponent estimate
$\approx 1.5\textit{--}2.1$ may depend on the spacing between electrodes in
the array [160]. Further, avalanches occur in short timescales of 1–100 ms
[180, 26], but electrographic long-range temporal correlations timescales are
in the range 1–1000 s, so it is plausible that the 3/2 exponent be restricted
to spreading of activity within (isofrequency) levels.

The densities of pulse activities and their avalanches (Fig. 7) are a priori
the variables most related to neural activity, but energy densities allow to
visualize the distribution of stored and relayed surprisal across levels (Fig.
16). Energy is also interesting because it is the only variable conserved in
PC, defines discharge rates, and has the physiological interpretation of
irreducible uncertainty. The scaling exponents of pulse (-1.95 for mixture,
-1.64 for sum density) and subthreshold energy densities (-.83 for mixture,
-.29 for sum density) are shown in Fig. 16. Both variables display power law
or DSI behavior within their range, with a right tail roll-off more
conspicuous for subthreshold activities caused by the thresholding.

Figure 16: Density estimates of pulse and subthreshold activity energies for
$n_{l}=7$ levels. First from above: pulse energy mixture density
$f_{\mathcal{E}_{mix}^{2}}$ (black), and individual levels
$f_{\mathcal{E}_{l}^{2}}$ (light grey) for all $l=1\ldots 7$. Second: sum of
pulse energies $f_{\Sigma\mathcal{E}^{2}}$. Third: subthreshold activity
mixture $f_{\mathcal{A}_{mix}^{2}}$ and individual level energy
$f_{A_{l}^{2}}$ densities. Fourth: sum of subthreshold activity energies
$f_{\Sigma\mathcal{A}^{2}}$. Pulse energy slopes (dashed, top -1.95, bottom
-1.64) were fitted in the least squares sense in a range that excluded finite-
size effects:
$[\langle\mathcal{E}_{1}^{2}\rangle,\langle\mathcal{E}_{7}^{2}\rangle]$ for
$f_{\mathcal{E}_{mix}^{2}}$ and
$[\langle\mathcal{E}_{1}^{2}\rangle,\sum_{i=1}^{7}\langle\mathcal{E}_{i}^{2}\rangle]$
for $f_{\Sigma\mathcal{E}^{2}}$. Subthreshold activity slopes (dashed, top
-.83, bottom -.29) were fitted in the range $[10^{-2},10^{2}]$.

###  VI.9 Autocorrelated and non-Gaussian inputs

For inferential systems that are good enough regulators [18, 19], the Hurst
exponent of the input stochastic process is $H\approx 0.5$, so $k\approx 1.5$
(Section V.2). This reflects that good regulators expect only unbiased,
temporally uncorrelated surprisal —which is why Brownian motion is a central
component of standard models of decision making in neuroscience [181].
Autocorrelated input may occur when the system is not a good enough model of
its environment or input. Then predictions may be biased, which in PC
manifests as error signal distribution being biased at all levels. For $H\neq
0.5$, there is a transition at $H\approx 0.75$, beyond which the subthreshold
density splits into two disjoint modes (section V.2).

The threshold adaptation rule of Eq. 4 allows PC variables to converge (if
there are enough levels) to the rescaled recurrent configuration state
(Section V; Fig. 9) if the input density is not too heavy-tailed: the tail
should have a finite first moment (mean). The reason is that iteratively
applying Eq. 4 converges to the mean (Section VI.5). Approximating the tail
with a Pareto distribution $\sim x^{-\alpha}$ (Appendix F), this is satisfied
if $\alpha>2$. The thresholds diverge if the input density mean is undefined
(e.g. Cauchy).

###  VI.10 Subthreshold bistability

Remarkably, the bimodal PC subthreshold activity densities (Fig. 2, 9, in
black, note the probability mass $g$ associated to zero is not shown) are
similar to in vivo intracellular spontaneous membrane potential shift
recordings (anesthetized rat corticostriatal, Fig. 5 of [182]; neo-striatal,
Fig. 3 of [183], reproduced in Fig. 17, left; cat neocortical pyramidal
neurons, Figs. 5A, 6A of [184], reproduced in Fig. 17, right) after equating
zero to resting membrane potential and $\langle\Theta\rangle$ to threshold
potential. Note that hyperpolarizing currents lead to near-unimodal densities,
whereas depolarizing currents lead to bimodal densities (Fig. 17), in a manner
that strongly resembles the PC densities induced by driving input of Hurst
parameter $H<0.5$ and $H>0.5$ respectively (Fig. 12). In agreement with [185],
for PC only the sensory unit lacks a bimodal distribution, and activity
densities with positively (negatively) autocorrelated exogenous input (Fig.
12) resemble membrane potential histograms under depolarizing
(hyperpolarizing) currents [183, 184].

Figure 17: Effect of polarizing current on membrane potential density
experimental recordings. Compare with Figs. 2 and 12. Left: In vivo
intracellular recordings of spontaneous activity of neo-striatal spiny cells
displaying Up-Down modes. -0.5 nA denotes a strong hyperpolarizing current, 0
nA no injected current and 0.75 nA a depolarizing current. Reused from [183].
Copyright 1996 Society for Neuroscience. Right: Cat deep pyramidal neuron
recordings under barbiturate anesthesia (top) and ketamine-xylazine anesthesia
(bottom) under hyperpolarizing and depolarizing currents. Reused with
permission from [184].

The membrane potential of cortical neurons responds to synaptic input in a
two-state manner [186]. Some neurons toggle between one very hyperpolarized
(Down state), and one more depolarized (Up state) membrane potential [187],
where fluctuations around the Up state are of higher amplitude and follow
stimulation [186], whereas the Down state is relatively free of noise [187].
Similarly, waiting times between avalanches have a bimodal distribution [116].
The possible coexistence of both SOC and both continuous and discontinuous
phase transitions via a cusp bifurcation has been demonstrated by modeling the
electrochemical dynamics of synaptic terminals with a few parameters [47],
which provides an explanation for Up-Down states. This behavior has been
construed as an effect of homeostatic regulation [188]. Bistability (in alpha
rhythms) has also been proposed to ensue from a (fifth-order branch)
subcritical bifurcation in corticothalamic activity dynamics, that crucially
is noise-driven [189]. The PC bimodal density of subthreshold activity in both
feedforward and rescaled recurrent loop configurations provides a simple
explanation for Up-Down states as the steady state expression of excitable
elements that accumulate synaptic input while adjusting the firing threshold
to the incoming pulse sizes, which in turn is a manifestation of the trade-off
between metabolic power consumption and performance.

##  VII Caveats and prospective research paths

The bare bones design of PC unavoidably brings about shortfalls. For instance,
we only examined one adaptive threshold rule that weights mean squared error,
but there are other biologically plausible options (Section VI.5). The simple
daisy chain linear topology of PC precludes modeling intra-level effects (e.g.
lateral connections between neurons or microcolumns in the same level). We
likened activity energy or surprisal to irreducible uncertainty, thereby
neglecting fluctuations derived from error signals that are successfully
suppressed, which are likely to account for a substantial fraction of activity
and surprisal in full-blown good regulators. Although simplifications warrant
caution in assessing the applicability of modeling results, they may also
indicate new appealing research paths.

Crucially, we left out top-down effects, which are a fundamental aspect of
hierarchical Bayesian generative models, where the causal asymmetry between
top-down modulatory effects and bottom-up driving effects is pivotal to
inference and learning [28, 30]. These lead to back and forth volleys between
contiguous levels, and the effect of bottom-up signals would be modulated by
supraordinate priors at every level. By shifting the reference or zero value
of activity, top-down effects by themselves could suppress discharges (by
reconfiguring internal states to accurately predict input, e.g. eureka moment)
or encourage them (e.g. under flawed perceptual inference). A top-down signal
could set off a surge of activity in subordinate levels, with its size being
exponentially related to the rank of the emitting level (Section VI.6). An
expedient course for future work on top-down effects could be to explicitly
model top-down signals that suppress uncertainty in subordinate levels and to
split the time step into the number of hierarchy levels traversed. Potentially
this could inform the presence of local peaks commonly superimposed on M/EEG
power spectra [124] and untangle avalanche clusters due to temporal
subsampling [151, 173] (Section VI.8). Finally, allowing PC to accept input
before stimulus-evoked cascades end would render it (more) subcritical
(Section VI.7).

In PC, energy (surprisal) is forced to traverse the whole chain before leaving
the system. One interesting accessory feature would be fixing the top level
threshold and considering the surprisal quenching at the top level as caused
by the action of the system on the environment (recall that top level
surprisal cannot be suppressed by supraordinate levels). Then, the PC
feedforward configuration would become a metaphor for passive perception under
strong priors and occasional stimulus-driven action (akin to mind-wandering).
In addition, difficult or surprising situations could be modeled by increasing
the input variance, or even its tail shape. With only bottom-up signals there
are few alternatives, but when top-down signals are included many plausible
variations appear. For instance, although PC’s behavior suggests that the
hierarchical coupling of brain levels is essential to its critical behavior,
the switching between transient states characteristic of cognitive dynamics
[95, 190, 191] is not explicitly modeled. Allowing to suppress surprisal by
acting on the environment (as opposed to only reconfiguring internal states)
would provide an explicit mechanism that displays the same sort of itinerant
dynamics arising in SOC, synergetics [91], and living creatures [36]. The
dichotomy between perception and action is a fundamental design principle of
living systems that dictates the architecture of hierarchical brain models
through their intimate relationship [112]. However, the implementation of
action effectors does not arise from the top level: the timescale of the
faster goal-directed motor planning is roughly $.1\text{--}.5$ Hz for humans,
which (not casually) corresponds to the time window of the stream of
consciousness [192, 193, 194]. Hence, the top level associated with action is
lower than the top level associated with perception. In motor control and
free-energy formulations of behavior, action is computed via an inverse model
[195, 196, 112] that maps sensory prediction errors to actions. Thus, a
hierarchy level operating roughly in the theta band, consistent with large
theta-cycle avalanches triggering smaller and faster ones [116], is likely to
encompass the top level of the inverse model, which would be akin to a
“sandbox” used to simulate future courses of action. A candidate for
implementing the highest modules of the inverse model is the default mode
network. The resting brain evinces large-amplitude spontaneous $<$0.1 Hz
fluctuations in neuromagnetic signals that are temporally correlated across
functionally related areas (chiefly medial prefrontal cortex, posterior
cingulate cortex, and angular gyrus [197]) that are posited to constitute a
default mode network subserving long term planning, recollection [198, 199,
200, 201] and mood [202]. These notions suggest research pathways for how
itinerant dynamics emerge from the exploration of the intertwined action-
perception space of possibilities and sensory inputs, and how this relates to
the trade-off between metabolic power saving and performance.

The analysis of autocorrelated input (Section VI.9) hints that multifractal
analysis [203] could be a useful tool in assessing potential time-varying
multiplicative cascading processes underlying dynamics of PC aggregate
activity. This is because multifractal analysis can be used to break down a
time series of multiplicative interactions that transport information or
energy across multiple temporal scales into overlapping contributions derived
from a time-varying scaling exponent. Its application to theoretical
neuroscience entails “a local threshold dynamic between adjacent neurons that
leads to a dissipation of exogenous or afferent stimuli across multiple
spatial scales of the nervous system … This intermittent spatial dissipation …
defines the sensitivity and adaptability of the cognitive system to exogenous
perturbations” [204]. This would allow examining the degree of incoming pulses
autocorrelation —which are likely to play an important role in the operation
of PC-like architectures— at both level-wise and whole system scales.

##  VIII Conclusion

The perfusive cascade is a minimal model of a good regulator [18] that is
wholly stripped of the salient behavior of its generative model, but displays
dynamics characteristic of highly optimized systems that model a complex,
multi-scale environment. Surprisal transport through its stylized hierarchy is
near self-similar [203, 204], and leads to scaling laws for all its derived
variables. The state units that embody its model of the world are self-
dissimilar [48] and are not explicitly modeled, although they admit the
assumption that activities at all levels reflect the mismatch between incoming
pulses (derived from sensory input) and expectations about the sensory input
[30]. Its essential features are no adjustable parameters, a universal scaling
function in the limit of infinite levels (under weak conditions; Section
VI.9), trade-off between metabolic power consumption and performance, SOC or
autovitiation [36], a hierarchy of timescales [35, 133] where each level
represents brain regions characterized by timescale, and bistable subthreshold
fluctuations closely resembling intracellular recordings (Section VI.10). The
subthreshold activities and pulses have Lorentzian and flat temporal power
spectra, respectively; their energy represents surprisal or irreducible
uncertainty that is accumulated and relayed across levels. With a biologically
plausible refractory period ($T_{0}=5$ms) and hierarchy depth ($n_{l}=7$
[134]), the resulting PC scaling range of 0.8–85 Hz roughly matches the
physiologically relevant brain rhythms. In contrast to other similar models
that predict a wide range of scaling exponents [39, 40], the spectrum scaling
exponent ensues from the weighted linear combination of spectra from each
level, which can yield any $\beta\in[0,2]$— as opposed parameter tuning. These
weights are proportional to the number of neurons in each level, which
decreases exponentially while ascending the hierarchy (Section VI.6). Thus,
the temporal scaling patterns found in electrophysiological recordings could
be explained as weighted linear combinations of subthreshold activities and
pulses from different hierarchy levels. Finally, this scheme suggests some
testable conjectures, pertinent to experimental neuroscience: (1) prediction
error neurons (presumably in supragranular layers) that exhibit Up-Down states
exist all over the brain; (2) power spectra of M/EEG traces can display any
scaling exponent in the range $[-2,0]$ (depending on the task and recorded
brain area composition) within the frequency band where circuit level neural
dynamics unfold (which corresponds roughly to 0–100 Hz for a refractory period
of 5 ms); and (3) temporal power spectra scaling exponents are steeper in rest
than under stimulation, because (stimuli-driven) error signals have flatter
spectra.

###### Acknowledgements.

The article was supported by the HSE University Basic Research Program.

##  Appendix A  Analytic expressions for stationary densities

Let $w\equiv w^{+}=w^{-}$, the exogenous driving input be
$I\sim\mathcal{N}(0,1)$, $f_{A_{l}}$ be the steady state subthreshold activity
density, and post-pulse activity be defined as
$\tilde{A}_{l}=A_{l}+\mathcal{E}_{l-1}$ (with $\mathcal{E}_{0}\equiv I$).
Then, for $l=1\ldots n_{l}$:

$\displaystyle f_{\tilde{A}_{l}}$ $\displaystyle=$ $\displaystyle
f_{A_{l}}*f_{\mathcal{E}_{l-1}},$ (58) $\displaystyle f_{A_{l}}$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{i=0}f_{A_{l}^{(i)}}q_{l}^{(i)}\prod^{i-1}_{j=0}(1-q_{l}^{(j)})$
$\displaystyle=$ $\displaystyle
f_{A_{l}^{(0)}}q_{l}^{(0)}+(1-q_{l}^{(0)})(f_{A_{l}^{(1)}}q_{l}^{(1)}+$
$\displaystyle(1-q_{l}^{(1)})((f_{A_{l}^{(2)}}q_{l}^{(2)}+...)),$

where $*$ denotes convolution product 666Alternatively to A, $f_{A_{l}}$ could
be calculated as
$f_{A_{l}}=B_{\theta_{l}}[f_{\tilde{A}_{l}}]q_{l}^{(0)}+(1-q_{l}^{(0)})(B_{\theta_{l}}[f_{\tilde{A}_{l}^{(1)}}]q_{l}^{(1)}+(1-q_{l}^{(1)})(\ldots))$,
with $f_{\tilde{A}_{l}^{(t)}}=B_{\theta_{l}}[f_{\tilde{A}_{l}}^{(t-1)}]*I$,
$f_{\tilde{A}_{l}^{(0)}}=f_{\mathcal{E}_{l-1}}$, and
$f_{\tilde{A}_{l}^{(0)}}=I$.. This is reminiscent of a iterative method
devised to study irreversible growth processes (Fixed Scale Transformation
[206]). $f_{A_{l}^{(t)}}$ is the stationary density of $A_{l}$ at iteration
$t$ after the last reset to zero, and
$q_{l}^{(t)}=1-\int_{-\theta_{l}}^{\theta_{l}}f_{\tilde{A}_{l}^{(t)}}$.
Explicitly,
$f_{A_{l}^{(t)}}=B_{\theta_{l}}[f_{A_{l}^{(t-1)}}*f_{\mathcal{E}_{l-1}}]$; in
particular, $f_{A_{l}^{(0)}}=\delta(0)$, where $\delta$ is the Dirac delta;
$f_{A_{l}^{(1)}}=B_{\theta_{l}}[f_{\mathcal{E}_{l-1}}]$,
$f_{A_{1}^{(1)}}=B_{\theta_{1}}[I]$. For example, for $l=1$:

$\displaystyle f_{A_{1}}$ $\displaystyle=$
$\displaystyle\delta(0)q_{1}^{(0)}+(1-q_{1}^{(0)})(B_{\theta_{1}}[I]q_{1}^{(1)}+$
$\displaystyle(1-q_{1}^{(1)})(B_{\theta_{1}}[B_{\theta_{1}}[I]*I]q_{1}^{(2)}+\dots)).$

$q_{l}^{(0)}$ is the probability of firing consecutively (twice or more), and
in general $q_{l}^{(t)}$ the probability of firing $t+1$ iterations after the
last firing time. In particular, note that
$q_{l}^{(0)}=1-\int_{-\theta_{l}}^{\theta_{l}}f_{\mathcal{E}_{l-1}}$ would be
on average $1/2$ if activities were reset after every iteration (which would
lead to $\langle
A_{l}\rangle=\langle\Theta_{l}\rangle=\langle\mathcal{E}_{l-1}\rangle$).

In the limit $t\to\infty$, the probability $P(\mathcal{E}_{l}\neq
0)=P(A_{l}=0)=P(|\tilde{A}_{l}|>\Theta_{l})$ of unit $l$ firing at any time
$t$ is

$\displaystyle P(\mathcal{E}_{l}\neq 0)$ $\displaystyle=$ $\displaystyle
1-\sum^{\infty}_{i=0}q_{l}^{(i)}\prod^{i}_{j=0}(1-q_{l}^{(j)})$
$\displaystyle=$
$\displaystyle\sum^{\infty}_{i=0}[q_{l}^{(i)}]^{2}\prod^{i-1}_{j=0}(1-q_{l}^{(j)})$
$\displaystyle=$ $\displaystyle[q_{l}^{(0)}]^{2}+(1-q_{l}^{(0)})\cdot$
$\displaystyle([q_{l}^{(1)}]^{2}+(1-q_{l}^{(1)})(...)).$

A coarse estimation of $\langle\Theta_{1}\rangle$ using Eqs. 61, 58 and A up
to first order approximation —one iteration (till $i=1$), ignoring
$O(A_{1}^{(2)})$ terms and above— yields $\langle\Theta_{1}\rangle\approx
0.8197$, in reasonable approximation to the numerical simulations (Table 2).
Plugging into Eq. A the numerical estimate for $\langle\Theta_{1}\rangle$ and
truncating the series at $i=3$, we get
$f_{\mathcal{E}_{1}}(\varepsilon_{1}\neq 0)\approx 0.4336$, which face-
validates Eq. A. Table 2 displays the average activities, thresholds, and
pulses or error signals obtained through numerical simulations with $L=7$
units and $n_{t}=9\cdot 10^{5}$ iterations and exogenous input distributed as
white Gaussian noise.

Table 2: Numerical estimates of mean subthreshold activities, thresholds, and pulse absolute values. In theory $P(A_{l}=0)=P(|A_{l}+\mathcal{E}_{l-1}|>\Theta_{l}\;|\;\mathcal{E}_{l-1}\neq 0)=P(\mathcal{E}_{l}\neq 0\;|\;\mathcal{E}_{l-1}\neq 0)\equiv g_{l}:\forall l>1$ should hold because, given a uncorrelated exogenous input, $\mathcal{E}_{l}^{(t)}$ and $\mathcal{E}_{l}^{(t+1)}$ must also be uncorrelated. The discrepancies are due to sampling error. Parameters as in Table 1.  Level $l$ | $\langle|A_{l}|\rangle$ | $\langle\Theta_{l}\rangle$ | $\langle|\mathcal{E}_{l}|\rangle$ | $P(A_{l}=0)$
---|---|---|---|---
1 | 0.2309(2) | 0.8464 | 1.4479(8) | .4251(5)
2 | 0.4782(5) | 1.4769 | 2.2214(14) | .4501(5)
3 | 0.7408(9) | 2.2470 | 3.2759(29) | .4605(5)
4 | 1.0686(13) | 3.3212(2) | 4.8099(60) | .4696(5)
5 | 1.5937(19) | 4.8576(2) | 7.0698(131) | .4620(5)
6 | 2.2914(28) | 7.1114(4) | 10.3020(282) | .4707(5)
7 | 3.4148(42) | 10.4729(6) | 15.2187(606) | .4648(5)

Solving for $f_{\Theta_{l}}$ would be laborious —it requires specifying $w$ to
solve simultaneously the interlocked pair of Eqs. A and 4. However, by the
central limit theorem, in the steady state $f_{\Theta_{l}}$ converges in
distribution to a normal function, for $l=1\ldots n_{l}$. A larger $w$ induces
flatter $f_{\Theta_{l}}$ and $f_{A_{l}}$ densities. By Eq. 4, the expectation
of the threshold equals the expectation of post-pulse activity
$\langle\Theta_{l}\rangle=\langle|\tilde{A}_{l}|\rangle$. Thus, if only the
expectation $\langle\Theta_{l}\rangle$ is sought, the calculation of joint
densities can be eschewed by calculating $\langle\tilde{A}_{l}\rangle$ or
equivalently by finding the value that on average matches the signs of
threshold fluctuations, which numerically can be estimated by minimizing the
functional

$\langle\Theta_{l}\rangle=\langle|\tilde{A}_{l}|\rangle\approx\arg\min_{\theta}\left|\langle\tilde{A}_{l}-\theta\rangle_{\tilde{A}_{l}>\theta}-\langle\theta-\tilde{A}_{l}\rangle_{\tilde{A}_{l}<\theta}\right|,$
(61)

where
$\langle\tilde{A}_{l}-\theta\rangle_{\tilde{A}_{l}>\theta}=\int_{\theta}^{\infty}(x-\theta)f_{\tilde{A}_{l}}(x)dx$,
and
$\langle\theta-\tilde{A}_{l}\rangle_{\tilde{A}_{l}<\theta}=\int_{0}^{\theta}(\theta-x)f_{\tilde{A}_{l}}(x)dx$
and $f_{\tilde{A}_{l}}$ is approximated with its numerical distribution. This
approximation becomes exact in the limit $w\to 0$.

##  Appendix B  Brownian neurons: firing behavior at fine-grained timescales

We have assumed that integrator units discharge only at discrete time steps.
This trait mirrors refractoriness, or the property of being unresponsive to
stimuli for a period of time (refractory period). Neurons are excitable
elements with refractory period, but it is also interesting to study the case
of non-refractory elements or “Brownian neurons”, which are continuously
responsive and may take in continuous input. These elements are a kind of
stochastic model of neuron membrane voltage, similar to Stein’s or Ornstein-
Uhlenbeck neuron models [207, 208]. Here we examine excitable behavior as the
refractory period $T_{0}$ approaches zero, which is the scaling limit where a
random walk converges to the Wiener process. With the adaptive threshold rule
of Eq. 4, a refractory period $T_{0}=0$ under a continuous-time driving
process entails a zero threshold, so we will fix the threshold to a positive
value.

###  B.1 Firing rate

First we derive bounds and approximations for the discharge rates of the
sensory unit and for units located far enough from the input. The firing rate
of the sensory unit is the probability that it fires at any time
$P(A_{1}=0)=P(\mathcal{E}_{1}\neq 0)$, given by Eq. A. As mentioned in
Appendix A, the probability of discharging at iteration $t$ is the exceedance
probability
$q_{1}^{(t)}=1-\int_{-\theta_{1}}^{\theta_{1}}f_{\tilde{A}_{1}^{(t)}}$. Since
$f_{\tilde{A}_{1}}$ falls exponentially-like, $\theta_{1}$ is guaranteed to
lie not far from the median of $f_{\tilde{A}_{1}}$, so the firing rate is
guaranteed to remain close to $1/2$, which would not occur if
$f_{\tilde{A}_{1}}$ decayed as a power law [209]. Because of diffusion, this
probability would be lower (higher) if the activity $\tilde{A}_{1}$ were
distributed as white Gaussian noise $\mathcal{N}(0,1)=W^{(1)}$ (Wiener process
$W^{(t)}$). A tighter upper bound can be achieved by approximating
$\tilde{A}_{1}^{(t)}$ with the sum of Gaussian input with the activity
attained in the limit of infinitely many iterations without firing, namely
$\tilde{A}_{1}^{(t)}\sim\mathcal{N}(0,1)+\tilde{A}_{1}^{(\infty)}$, where
$\tilde{A}_{1}^{(\infty)}\sim\mathcal{U}_{[-\theta_{1},\theta_{1}]}$ in
$[-\theta,\theta]$ as a result of iteratively applying convolution and
truncation. For fixed $t$, this random variable has lower variance than
Brownian motion, but higher than $\tilde{A}_{1}$. In particular, Eq. A yields
for the sensory unit $P_{\mathcal{N}(0,1)}(\mathcal{E}_{1}\neq
0)=.3973<P_{\tilde{A}_{1}^{(t)}}(\mathcal{E}_{1}\neq
0)=.4251<P_{\tilde{A}_{1}^{(\infty)}+\mathcal{N}(0,1)}(\mathcal{E}_{1}\neq
0)=.4608<P_{B^{(t)}}(\mathcal{E}_{1}\neq 0)=.5157$.

###  B.2 Hitting time density

The hitting time is the random time when a stochastic process first reaches
some target subset of the codomain. For Brownian motion on the real line
starting at zero, the hitting time for the interval $[\theta,\infty)$ follows
[210] a Lévy distribution with scale parameter $\theta^{2}$, as
$\tau_{\theta_{1}}=\frac{\theta}{\sqrt{2\pi}}e^{-\theta^{2}/(2t)}t^{-3/2}$.
Its heavy tail falls off as $\sim t^{-3/2}$, and its mean is undefined. For
fractional Brownian motion with Hurst parameter $H$, it can be shown that it
follows $\sim t^{2-H}$ [74]. However, the problem is essentially changed if
the target subset is defined by two flanking thresholds as
$\mathbb{R}\setminus[-\theta,\theta]$. Then $\tau_{\theta}=\min{\\{t\geq
0\,:\,W^{(t)}=\theta\\}}$, which is the same as the first exit time from
$[-\theta,\theta]$. If we assume no reset, the first exit time of $A_{1}$ from
$[-\theta_{1},\theta_{1}]$ (or first hit time of its complement interval)
would be the same as the hitting time of $\theta_{1}$ for reflected Brownian
motion $|W|$ (RBM). The expected value of RBM is
$\langle|W^{(t)}|\rangle=\sqrt{2t/\pi}$ (from the Wiener process transition
density function with origin at zero $f_{W^{(t)}}(x)=\frac{1}{\sqrt{2\pi
t}}e^{-\frac{x^{2}}{2t}}$ [210]), so the threshold
$\langle|W^{(t_{1})}|\rangle=\theta_{1}$ should be attained at time
$t_{1}=\frac{\pi}{2}\theta_{1}^{2}$. Thus, if $\tilde{A}_{1}\sim|W|$, the
threshold $\theta_{1}=\sqrt{2/\pi}$ would be on average reached at $t_{1}=1$.
However, the actual values are $t_{1}=1.125$ and $\theta_{1}=.8464$ because
leftover activity from quiescent trials and the threshold updating rule (Eq.
16) push the threshold outwards.

The first exit time of $\tilde{A}_{1}$ can be calculated as the weighted sum
of each time step by its firing probability. Simulations yield
$\langle\tau_{\pm\theta_{1}}\rangle=1.655$, but the expected hitting time or
residence time for RBM is

$\langle\tau_{\pm\theta}\rangle=\theta^{2},$ (62)

which can be derived from the knowledge that $(W^{(t)})^{2}-t$ is a martingale
[70]. $\langle\tau_{\pm\theta_{1}}\rangle=.7164$ if $\theta_{1}=.8464$. This
discrepancy (a factor of $\approx 2.31$) is due to the discretization of
allowed firing times imposed by the existence of a refractory period
$T_{0}=1>\langle\tau_{\pm\theta_{1}}\rangle$. The discrepancy is roughly the
fraction of energy “skipping the threshold during the refractory period”,
which coincides with the inverse of the gain $g_{1}^{-1}\approx 2.35$ (see
Table 1 to look up $g_{1}$).

We could stop here, but in pursuance of getting a comprehensive picture of the
dynamics of integrator units in the limit of non-refractoriness, let us derive
the RBM hitting time density $f_{\tau_{\pm\theta}}$. The reflection principle
[70], which stands on the observation that the distribution of future paths is
symmetric at any given time, relates the running maximum of a Wiener process
up to time $t$ to its tail distribution at time $t$ 777Notice this implies
that for fixed $t$, the distribution of RBM coincides with the running maximum
distribution for Brownian motion $f_{|B^{(t)}|}(x)=f_{\max_{0\leq u\leq
t}W^{(u)}}(x)$. In turn, the running maximum is equal to the cumulative
distribution of the hitting time
$F_{\tau_{\pm\theta}}(t)=P(\tau_{\pm\theta}\leq t)$ [210], so

$P(\tau_{\pm\theta}\leq t)=P(\max_{0\leq u\leq
t}W^{(u)}\geq\theta)=2P(W^{(t)}>\theta).$

This can also be applied to calculate the hitting time of RBM, but some
caution is called for. The idea is applying the reflection principle
iteratively to the thresholds $-\theta$ and $\theta$. The RBM hitting time is
the first exit time from the $[-\theta,\theta]$ stripe:
$\tau_{\pm\theta}=\min{\\{t\,:\,W^{(t)}\notin[-\theta,\theta]\\}}$. As in
Brownian motion, when a particle first reaches the threshold $\theta$, it
carries on outwards or reverts with the same probability (by symmetry, the
same reasoning applies to particles attaining first $-\theta$). Hence,
$P(\max_{0\leq u\leq t}|B^{(u)}|\geq\theta)=2P(|W^{(t)}|>\theta)$ as long as
paths do not travel further than $2\theta$ from the first hit, i.e., if
$|B^{(t)}|<3\theta$. Beyond this point, some of the reverting paths will
attain the opposite threshold $-\theta$. To exclude the contribution from
virtual paths that hit thresholds more than once, we have to remove the path
densities located in $[(4k-1)\theta,(4k+1)\theta]$ for $k\in\mathbb{Z}$ —these
are the odd-numbered interleaved stripes, resulting from iteratively
reflecting the $[-\theta,\theta]$ stripe over the thresholds 888This is a
common procedure in finance for pricing double barrier options. This ensures
we count only paths that hit at least one of the thresholds, and not paths
hitting multiple times. To sum up,

$\displaystyle F_{\tau_{\pm\theta}}(t)$ $\displaystyle=$ $\displaystyle
P(\tau_{\pm\theta}\leq t)=P(\max_{0\leq u\leq t}|W^{(u)}|\geq\theta)$
$\displaystyle=$ $\displaystyle\int_{-\theta}^{\theta}P(|W^{(t)}|>\theta)dy$
$\displaystyle=$
$\displaystyle\int_{-\theta}^{\theta}\sum_{k\in\mathbb{Z}}2f_{W^{(t)}}(y+(4k+2)\theta)dy$
$\displaystyle=$
$\displaystyle\sum_{k\in\mathbb{Z}}2\int_{-\theta}^{\theta}\frac{1}{\sqrt{2\pi
t}}e^{-\frac{(y+(4k+2)\theta)^{2}}{2t}}dy$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2}{\pi
t}}\sum_{k\in\mathbb{Z}}\int_{(4k+1)\theta}^{(4k+3)\theta}e^{-\frac{x^{2}}{2t}}dx$
$\displaystyle=$
$\displaystyle\sum_{k\in\mathbb{Z}}\operatorname{erf}(\frac{(4k+3)\theta}{\sqrt{2t}})-\operatorname{erf}(\frac{(4k+1)\theta}{\sqrt{2t}}),$

where we performed the change of variables $x=y+(4k+2)\theta$, and the error
function $\operatorname{erf}$ takes $\frac{x}{\sqrt{2t}}$ as variable. This is
equivalent to the generalized reflection principle formula [213]:
$P(\max_{0\leq u\leq
t}|W^{(u)}|\leq\theta)=\sum_{k\in\mathbb{Z}}f_{W^{(t)}}(y+4k\theta)-f_{W^{(t)}}(y+(4k+2)\theta)=1-P(\max_{0\leq
u\leq t}|W^{(u)}|>\theta)$. Finally, differentiating $P(\tau_{\pm\theta}\leq
t)$ returns the hitting time density:

$\displaystyle f_{\tau_{\pm\theta}}(t)$ $\displaystyle=$
$\displaystyle\frac{d}{dt}P(\tau_{\pm\theta}\leq t)$ (64) $\displaystyle=$
$\displaystyle\sum_{k\in\mathbb{Z}}\frac{d}{dt}(\operatorname{erf}(\frac{(4k+3)\theta}{\sqrt{2t}})-\operatorname{erf}(\frac{(4k+1)\theta}{\sqrt{2t}}))$
$\displaystyle=$ $\displaystyle\frac{\theta}{\sqrt{2\pi
t^{3}}}\sum_{k\in\mathbb{Z}}(4k+1)e^{-\frac{\theta^{2}(4k+1)^{2}}{2t}}-(4k-1)e^{-\frac{\theta^{2}(4k-1)^{2}}{2t}}$
$\displaystyle=$ $\displaystyle\frac{\theta}{\sqrt{2\pi
t^{3}}}\sum_{k\in\mathbb{Z}}(-1)^{k}(2k+1)e^{-\frac{\theta^{2}(2k+1)^{2}}{2t}}$
$\displaystyle=$ $\displaystyle\frac{\theta\sqrt{2}}{\sqrt{\pi
t^{3}}}\sum_{k=0}^{\infty}(-1)^{k}(2k+1)e^{-\frac{\theta^{2}(2k+1)^{2}}{2t}}$

Alternatively to Eq. 64, $f_{\tau_{\pm\theta}}$ can be calculated from the
moment-generating function of $f_{W^{(t)}}$, which can be derived from the
properties of martingales. $e^{\sigma W^{(t)}-\frac{\sigma^{2}}{2}t}$ is a
martingale [70], so for a Wiener process starting at the origin, it holds that
$\langle e^{0}\rangle=\langle\exp{(\sigma
W^{(\min(\tau_{+\theta},\tau_{-\theta}))}-\frac{\sigma^{2}}{2}\min(\tau_{+\theta},\tau_{-\theta}}))\rangle$,
so
$1=\frac{1}{2}\langle\exp{(\sigma\theta-\frac{\sigma^{2}}{2}\tau_{\pm\theta})}+\exp{(-\sigma\theta-\frac{\sigma^{2}}{2}\tau_{\pm\theta})}\rangle$,
and rearranging $\langle
e^{-\frac{\sigma^{2}}{2}\tau_{\pm\theta}}\rangle\frac{e^{\sigma\theta}+e^{-\sigma\theta}}{2}=1$.
After the change of variables $\sigma=\sqrt{2\theta}$, we get an expression
for the Laplace transform of the hitting time density [70]:

$\langle
e^{-s\,\tau_{\pm\theta}}\rangle=\operatorname{sech}{\theta\sqrt{2s}}.$

Its inverse Laplace transform
$\mathcal{L}^{-1}\\{\operatorname{sech}{\theta\sqrt{2s}}\\}(t)=\frac{1}{2\pi
i}\lim_{R\to\infty}\int_{c-iR}^{c+iR}\frac{e^{st}ds}{\cosh{\theta\sqrt{2s}}}$
is the hitting time density $f_{\tau_{\pm\theta}}(t)$, with $c\in\mathbb{R}$
such that the real part of any of the singularities of the integrand is less
than $c$. Cauchy’s residue theorem allows to compute the line integral of an
analytic function $g$ over a closed rectifiable curve $\Gamma$ in terms of the
function’s residues [214]. Let $g(s)=\frac{e^{st}}{\cosh{\theta\sqrt{2s}}}$ be
the integrand, and $\Gamma=C\cup\gamma$, where $\gamma=[c-iR,c+iR]$ is the
line integral we wish to calculate, and
$C=c+Re^{i\varphi}:\varphi\in[\frac{\pi}{2},\frac{3\pi}{2}]$ is the arc of the
semicircle with center $c$ and radius $R$. Then, the residue theorem states
that $\oint_{\Gamma}g=\int_{\gamma}g+\int_{C}g=2\pi
i\sum_{k}\operatorname{Res}{\left(g,s_{k}\right)}$, where $s_{k}$ are the
poles of $g$ located inside $\Gamma$. Hence, letting the radius $R$ of the
semicircle $\Gamma$ grow to infinity,

$\frac{1}{2\pi
i}\lim_{R\to\infty}\int_{c-iR}^{c+iR}g=\sum_{k}\operatorname{Res}{\left(g,s_{k}\right)}-\frac{1}{2\pi
i}\lim_{R\to\infty}\int_{C}g.$ (65)

Since $e^{st}$ is an entire function, the poles $s_{k}$ of $g$ are the zeroes
of $\cosh{\theta\sqrt{2s}}$, which are defined by
$\cosh{\theta\sqrt{2s_{k}}}=0$; this yields the distinct poles
$s_{k}=-\frac{\pi^{2}}{8\theta^{2}}(2k+1)^{2}$, with $k\in\mathbb{N}_{0}$;
seeing that all poles satisfy $s_{k}<0$, we can take $c=0$. We can evaluate
the second term of the r.h.s. of Eq. 65 through an approach akin to that of
Jordan’s lemma [214], i.e. by finding an upper bound to the contour integral
modulus

$\displaystyle\left|\int_{C}g\right|$ $\displaystyle\leq$
$\displaystyle\int_{C}|g|=\int_{\frac{\pi}{2}}^{\frac{3\pi}{2}}\left|\frac{e^{tR(\cos{\varphi}+i\sin{\varphi})}Rie^{i\varphi}}{\cosh{(\theta\sqrt{2R}e^{i\varphi/2}})}\right|d\varphi$
$\displaystyle=$ $\displaystyle
R\int_{\frac{\pi}{2}}^{\frac{3\pi}{2}}\frac{e^{tR\cos{\varphi}}}{|\cosh{(\theta\sqrt{2R}e^{i\varphi/2}})|}d\varphi$
$\displaystyle\leq$ $\displaystyle
2R\int_{\frac{\pi}{2}}^{\frac{3\pi}{2}}e^{tR\cos{\varphi}}d\varphi,$

which has limit zero as $R\to\infty$ because $\cos{\varphi}\leq 0$ for
$\varphi\in[\frac{\pi}{2},\frac{3\pi}{2}]$. The last step requires calculating
$\max_{\varphi\in[\frac{\pi}{2},\frac{3\pi}{2}]}|\cosh{(\theta\sqrt{2R}e^{i\varphi/2}})|^{-1}=\max_{\varphi\in[\frac{\pi}{2},\frac{3\pi}{2}]}2|\cosh{(2\theta\sqrt{2R}\cos{\frac{\varphi}{2}})}+\cos{(2\theta\sqrt{2R}\cos{\frac{\varphi}{2}})}|^{-2}=2$,
which is achieved for $\varphi=\pi$. Consequently,

$f_{\tau_{\pm\theta}}=\frac{1}{2\pi
i}\lim_{R\to\infty}\int_{-iR}^{iR}g=\sum_{k}\operatorname{Res}{\left(g,s_{k}\right)}.$

Since
$\frac{d}{ds}\cosh{\theta\sqrt{2s}}\rvert_{s=s_{k}}=\frac{\theta}{\sqrt{2s_{k}}}\sinh{\theta\sqrt{2s_{k}}}=\frac{\theta^{2}}{\frac{\pi}{2}(2k+1)}\sin{\frac{\pi}{2}(2k+1)}\neq
0$, all poles ($\forall k\in\mathbb{N}^{0}$) are simple and isolated. Hence,
$\frac{e^{st}}{\cosh{\theta\sqrt{2s}}}$ is a meromorphic function (except at
infinity), and it can be expanded (by Mittag-Leffler’s theorem) around each of
its poles $s_{k}$ as a sum of partial fractions, where each partial fraction
expansion is the principal part of $\operatorname{sech}{\theta\sqrt{2s}}$ at
the pole $s_{k}$ [214]. The residue at pole $s_{k}$ is thus

$\displaystyle\operatorname{Res}{\left(g,s_{k}\right)}$ $\displaystyle=$
$\displaystyle\lim_{s\to s_{k}}(s-s_{k})g(s)$ $\displaystyle=$
$\displaystyle\frac{e^{s_{k}t}}{\frac{\theta^{2}}{\frac{\pi}{2}(2k+1)}\sinh{\theta\sqrt{2s_{k}}}}$
$\displaystyle=$
$\displaystyle\frac{\frac{\pi}{2}(2k+1)}{\theta^{2}}\frac{e^{-\frac{\pi^{2}}{8\theta^{2}}(2k+1)^{2}t}}{\sin{\frac{\pi}{2}(2k+1)}}$
$\displaystyle=$
$\displaystyle\frac{\pi}{2\theta^{2}}(2k+1)(-1)^{k}e^{-\frac{\pi^{2}}{8\theta^{2}}(2k+1)^{2}t}.$

Finally, by the residue theorem

$\displaystyle f_{\tau_{\pm\theta}}(t)$ $\displaystyle=$
$\displaystyle\frac{\pi}{2\theta^{2}}\sum_{k=0}^{\infty}(-1)^{k}(2k+1)e^{-\frac{\pi^{2}}{8\theta^{2}}(2k+1)^{2}t}.$
(66)

Numerical simulations show that Eq. 64 converges much faster than Eq. 66 for
$t\ll 1$, and vice versa for $t\gg 1$. We can use Eq. 66 to calculate the
expected hitting time $\tau_{W}$ of a Brownian neuron, with threshold
$\theta_{1}$, where firing of sensory units is possible only at specific times
separated by $T_{0}=1$, as
$\langle\tau_{W}\rangle=\sum_{k\in\mathbb{Z}^{+}}k\frac{f_{\tau_{\pm\theta}}(k)}{\sum_{i\in\mathbb{Z}^{+}}f_{\tau_{\pm\theta}}(i)}=1.21756$,
which lies between the values of $\langle\tau_{\theta_{1}}\rangle$ for RBM and
the sensory unit. Importantly, this shows that the hitting time density of
Brownian neurons has an exponentially decaying tail.

###  B.3 Brownian neurons mimicking refractory units: reduced variance and
noisy relaxation processes

How to tinker with Brownian neurons so they reproduce the firing rates and
mean hitting times of discrete (sensory level) neurons? For mean hitting
times, perhaps the simplest way is reducing the variance of the Brownian
neuron. Reusing Eq. B.2 gives for the standard deviation of the Brownian
neuron that equates its mean hitting time to that of the sensory unit
$\sigma_{\text{red}}=\theta_{1}/\sqrt{\langle\tau_{1}\rangle}=.6579$. For
firing rates, we can counterpoise the discrepancy between Brownian and sensory
neurons by adding a zero-reverting term to the Brownian neuron transition
density, thus rendering it a noisy relaxation or (driftless) “Ornstein-
Uhlenbeck neuron”. This mimics the activity resetting effect of discharges
(Eq. 2). The problem now is to find the mean-reverting coefficient $\rho$ of
the stochastic differential equation $dx=-\rho xdt+dW$ obeyed by the Ornstein-
Uhlenbeck neuron that matches the firing rate of the sensory unit
$P_{\tilde{A}_{1}^{(t)}}(\mathcal{E}_{1}\neq 0)$. Numerical optimization of
$\rho$ with respect to its transition density [210] yields $\rho=.2819$, which
compared to $0$ shows the influence of discharge-induced suppression on the
dynamics of $\tilde{A}_{1}$.

##  Appendix C  Dynamics of thresholds, pulses, and activities

Threshold dynamics can be described with an AR(1) model. Recalling Eq. 4, and
assuming $w=w^{+}=w^{-}$:

$\displaystyle\theta^{(t+1)}$ $\displaystyle=$
$\displaystyle(1-w)\theta^{(t)}+w|\tilde{a}^{(t)}|$ (67) $\displaystyle=$
$\displaystyle w\langle\Theta\rangle+(1-w)\theta^{(t)}+\tilde{\epsilon}^{(t)}$

where $w\langle\Theta\rangle$ is the average, $1-w$ is the regressor
parameter, and
$\tilde{\epsilon}^{(t)}=w(|\tilde{a}^{(t)}|-\langle\Theta\rangle)$ is a white
noise process with zero mean and variance
$\sigma_{\tilde{\epsilon}}^{2}=w^{2}\operatorname{var}{|\tilde{A}|}$. The $l$
subscripts have been omitted for clarity. For $1-w\approx 1$, this constitutes
a wide-sense stationary process and the central limit theorem ensures that (in
the steady state) $\Theta$ is approximately Gaussian-distributed, with mean
$\langle\Theta\rangle$ and covariance
$C_{\theta,\theta}(n)\equiv\operatorname{cov}{[\theta^{(t+n)]},\theta^{(t)}]}=\frac{\sigma_{\tilde{\epsilon}}^{2}}{1-(1-w)^{2}}(1-w)^{|n|}=\frac{w^{2}\operatorname{var}{|\tilde{A}|}}{1-(1-w)^{2}}(1-w)^{|n|}\approx\frac{\operatorname{var}{|\tilde{A}|}}{2}w(1-w)^{|n|}=\frac{\sigma_{\tilde{\epsilon}}^{2}}{2w}(1-w)^{|n|}$,
whence

$C_{\theta,\theta}(n)\approx\frac{\sigma_{\tilde{\epsilon}}^{2}}{2w}(1-w)^{|n|},$
(68)

where the approximation is valid for small $w$ [215]. $C_{\theta,\theta}(n)$
decays with a time constant $-T_{s}/\log{(1-w)}\approx T_{s}/w$ ($T_{s}=T_{l}$
sampling period) 999The expected mean
$\langle\theta^{(t+n)}|\theta^{(t)}\rangle$ and variance
$\operatorname{var}{(\theta^{(t+n)}|\theta^{(t)})}$ after $n$ iterations can
be derived explicitly by iterated application of Eq. 67
$\langle\theta^{(t+n)}|\theta^{(t)}\rangle=\langle\Theta\rangle(1-(1-w)^{n})+\theta^{(t)}(1-w)^{n}$
and
$\operatorname{var}{(\theta^{(t+n)}|\theta^{(t)})}\approx(1-(1-w)^{2n})\frac{w}{2}\operatorname{var}{|\tilde{a}^{(t)}|}$.
Notice that we set different sampling periods $T_{s}=T_{l}=\nu_{l-1}^{-1}$
corresponding to each level $l=1\ldots 7$ (with $\nu_{0}^{-1}=T_{0}$), because
thresholds can be updated only upon arrival of pulses, so the firing rate of
the preceding unit determines the threshold updating rate at each unit. Since
$\tau\gg T_{l}$ for small $w$, we can calculate the power spectrum through the
Fourier transform 101010Specifically, the non-unitary Fourier transform, in
ordinary frequency units. of the (continuum approximation) covariance:

$\displaystyle S_{\theta,\theta}(\nu|l)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{\tilde{\epsilon}}^{2}}{1-(1-w)^{2}}\frac{2\tau^{-1}}{\tau^{-2}+(2\pi\nu)^{2}}$
(69) $\displaystyle\approx$
$\displaystyle\frac{w^{2}\operatorname{var}{|\tilde{A}|}}{2w}\frac{2w/T_{l}}{w^{2}/T_{l}^{2}+(2\pi\nu)^{2}}$
$\displaystyle=$
$\displaystyle\operatorname{var}{|\tilde{A}|}\frac{T_{l}}{1+\left(\frac{T_{l}}{w}2\pi\nu\right)^{2}},$

which is a Lorentzian function, with maximum
$S_{\theta,\theta}(0)=T_{l}\operatorname{var}{|\tilde{A}|}$ at $\nu=0$ and
corner frequency $\gamma_{\theta}=\frac{w}{2\pi T_{l}}$. The corner
frequencies for each unit firing period are shown in Table 3. Note that the
time resolution a each level $T_{s}=T_{l}$ decreases with the level number,
reflecting the slower dynamics of later units.

The continuous-time analogue of the AR(1) process of Eq. 67 is a Ornstein-
Uhlenbeck process (OUP) with mean $\langle\Theta_{l}\rangle$, mean-reversal
term $-w$, and noise variance $\sigma_{\tilde{\epsilon}}$. This can be deduced
by rearranging Eq. 67 as
$\theta_{l}^{(t)}=\theta_{l}^{(t-1)}-w(\theta_{l}^{(t-1)}-\langle\Theta_{l}\rangle)+\tilde{\epsilon}^{(t)}$
and taking the limit as the step size approaches zero:

$d\theta=-w(\theta-\langle\Theta\rangle)dt+\sigma_{\tilde{\epsilon}}dW.$ (70)

Donsker’s invariance principle guarantees that the random walk series of
$\tilde{\epsilon}^{(t)}$ increments
$\lim_{n\to\infty}\frac{1}{\sqrt{n}}\sum_{i=1}^{nt}\tilde{\epsilon}^{(i)}$
converges in distribution to $\sigma_{\tilde{\epsilon}_{l}}W(t)$ for
$t\in[0,1]$, where $W(t)$ is the standard Wiener process. This corresponds to
the scaling limit $T_{s}\to 0$ of the discrete AR(1) process in Eq. 67. When
$T_{s}$ is not negligible, Eq. 70 is no longer a good approximation and it has
to be corrected with additional factors that depend on the step size. In the
steady state, the OUP covariance function [210] is
$\operatorname{cov}{[\theta^{(t+u)},\theta^{(t)}]}=\frac{\operatorname{var}{|\tilde{A}|}}{2}we^{-wu}=\frac{\sigma_{\tilde{\epsilon}}^{2}}{2w}e^{-wu}$,
which for $w\to 0$ approaches the covariance function of the discrete AR(1)
process (Eq. 68). Numerical simulations with $w=0.01$ yielded
$\sigma_{|\tilde{A}|}\equiv\sqrt{\operatorname{var}{|\tilde{A}|}}\approx
0.53\langle\Theta_{l}\rangle$ for $l\geq 3$ (Table 3) 111111cf. variance of
pulses $\sigma_{\mathcal{E}}=0.35\langle\Theta_{l}\rangle$.. This gives the
following reasonable continuum approximation to the threshold dynamics for
$l\geq 3$:

$\frac{1}{w}d\theta\approx-(\theta-\langle\Theta\rangle)dt+0.53\langle\Theta\rangle
dW.$ (71)

The solid curves in Fig. 3 are derived from Eq. 71, except the noise variances
from Table 3 were used instead of 0.53.

Table 3: Numerical estimates of threshold noise variance components $\sigma_{|\tilde{A}_{l}|}$ normalized by the average threshold, and corner frequencies $\gamma_{\theta,l}$ and $\gamma_{a,l}$ (time series of length $n_{t}=9\cdot 10^{6}$). Level $l$ | $\frac{\sigma_{|\tilde{A}_{l}|}}{\langle\Theta_{l}\rangle}$ | $\frac{\sigma_{|\mathcal{E}_{l}|}}{\langle|\mathcal{E}_{l}|\rangle}$ |  $\gamma_{\theta}(l)$ (Hz) |  $\gamma_{a}(l)$ (Hz)
---|---|---|---|---
1 | .7557 | .3493 | .318 | 59.2
2 | .5664 | .2576 | .135 | 19.5
3 | .5393 | .2491 | .061 | 8.41
4 | .5301 | .2421 | .028 | 3.85
5 | .5395 | .2419 | .013 | 1.80
6 | .5253 | .2425 | .006 | 0.82
7 | .5452 | .2411 | .003 | 0.38

Computing the power spectral density (PSD) using the discrete Fourier
transform on windowed finite time series entails unavoidably an estimation
bias. This is because using a finite-length signal is equivalent to
multiplying the original infinite-length signal by a window function, and the
frequency response of a finite window (Fejér’s kernel) is not just a Dirac
delta, but a distribution with a finite-sized mainlobe and many tailing off
sidelobes that spread out to cover the whole frequency range [219]. As a
result, the estimate of a windowed signal has power leaked into all
frequencies around any true signal frequency. This effect is called spectral
leakage, and it affects the periodogram and its derivatives such as Welch’s
method, typically biasing PSD estimates upwards [219] (but flat spectrum or
white noise signals are not affected). In particular, the high frequency power
of AR(1) models does not decay as a power law —as it does for the continuous
OUP— because these frequencies correspond to trial-by-trial changes, that
cannot vanish in discrete implementations. In general, on top of estimation
bias, there is always a trade-off between variance reduction and resolution
when estimating PSDs. Using Welch’s method with a (rectangular) window length
$n_{w}=10^{5}$, the squared Fejér’s kernel has a peak mainlobe of $10^{10}$
and its largest sidelobes peak at about 27dB below it (at a distance of
$f_{s}/n_{w}=\approx 10^{-3}$Hz, the lowest resolvable frequency). In Fig. 3,
observe that spectral leakage inflates the PSD at the lowest and highest
frequencies; the PSD cannot be estimated at the lowest frequencies ($\approx
10^{-3}$Hz) due to the loss of resolution attending the improvement in
variance afforded by Welch’s method. To verify that the bias is due to
spectral leakage, we also estimated PSDs using the computationally expensive
Thompson’s multitaper method (not shown), which uses a bank of optimal
bandpass filters called Slepian sequences (instead of a bank of rectangular
bandpass filters as in and Welch’s method) that reduces Fejér’s kernel’s
sidelobes. For a time-bandwidth parameter $3$, the bias is substantially
reduced at both lowest and highest frequencies, although at the expense of
higher variance.

The PSD of a sum of signals is the sum of the signals’ PSDs only if the
signals are uncorrelated. Otherwise,
$S_{\theta_{i}+\theta_{j}}(\nu)=S_{\theta_{i},\theta_{i}}(\nu)+S_{\theta_{j},\theta_{j}}(\nu)+2\operatorname{Re}({S_{\theta_{i},\theta_{j}}(\nu)})$,
where $S_{\theta_{i},\theta_{j}}(\nu)$ is the cross power spectral density
(CPSD). Hence, estimating the PSD of $\sum_{i=1}^{L}\theta_{l}^{(t)}$ in
general requires computing numerically $L(L-1)/2=21$ pairwise CPSD functions
for $L=7$, with no closed-form expression. However, a reasonable approximation
can be achieved computing only the CPSD function that contributes the highest
power:
$S_{\sum_{i=1}^{7}\theta_{l}}(\nu)\approx\sum_{i=1}^{7}S_{\theta_{i},\theta_{i}}(\nu)+2\operatorname{Re}({S_{\theta_{6},\theta_{7}}(\nu)})$.
Inspection of the PSD of the sum of variables in Figs. 3 and 4 reveals that
pulses $\varepsilon_{l}$ have pairwise CPSD components that contribute to the
CPSD of the sum only at low frequencies. This is confirmed by their close to 1
coherence estimates and zero-lag cross-phase (Fig. 18).

Figure 18: Cross-spectral phase (cells above the main diagonal) and magnitude-
squared coherence (below the main diagonal) pair-wise estimates between pulses
$\varepsilon_{l}^{(t)}$ for $l=1,3,5,7$. The main diagonal is populated by PSD
plots, as in Fig. 6. The cross-spectral phase was defined as the argument of
the cross-spectral density $S_{\varepsilon_{i},\varepsilon_{j}}$; and the
magnitude-square coherence was defined as
$\frac{|S_{\varepsilon_{i},\varepsilon_{j}}|^{2}}{S_{\varepsilon_{i}}S_{\varepsilon_{j}}}$,
with PSD estimation parameters set as in Fig. 3.

The activities $a_{l}^{(t)}$ follow a random walk with stationary and
independent increments distributed as $\varepsilon_{l}^{(t)}$ that resets to
zero upon exiting $[-\theta_{l},\theta_{l}]$. Rewriting Eq. 2:

$a_{l}^{(t)}=(a_{l}^{(t-1)}+\varepsilon_{l-1}^{(t)})\mathbb{1}_{a_{l}^{(t-1)}+\varepsilon_{l-1}^{(t)}\in[-\theta_{l},\theta_{l}]}\,,$
(72)

where $\mathbb{1}_{x}$ is the indicator function, which equals zero when its
argument $x$ is false and one when true, and $\theta_{l}=\theta_{l}^{(t)}$ is
also a function of $t$. In general, subthreshold activity can be modeled as a
Markov renewal process with discrete discharge times $t$, inter-event times
$\nu_{a,l-1}^{-1}$ that equal the inverse of the firing rate of the preceding
unit, and activity as a continuously valued state $A_{l}$, where the jumps are
the incoming pulses from the preceding unit $\mathcal{E}_{l-1}$.
Alternatively, it can be construed as a discrete version of Brownian motion
with stochastic resetting [220]: since $\varepsilon^{(t)}$ has a well defined
variance, in the scaling limit $a_{l}^{(t)}$ behaves as a Brownian particle
that resets to zero as soon as it exits $[-\theta_{l},\theta_{l}]$, so Eq. 72
becomes:

$da=\sigma_{\mathcal{E}}dW\quad\text{if}\;a\in[-\theta,\theta]$ (73)

and $a=0$ otherwise, with $\sigma_{\mathcal{E}}$ the standard deviation of
$\mathcal{E}$.

The covariance function $C_{a,a}(t_{1},t_{2})$ of the process $a^{(t)}$ with
stochastic reset can be expressed [221] in terms of the covariance function of
the same process without reset —which is a Wiener process scaled by
$\sigma_{\mathcal{E}_{0}}$—
$\sigma_{\mathcal{E}_{0}}^{2}C_{W,W}(t_{1},t_{2})$, with $t_{2}>t_{1}$, as
$C_{a,a}(t_{1},t_{2})=P(\tau>t_{2}-t_{1})\left[\int_{0}^{t_{1}}P(\tau<t^{\prime})P(t^{\prime}<\tau<t^{\prime}+dt^{\prime})\sigma_{\mathcal{E}_{0}}^{2}C_{W,W}(t^{\prime},t_{2}-t_{1}+t^{\prime})+P(\tau>t_{1})\sigma_{\mathcal{E}_{0}}^{2}C_{W,W}(t_{1},t_{2})\right]$,
where $\tau$ is used to denote $\tau_{\pm\theta}$ for clarity. Since the
activity process without reset becomes a Wiener process in the scaling limit,
$C_{W,W}(t_{1},t_{2})=\min(t_{1},t_{2})=t_{1}$, and in the steady state
$P(\tau>t_{1})C_{W,W}(t_{1},t_{2})=\bar{F}_{\tau}(t_{1})t_{1}\to 0$ as
$t_{1}\to\infty$, so we obtain

$C_{a,a}(t_{1},t_{2})=\bar{F}_{\tau}(t_{2}-t_{1})\int_{0}^{t_{1}}F_{\tau}(t^{\prime})f_{\tau}(t^{\prime})\sigma_{\mathcal{E}_{0}}^{2}t^{\prime}dt^{\prime},$
(74)

where $\bar{F}_{\tau}(t)=1-F_{\tau}(t)$ is the tail distribution of $\tau$,
and $\sigma_{\mathcal{E}_{0}}^{2}=2D$ is just twice the diffusivity parameter
$D$ used in diffusion problems. There is no closed-form expression for this
equation that we know of, so we shall approximate it.

Repeated firing or discharging can be modeled as a renewal process, which
involves ignoring the activity values and focusing only on the accumulated
number of discharges $N(t)$ at each iteration $t$, where the jumps are resets
to zero. The inter-event times (the duration between successive discharges)
are distributed as the first exit time $f_{\tau}$ (Eq. 64 or 66). Although the
probability of firing in a given time interval depends on the last discharge
time, this dependency vanishes when averaging over long time spans. By the
elementary renewal theorem [210], the number of discharges is asymptotically
linear, i.e. in the long run $t\to\infty$ the averaged firing frequency
approaches the inverse of the hitting time mean

$\lim_{t\to\infty}\langle\nu\rangle\equiv\lim_{t\to\infty}\frac{1}{t}\langle
N(t)\rangle=\langle\tau\rangle^{-1}=\theta^{-2},$ (75)

where the last equality follows from Eq. 62. Since the inverse is a convex
function, by Jensen’s inequality the expectation of the inverse is larger or
equal than the inverse of the expectation:
$\langle\nu\rangle=\langle\tau^{-1}\rangle\geq\langle\tau\rangle^{-1}$, which
indicates that $\theta^{-2}$ is also a lower bound of the average firing rate.
For large $t$ the asymptotic variance of $N(t)$ behaves [210] according to

$\lim_{t\to\
\infty}\frac{\operatorname{var}{[N(t)]}}{t}=\frac{\operatorname{var}{[\tau]}}{\langle\tau\rangle^{3}}=\frac{2}{3}\theta^{-2},$
(76)

where in the rightmost equality we used
$\operatorname{var}{[\tau_{\pm\theta}]}=\frac{2}{3}\theta^{4}$, which can be
derived from the properties of martingales [70] in a similar manner to Eq. 62.
Together with Eq. 75, this determines the magnitude of the fluctuations of
inter-firing times about their expectation $\theta^{-2}$; these fluctuations
are approximately normally distributed [210].

From the first exit time density, we can deduce the probability that a unit
fires with rate $r=t^{-1}$ from $f_{\tau}(t)dt=f_{\nu}(r)dr$. Thus, the
sequence of firing times can be interpreted as a inhomogeneous Poisson point
process over $t$ with rate intensity function $rf_{\tau}(t)$ (cf. Eq. 66).
Hence the probability that the activity fires or reaches the threshold $n$
times within a segment $[0,t]$ starting from the previous reset at zero time
is

$P(\text{firing }n\text{ times})=\frac{[\Lambda(t)]^{n}}{n!}e^{-\Lambda(t)},$

where $\Lambda(t)=\int_{0}^{t}\tau^{-1}f_{\tau}(\tau)d\tau=\langle
t^{-1}\rangle$ before the first discharge, but as $t\to\infty$, by the
elementary renewal theorem $\Lambda(t)\to\langle\tau\rangle^{-1}t=rt$.

Thus, by neglecting short time correlations, the elementary renewal theorem
furnishes a simple approximation to the firing rate. This effectively amounts
to construing repeated firing as a Poisson point process with rate
$r=\theta^{-2}$. By considering the inter-event times of $a^{(t)}$ as a
(homogeneous) Poisson process with constant resetting rate, we can derive
again (Eq. 74) the covariance function of $a^{(t)}$ in terms of the covariance
function of Brownian motion [221] by substituting $F_{\tau}$ for its
homogeneous counterpart as

$\displaystyle C_{a,a}(t_{1},t_{2})$ $\displaystyle\approx$ $\displaystyle
e^{-r(t_{2}-t_{1})}\int_{0}^{t_{1}}re^{-rt^{\prime}}\sigma_{\mathcal{E}}^{2}t^{\prime}dt^{\prime}$
(77) $\displaystyle=$ $\displaystyle
r\sigma_{\mathcal{E}_{0}}^{2}e^{-r(t_{2}-t_{1})}\int_{0}^{t_{1}}t^{\prime}e^{-rt^{\prime}}dt^{\prime}$
$\displaystyle=$
$\displaystyle\frac{\sigma_{\mathcal{E}_{0}}^{2}}{r}e^{-r(t_{2}-t_{1})}[1-e^{-rt_{1}}(rt_{1}+1)]$
$\displaystyle\underset{t_{1}\to\infty}{=}$
$\displaystyle\frac{\sigma_{\mathcal{E}_{0}}^{2}}{r}e^{-r(t_{2}-t_{1})},$

where, in the first line, the first factor is the probability that there is no
reset between $t_{1}$ and $t_{2}$ (Eq. C with $n=0$); and the integrand
denotes the probability that the last reset before $t_{1}$ occurs between
$t^{\prime}$ and $t^{\prime}+dt^{\prime}$ [221]. Note that this approximated
covariance function is stationary only if fluctuations around the expectation
are ignored (from the elementary renewal theorem). These fluctuations
introduce a bias, but we can alleviate it by making use of Eq. 76 and taking
the expectation of Eq. 77 with respect to the approximately normally
distributed $r$:

$C_{a,a}(t)\approx\sigma_{\mathcal{E}_{0}}^{2}\langle\frac{e^{-r|t|}}{r}\rangle\approx\theta^{2}\sigma_{\mathcal{E}_{0}}^{2}e^{-\frac{4}{3}\theta^{-2}|t|},$
(78)

where $t=t_{2}-t_{1}$ (when $t_{2}>t_{1}$); we approximated $\langle
r^{-1}\rangle\approx\langle r\rangle^{-1}$; we used
$r=\langle\nu\rangle\approx\mathcal{N}(\theta^{-2},\frac{2}{3}\theta^{-2})$
from Eqs. 75 and 76; and $\langle
e^{\mathcal{N}(\theta^{-2},\frac{2}{3}\theta^{-2})}\rangle=e^{\frac{4}{3}\theta^{-2}}$
by the properties of the lognormal distribution. $C_{a,a}$ decays with time
constant $\frac{3}{4}\theta^{2}$. Finally, the PSD of $a^{(t)}$ can be
computed by Fourier-transforming Eq. 78:

$S_{a,a}(\nu|l)=\theta_{l}^{2}\sigma_{\mathcal{E}_{0}}^{2}\frac{2\frac{4}{3}\theta_{l}^{-2}T_{0}^{-1}}{(2\pi\nu)^{2}+(\frac{4}{3}\theta_{l}^{-2}T_{0}^{-1})^{2}}=\frac{\sigma_{\mathcal{E}_{0}}^{2}\theta_{l}^{2}\gamma_{a}}{\pi(\nu^{2}+\gamma_{a}^{2})},$
(79)

which is again a Lorentzian function —like the PSD of $\theta^{(t)}$— with
corner frequency $\gamma_{a}=(2\pi\frac{3}{4}\theta_{l}^{2}T_{0})^{-1}$ and
maximum
$S_{a,a}(0)=2\sigma_{\mathcal{E}_{0}}^{2}T_{0}(\frac{3}{4})^{2}\theta_{l}^{4}$.
Notice that, unlike for $S_{\theta,\theta}(\nu|l)$, the “sampling rate” is
fixed to $T_{0}$ across all levels. This is because the Brownian motion
approximation to activities (Eq. 73) and its subsequent assumptions are not
valid at the large timescales of higher levels: they incorporate large and
sparse fluctuations induced by $\varepsilon^{(t)}_{l}$ that can only be
modeled as Brownian motion by breaking up the $\varepsilon^{(t)}_{l}$ into
small fluctuations across many small $T_{0}$ iterations.

This still does not yield a good estimate of the numerically computed PSD
because it overestimates power across the whole frequency range with respect
to the numerical simulation output. The reason lies in its discretization and
recording schemes: since activities are reset at every iteration when a
discharge occurs, $P(A_{l}=0)=g_{l}$, so the amplitude of the simulated time
series of activities is a fraction $1-g_{l}$ of the amplitude of the
continuous-time analog, upon which we derived the PSD estimate. With this last
modification, the PSD estimate of the discretized time series becomes

$S_{a,a}(\nu|l)\approx\frac{\sigma_{\mathcal{E}_{0}}^{2}\theta_{l}^{2}(1-g_{l})^{2}\gamma_{a}}{\pi(\nu^{2}+\gamma_{a}^{2})}.$
(80)

Unlike pulses, activities are pair-wise incoherent over the whole frequency
range (Fig. 19), so a specific frequency power in a given level is not
transmuted into power in the subsequent level. As a corollary, the PSD of the
sum of activities can be approximated by the sum of their PSDs (Fig. 5).

Figure 19: Cross-spectral phase and magnitude-squared coherence pair-wise
estimates of subthreshold activities $a_{l}^{(t)}$ for $l=1,3,5,7$. Layout and
estimation parameters as in Fig. 18.

Regardless of the specific parameters, the takeaway is that activity dynamics,
characterized by interleaved periods of diffusion and resetting forces, are
well described by the OUP, a similar stochastic process where mean-reverting
forces replace jump resets.

##  Appendix D  Some relationships among activity variables’ statistics

From the definition of variance:

$\langle|\tilde{A}_{l}|^{2}\rangle=\langle\tilde{A}_{l}^{2}\rangle=\operatorname{var}{|\tilde{A}_{l}|}+\langle|\tilde{A}_{l}|\rangle^{2}=\operatorname{var}{|\tilde{A}_{l}|}+\langle\Theta_{l}\rangle^{2}=\langle
A_{l}^{2}\rangle+\langle\mathcal{E}_{l-1}^{2}\rangle=\langle
A_{l}^{2}\rangle+\langle\mathcal{E}_{l}^{2}\rangle g_{l}$

$\langle|\tilde{A}_{l}|^{2}\rangle=\langle|A_{l}+\mathcal{E}_{l-1}|^{2}\rangle=\langle(A_{l}+\mathcal{E}_{l-1})^{2}\rangle=\langle
A_{l}^{2}+\mathcal{E}_{l-1}^{2}+2A_{l}\mathcal{E}_{l-1}\rangle=\langle
A_{l}^{2}\rangle+\langle\mathcal{E}_{l-1}^{2}\rangle+2\langle
A_{l}\mathcal{E}_{l-1}\rangle=\langle|A_{l}|^{2}\rangle+\langle|\mathcal{E}_{l-1}|^{2}\rangle$

From the law of total variance applied to $|A|$, $|A^{+}|$, and $|\tilde{A}|$:

$\operatorname{var}{|\tilde{A}_{l}|}=(1-g_{l})\operatorname{var}{|A_{l}^{+}|}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})(\langle|\mathcal{E}_{l}|\rangle^{2}+\langle|A_{l}^{+}|\rangle^{2})-2g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}^{+}|\rangle$

$\operatorname{var}{|A_{l}|}=\operatorname{var}{A_{l}^{+}}(1-g_{l})-\langle|A_{l}^{+}|\rangle^{2}(1-g_{l})^{2}$

$\operatorname{var}{|A_{l}^{+}|}=\operatorname{var}{A_{l}}(1-g_{l})^{-1}-\langle|A_{l}|\rangle^{2}(1-g_{l})^{-2}$

From the definition of variance:

$\operatorname{var}{|\tilde{A}_{l}|}=(1-g_{l})\operatorname{var}{|A_{l}^{+}|}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})(\langle|\mathcal{E}_{l}|\rangle-\langle|A_{l}^{+}|\rangle)^{2}=\langle
A_{l}^{2}\rangle-\langle|A_{l}|\rangle^{2}(1-g_{l})^{-1}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle^{2}+g_{l}(1-g_{l})^{-1}\langle|A_{l}|\rangle^{2}-2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle=\langle
A_{l}^{2}\rangle+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle^{2}-\langle|A_{l}|\rangle^{2}-2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle=\operatorname{var}{|A_{l}|}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle^{2}-2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle$

$\langle|\tilde{A}_{l}|^{2}\rangle=\operatorname{var}{|\tilde{A}_{l}|}+\langle\Theta_{l}\rangle^{2}=\operatorname{var}{|A_{l}|}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle^{2}-2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle+g_{l}^{2}\langle|\mathcal{E}_{l}|\rangle^{2}+\langle|A_{l}|\rangle^{2}+2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle=\operatorname{var}{|A_{l}|}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle^{2}+g_{l}^{2}\langle|\mathcal{E}_{l}|\rangle^{2}+\langle|A_{l}|\rangle^{2}=\operatorname{var}{|A_{l}|}+g_{l}\operatorname{var}{|\mathcal{E}_{l}|}+g_{l}\langle|\mathcal{E}_{l}|\rangle^{2}+\langle|A_{l}|\rangle^{2}=\langle
A_{l}^{2}\rangle+g_{l}\langle\mathcal{E}_{l}^{2}\rangle$

$\langle\Theta_{l}\rangle^{2}=\langle
A_{l}^{2}\rangle+\langle\mathcal{E}_{l}^{2}\rangle
g_{l}-\operatorname{var}{|A_{l}|}-g_{l}\operatorname{var}{|\mathcal{E}_{l}|}-g_{l}(1-g_{l})\langle|\mathcal{E}_{l}|\rangle^{2}+2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle=\langle|A_{l}|\rangle^{2}+g_{l}^{2}\langle|\mathcal{E}_{l}|\rangle^{2}+2g_{l}\langle|\mathcal{E}_{l}|\rangle\langle|A_{l}|\rangle$

The relation $\langle
A^{2}\rangle\approx\langle(\theta_{*}-|A^{+}|)^{2}\rangle$ was verified to
reasonably hold in numerical simulations. It equates the average subthreshold
activity energy to the average post-pulse subthreshold energy deficit below
the threshold energy in quiescent trials, which is a manifestation of energy
conservation.

Figure 20: Surface (top) and contour (bottom) plots of the joint densities
$P(a_{t+1}/k,a_{t})$ in left, $P(\varepsilon_{t+1},k\varepsilon_{t})$ in
middle left, $P(a_{t+1}/k,k\varepsilon_{t})$ in middle right, and
$P(\varepsilon_{t+1},a_{t})$ in right. Variables have been scaled by $k$ or
$1/k$ to facilitate comparison. The color bar at the right of the top right
plot is the color legend of the bottom row contour plots, it represents
probability.

##  Appendix E  Numerical analysis of the rescaled recurrent loop
configuration

The recurrence relation specified by Eqs. 34–39 in Section V implements a
hybrid of bootstrap resampling and discrete-time dynamical system that enables
estimating the densities of PC variables, in rescaled recurrent loop
configuration, a sort of “perfusive rescaled loop” that bites its own tail.
This can be programmed as an iterative loop.

One pool of $N=10^{5}$ data points was defined for each variable. Each pool
was initialized with values randomly picked out from the time series ($9\cdot
10^{6}$ points) of the homologous variable at level $l=1$ (see Fig. 2, second
row of left column) in the PC feedforward configuration described in Section
III.1. At each iteration $t$, pulses $\varepsilon_{i}$ were resampled randomly
from the pulse pool $\varepsilon_{1:N}$, whereas the gains $k_{t}$, thresholds
$\theta_{t}$ and subthreshold activities $a_{t}$ took the values they had in
the previous iteration. Then each variable was allocated a new value at
iteration $t+1$ according to Eqs. 34–39. This routine was looped $10^{6}$
times. Instead of using the last trial value $a_{t}$, subthreshold activities
could also be resampled randomly from the pool $a_{1:N}$ (just as pulses were
resampled from $\varepsilon_{1:N}$) to achieve the same results, because same-
trial subthreshold activities and incoming pulses are independent. The same
holds for $\theta_{t}$ and $k_{t}$ (in Eq. 38). This was verified in
simulations. This scheme generated the data underlying Figs. 9, 10, 20. Some
histogram functions (for plotting transition matrices) are susceptible to
binning artifacts; this can be often averted by adding a randomly signed
machine epsilon to “jitter” the samples.

Unlike in the feedforward configuration, here we need to estimate $k_{t}$ to
iterate the map. Since $k_{t}$ is a function of the pool of pulses
$\varepsilon_{1:N}$ (Eq. 38), it is susceptible to estimation errors (through
pulses), which in turn propagate to the other variables at each iteration.
This feature compounds the estimation uncertainty of all variables, and
precludes simulating the recurrent loop with fixed threshold and scaling
factor (i.e. $\theta_{*}$ and $k_{*}$ in FTA) —except in the limit of
infinitely large bootstrapping pools. Another caveat is that a recurrent loop
embodies a positive feedback that tends to magnify skewness: if left
unchecked, the densities will become increasingly lopsided, to the point that
the resulting densities will correspond to the case of perfectly correlated
driving input ($H=1$, see below). To avoid this effect, every few iterations
the pools must be re-symmetrized, by e.g. randomly flipping the sign of all
pulse samples. By definition, FTA is achieved in the limit of vanishing
$\theta_{t}$ fluctuations: by calculating $\theta$ with a pool of subthreshold
activity samples $\tilde{a}_{1:N}$, i.e. using Eq. 40 instead of Eq. 38, and
setting a large enough $N$ —by the law of large numbers, threshold
fluctuations wane along with $N$. However, the larger is the pool, the more
“inertia” it has, so in practice it is expedient to start with a small
(“light”) pool, and only when it reaches a steady state fill it up with new
samples; a large biased pool will remain so for many iterations, which is
computationally wasteful.

The above scheme can accommodate autocorrelated driving input with little
modification. First, the pulse $\varepsilon_{1:N}$ and subthreshold activity
$a_{1:N}$ pools should take in samples after applying the absolute value
function, so they contain only positive or zero values; since all densities
are symmetrical about the origin, no information is thereby lost. Second, the
incoming pulse magnitude is defined by the pool sample, whereas its sign is
pseudo-randomly generated as $\operatorname{sgn}{(\mathcal{U}_{[0,1]}<H)}$,
where $\mathcal{U}_{[0,1]}$ the uniform density with support $[0,1]$, and
$H\in[0,1]$ is a parameter describing the autocorrelation of driving input,
which for fractional Gaussian noise is the Hurst exponent. Thus, the larger is
$H$, the higher is the probability that two consecutive pulses have the same
sign, with $H=0.5$ characterizing white noise. The duplets used to estimate
transition matrices in Figs. 10, 20 were stored at every iteration in separate
variables.

A numerically simulated trajectory of the triplet
$[\sigma_{A}/\sigma_{A_{*}},\langle\Theta\rangle/\theta_{*},\langle|\mathcal{E}|\rangle/\langle|\mathcal{E}_{*}|\rangle]$
in phase space for a $10^{5}$ iterations fragment is shown in Fig. 21.

Figure 21: Simulated trajectory of states for $10^{5}$ iterations in rescaled
recurrent loop configuration (Eqs. 34, 37, 40, 39), embedded in the three-
dimensional manifold of the normalized variables $\sigma_{A}/\sigma_{A_{*}}$,
$\langle\Theta\rangle/\theta_{*}$ and
$\langle|\mathcal{E}|\rangle/\langle|\mathcal{E}_{*}|\rangle$. The dot
indicates the fixed threshold approximation (FTA) estimate location. This and
the trajectory are also shown projected (in gray) onto the three planes
defined pair-wise by the FTA estimates. Parameters: Pool size $N=10^{5}$,
Gaussian noise input Hurst exponent $H=0.5$, $10^{6}$-iteration burn-in.

The system was let settle in the steady state beforehand for $10^{6}$
iterations, so the trajectory wiggles in an attractor located nearby the FTA
point (although visually it is similar to a 3-dimensional Brownian particle)
while displaying a small bias respect to the FTA estimate (Fig. 11). In
general, reducing the pool size $N$ leads to larger fluctuations, but
increasing it too much hinders the trajectory from returning to the vicinity
of the expected value after a (random) fluctuation. Increasing the learning
rate while using the delta learning rule (Eq. 38 instead of Eq. 40) results in
larger oscillations in the $\langle\Theta\rangle/\theta_{*}$ dimension (not
shown).

Fig. 22 displays variable estimates using the delta learning rule for
thresholds in recurrent configuration as a function of input Hurst parameter.
Although it conveys similar information with respect to its homologous Fig. 13
with explicit threshold averaging, we can observe that with the smallest
learning rate of $w=10^{-4}$ (dashed line), the hill peak at $H\approx.74$ is
shifted to $H\approx.8$. The reason seems to lie in that under superdiffusion,
this configuration greatly overestimates the discharge rate $g$ and
underestimates the threshold with respect to the fixed threshold
approximation, as can be seen in Fig. 22. Finally, note that in the
feedforward configuration simulations in Section III, we used a learning rate
$w=0.01$, which corresponds to the dashed-dotted line in Fig. 22.

Figure 22: Perfusive cascade variables’ (see legend) estimates with delta
learning rule, as a function of driving input autocorrelation, in recurrent
rescaled loop configuration. $A$, $\mathcal{E}$ are scaled consistently with
$\langle\mathcal{E}^{2}\rangle=1$. Line style indicates the estimation
procedure: closed form expression (Eqs. 48–51) derived through the fixed
threshold approximation (solid thin); simulation with learning rate
$w=10^{-1}$ (dotted), $N=10^{-2}$ (dash-dotted), $N=10^{-3}$ (solid thick),
$N=2\cdot 10^{-4}$ (dashed). Dashed lines largely overlap with solid thick
lines.

##  Appendix F  Modeling of pulse size densities with the generalized Pareto
distribution

We suggested in Section III.1 that pulse density profiles follow approximately
a generalized Pareto distribution (GPD), which is a family of densities, often
used to model the shape of the tail of other distributions [222]. Here we use
the GPD as a parametric approximation to the tails of $f_{\tilde{A}}$, or
equivalently to the pulse sizes $f_{|\mathcal{E}|}$. The GPD density is
defined as

$f_{GPD}(x|\mu,\sigma,\xi)=\frac{1}{\sigma}\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-1-\frac{1}{\xi}},$
(81)

with parameters location $\mu$ (which defines the minimum of $\mathcal{E}$),
scale $\sigma$ (note this runs counter to its meaning of standard deviation in
other sections), and shape $\xi$ (which characterizes how steeply the tail
falls off). Its support is $x\geq\mu$ when $\xi\geq 0$ and $\mu\leq
x\leq\mu-\sigma/\xi\equiv M$ when $\xi<0$. Its mean is
$\mu+\frac{\sigma}{1-\xi}$ when $\xi<1$, its median is
$\mu+\sigma\frac{2^{\xi}-1}{\xi}$, and its mode is $\mu$.

By definition, $f_{|\mathcal{E}_{l}|}$ is the tail distribution of post-pulse
activity absolute value $|\tilde{A}_{l}|$. Its cumulative distribution
function $F_{|\mathcal{E}_{l}|}(x)\equiv P(|\mathcal{E}_{l}|\leq x)$ is also
the excess distribution of post-pulse activities above the firing threshold
$P(|\tilde{A}_{l}|-\theta_{l}\leq
x\;|\;|\tilde{A}_{l}|>\theta_{l})=\frac{F_{|\mathcal{E}_{l}|}(\theta_{l}+x)-F_{|\mathcal{E}_{l}|}(\theta_{l})}{1-F_{|\mathcal{E}_{l}|}(\theta_{l})}$.
For larger enough $l$, this can be interpreted as the tail of the cumulative
distribution of a sequence of identically distributed and independent random
variables conditioned to be larger than some threshold: each $\varepsilon_{l}$
can be decomposed as the sum of inputs $\varepsilon_{0}$’s by recursively
substituting each $\varepsilon_{l}$ by its summands
$\varepsilon_{l-1}+a_{l-1}$, which by construction are already conditioned on
their absolute value exceeding $\theta_{l-1}$. This supplies the conditions of
applicability of the Gnedenko-Pickands-Balkema-de Haan theorem, which states
that the tail of the excess distribution is asymptotically given by a GPD as
the variable approaches the end of the tail of the distribution [65]. This
implies that in the steady state, and for large enough $l$, the tail
distribution $F_{\mathcal{E}_{l}}$ can be asymptotically characterized by a
GPD.

This justifies using the GPD to model $f_{|\mathcal{E}_{l}|}$ for large enough
$l$. Maximum likelihood fits (through MATLAB’s fitdist function) of the GDP
parameters to the numerically obtained histograms of $f_{|\mathcal{E}_{l}|}$
yielded the estimates displayed in Table 4.

Table 4: Maximum likelihood estimates of GPD parameters fit to the numerical pulse size distributions of the perfusive cascade levels $l=1\ldots 7$. Location parameters $\mu_{l}$ were fixed to their corresponding level threshold means (cf. Table 2); scale parameters are scaled by their corresponding location parameters. Only significant digits are shown. Numerical simulation parameters are as in Fig. 2. Level $l$ | $\mu_{l}$ | $\sigma_{l}/\mu_{l}$ | $\xi_{l}$
---|---|---|---
1 | $\langle\Theta_{1}\rangle$ | .832 | -.151
2 | $\langle\Theta_{2}\rangle$ | .596 | -.169
3 | $\langle\Theta_{3}\rangle$ | .571 | -.214
4 | $\langle\Theta_{4}\rangle$ | .564 | -.239
5 | $\langle\Theta_{5}\rangle$ | .606 | -.304
6 | $\langle\Theta_{6}\rangle$ | .606 | -.319
7 | $\langle\Theta_{7}\rangle$ | .609 | -.315

The GPD parameters corresponding to the fixed threshold approximation can be
obtained through the estimates for $\theta$ and $\sigma_{|\mathcal{E}|}$ given
at the end of Section IV.4; thereby we can establish the equalities
$\mu_{*}=\theta_{*}$ and
$\mu_{*}+\frac{\sigma_{*}}{1-\xi_{*}}=\sigma_{|\mathcal{E}_{*}|}$. This yields

$\displaystyle\frac{\sigma_{*}}{\mu_{*}}$ $\displaystyle\approx$
$\displaystyle 0.6579$ $\displaystyle\xi_{*}$ $\displaystyle\approx$
$\displaystyle-0.4274.$

It follows that $-1-\frac{1}{\xi_{*}}\approx 1.3396$ is the exponent of
$f_{GPD}$, the mean is $\frac{\sigma_{|\mathcal{E}_{*}|}}{\mu_{*}}\approx
1.4609$ (which entails a gain of $.4349$), the median
$1+\frac{\sigma_{*}}{\mu_{*}}\frac{2^{\xi}_{*}-1}{\xi_{*}}\approx 1.3947$
(gain $.5$), and the maximum of the support $\frac{M_{*}}{\mu_{*}}\approx
2.5393$ (notice all these values have been scaled by the threshold
$\mu_{*}=\theta_{*}$). Finally, the derivative of the density at the threshold
is $f_{GPD}^{\prime}(\mu_{*})=-\frac{1+\xi_{*}}{\sigma_{*}}\approx-0.8703$.

## References

  * Prigogine and Nicolis [1985] I. Prigogine and G. Nicolis, Self-Organisation in Nonequilibrium Systems: Towards A Dynamics of Complexity, in  _Bifurcation Analysis_ (Springer Netherlands, 1985) pp. 3–12.
  * Vespignani and Zapperi [1998] A. Vespignani and S. Zapperi, How self-organized criticality works: A unified mean-field picture, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 57, 6345 (1998), arXiv:9709192 [cond-mat] .
  * Bonachela and Muñoz [2009] J. A. Bonachela and M. A. Muñoz, Self-organization without conservation: True or just apparent scale-invariance?, Journal of Statistical Mechanics: Theory and Experiment , P09009 (2009), arXiv:0905.1799 .
  * Bak  _et al._ [1987] P. Bak, C. Tang, and K. Wiesenfeld, Self-organized criticality: An explanation of the 1/f noise, Physical Review Letters 59, 381 (1987).
  * Gutenberg and Richter [1956] B. Gutenberg and C. F. Richter, Earthquake magnitude, intensity, energy, and acceleration: (Second paper), Bulletin of the Seismological Society of America 46, 105 (1956).
  * Feder and Feder [1991] H. J. S. Feder and J. Feder, Self-organized criticality in a stick-slip process, Physical Review Letters 66, 2669 (1991).
  * Olami  _et al._ [1992] Z. Olami, H. J. S. Feder, and K. Christensen, Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes, Physical Review Letters 68, 1244 (1992).
  * Nagel and Paczuski [1995] K. Nagel and M. Paczuski, Emergent traffic jams, Physical Review E 51, 2909 (1995).
  * Drossel and Schwabl [1992] B. Drossel and F. Schwabl, Self-organized critical forest-fire model, Physical Review Letters 69, 1629 (1992).
  * Rhodes  _et al._ [1997] C. J. Rhodes, H. J. Jensen, and R. M. Anderson, On the critical behaviour of simple epidemics, Proceedings of the Royal Society B: Biological Sciences 264, 1639 (1997).
  * Corral  _et al._ [1995] A. Corral, C. J. Pérez, A. Díaz-Guilera, and A. Arenas, Self-organized criticality and synchronization in a lattice model of integrate-and-fire oscillators, Physical Review Letters 74, 118 (1995).
  * Herz and Hopfield [1995] A. V. Herz and J. J. Hopfield, Earthquake cycles and neural reverberations: Collective oscillations in systems with pulse-coupled threshold elements., Physical review letters 75, 1222 (1995).
  * Chialvo and Bak [1999] D. Chialvo and P. Bak, Learning from mistakes, Neuroscience 90, 1137 (1999).
  * Lu  _et al._ [1993] E. T. Lu, R. J. Hamilton, J. M. McTiernan, and K. R. Bromund, Solar flares and avalanches in driven dissipative systems, The Astrophysical Journal 412, 841 (1993).
  * Beggs and Plenz [2003] J. M. Beggs and D. Plenz, Neuronal Avalanches in Neocortical Circuits, Journal of Neuroscience 23, 11167 (2003).
  * Langton [1990] C. G. Langton, Computation at the edge of chaos: Phase transitions and emergent computation, Physica D: Nonlinear Phenomena 42, 12 (1990).
  * Shew and Plenz [2013] W. L. Shew and D. Plenz, The functional benefits of criticality in the cortex., The Neuroscientist : a review journal bringing neurobiology, neurology and psychiatry 19, 88 (2013).
  * Conant and Ashby [1970] R. C. Conant and W. R. Ashby, Every good regulator of a system must be a model of that system, International Journal of Systems Science 1, 89 (1970).
  * Friston  _et al._ [2006] K. Friston, J. Kilner, and L. Harrison, A free energy principle for the brain, Journal of Physiology Paris 100, 70 (2006).
  * Chialvo [2006] D. R. Chialvo, Psychophysics: Are our senses critical?, Nature Physics 2, 301 (2006).
  * Hasenstaub  _et al._ [2010] A. Hasenstaub, S. Otte, E. Callaway, and T. J. Sejnowski, Metabolic cost as a unifying principle governing neuronal biophysics, Proceedings of the National Academy of Sciences of the United States of America 107, 12329 (2010).
  * Roberts  _et al._ [2014] J. A. Roberts, K. K. Iyer, S. Vanhatalo, and M. Breakspear, Critical role for resource constraints in neural models, Frontiers in Systems Neuroscience 8, 154 (2014).
  * Novikov  _et al._ [1997] E. Novikov, A. Novikov, D. Shannahoff-Khalsa, B. Schwartz, and J. Wright, Scale-similar activity in the brain, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 56, R2387 (1997).
  * Linkenkaer-Hansen  _et al._ [2001] K. Linkenkaer-Hansen, V. V. Nikouline, J. M. Palva, and R. J. Ilmoniemi, Long-Range Temporal Correlations and Scaling Behavior in Human Brain Oscillations, The Journal of Neuroscience 21, 1370 (2001).
  * Kello  _et al._ [2010] C. T. Kello, G. D. A. Brown, R. Ferrer-i Cancho, J. G. Holden, K. Linkenkaer-Hansen, T. Rhodes, and G. C. Van Orden, Scaling laws in cognitive sciences, Trends in Cognitive Sciences 14, 223 (2010).
  * Palva  _et al._ [2013] J. M. Palva, A. Zhigalov, J. Hirvonen, O. Korhonen, K. Linkenkaer-Hansen, and S. Palva, Neuronal long-range temporal correlations and avalanche dynamics are correlated with behavioral scaling laws., Proceedings of the National Academy of Sciences of the United States of America 110, 3585 (2013).
  * Zhigalov  _et al._ [2015] A. Zhigalov, G. Arnulfo, L. Nobili, S. Palva, and J. Matias Palva, Relationship of fast- and slow-timescale neuronal dynamics in human MEG and SEEG, Journal of Neuroscience 35, 5385 (2015).
  * Dayan  _et al._ [1995] P. Dayan, G. E. Hinton, R. M. Neal, and R. S. Zemel, The Helmholtz Machine, Neural Computation 7, 889 (1995).
  * Mumford [1992] D. Mumford, On the computational architecture of the neocortex - II The role of cortico-cortical loops, Biological Cybernetics 66, 241 (1992).
  * Friston [2003] K. Friston, Learning and inference in the brain, Neural Networks 16, 1325 (2003).
  * Lee and Mumford [2003] T. S. Lee and D. Mumford, Hierarchical Bayesian inference in the visual cortex, Journal of the Optical Society of America A 20, 1434 (2003).
  * Friston [2005] K. Friston, A theory of cortical responses, Philosophical Transactions of the Royal Society B: Biological Sciences 360, 815 (2005).
  * Goldman-Rakic and Rakic [1991] P. S. Goldman-Rakic and P. Rakic, Preface: Cerebral cortex has come of age, Cerebral Cortex 1, 1 (1991).
  * Phillips and Singer [1997] W. A. Phillips and W. Singer, In search of common foundations for cortical computation (1997).
  * Kiebel  _et al._ [2008] S. J. Kiebel, J. Daunizeau, and K. J. Friston, A Hierarchy of Time-Scales and the Brain, PLoS Computational Biology 4, e1000209+ (2008).
  * Friston  _et al._ [2012] K. Friston, M. Breakspear, and G. Deco, Perception and self-organized instability, Frontiers in Computational Neuroscience 6, 10.3389/fncom.2012.00044 (2012).
  * Eurich  _et al._ [2002] C. W. Eurich, J. M. Herrmann, and U. A. Ernst, Finite-size effects of avalanche dynamics, Physical Review E 66, 15 (2002).
  * Levina  _et al._ [2007] A. Levina, J. M. Herrmann, and T. Geisel, Dynamical synapses causing self-organized criticality in neural networks, Nature Physics 3, 857 (2007).
  * Poil  _et al._ [2012] S.-S. S. Poil, R. Hardstone, H. D. Mansvelder, and K. Linkenkaer-Hansen, Critical-state dynamics of avalanches and oscillations jointly emerge from balanced excitation/inhibition in neuronal networks., Journal of Neuroscience 32, 9817 (2012).
  * De Arcangelis  _et al._ [2006] L. De Arcangelis, C. Perrone-Capano, and H. J. Herrmann, Self-organized criticality model for brain plasticity, Physical Review Letters 96, 10.1103/PhysRevLett.96.028107 (2006).
  * Bertschinger and Natschläger [2004] N. Bertschinger and T. Natschläger, Real-Time Computation at the Edge of Chaos in Recurrent Neural Networks, Neural Computation 16, 1413 (2004).
  * Friston  _et al._ [2007] K. J. Friston, J. Mattout, N. Trujillo-Barreto, J. Ashburner, and W. Penny, Variational free energy and the Laplace approximation, NeuroImage 34, 220 (2007).
  * Loreto  _et al._ [1996] V. Loreto, A. Vespignani, and S. Zapperi, Renormalization scheme for forest-fire models, Journal of Physics A: Mathematical and General 29, 2981 (1996).
  * Izhikevich [2007] E. M. Izhikevich,  _Dynamical systems in neuroscience : the geometry of excitability and bursting_, Computational neuroscience (MIT Press, 2007).
  * Note [1] The PC system is similar to a sequence of sluices or a stack of sōzus. A sōzu is a water-powered rocking device —a bamboo pipe pivoted to one side of its center of gravity. At rest its mouth is tilted upwards while its bottom lies on a support. A trickle of water fills it till its tipping point (the center of gravity is displaced past the pivot), when the pipe rotates around its pivot and dumps out _all_ the water it contained. A sōzu operates intermittently in the seesaw fashion of a relaxation oscillator [223] as long as the water flow persists. Thus, in a hierarchy of sōzus arranged so that each lies below the preceding, outflows and inflows are coupled sequentially, and the ensemble constitutes a cascade with level-wise intermittent falls. However, there are two important differences: water can only steadily build up, until it is discharged (there is no anti-water equivalent) and thresholds are fixed. Threshold adaptation could be accomplished by fastening the pivot shaft to the pipe, filling it with a dilatant (shear-thickening fluid), and making it rotate about another internal coaxial shaft in turn attached to the support. If the dilatant viscosity change had a long enough time constant, it would mimic approximately the behavior of adaptive thresholds.
  * De Los Rios and Zhang [1999] P. De Los Rios and Y. C. Zhang, Universal 1/f noise from dissipative SOC models, Physical Review Letters 82, 472 (1999).
  * Levina  _et al._ [2009] A. Levina, J. M. Herrmann, and T. Geisel, Phase transitions towards criticality in a neural system with adaptive interactions., Physical review letters 102, 118110 (2009).
  * Carlson and Doyle [1999] J. M. Carlson and J. Doyle, Highly optimized tolerance: A mechanism for power laws in designed systems, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 60, 1412 (1999).
  * Chialvo [2008] D. R. Chialvo, Emergent complexity: What uphill analysis or downhill invention cannot do, New Ideas in Psychology 26, 158 (2008), arXiv:0804.0008 .
  * Fairhall  _et al._ [2001] A. L. Fairhall, G. D. Lewen, W. Bialek, and D. Ruyter, Efficiency and ambiguity in an adaptive neural code, Nature 412, 787 (2001).
  * Kohn [2007] A. Kohn, Visual adaptation: Physiology, mechanisms, and functional benefits (2007).
  * Nessler  _et al._ [2013] B. Nessler, M. Pfeiffer, L. Buesing, and W. Maass, Bayesian Computation Emerges in Generic Cortical Microcircuits through Spike-Timing-Dependent Plasticity, PLoS Computational Biology 9, e1003037 (2013).
  * Bernacchia [2014] A. Bernacchia, The interplay of plasticity and adaptation in neural circuits: A generative model, Frontiers in Synaptic Neuroscience 6, 26 (2014).
  * Kauffman [1992] S. A. Kauffman, Origins of Order in Evolution: Self-Organization and Selection, in  _Understanding Origins_ (Springer, Dordrecht, 1992) pp. 153–181.
  * Platt  _et al._ [1993] N. Platt, E. A. Spiegel, and C. Tresser, On-off intermittency: A mechanism for bursting, Physical Review Letters 70, 279 (1993).
  * Evans and Majumdar [2011a] M. R. Evans and S. N. Majumdar, Diffusion with optimal resetting, Journal of Physics A: Mathematical and Theoretical 44, 10.1088/1751-8113/44/43/435001 (2011a), arXiv:1107.4225 .
  * Holland [1998] J. H. J. H. Holland, Emergence: from chaos to order, Choice Reviews Online 35, 35 (1998).
  * Dayan and Abbott [2001] P. Dayan and L. F. Abbott,  _Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems_, 1st ed. (The MIT Press, 2001).
  * Dehaene and Changeux [2005] S. Dehaene and J. P. Changeux, Ongoing spontaneous activity controls access to consciousness: A neuronal model for inattentional blindness, PLoS Biology 3, 0910 (2005).
  * Colaiori  _et al._ [2004] F. Colaiori, A. Baldassarri, and C. Castellano, Average trajectory of returning walks, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 69, 11 (2004), arXiv:0402285 [cond-mat] .
  * Schottky [1926] W. Schottky, Small-shot effect and flicker effect, Physical Review 28, 74 (1926).
  * Lowen and Teich [1993] S. B. Lowen and M. C. Teich, Fractal renewal processes generate 1/f noise, Physical Review E 47, 992 (1993).
  * Barnes and Allan [1966] J. A. Barnes and D. W. Allan, A Statistical Model of Flicker Noise, Proceedings of the IEEE 54, 176 (1966).
  * Mandelbrot and Van Ness [1968] B. B. Mandelbrot and J. W. Van Ness, Fractional Brownian Motions, Fractional Noises and Applications, SIAM Review 10, 422 (1968).
  * Sornette [2000] D. Sornette,  _Critical Phenomena in Natural Sciences_, Springer Series in Synergetics (Springer, Berlin, Heidelberg, 2000).
  * Zhang [1989] Y. C. Zhang, Scaling theory of self-organized criticality, Physical Review Letters 63, 470 (1989).
  * Straatsma  _et al._ [1986] T. P. Straatsma, H. J. Berendsen, and A. J. Stam, Estimation of statistical errors in molecular simulation calculations, Molecular Physics 57, 89 (1986).
  * Kass and Steffey [1989] R. E. Kass and D. Steffey, Approximate bayesian inference in conditionally independent hierarchical models (Parametric empirical bayes models), Journal of the American Statistical Association 84, 717 (1989).
  * Sethna  _et al._ [2001] J. P. Sethna, K. A. Dahmen, and C. R. Myers, Crackling noise (2001).
  * Durrett [2019] R. Durrett,  _Probability: Theory and Examples_, 5th ed. (Cambridge University Press, 2019).
  * Chialvo  _et al._ [2000] D. R. Chialvo, G. A. Cecchi, and M. O. Magnasco, Noise-induced memory in extended excitable systems, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 61, 5654 (2000).
  * Strogatz [2000] S. H. Strogatz, Nonlinear dynamics and chaos : with applications to physics, biology, chemistry, and engineering, Paperback (2000).
  * Geisel and Thomae [1984] T. Geisel and S. Thomae, Anomalous diffusion in intermittent chaotic systems, Physical Review Letters 52, 1936 (1984).
  * Ding and Yang [1995] M. Ding and W. Yang, Distribution of the first return time in fractional Brownian motion and its application to the study of on-off intermittency, Physical Review E 52, 207 (1995).
  * Note [2] We use the term metabolic energy (or metabolic power, for metabolic energy rate) when referring to the energy required by living organisms to maintain homeostasis, to distinguish it from squared activity.
  * Rao and Ballard [1999] R. P. N. Rao and D. H. Ballard, Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects, Nature Neuroscience 2, 79 (1999).
  * Hochstein and Ahissar [2002] S. Hochstein and M. Ahissar, View from the top: Hierarchies and reverse hierarchies in the visual system, Neuron 36, 791 (2002).
  * Friston  _et al._ [2008] K. J. Friston, N. Trujillo-Barreto, and J. Daunizeau, DEM: A variational treatment of dynamic systems, NeuroImage 41, 849 (2008).
  * Note [3] Bernstein-Von Mises theorem.
  * Beal and Ghahramani [2003] M. J. Beal and Z. Ghahramani, The Variational Bayesian EM Algorithm for Incomplete Data: with Application to Scoring Graphical Model Structures, in  _Bayesian Statistics_ , edited by J. M. Bernardo, M. J. Bayarri, J. O. Berger, A. P. Dawid, D. Heckerman, A. F. M. Smith, and M. West (Oxford University Press, 2003) Chap. 7.
  * Kauffman and Johnsen [1991] S. A. Kauffman and S. Johnsen, Coevolution to the edge of chaos: Coupled fitness landscapes, poised states, and coevolutionary avalanches, Journal of Theoretical Biology 149, 467 (1991).
  * Taleb [2010] N. N. Taleb,  _The Black Swan: The Impact of the Highly Improbable_ (Random House Trade Paperbacks, 2010) p. 480.
  * Raup [1986] D. M. Raup, Biological extinction in earth history (1986).
  * McKinney [1997] M. L. McKinney, Extinction vulnerability and selectivity: Combining ecological and paleontological views (1997).
  * Newman [2000] M. Newman, The power of design (2000).
  * Beggs [2019] J. M. Beggs, The Critically Tuned Cortex (2019).
  * Lennie [2003] P. Lennie, The cost of cortical computation, Current Biology 13, 493 (2003).
  * Olshausen and Field [1996] B. A. Olshausen and D. J. Field, Emergence of simple-cell receptive field properties by learning a sparse code for natural images., Nature 381, 607 (1996).
  * Stevens and Zador [1996] C. F. Stevens and A. Zador, Information through a spiking neuron, Advances in neural information processing systems , 75 (1996).
  * Note [4] Although information content and surprisal are usually synonyms to Shannon information, here information content denotes the entropy of the living system’s internal model, whereas surprisal denotes the fraction of that entropy attributable to surprisal, corresponding to mismatches between the sensory input and the expected sensory input.
  * Haken [1983] H. Haken,  _Synergetics_, Springer Series in Synergetics (Springer Berlin Heidelberg, Berlin, Heidelberg, 1983).
  * Kleidon [2009] A. Kleidon, Nonequilibrium thermodynamics and maximum entropy production in the Earth system: applications and implications. (2009).
  * Gamble [2012] C. Gamble,  _Settling the Earth: The Archaeology of Deep Human History_ (Cambridge University Press, 2012) pp. 1–377.
  * Bromer  _et al._ [2018] C. Bromer, T. M. Bartol, J. B. Bowden, D. D. Hubbard, D. C. Hanka, P. V. Gonzalez, M. Kuwajima, J. M. Mendenhall, P. H. Parker, W. C. Abraham, T. J. Sejnowski, and K. M. Harris, Long-term potentiation expands information content of hippocampal dentate gyrus synapses, Proceedings of the National Academy of Sciences of the United States of America 115, E2410 (2018).
  * Friston [1997a] K. J. Friston, Transients, metastability, and neuronal dynamics., NeuroImage 5, 164 (1997a).
  * Evans and Searles [2002] D. J. Evans and D. J. Searles, The fluctuation theorem, Advances in Physics 51, 1529 (2002).
  * Linsker [1990] R. Linsker, Perceptual neural organization: Some approaches based on network models and information theory, Annual Review of Neuroscience 13, 257 (1990).
  * Verghese and Pelli [1992] P. Verghese and D. G. Pelli, The information capacity of visual attention, Vision Research 32, 983 (1992).
  * Wu  _et al._ [2016] T. Wu, A. J. Dufford, M. A. Mackie, L. J. Egan, and J. Fan, The Capacity of Cognitive Control Estimated from a Perceptual Decision Making Task, Scientific Reports 6, 10.1038/srep34025 (2016).
  * Note [5] Note that the width of the conscious integration kernel (the interval over which events are integrated into a single percept) is $\approx 100\text{--}400$ ms [224, 225, 194], which suggests an information content for each “conscious snapshot” of $\approx 10$ b.
  * Hwa and Kardar [1989] T. Hwa and M. Kardar, Dissipative transport in open systems: An investigation of self-organized criticality, Physical Review Letters 62, 1813 (1989).
  * Manna  _et al._ [1990] S. S. Manna, L. B. Kiss, and J. Kertész, Cascades and self-organized criticality, Journal of Statistical Physics 61, 923 (1990).
  * Dickman  _et al._ [1998] R. Dickman, A. Vespignani, and S. Zapperi, Self-organized criticality as an absorbing-state phase transition, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 57, 5095 (1998), arXiv:9712115 [cond-mat] .
  * Efron and Morris [1973] B. Efron and C. Morris, Stein’s estimation rule and its competitors—an empirical bayes approach, Journal of the American Statistical Association 68, 117 (1973).
  * Fujimoto and Kaneko [2003] K. Fujimoto and K. Kaneko, How fast elements can affect slow dynamics, Physica D: Nonlinear Phenomena 180, 1 (2003).
  * Jirsa  _et al._ [1998] V. K. Jirsa, A. Fuchs, and J. A. Kelso, Connecting Cortical and Behavioral Dynamics: Bimanual Coordination, Neural Computation 10, 2019 (1998).
  * Zorick and Mandelkern [2015] T. Zorick and M. A. Mandelkern, Discrete scale invariance of human large eeg voltage deflections is more prominent in waking than sleep stage 2, Frontiers in Human Neuroscience 9, 638 (2015).
  * Jürgens  _et al._ [1995] E. Jürgens, F. Rösier, E. Hennighausen, and M. Heil, Stimulus-induced gamma oscillations: harmonics of alpha activity?, NeuroReport 6, 813 (1995).
  * Penttonen and Buzsáki [2003] M. Penttonen and G. Buzsáki, Natural logarithmic relationship between brain oscillators, Thalamus and Related Systems 2, 145 (2003).
  * Erland and Greenwood [2007] S. Erland and P. E. Greenwood, Constructing 1 $\omega$$\alpha$ noise from reversible Markov chains, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 76, 031114 (2007).
  * van Vreeswijk and Sompolinsky [1996] C. van Vreeswijk and H. Sompolinsky, Chaos in neuronal networks with balanced excitatory and inhibitory activity, Science 274, 1724 (1996).
  * Friston  _et al._ [2010] K. J. Friston, J. Daunizeau, J. Kilner, and S. J. Kiebel, Action and behavior: A free-energy formulation, Biological Cybernetics 102, 227 (2010).
  * Blakemore  _et al._ [2002] S. J. Blakemore, D. M. Wolpert, and C. D. Frith, Abnormalities in the awareness of action, Trends in Cognitive Sciences 6, 237 (2002).
  * Friston [2018] K. J. Friston, Am I Self-Conscious? (Or Does Self-Organization Entail Self-Consciousness?), Frontiers in Psychology 9, 579+ (2018).
  * Martinez-Saito [2021] M. Martinez-Saito, A Skeptical View on the Physics-Consciousness Explanatory Gap, Axiomathes , 1 (2021).
  * Lombardi  _et al._ [2014] F. Lombardi, H. J. Herrmann, D. Plenz, and L. De Arcangelis, On the temporal organization of neuronal avalanches, Frontiers in Systems Neuroscience 8, 204 (2014).
  * Srinivasan  _et al._ [1996] R. Srinivasan, P. L. Nunez, D. M. Tucker, R. B. Silberstein, and P. J. Cadusch, Spatial sampling and filtering of EEG with spline laplacians to estimate cortical potentials, Brain Topography 8, 355 (1996).
  * Freeman  _et al._ [2000] W. J. Freeman, L. J. Rogers, M. D. Holmes, and D. L. Silbergeld, Spatial spectral analysis of human electrocorticograms including the alpha and gamma bands, Journal of Neuroscience Methods 95, 111 (2000).
  * Bédard  _et al._ [2006] C. Bédard, H. Kröger, and A. Destexhe, Does the 1/f frequency scaling of brain signals reflect self-organized critical states?, Physical Review Letters 97, 118102 (2006), arXiv:0608026 [q-bio] .
  * Miller  _et al._ [2009] K. J. Miller, L. B. Sorensen, J. G. Ojemann, and M. Den Nijs, Power-law scaling in the brain surface electric potential, PLoS Computational Biology 5, e1000609 (2009).
  * Smit  _et al._ [2011] D. J. Smit, E. J. de Geus, M. E. van de Nieuwenhuijzen, C. E. van Beijsterveldt, G. C. M. van Baal, H. D. Mansvelder, D. I. Boomsma, and K. Linkenkaer-Hansen, Scale-free modulation of resting-state neuronal oscillations reflects prolonged brain maturation in humans, Journal of Neuroscience 31, 13128 (2011).
  * Peng  _et al._ [1994] C. K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger, Mosaic organization of DNA nucleotides, Physical Review E 49, 1685 (1994).
  * Engel  _et al._ [2001] A. K. Engel, P. Fries, and W. Singer, Dynamic predictions: Oscillations and synchrony in top–down processing, Nature Reviews Neuroscience 2, 704 (2001).
  * Buzsáki and Draguhn [2004] G. Buzsáki and A. Draguhn, Neuronal oscillations in cortical networks., Science 304, 1926 (2004).
  * Kopell  _et al._ [2000] N. Kopell, G. B. Ermentrout, M. A. Whittington, and R. D. Traub, Gamma rhythms and beta rhythms have different synchronization properties, Proceedings of the National Academy of Sciences of the United States of America 97, 1867 (2000).
  * Von Stein  _et al._ [1999] A. Von Stein, P. Rappelsberger, J. Sarnthein, and H. Petsche, Synchronization between temporal and parietal cortex during multimodal object processing in man, Cerebral Cortex 9, 137 (1999).
  * Varela  _et al._ [2001] F. Varela, J.-P. P. Lachaux, E. Rodriguez, and J. Martinerie, The brainweb: Phase synchronization and large-scale integration, Nature Reviews Neuroscience 2, 229 (2001).
  * Hahn  _et al._ [2006] T. T. Hahn, B. Sakmann, and M. R. Mehta, Phase-locking of hippocampal interneurons’ membrane potential to neocortical up-down states, Nature Neuroscience 9, 1359 (2006).
  * Gireesh and Plenz [2008] E. D. Gireesh and D. Plenz, Neuronal avalanches organize as nested theta- and beta/gamma-oscillations during development of cortical layer 2/3, Proceedings of the National Academy of Sciences of the United States of America 105, 7576 (2008).
  * He  _et al._ [2010] B. J. He, J. M. Zempel, A. Z. Snyder, and M. E. Raichle, The temporal structures and functional significance of scale-free brain activity, Neuron 66, 353 (2010).
  * Yu and Margoliash [1996] A. C. Yu and D. Margoliash, Temporal hierarchical control of singing in birds, Science 273, 1871 (1996).
  * Hasson  _et al._ [2008] U. Hasson, E. Yang, I. Vallines, D. J. Heeger, and N. Rubin, A Hierarchy of Temporal Receptive Windows in Human Cortex, Journal of Neuroscience 28, 2539 (2008).
  * Hasson  _et al._ [2015] U. Hasson, J. Chen, and C. J. Honey, Hierarchical process memory: Memory as an integral component of information processing (2015).
  * Mesulam [1998] M.-M. Mesulam, From sensation to cognition, Brain 121, 1013 (1998).
  * Sherman and Guillery [1998] S. M. Sherman and R. W. Guillery, On the actions that one nerve cell can have on another: Distinguishing ”drivers” from ”modulators”, Proceedings of the National Academy of Sciences of the United States of America 95, 7121 (1998).
  * Huk and Shadlen [2005] A. C. Huk and M. N. Shadlen, Neural Activity in Macaque Parietal Cortex Reflects Temporal Integration of Visual Motion Signals during Perceptual Decision Making, Journal of Neuroscience 25, 10420 (2005).
  * Poeppel  _et al._ [2008] D. Poeppel, W. J. Idsardi, and V. Van Wassenhove, Speech perception at the interface of neurobiology and linguistics (2008).
  * Quian Quiroga  _et al._ [2005] R. Quian Quiroga, L. Reddy, G. Kreiman, C. Koch, I. Fried, R. Q. Quiroga, L. Reddy, G. Kreiman, C. Koch, and I. Fried, Invariant visual representation by single neurons in the human brain, Nature 435, 1102 (2005).
  * Friston [2010] K. J. Friston, The free-energy principle: a unified brain theory?, Nature Reviews Neuroscience 11, 127 (2010).
  * West and Gundersen [1990] M. J. West and H. J. Gundersen, Unbiased stereological estimation of the number of neurons in the human hippocampus, Journal of Comparative Neurology 296, 1 (1990).
  * Buzsáki [1996] G. Buzsáki, The hippocampo-neocortical dialogue (1996).
  * Friston and Buzsáki [2016] K. Friston and G. Buzsáki, The Functional Anatomy of Time: What and When in the Brain (2016).
  * Thompson and Kim [1996] R. F. Thompson and J. J. Kim, Memory systems in the brain and localization of a memory, Proceedings of the National Academy of Sciences of the United States of America 93, 13438 (1996).
  * Honey  _et al._ [2012] C. J. Honey, T. Thesen, T. H. Donner, L. J. Silbert, C. E. Carlson, O. Devinsky, W. K. Doyle, N. Rubin, D. J. Heeger, and U. Hasson, Slow Cortical Dynamics and the Accumulation of Information over Long Timescales, Neuron 76, 423 (2012).
  * Murray  _et al._ [2014] J. D. Murray, A. Bernacchia, D. J. Freedman, R. Romo, J. D. Wallis, X. Cai, C. Padoa-Schioppa, T. Pasternak, H. Seo, D. Lee, and X. J. Wang, A hierarchy of intrinsic timescales across primate cortex, Nature Neuroscience 17, 1661 (2014).
  * Paczuski  _et al._ [1996] M. Paczuski, S. Maslov, and P. Bak, Avalanche dynamics in evolution, growth, and depinning models, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 53, 414 (1996), arXiv:9510002 [adap-org] .
  * Dong and Atick [1995] D. W. Dong and J. J. Atick, Statistics of natural time-varying images, Network: Computation in Neural Systems 6, 345 (1995).
  * Ruderman [1997] D. L. Ruderman, Origins of scaling in natural images, Vision Research 37, 3385 (1997).
  * Billock  _et al._ [2001] V. A. Billock, G. C. De Guzman, and J. A. Scott Kelso, Fractal time and 1/f spectra in dynamic images and human vision, Physica D: Nonlinear Phenomena 148, 136 (2001).
  * Kitzbichler  _et al._ [2009] M. G. Kitzbichler, M. L. Smith, S. R. Christensen, and E. Bullmore, Broadband criticality of human brain network synchronization, PLoS Computational Biology 5, e1000314 (2009).
  * Priesemann  _et al._ [2014] V. Priesemann, M. Wibral, M. Valderrama, R. Pröpper, M. Le Van Quyen, T. Geisel, J. Triesch, D. Nikolić, and M. H. Munk, Spike avalanches in vivo suggest a driven, slightly subcritical brain state, Frontiers in Systems Neuroscience 8, 108 (2014).
  * Gollo [2017] L. L. Gollo, Coexistence of critical sensitivity and subcritical specificity can yield optimal population coding, Journal of the Royal Society Interface 14, 10.1098/rsif.2017.0207 (2017), arXiv:1707.04484 .
  * Freeman [1987] W. J. Freeman, Simulation of chaotic EEG patterns with a dynamic model of the olfactory system, Biological Cybernetics 56, 139 (1987).
  * Bonachela  _et al._ [2010] J. A. Bonachela, S. De Franciscis, J. J. Torres, and M. A. Munoz, Self-organization without conservation: Are neuronal avalanches generically critical?, Journal of Statistical Mechanics: Theory and Experiment 2010, P02015 (2010), arXiv:1001.3256 .
  * Shew  _et al._ [2015] W. L. Shew, W. P. Clawson, J. Pobst, Y. Karimipanah, N. C. Wright, and R. Wessel, Adaptation to sensory input tunes visual cortex to criticality, Nature Physics 11, 659 (2015).
  * Moretti and Muñoz [2013] P. Moretti and M. A. Muñoz, Griffiths phases and the stretching of criticality in brain networks, Nature Communications 4, 1 (2013).
  * Larremore  _et al._ [2014] D. B. Larremore, W. L. Shew, E. Ott, F. Sorrentino, and J. G. Restrepo, Inhibition causes ceaseless dynamics in networks of excitable nodes, Physical Review Letters 112, 138103 (2014).
  * Levina  _et al._ [2014] A. Levina, J. M. Herrmann, and T. Geisel, Theoretical Neuroscience of Self-Organized Criticality: From Formal Approaches to Realistic Models, in  _Criticality in Neural Systems_, Vol. 9783527411 (Wiley Blackwell, 2014) pp. 417–436.
  * De Andrade Costa  _et al._ [2015] A. De Andrade Costa, M. Copelli, and O. Kinouchi, Can dynamical synapses produce true self-organized criticality?, Journal of Statistical Mechanics: Theory and Experiment 2015, P06004 (2015).
  * Beggs [2008] J. M. Beggs, The criticality hypothesis: How local cortical networks might optimize information processing, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 366, 329 (2008).
  * Beggs and Timme [2012] J. M. Beggs and N. Timme, Being critical of criticality in the brain, Frontiers in Physiology 3 JUN, 163 (2012).
  * Cocchi  _et al._ [2017] L. Cocchi, L. L. Gollo, A. Zalesky, and M. Breakspear, Criticality in the brain: A synthesis of neurobiology, models and cognition (2017), arXiv:1707.05952 .
  * Shriki  _et al._ [2013] O. Shriki, J. Alstott, F. Carver, T. Holroyd, R. N. A. Henson, M. L. Smith, R. Coppola, E. Bullmore, and D. Plenz, Neuronal Avalanches in the Resting MEG of the Human Brain, Journal of Neuroscience 33, 7079 (2013).
  * Scott  _et al._ [2014] G. Scott, E. D. Fagerholm, H. Mutoh, R. Leech, D. J. Sharp, W. L. Shew, and T. Knöpfel, Voltage imaging of waking mouse cortex reveals emergence of critical neuronal dynamics, Journal of Neuroscience 34, 16611 (2014).
  * Tagliazucchi  _et al._ [2012] E. Tagliazucchi, P. Balenzuela, D. Fraiman, and D. R. Chialvo, Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis, Frontiers in Physiology 3 FEB, 15 (2012).
  * Petermann  _et al._ [2009] T. Petermann, T. C. Thiagarajan, M. A. Lebedev, M. A. Nicolelis, D. R. Chialvo, and D. Plenz, Spontaneous cortical activity in awake monkeys composed of neuronal avalanches., Proceedings of the National Academy of Sciences of the United States of America 106, 15921 (2009).
  * McWhorter [1957] A. McWhorter, 1/f noise and related surface effects in germanium., in  _Semiconductor Surface Physics_ (MIT Lincoln Laboratory, 1957) pp. 207–228.
  * Bastos  _et al._ [2012] A. M. Bastos, W. M. Usrey, R. A. Adams, G. R. Mangun, P. Fries, and K. J. Friston, Canonical Microcircuits for Predictive Coding (2012).
  * Stewart and Plenz [2006] C. V. Stewart and D. Plenz, Inverted-U profile of dopamine-NMDA-mediated spontaneous avalanche recurrence in superficial layers of rat prefrontal cortex, Journal of Neuroscience 26, 8148 (2006).
  * Fagerholm  _et al._ [2016] E. D. Fagerholm, G. Scott, W. L. Shew, C. Song, R. Leech, T. Knöpfel, and D. J. Sharp, Cortical Entropy, Mutual Information and Scale-Free Dynamics in Waking Mice, Cerebral Cortex 26, 3945 (2016).
  * Destexhe  _et al._ [1999] A. Destexhe, D. Contreras, and M. Steriade, Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states, Journal of Neuroscience 19, 4595 (1999).
  * Srinivasan  _et al._ [1999] R. Srinivasan, D. P. Russell, G. M. Edelman, and G. Tononi, Increased synchronization of neuromagnetic responses during conscious perception, Journal of Neuroscience 19, 5435 (1999).
  * Blythe and Nikulin [2017] D. A. Blythe and V. V. Nikulin, Long-range temporal correlations in neural narrowband time-series arise due to critical dynamics, PLoS ONE 12, e0175628 (2017).
  * Zapperi  _et al._ [1995] S. Zapperi, K. B. Lauritsen, and H. E. Stanley, Self-organized branching processes: Mean-field theory for avalanches, Physical Review Letters 75, 4071 (1995), arXiv:9509161 [cond-mat] .
  * Kuntz and Sethna [2000] M. C. Kuntz and J. P. Sethna, Noise in disordered systems: The power spectrum and dynamic exponents in avalanche models, Physical Review B - Condensed Matter and Materials Physics 62, 11699 (2000), arXiv:9911207 [cond-mat] .
  * Muñoz  _et al._ [1999] M. A. Muñoz, R. Dickman, A. Vespignani, and S. Zapperi, Avalanche and spreading exponents in systems with absorbing states, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 59, 6175 (1999), arXiv:9811287 [cond-mat] .
  * Dickman  _et al._ [2000] R. Dickman, M. A. Muñoz, A. Vespignani, and S. Zapperi, Paths to self-organized criticality, Brazilian Journal of Physics 30, 27 (2000), arXiv:9910454 [cond-mat] .
  * Lübeck [2004] S. Lübeck, Universal scaling behavior of non-equilibrium phase transitions (2004), arXiv:0501259 [cond-mat] .
  * Michiels Van Kessenich  _et al._ [2016] L. Michiels Van Kessenich, L. De Arcangelis, and H. J. Herrmann, Synaptic plasticity and neuronal refractory time cause scaling behaviour of neuronal avalanches, Scientific Reports 6, 10.1038/srep32071 (2016).
  * Plenz and Thiagarajan [2007] D. Plenz and T. C. Thiagarajan, The organizing principles of neuronal avalanches: cell assemblies in the cortex?, Trends in neurosciences 30, 101 (2007).
  * Gold and Shadlen [2007] J. I. Gold and M. N. Shadlen, The Neural Basis of Decision Making, Annual Review of Neuroscience 30, 535 (2007).
  * Cowan and Wilson [1994] R. L. Cowan and C. J. Wilson, Spontaneous firing patterns and axonal projections of single corticostriatal neurons in the rat medial agranular cortex, Journal of Neurophysiology 71, 17 (1994).
  * Wilson and Kawaguchi [1996] C. J. Wilson and Y. Kawaguchi, The origins of two-state spontaneous membrane potential fluctuations of neostriatal spiny neurons, Journal of Neuroscience 16, 2397 (1996).
  * Paré  _et al._ [1998] D. Paré, E. Shink, H. Gaudreau, A. Destexhe, and E. J. Lang, Impact of spontaneous synaptic activity on the resting properties of cat neocortical pyramidal neurons in vivo, Journal of Neurophysiology 79, 1450 (1998).
  * Mihalas  _et al._ [2014] S. Mihalas, D. Millman, R. Iyer, A. Kirkwood, and E. Niebur, Nonconservative Neuronal Networks During Up-States Self-Organize Near Critical Points, in  _Criticality in Neural Systems_, Vol. 9783527411 (Wiley-Blackwell, 2014) pp. 437–464.
  * Cossart  _et al._ [2003] R. Cossart, D. Aronov, and R. Yuste, Attractor dynamics of network UP states in the neocortex, Nature 423, 283 (2003).
  * Wilson [2008] C. Wilson, Up and down states, Scholarpedia 3, 1410 (2008).
  * Lombardi  _et al._ [2012] F. Lombardi, H. J. Herrmann, C. Perrone-Capano, D. Plenz, and L. De Arcangelis, Balance between excitation and inhibition controls the temporal organization of neuronal avalanches, Physical Review Letters 108, 228703 (2012).
  * Freyer  _et al._ [2011] F. Freyer, J. A. Roberts, R. Becker, P. A. Robinson, P. Ritter, and M. Breakspear, Biophysical mechanisms of multistability in resting-state cortical rhythms, Journal of Neuroscience 31, 6353 (2011).
  * Friston [1997b] K. J. Friston, Another neural code?, NeuroImage 5, 213 (1997b).
  * Rabinovich  _et al._ [2008] M. Rabinovich, R. Huerta, and G. Laurent, Neuroscience: Transient dynamics for neural processing, Science 321, 48 (2008).
  * Varela [1999] F. J. Varela, The specious present: A neurophenomenology of time consciousness, Naturalizing phenomenology: Issues in contemporary phenomenology and cognitive science LXIV, 266 (1999).
  * Eagleman and Sejnowski [2000] D. M. Eagleman and T. J. Sejnowski, Motion Integration and Postdiction in Visual Awareness, Science 287, 2036 (2000).
  * Herzog  _et al._ [2016] M. H. Herzog, T. Kammer, and F. Scharnowski, Time Slices: What Is the Duration of a Percept?, PLoS Biol 14, e1002433+ (2016).
  * Wolpert and Kawato [1998] D. M. Wolpert and M. Kawato, Multiple paired forward and inverse models for motor control., Neural Networks 11, 1317 (1998).
  * Wolpert  _et al._ [1998] D. M. Wolpert, S. J. Goodbody, and M. Husain, Maintaining internal representations: the role of the human superior parietal lobe, Nature Neuroscience 1, 529 (1998).
  * Leech  _et al._ [2012] R. Leech, R. Braga, and D. J. Sharp, Echoes of the brain within the posterior cingulate cortex, Journal of Neuroscience 32, 215 (2012).
  * Biswal  _et al._ [1995] B. Biswal, F. Zerrin Yetkin, V. M. Haughton, and J. S. Hyde, Functional connectivity in the motor cortex of resting human brain using echo‐planar mri, Magnetic Resonance in Medicine 34, 537 (1995).
  * Achard  _et al._ [2006] S. Achard, R. Salvador, B. Whitcher, J. Suckling, and E. Bullmore, A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs., The Journal of neuroscience : the official journal of the Society for Neuroscience 26, 63 (2006).
  * Fox and Raichle [2007] M. D. Fox and M. E. Raichle, Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging, Nat Rev Neurosci 8, 700 (2007).
  * Haimovici  _et al._ [2013] A. Haimovici, E. Tagliazucchi, P. Balenzuela, and D. R. Chialvo, Brain organization into resting state networks emerges at criticality on a model of the human connectome, Physical Review Letters 110, 178101 (2013), arXiv:1209.5353 .
  * Gollo  _et al._ [2015] L. L. Gollo, A. Zalesky, R. Matthew Hutchison, M. Van Den Heuvel, and M. Breakspear, Dwelling quietly in the rich club: Brain network determinants of slow cortical fluctuations, Philosophical Transactions of the Royal Society B: Biological Sciences 370, 10.1098/rstb.2014.0165 (2015).
  * Mandelbrot [1974] B. B. Mandelbrot, Intermittent turbulence in self-similar cascades: Divergence of high moments and dimension of the carrier, Journal of Fluid Mechanics 62, 331 (1974).
  * Ihlen and Vereijken [2010] E. A. Ihlen and B. Vereijken, Interaction-dominant dynamics in human cognition: Beyond 1/f $\alpha$fluctuation, Journal of Experimental Psychology: General 139, 436 (2010).
  * Note [6] Alternatively to A, $f_{A_{l}}$ could be calculated as $f_{A_{l}}=B_{\theta_{l}}[f_{\tilde{A}_{l}}]q_{l}^{(0)}+(1-q_{l}^{(0)})(B_{\theta_{l}}[f_{\tilde{A}_{l}^{(1)}}]q_{l}^{(1)}+(1-q_{l}^{(1)})(\ldots))$, with $f_{\tilde{A}_{l}^{(t)}}=B_{\theta_{l}}[f_{\tilde{A}_{l}}^{(t-1)}]*I$, $f_{\tilde{A}_{l}^{(0)}}=f_{\mathcal{E}_{l-1}}$, and $f_{\tilde{A}_{l}^{(0)}}=I$.
  * Pietronero  _et al._ [1988] L. Pietronero, A. Erzan, and C. Evertsz, Theory of Laplacian fractals: Diffusion limited aggregation and dielectric breakdown model, Physica A: Statistical Mechanics and its Applications 151, 207 (1988).
  * Stein [1965] R. B. Stein, A Theoretical Analysis of Neuronal Variability, Biophysical Journal 5, 173 (1965).
  * Ditlevsen and Lansky [2005] S. Ditlevsen and P. Lansky, Estimation of the input parameters in the Ornstein-Uhlenbeck neuronal model, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 71 (2005).
  * Sornette  _et al._ [2004] D. Sornette, F. Deschâtres, T. Gilbert, and Y. Ageon, Endogenous Versus Exogenous Shocks in Complex Networks: An Empirical Test Using Book Sale Rankings, Physical Review Letters 93, 228701 (2004).
  * Taylor and Karlin [1998] H. M. Taylor and S. Karlin,  _An Introduction to Stochastic Modeling, 3rd ed._ (Academic Press Limited, 1998).
  * Note [7] Notice this implies that for fixed $t$, the distribution of RBM coincides with the running maximum distribution for Brownian motion $f_{|B^{(t)}|}(x)=f_{\mathop{max}\displaylimits_{0\leq u\leq t}W^{(u)}}(x)$.
  * Note [8] This is a common procedure in finance for pricing double barrier options.
  * Faulhaber [2002] O. Faulhaber,  _Analytic methods for pricing double barrier options in the presence of stochastic volatility_ , Ph.D. thesis, University of Kaiserslautern, Germany (2002).
  * Markushevich [1950] A. I. Markushevich,  _Theory of analytic functions_ (1950).
  * Shumway and Stoffer [2006] R. H. Shumway and D. S. Stoffer,  _Time Series Analysis and Its Applications: With R Examples (Springer Texts in Statistics)_, 2nd ed. (Springer, 2006).
  * Note [9] The expected mean $\langle\theta^{(t+n)}|\theta^{(t)}\rangle$ and variance $\operatorname{var}{(\theta^{(t+n)}|\theta^{(t)})}$ after $n$ iterations can be derived explicitly by iterated application of Eq. 67 $\langle\theta^{(t+n)}|\theta^{(t)}\rangle=\langle\Theta\rangle(1-(1-w)^{n})+\theta^{(t)}(1-w)^{n}$ and $\operatorname{var}{(\theta^{(t+n)}|\theta^{(t)})}\approx(1-(1-w)^{2n})\frac{w}{2}\operatorname{var}{|\tilde{a}^{(t)}|}$.
  * Note [10] Specifically, the non-unitary Fourier transform, in ordinary frequency units.
  * Note [11] Cf. variance of pulses $\sigma_{\mathcal{E}}=0.35\langle\Theta_{l}\rangle$.
  * Percival [1994] D. B. Percival, Spectral Analysis of Univariate and Bivariate Time Series, in  _Methods in Experimental Physics_, Vol. 28 (Academic Press, 1994) pp. 313–348.
  * Evans and Majumdar [2011b] M. R. Evans and S. N. Majumdar, Diffusion with stochastic resetting, Physical Review Letters 106, 160601 (2011b), arXiv:1102.2704 .
  * Majumdar and Oshanin [2018] S. N. Majumdar and G. Oshanin, Spectral content of fractional Brownian motion with stochastic reset, Journal of Physics A: Mathematical and Theoretical 51, 435001 (2018).
  * Coles [2001] S. Coles,  _An Introduction to Statistical Modeling of Extreme Values_, Springer Series in Statistics (Springer London, London, 2001).
  * Ginoux and Letellier [2012] J. M. Ginoux and C. Letellier, Van der Pol and the history of relaxation oscillations: Toward the emergence of a concept, Chaos 22, 023120 (2012), arXiv:1408.4890 .
  * Efron [1970] R. Efron, The relationship between the duration of a stimulus and the duration of a perception, Neuropsychologia 8, 37 (1970).
  * Burr [1980] D. Burr, Motion smear, Nature 284, 164 (1980).