PHYSICAL REVIEW LETTERS 125, 028101 (2020)

Edge of Chaos and Avalanches in Neural Networks with Heavy-Tailed
Synaptic Weight Distribution1Łukasz Kuśmierz ,1,* Shun Ogawa ,1 and Taro Toyoizumi

1,2,†

Laboratory for Neural Computation and Adaptation, RIKEN Center for Brain Science, 2-1 Hirosawa,
Wako, Saitama 351-0198, Japan2
Department of Mathematical Informatics, Graduate School of Information Science and Technology,
The University of Tokyo, Tokyo 113-8656, Japan(Received 12 November 2019; revised 3 March 2020; accepted 26 May 2020; published 7 July 2020)

We propose an analytically tractable neural connectivity model with power-law distributed synaptic
strengths. When threshold neurons with biologically plausible number of incoming connections are
considered, our model features a continuous transition to chaos and can reproduce biologically relevant low
activity levels and scale-free avalanches, i.e., bursts of activity with power-law distributions of sizes andlifetimes. In contrast, the Gaussian counterpart exhibits a discontinuous transition to chaos and thus cannot
be poised near the edge of chaos. We validate our predictions in simulations of networks of binary as well as
leaky integrate-and-fire neurons. Our results suggest that heavy-tailed synaptic distribution may form aweakly informative sparse-connectivity prior that can be useful in biological and artificial adaptive systems.
DOI: 10.1103/PhysRevLett.125.028101

Scale-free neuronal avalanches, commonly associatedwith criticality, have been observed in cortical networks invarious settings, including cultured and acute slices fromrat somatosensory cortex [1], eye-attached ex vivo preparation of turtle visual cortex [2], visual cortex in anesthetized rats [3], primary visual cortex in anesthetizedmonkeys [3], and premotor, motor, and somatosensory
cortex in awake monkeys [4]. Criticality implies theexistence of a continuous transition between two distinct
collective phases. In the context of neuronal avalanches,most commonly studied transitions are between quiescentand active states [1,5,6] or synchronuous and asynchronuous states [3,7]. In addition to providing a plausible
generating mechanism for the neuronal avalanches, theexistence of a continuous transition would have importantfunctional implications, as it has been shown that computation is most efficient around a critical point [6,8,9], often
associated with the edge of chaos [10–14]. However, therelation between neuronal avalanches, criticality, and theedge of chaos is not fully understood [9,15,16].
Different scenarios of the transition to chaos in randomly
connected neural networks were extensively studied overthe last 30 years [14,17–27]. According to the prevailingassumption rooted in the central limit theorem, the total
synaptic input current of each neuron can be modeled as a
Gaussian random variable (Gaussian assumption). Here we
argue that the Gaussian assumption cannot account forsome of the experimentally observed features of neuronal
circuits.
In particular, the continuous nature of the phase transition observed in the conventional models is sensitive totheoretical assumptions that are not biologically grounded.
0031-9007=20=125(2)=028101(6)Most works that study transition to chaos employ rate
models with continuous nonthresholded activation functions, often of a sigmoidal shape [14,17,18,21,23,24,27].
Sometimes thresholds are introduced, but often the analysis
is restricted to the suprathreshold regime [25,28]. But mostneurons in the brain spike only when driven by strong
enough excitatory synaptic input above a threshold [29–31].
Thus, we model a self-sustained (autonomous) activity in a
network of individually subthreshold neurons. Other models
exhibiting continuous transition to chaos [11] or neuronalavalanches [6,32] rely on extremely sparse (∼10 connections
per neuron) networks. However, many neurons in thevertebrate brain receive a large number of inputs from otherneurons (∼104 ) [33]. We observed that the transition to chaos
becomes discontinuous when densely connected subthreshold units are used in tandem with the Gaussian assumption(Fig. 3) [20,34]. This discontinuous character of the transition makes it hard for the network to robustly exhibit theedge of chaos, low activity levels, or avalanches. Although a
discontinuous transition can be smoothed by noise, leading
again to critical behavior away from the edge of chaos if thenoise level is appropriate [35], such noise-induced criticality
requires extra fine-tuning. We explore instead the possibility
that an autonomous network exhibits critical behavior at theedge of chaos.
To fix this issue, we draw on the experimental works
reporting heavy-tailed distributions of synaptic weights invarious areas of the brain [36–41]. Multiple theoreticalmechanisms have been suggested to realize such distributions, e.g., modified spike-timing-dependent plasticity(STDP) rule [42] or STDP combined with homeostatic
plasticity [43]. Notably, recent studies have suggested that028101-1© 2020 American Physical Society

PHYSICAL REVIEW LETTERS 125, 028101 (2020)
experimentally observed activity-independent intrinsic spine
dynamics can straightforwardly explain the heavy-tailed
distributions of synaptic weights [44–48].
Although extensively studied, the computational roleof synaptic heavy tails is still not fully understood. A lognormal distribution is often assumed and the results are
obtained by means of computer simulations [40,49–52].
Unfortunately, the log-normal distribution is not a stable
distribution [53]. Consequently, the corresponding distribution of the membrane potential depends in a nontrivial
way on the details of the connectivity, including number ofincoming connections. Moreover, if the number of incoming
connections is scaled linearly with the number of neurons,
the Gaussian assumption is recovered in the thermodynamic
limit. This hinders theoretical approaches to study the effectsof heavy tails in these models. Therefore, a simple modelthat robustly predicts the effects of synaptic heavy tails isneeded. We fill this gap by assuming random, power-law
distributed synaptic weights.
Our aim is to inspect how the distribution of synaptic
efficacies, modulated by the activation function, affects thetransition to chaos and the associated avalanches. To thisend, in our calculations we focus on the network effects andhence simplify the dynamics of individual neurons byconsidering the following discrete-time network dynamics
xi ðt þ 1Þ ¼

NX

Jij ϕ(xj ðtÞ);

ð1Þ

j¼1

where ϕðxÞ is the activation function, assumed to be
identical across the network, and J is the connectivity
matrix. The network is fully connected and the synaptic
weights are independently drawn from the common
Cauchy distribution (Fig. 1)
ρðJ ij Þ ¼

1g=N;
π ðg=NÞ2 þ J2ijð2Þwith the characteristic functionΦJ ðkÞ ¼ e−γjkj ;ð3Þwhere γ ¼ g=N defines the width of the distribution.
We refer to the model prescribed by (1) and (2) as the
Cauchy network.
Due to the generalized central limit theorem [53], in thethermodynamic limit of N → ∞ results obtained for the
Cauchy model are applicable to networks with connections
drawn independently from any symmetric distribution with
1=x2 tails that are scaled with the number of neurons as
1=N. In contrast, in the more commonly used Gaussian
networks, the synaptic weights are independently drawn
from the normal distribution Jij ∼ N ð0; g2 =NÞ. In thethermodynamic limit this corresponds to connectivity
matrices with entries independently drawn from any distribution with zero mean andpﬃﬃﬃﬃa finite variance, as long as theweights are scaled as 1= N .

100
Prob.(|Ji j | > x)

Densities

0.3
0.2
0.1
0.0

−6

−3

0
Ji j

3

6

CauchyGaussian

10−1
10−2
10−3
10−4

100

101

102

103

x

FIG. 1. (Top) Visualizations of neural networks with Gaussian(cyan) and Cauchy (orange) distribution of weights. Thickness
and color saturation of edges correspond to the (nonlinearlytransformed) strengths of the connections. (Middle) Probability
density functions (left) and cumulative distribution functions(right) of Cauchy and Gaussian random variables. The Cauchy
distribution features much thicker tails than the Gaussian distribution. (Bottom) Sample realizations of the Poisson criticalbranching process with the duration T ¼ 11, size S ¼ 18 (left),and T ¼ 35, S ¼ 367 (right). The initial seeds are marked withthe green color. In this Letter, we show that activity of a fully
connected Cauchy (but not Gaussian) network around the criticalpoint can be mapped to the critical branching process.

The natural order parameter in the system at hand is themean network activity, defined as
mt ¼

N
1X
jϕ(xi ðtÞ)j:
N i¼1

ð4Þ

The state at time t þ 1 depends on J and xðtÞ. We fix theactivity vector at time t and treat xðt þ 1Þ as a function of J,
which allows us to characterize the distribution of xi ðt þ 1Þusing Φxi ðtþ1Þ ðkÞ asNX
ikxi ðtþ1Þ
−1
iJ ¼ exp −gjkjNjϕ(xj ðtÞ)jhej¼1¼ expð−gmt jkjÞ:ð5Þ
The activity of a neuron at time t þ 1, as a function ofsynaptic weights, is a Cauchy random variable whose

028101-2

PHYSICAL REVIEW LETTERS 125, 028101 (2020)
width depends on the activity at time t only through its
mean value.
To proceed we assume self-averaging, i.e., that the mean
activity is the same for each realization of the network.
Since in our model synaptic weights are statistically thesame for all neurons, in the limit of N → ∞ the mean
activity can alternatively be expressed as
mt ¼ hjϕðxi ðtÞÞjiJ

ð ∀ iÞ:ð6Þ

We use the result of (5), i.e., that xi ðtÞ averaged over J is a
Cauchy variable with γ ¼ gmt , together with (6) and arrive at
the evolution of the mean activity in a simple integral form
Z ∞mtþ1 ¼
Dzjϕðgmt zÞj;ð7Þ−∞

where Dz ¼ π −1 dz=ð1 þ z2 Þ denotes that the integral iscalculated with respect to the standard Cauchy measure.
The steady-state mean activity can be obtained from (7) in aself-consistent manner.
We are now in the position to analyze the dependence ofthe dynamics of the Cauchy network on the activation
function. For ϕðxÞ ¼ x, the integral on the right-hand sideof (7) diverges, suggesting that the network is unstable.
Indeed, it is easy to understand why this is the case. For
linear networks the dynamics is fully determined by the
eigenvalues of the connectivity matrix J. It is known that,in contrast to random matrices with Gaussian entries, aCauchy random matrix features an unbounded supportof the eigenvalues density, even in the limit of N → ∞[54–56]. Thus, we can conclude that, regardless of thevalue of the g, the dynamics of a Cauchy neural network isin this case divergent. For the same reason, any ϕðxÞ that islinear around x ≈ 0 and grows sufficiently slow for large x
leads to a self-sustained, active dynamics for any g [57].
However, in the biologically relevant regime neurons
exhibit saturation and thresholding at, respectively, large
and low values of total synaptic input. The corresponding
Cauchy network generically exhibits two phases: quiescentand active, and an associated transition between them [57].
In general the nature of the active phase will depend on thedetails of the activation function.
To further simplify the calculations and to simultaneously model edge of chaos and avalanches, in the following
we focus our attention on the binary activation functionϕðxÞ ¼ Θðx − θÞ, where Θ denotes the Heaviside functionand θ denotes the threshold. In this case the mean-field
equation (7) simplifies to1
mtþ1 ¼ arctan ðmt g=θÞ:πð8ÞThe stability of the trivial fixed point can be checked byexpanding the rhs of (8) around mt ¼ 0: mtþ1 ¼
ðg=πθÞmt þ Oðm3t Þ. The fixed point at mt ¼ 0, corresponding to the quiescent phase, is unstable for g > πθ. Since

arctanðxÞ=π is saturating and concave for all x > 0, anotherstable fixed point m close to 0 appears, through thesupercritical pitchfork bifurcation, exactly when
the trivial
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
pﬃﬃﬃ

−3=2
fixed point loses its stability [m ≈ 3ðg=θÞ
ðg=θÞ − πnear the transition point]. Due to the quenched, asymmetric
disorder of the connectivity matrix we can expect this fixed
point to represent a chaotic attractor of the network [61],with a large sensitivity to small perturbations. Our computer simulations confirm this prediction [57].
The transition from the quiescent to the chaotic phasecan be understood from the underlying structure of connections. Due to the power-law connectivity density, we
can expect that only a small fraction of the connections
contribute to the activity profile of the network. Indeed, as
we show in the following, the transition to chaos is driven
by the percolation transition of autocrat connections forwhich Jij > θ; i.e., an active presynaptic neuron will
activate the postsynaptic neuron in the absence of other
inputs. Around the critical point the mean activity of the
network is infinitesimal and thus the higher order interaction events (e.g., two neurons activating another neuron)are negligible. In other words, to a good approximation, aneuron can only be activated by another single neuron
through an autocrat connection, independently from otherneurons. This suggests that the transition to chaos in theneural network model is related to the critical branchingprocesses [64] (Fig. 1).
In the Cauchy case the probability that a given connection is an autocrat reads 
1gProbðJij > θÞ ¼ arctan
:
ð9Þ
π
Nθ
For a given neuron, the number of outgoing (or incoming)
autocrat connections is a binomial random variable with N
trials and the probability of success given by (9). In thelimit of N → ∞ it converges to the Poisson randomvariable with intensity
 
N
g
g
λ ¼ lim arctan
¼ :
ð10Þ
N→∞ π
Nθ
θπ
Now, let the initial state of the network be such that only asingle neuron (seed) is active. The number of active
neurons (descendants) in the next step is given by the
Poisson distribution and the mean number of active neurons
is given by (10). The theory of branching processes predicts
that the population will eventually die out almost surelyfor λ ≤ 1 and has a finite survival probability for λ > 1.
At λ ¼ 1 the process is critical and features scale-freeavalanches. The critical point predicted by the branchingprocess formulation of the network dynamics, g ¼ πθ, isthe same as the mean-field critical point predicted by (8).
The mapping to the branching process explains many
features of the Cauchy neural network around the critical

028101-3

PHYSICAL REVIEW LETTERS 125, 028101 (2020)

FIG. 2. Avalanche size and lifetime distributions in the networks of binary units. As expected from our theoretical predictions, at the transition point these distributions are describedby power laws, and the critical exponents match those of thecritical branching process.

point. Below the critical point the steady state is quiescent,and a bit-flip perturbation corresponds to a single neuron(seed) being activated. The local expansion rate ofsuch perturbation is given by λ. Above the critical point(g > πθ) each bit flip contributes in the same manner as asingle seed and, additionally, interacts with other active
neurons to activate and deactivate other descendants.
Thus, in the vicinity of the transition point λ gives alower bound on the local expansion rate of a perturbationin the steady state, and for λ > 1 the network is expected
to be chaotic in the thermodynamic limit [57]. Moreover,
the transition to chaos belongs to the mean-field directed
percolation universality class [65–67]. The propagationof the corresponding avalanches is characterized by [57]
power-law distributed sizes S: ProbðS > sÞ ∼ s−1=2 , andpower-law distributed lifetimes T: ProbðT > tÞ ∼ t−1 .
These theoretical predictions were corroborated by our
computer simulations of the Cauchy network, as shownin Fig. 2.

For a comparison, we have also studied Gaussian networks of threshold units with a fixed number of connections per neuron K [57,68]. While extremely sparsely
connected Gaussian networks (K ⪅ 12) behave qualitatively similar to the Cauchy network, the transition to chaos
becomes discontinuous in the biologically relevant regimeof K⪆13. With a biologically realistic K and finite N, thenetwork activity jumps between two metastable states near
the transition point, and cannot be robustly posed at theedge of chaos. The discontinuous property is due to theemergence of a metastable active state by the saddle-nodebifurcation as g=θ increases.
Importantly, although our theoretical predictions were
derived assuming simplistic threshold neural units, theytranslate directly to networks of more biologically plausible
leaky integrate-and-fire (LIF) neurons. The difference ofcontinuous and discontinuous transition is confirmed by thepresence or absence of a hysteresis loop in more realistic
networks of LIF neurons (Fig. 3). Hence, unlike theGaussian networks with realistic K, Cauchy networks
demonstrate critical phenomena and can reproduce experimentally observed scale-free avalanches at the criticalpoint. Moreover, a large Cauchy network can exhibit
arbitrarily low, self-sustained activity levels. In contrast,
the lowest possible activity level that can be achieved by theGaussian network with realistic K is about 11% in thebinary case and 40 Hz in the LIF case (Fig. 3).
For clarity we chose to limit our presentation to the Cauchy
distribution of Jij , but our results naturally extend to otherpower-law distributions. Indeed, let the synaptic efficacy
density asymptotically behave like a power law ρðJij Þ ∼Cαgα N −1 jJij j−1−α [69]. We then have ProbðJ ij > θÞ ¼
CN −1 ðg=θÞα , which holds for large enough N. The branching parameter is calculated as in (10) and reads λ ¼ Cðg=θÞα .
A continuous transition takes place at λ ¼ 1 and its features

FIG. 3. Continuous vs discontinuous transition in networks of leaky integrate-and-fire neurons. A slowly changing current was
injected and an average firing rate of the network was recorded as a function of time and the injected current amplitude. As predicted byour theory, a network with Gaussian weights exhibits a discontinuous transition between active and inactive states, which generates acharacteristic hysteresis loop. In contrast, the Cauchy network exhibits a continuous transition and thus shows no signs of the hysteresis
loop.

028101-4

PHYSICAL REVIEW LETTERS 125, 028101 (2020)
are, as before, described by the directed percolation universality class.
The connectivity in the current model is unstructured.
It would be interesting to combine power-law synaptic
weight distributions with structured networks, for example,
with hierarchical modules [70] or oscillations [3,7,71].
Incidentally, the same Cauchy distribution of membranepotential was found in the quadratic integrate and fire neuronmodel [72] due to single-neuron dynamics. Nontrivial effects
may arise from a combination of single-neuron and networkdriven heavy-tail statistics.
Power-law distributions of synaptic weights feature many
very weak synapses, that do not directly contribute to thecomputation. Even though this may seem wasteful, we thinkthat such architectures are not only biologically plausible[73] but may be beneficial. One possibility is that even weakconnections can activate a neuron once contextual input fromanother part of the brain increases the baseline membranepotential close to its spiking threshold. Such contextual input
can also raise the spiking probability of nearby neurons sothat synchronous activation of weak connections is more
likely. Weak synapses have also been reported to play a rolein unsupervised features extraction [74]. In the context ofreservoir computing [75,76] high computational capabilities
were achieved using nonbiologically sparse connectivity[77]. Our model provides a more biologically plausible
solution that the connectivity can be anatomically dense
but effectively sparse due to heavy-tailed synaptic weight
distribution. More generally, the optimal degree of sparsity
depends on the role of a given brain structure and the type ofthe employed plasticity [78]. Power-law distributed synaptic
weights may in this context provide a weakly informative[79,80] sparse connectivity prior, with weak and effectively
silent synapses providing a pool of potential connections thatcan be recruited when and if needed, as observed in the brain
during development [81,82].
Our results demonstrate that the shape of the synaptic
weight distribution can dramatically affect dynamics ofneural networks. A biological distribution of synaptic
weights can give distinct predictions from the frequently
assumed Gaussian distribution. The proposed mathematicalframework with power-law synaptic weights can easily be
adapted to other scenarios in future studies.
We thank Francesco Fumarola and anonymous referees
for invaluable discussions and comments on the manuscript.
Supported by RIKEN Center for Brain Science, Brain/
MINDS from AMED under Grant No. JP20dm020700,
and JSPS KAKENHI Grant No. JP18H05432.

*

nalewkoz@gmail.com
taro.toyoizumi@riken.jp[1] J. M. Beggs and D. Plenz, J. Neurosci. 23, 11167 (2003).
[2] W. L. Shew, W. P. Clawson, J. Pobst, Y. Karimipanah, N. C.
Wright, and R. Wessel, Nat. Phys. 11, 659 (2015).
†[3] A. J. Fontenele, N. A. P. de Vasconcelos, T. Feliciano,
L. A. A. Aguiar, C. Soares-Cunha, B. Coimbra, L. Dalla
Porta, S. Ribeiro, A. J. Rodrigues, N. Sousa et al., Phys.
Rev. Lett. 122, 208101 (2019).
[4] T. Petermann, T. C. Thiagarajan, M. A. Lebedev, M. A.
Nicolelis, D. R. Chialvo, and D. Plenz, Proc. Natl. Acad.
Sci. U.S.A. 106, 15921 (2009).
[5] C. Haldeman and J. M. Beggs, Phys. Rev. Lett. 94, 058101(2005).
[6] O. Kinouchi and M. Copelli, Nat. Phys. 2, 348 (2006).
[7] S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz, Proc.
Natl. Acad. Sci. U.S.A. 115, E1356 (2018).
[8] R. Stoop and F. Gomez, Phys. Rev. Lett. 117, 038102 (2016).
[9] M. A. Munoz, Rev. Mod. Phys. 90, 031001 (2018).
[10] C. G. Langton, Physica (Amsterdam) 42D, 12 (1990).
[11] N. Bertschinger and T. Natschläger, Neural Comput. 16,1413 (2004).
[12] R. Legenstein and W. Maass, in New Directions in Statistical Signal Processing: From Systems to Brain, edited by
S. Haykin, J. C. Principe, T. Sejnowski, and J. McWhirter(MIT Press, Cambridge, 2007), pp. 127–154.
[13] T. Toyoizumi and L. F. Abbott, Phys. Rev. E 84, 051908(2011).
[14] J. Schuecker, S. Goedeke, and M. Helias, Phys. Rev. X 8,041029 (2018).
[15] M. Benayoun, J. D. Cowan, W. van Drongelen, and E.
Wallace, PLoS Comput. Biol. 6, e1000846 (2010).
[16] J. M. Beggs and N. Timme, Front. Physiol. 3, 163 (2012).
[17] H. Sompolinsky, A. Crisanti, and H.-J. Sommers, Phys. Rev.
Lett. 61, 259 (1988).
[18] L. Molgedey, J. Schuchhardt, and H. G. Schuster, Phys. Rev.
Lett. 69, 3717 (1992).
[19] D. J. Amit and N. Brunel, Cereb. Cortex 7, 237 (1997).
[20] N. Brunel, J. Comput. Neurosci. 8, 183 (2000).
[21] K. Rajan, L. F. Abbott, and H. Sompolinsky, Phys. Rev. E
82, 011903 (2010).
[22] M. Massar and S. Massar, Phys. Rev. E 87, 042809 (2013).
[23] M. Stern, H. Sompolinsky, and L. F. Abbott, Phys. Rev. E
90, 062710 (2014).
[24] J. Aljadeff, M. Stern, and T. Sharpee, Phys. Rev. Lett. 114,088101 (2015).
[25] J. Kadmon and H. Sompolinsky, Phys. Rev. X 5, 041030(2015).
[26] G. Wainrib and M. N. Galtier, Neural Netw. 76, 39 (2016).
[27] A. Crisanti and H. Sompolinsky, Phys. Rev. E 98, 062120(2018).
[28] O. Harish and D. Hansel, PLoS Comput. Biol. 11, e1004266(2015).
[29] P. Dayan and L. F. Abbott, Theoretical Neuroscience (TheMIT Press, Cambridge, 2001).
[30] M. Häusser and P. Monsivais, Neuron 40, 449 (2003).
[31] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski,
Neuronal Dynamics: From Single Neurons to Networks andModels of Cognition (Cambridge University Press, Cambridge, England, 2014).
[32] D. Millman, S. Mihalas, A. Kirkwood, and E. Niebur, Nat.
Phys. 6, 801 (2010).
[33] P. Dayan and L. F. Abbott, Theoretical Neuroscience:
Computational and Mathematical Modeling of Neural
Systems (The MIT Press, 2005).

028101-5

PHYSICAL REVIEW LETTERS 125, 028101 (2020)[34] J. H. Marshel, Y. S. Kim, T. A. Machado, S. Quirin, B. Benson,
J. Kadmon, C. Raja, A. Chibukhchyan, C. Ramakrishnan,
M. Inoue et al., Science 365, eaaw5202 (2019).
[35] S. Scarpetta, I. Apicella, L. Minati, and A. de Candia, Phys.
Rev. E 97, 062305 (2018).
[36] R. Sayer, M. Friedlander, and S. Redman, J. Neurosci. 10,826 (1990).
[37] D. Feldmeyer, V. Egger, J. Lübke, and B. Sakmann, J.
Physiol. 521, 169 (1999).
[38] S. Song, P. J. Sjöström, M. Reigl, S. Nelson, and D. B.
Chklovskii, PLoS Biol. 3, e68 (2005).
[39] S. Lefort, C. Tomm, J.-C. F. Sarria, and C. C. Petersen,
Neuron 61, 301 (2009).
[40] Y. Ikegaya, T. Sasaki, D. Ishikawa, N. Honma, K. Tao, N.
Takahashi, G. Minamisawa, S. Ujita, and N. Matsuki, Cereb.
Cortex 23, 293 (2013).
[41] Y. Loewenstein, A. Kuras, and S. Rumpel, J. Neurosci. 31,
9481 (2011).
[42] M. Gilson and T. Fukai, PLoS One 6, e25339 (2011).
[43] P. Zheng, C. Dimitrakakis, and J. Triesch, PLoS Comput.
Biol. 9, e1002848 (2013).
[44] N. Yasumatsu, M. Matsuzaki, T. Miyazaki, J. Noguchi, and
H. Kasai, J. Neurosci. 28, 13592 (2008).
[45] A. Nagaoka, H. Takehara, A. Hayashi-Takagi, J. Noguchi,
K. Ishii, F. Shirai, S. Yagishita, T. Akagi, T. Ichiki, and H.
Kasai, Sci. Rep. 6, 26651 (2016).
[46] J. Humble, K. Hiratsuka, H. Kasai, and T. Toyoizumi, Front.
Comput. Neurosci. 13, 38 (2019).
[47] K. Ishii, A. Nagaoka, Y. Kishida, H. Okazaki, S. Yagishita,
H. Ucar, N. Takahashi, N. Saito, and H. Kasai, eNeuro 5,0282 (2018).
[48] H. Okazaki, A. Hayashi-Takagi, A. Nagaoka, M. Negishi,
H. Ucar, S. Yagishita, K. Ishii, T. Toyoizumi, K. Fox, and H.
Kasai, Neurosci. Lett. 671, 99 (2018).
[49] J.-n. Teramae, Y. Tsubo, and T. Fukai, Sci. Rep. 2, 485 (2012).
[50] G. Buzsáki and K. Mizuseki, Nat. Rev. Neurosci. 15, 264(2014).
[51] J.-n. Teramae and T. Fukai, Proc. IEEE 102, 500 (2014).
[52] Y. Omura, M. M. Carvalho, K. Inokuchi, and T. Fukai, J.
Neurosci. 35, 14585 (2015).
[53] W. Feller, An Introduction to Probability and Its Applications (Wiley, New York, 1971), Vol. II.
[54] P. Cizeau and J.-P. Bouchaud, Phys. Rev. E 50, 1810 (1994).
[55] Z. Burda, R. A. Janik, J. Jurkiewicz, M. A. Nowak, G. Papp,
and I. Zahed, Phys. Rev. E 65, 021106 (2002).
[56] E. Gudowska-Nowak, M. A. Nowak, D. R. Chialvo, J. K.
Ochab, and W. Tarnowski, Neural Comput. 32, 395 (2020).
[57] See Supplemental Material at http://link.aps.org/
supplemental/10.1103/PhysRevLett.125.028101 for details,
which includes Refs. [58–60].
[58] E. Goles, Discrete Appl. Math. 13, 97 (1986).
[59] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah,
SIAM Rev. 59, 65 (2017).

[60] C. Linssen et al., Nest 2.16.0, 2018, https://doi.org/10.5281/
zenodo.1400175.
[61] In the binary case the system has a finite number of states forany finite N and the attractor has to be periodic. Thus,
formally, irregular aperiodic behavior can only be observedin the limit of N → ∞. Chaoslike signatures can nonethelessbe observed for finite N; e.g., the typical lengths oftransients and cycles rapidly change around the predictedtransition point, and grow exponentially with N in the“chaotic” phase [62,63].
[62] C. van Vreeswijk and H. Sompolinsky, Neural Comput. 10,1321 (1998).
[63] B. Luque and R. V. Solé, Physica (Amsterdam) 284A, 33(2000).
[64] T. E. Harris, The Theory of Branching Processes (Springer,
Berlin, 1963).
[65] P. Alstrøm, Phys. Rev. A 38, 4905 (1988).
[66] M. A. Munoz, R. Dickman, A. Vespignani, and S. Zapperi,
Phys. Rev. E 59, 6175 (1999).
[67] G. Ódor, Rev. Mod. Phys. 76, 663 (2004).
[68] B. Derrida, E. Gardner, and A. Zippelius, Europhys. Lett. 4,167 (1987).
[69] Note that we scale Jij with the number of neurons as N −1=α ,which assures the existence of a nontrivial limit N → ∞.
For α > 2 another choice is to scale the synaptic strengths
as N −1=2 , which in the limit of N → ∞ corresponds to theGaussian network.
[70] E. J. Friedman and A. S. Landsberg, Chaos 23, 013135(2013).
[71] S.-S. Poil, R. Hardstone, H. D. Mansvelder, and K.
Linkenkaer-Hansen, J. Neurosci. 32, 9817 (2012).
[72] E. Montbrió, D. Pazó, and A. Roxin, Phys. Rev. X 5,021028 (2015).
[73] L. Cossell, M. F. Iacaruso, D. R. Muir, R. Houlton, E. N.
Sader, H. Ko, S. B. Hofer, and T. D. Mrsic-Flogel, Nature(London) 518, 399 (2015).
[74] H. Huang, J. Phys. A 51, 08LT01 (2018).
[75] H. Jaeger, German National Research Center for Information Technology, Bonn, Germany, GMD Technical Report
No. 148, p. 13 2001.
[76] W. Maass, T. Natschläger, and H. Markram, Neural Comput.
14, 2531 (2002).
[77] L. Büsing, B. Schrauwen, and R. Legenstein, NeuralComput. 22, 1272 (2010).
[78] A. Litwin-Kumar, K. D. Harris, R. Axel, H. Sompolinsky,
and L. Abbott, Neuron 93, 1153 (2017).
[79] A. Gelman, A. Jakulin, M. G. Pittau, Y.-S. Su, Ann. Appl.
Stat. 2, 1360 (2008).
[80] S. Van Dongen, J. Theor. Biol. 242, 90 (2006).
[81] D. Liao, N. A. Hessler, and R. Malinow, Nature (London)
375, 400 (1995).
[82] G. A. Kerchner and R. A. Nicoll, Nat. Rev. Neurosci. 9, 813(2008).

028101-6