Supplementary Material: Edge of chaos and avalanches in neural networks with heavy-tailed
synaptic weight distributionŁukasz Kuśmierz,1, ∗ Shun Ogawa,1 and Taro Toyoizumi1, 2
1

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, JapanCONTENTS

I. Binary Gaussian network
A. Fully connected network
B. Sparse network

1
1
2

II. Mean field equations correctly predict the mean activity in the thermodynamic limit

3III. Active phase is chaotic in random networks of binary neurons

3

IV. Details of the simulations

5

V. Avalanche statistics
A. Avalanche size distribution
B. Avalanche lifetime distribution

5
5
6

VI. Existence of a continuous transition in Cauchy networks with a positive threshold

6

References

7
I.

BINARY GAUSSIAN NETWORKA.

Fully connected network

For a general activation function it is convenient to describe the behavior of the Gaussian networks in terms of another orderparameter,
1 X
q0 (t) =φ(xi (t))2 .
(S.1)
N i
However, in the binary case analyzed here it is equivalent to m(t). The corresponding dynamical mean-field equation [1] reads"!#1θm(t + 1) =1 − erf √,(S.2)
22m(t)gRx
where the error function is given by erf(x) = 2π−1/2 0 exp(−z2 )dz. Expanding the RHS of (S.2) around m(t) ≈ 0 givesr!
!!
m(t)g2m(t)g2θ2m(t + 1) =
exp −
1+O
.
(S.3)
2πθ2
2m(t)g2
θ2
Due to the exponential factor m(t + 1) < m(t) for small enough m(t), which proves that the fixed point m(t) = 0 is locally stable
under the evolution (S.2). Since the quiescent state is always stable, the transition to chaos, if present, must be discontinuous.
This is confirmed by the graphical inspection of (S.2) (Fig. S.1) and in computer simulations (Fig. S.2).

∗ nalewkoz@gmail.com

2

0.3
1.0

g = 2.0

g = 2.5

g = 3.2

g = 4.00.2
0.8m(t + 1)

0.1
0.6
0.0
0.3
0.4
0.2
Dense Gauss
Dense Cauchy
Sparse K=6
Sparse K=200

0.2
0.1

0.0
0.00
0.0

0.05

0.10 0.20.15

0.20

0.25
0.4 0.30
0.00 0.6
0.05
m(t)

0.10

0.150.8 0.20

0.25

0.30
1.0

FIG. S.1. Mean field m(t) 7→ m(t + 1) mapping for fully connected (dense) Gaussian and Cauchy networks and sparse Gaussian network witha fixed number of incoming connections per neuron K. For large values of K, the MF equation looks similar to the fully connected Gaussiancase. In particular, we observe a discontinuous transition at g ≈ 2.5. (note that this similarity only holds if g is not too large, since for any
finite K there is a continuous transition at which the trivial fixed point loses its stability). For an intermediate sparsity (2 < K / 12), only asecond order transition is observed and the dynamics looks qualitatively similar to the dense Cauchy case. For K ≤ 2 no transition to chaos isobserved and the trivial fixed point is always stable (results not shown).

The mapping to a branching process offers a simple way of understanding this result. The probability that a given synapse isan autocrat reads

 r


1 
N θ 

Prob(Ji j > θ) = 1 − erf −
(S.4)

2
2 g 
In the thermodynamic limit the average number of autocrat connections per neuron can be calculated as before as
λ = limN→∞ NP(Ji j > θ), which in our case leads to
r
!

Nθ2 
Ng
exp − 2 1 + O(1/N 2 ) = 0,
(S.5)
λ = limN→∞2θ
2gi.e. an activity starting from a single seed almost surely dies out.
Interestingly, the discontinuous character of the transition between quiescent and active phases have already been reported invarious models that employ Gaussian weights [2, 3], which confirms the generic nature of this result.
B.

Sparse network

Let each neuron receive exactly K incoming connections, randomly chosen from the network, and let Ji j ∼ N(0, g2 /K). If themean activity of the network at time t is m(t), the probability that exactly n incoming neurons are active is given by the binomial
distribution!
K
P(n; m(t)) =m(t)n (1 − m(t))K−n .
(S.6)
nThe membrane potential of a neuron x, conditioned on n incoming connections being active, is a normal random variable with
µ = 0 and σ2 = ng2 /K. The mean activity of the network in the next step is equal to the probability that the membrane potentialof any given neuron crosses the threshold, and so it reads √ 
!
KX θ K 
K
−1
n
K−n 
m(t + 1) = 2
m(t) (1 − m(t))(S.7)
1 − erf  √ .
n
2ng
n=1

3
Cauchy dense

Gauss dense

0.3

0.3

0.2

0.2
Theory

0.1

0.1

0.0

0.0
2.0

2.5

3.0

3.5

4.0

4.5

5.0

1.5

2.0

2.5

Gauss K=7

3.0

3.5

4.0

4.5

3.5

4.0

4.5

Gauss K=20

0.3

0.3

0.2

0.2

0.1

0.1

0.0

0.0
1.5

2.0

2.5

3.0

3.5

4.0

4.5

1.5

2.0

2.5

3.0

FIG. S.2. Steady-state mean activity as a function of the control parameter g/θ. Lines were obtained by solving our mean field equations(8), (S.2), and (S.7) self-consistently. The points were obtained from computer simulations of M = 105 /N independent realizations of J. For
each realization of the weight matrix, the network was first evolved for 400 steps from an active state (m(0) ≈ 0.5). The following 200 steps
were used in calculating the average, which was performed over the steps, neurons, and realizations of J. Error bars denote ±3σ confidence
intervals.

Our computer simulations corroborate the validity of (S.7), see Fig. S.2. Note, additionally, that it is not difficult to show that
(S.7) simplifies to the dense case (S.2) in the limit of K → ∞. Analogous calculations in the case of stochastic units were
presented in [4], where the range of validity of the annealed approximation is also discussed.

II.

MEAN FIELD EQUATIONS CORRECTLY PREDICT THE MEAN ACTIVITY IN THE THERMODYNAMIC LIMIT

We have performed computer simulations of Cauchy, fully connected Gaussian, and sparse Gaussian networks with binary
activation functions (Fig. S.2). The good match between our theoretical predictions and the simulation results strongly suggeststhat the mean field approach is exact in the thermodynamic limit.

III.

ACTIVE PHASE IS CHAOTIC IN RANDOM NETWORKS OF BINARY NEURONS

Here we show, with a simple mathematical argument and computer simulations, that whenever the mean network activity m is
non-zero in the steady-state, the corresponding attractor is chaotic. This is true for both Gaussian and Cauchy networks with thebinary activation function, and is expected to hold in more general settings, unless some non-trivial structure of the connectivity
is present, e.g. symmetric or antisymmetric connectivity (the latter with a symmetric version of the binary activation function,i.e. with θ = 0 and taking values ±1) leads to non-chaotic, periodic attractors [5].
It is easy to prove our result in the Gaussian case. In this case, with N → ∞, the maximal Lyapunov exponent is given by[1, 6, 7]1
gauss
λmax = ln g2

Z∞

2


√ 2
Dx φ0 g q0 x

(S.8)

−∞

where Dx denotes the standard Gaussian measure, q0 is the average squared activity (S.1). It is clear that whenever φ features√gauss
any discontinuity and the phase is active (g q0 , 0), formula (S.8) predicts λmax = ∞. In particular, the binary activation
function has a discontinuity at x = θ and thus the corresponding active phase is always chaotic. The analysis for the binary
Cauchy network is more involved and will be published elsewhere. However, the intuition is the same: small perturbations in
the active phase are expanded due to the discontinuity of the activation function. Note that these results are in line with previous

4Cauchy (replica 1)

C

100

E

80
Neuron #Neuron #80

Activity difference100

ln d(t) (log-distance)

A

60
40
20

60
40

θ = 0.5
θ = 0.9
θ = 0.95
θ = 1.05
θ = 1.1
θ = 1.5

6

4

2

20
0

B

100
time

150

200

Cauchy (replica 2)

50

100
time

80
m(t)

0.15

0.13

20

10

15

20

λ (unstable)λ (stable)

4
Cauchy (replica 1)
Cauchy (replica 2)
Perturbation time

0.14

40

5

t (from the perturbation time)

F
0.16

60

0

200

D

100

0.17

Neuron #150

Average expansion rate

50

Quiescent state
Stable steady state

3

2

1
0.12
50

100
time

150

200

0

50

100
t

150

200

0.6

0.8

1.0

1.2

1.4

θ

FIG. S.3. Numerical perturbation analysis of the binary Cauchy network. (A-D): Two replicas were evolved from identical initial conditions(θ = 0.95). A single-neuron perturbation, introduced at T 0 = 100, is quickly expanded, which can be seen by comparing individual neuron
activities (C) or mean network activities (D). (E): The average Hamming distances between perturbed and unperturbed states as functions of
the number of steps T after introducing a single-neuron state flip to a steady-state trajectory. The network is sensitive to the perturbations forθ < 1, corroborating the existence of a transition to chaos at the predicted value of θ = 1. (F): The average expansion rate as a function of θ. As
expected, the average expansion rate from the quiescent initial state is equal to λ = g/(πθ). The average expansion rate from the steady state isdifferent from λ for θ < 1; in this regime the quiescent state does not correspond to the stable steady state (the quiescent state is unstable). In(E) and (F) the distances are averaged over all possible neurons and over 10 realizations of J. See text for more details.

studies, e.g., it has been shown that in random Boolean networks Lyapunov exponents scale as ln K with the number of incoming
connections [8]. The issue of infinite Lyapunov exponents in binary networks is also discussed in detail in [9].
To support our claim we simulated Cauchy networks with binary activation functions. Two replicas that shared J were evolved
from identical initial conditions. At a fixed time T 0 = 100 a perturbation was introduced by flipping the state of a single neuron
in one of the replicas. We then calculate the activity differenceδi (t) = φ xi1 (t) − φ xi2 (t) ,(S.9)
the Hamming distanced(t) =NX

δi (t),(S.10)d(T 0 + 1),
d(T 0 )

(S.11)

i=1

and the average expansion rater=

where d(T 0 ) = 1. Similarly, we calculated the average expansion rate from the quiescent initial state λ. In this case we evolved
the (unperturbed) network from a quiescent initial state, which corresponds to the stable steady state for θ > θ∗ = 1 (g = π),
whereas for θ < 1 is unstable. The theoretically calculated branching parameter agrees with λ obtained from the simulations(Fig. S.3F). The evolution of the network is sensitive to the perturbation for θ < 1, which can be observed at the level of individual
neuron’s activities (Fig. S.3A-C) and, due to the finite size of the network, the mean activity (Fig. S.3D). In the chaotic phase,the Hamming distance initially grows exponentially and then saturates at d(∞) ∼ N. The chaotic phase is signaled by r > 1,
which is directly related to the maximal Lyapunov exponent being positive [8–10]. Although r in general is not equal to λ, bothquantities cross the value 1 at the critical point θ = 1. The discrepancy observed around the transition point (the slowly growing
distance at θ = 1.05 in Fig. S.3E and the corresponding r ≈ 1.04 > 1 (Fig. S.3F)) is expected to be a finite size effect.

5
IV.

DETAILS OF THE SIMULATIONS

Computer simulations of binary networks were performed using custom-written codes in MATLAB and Julia [11]. The results
for spiking neural networks were obtained with NEST Simulator [12] through a custom-written code in Python. Details of thesimulations:Fig. 2: Networks had size N = 104 and results were averaged over 10 independent realizations of J. For a given realizationof J, we run 104 simulations, each starting from different seed neuron (i.e. one neuron active, all other neurons inactive). g = π.
Fig. 3: Both Gaussian and Cauchy networks were fully connected with N = 104 . The injected current changed between −400
pA and 400 pA in small increments every 5 ms. All neurons were of type iaf psc alpha, which denotes a leaky integrate-andfire neuron with alpha-shaped postsynaptic currents. Default parameters of the model neuron were used, i.e.: resting potential
E L = −70 mV, capacity of the membrane Cm = 250 pF, membrane time constant τm = 10 ms, refractory period tre f = 2 ms,
spike threshold Vth = −55 mV, reset potential Vreset = −70 mV, rise time of the excitatory and inhibitory synaptic alpha function
τ syn = 2 ms. Static synapses were used with the default delay of 1 ms, and weights were randomly drawn from symmetric (a)√Gaussian distribution with σ = 2.4 × 103 / N pA, and (b) Cauchy distribution with γ = 1.92 × 103 /N pA. Poisson noise wasinjected randomly into the network, activating each neuron approximately twice every second. Note that without any external
input the analyzed networks cannot recover from the quiescent state, since the model neurons are never spontaneously active.
Sub-sampling: For the sake of clarity and drawing efficiency, in the raster plots only 10% of spikes of 100 randomly chosenneurons were drawn. Activity histograms were created using all data.
Fig. S.3: The averaging was performed over all N = 104 neurons (i.e., including the unperturbed network, N + 1 replicas were
simulated) and over 10 realizations of J. At time T = 0 the unperturbed network was prepared in the (a) steady state by evolving
it for 100 steps before introducing perturbations or (b) quiescent state. g = π. Sub-sampling: For the sake of clarity and drawing
efficiency, activities of a randomly chosen subpopulation of neurons (N 0 = 100) is shown.

V.

AVALANCHE STATISTICS

The mapping to the branching process together with the general results known for critical branching processes provide thecritical exponents [13–15]. Here, for completeness, we calculate two critical exponents in our specific case.

A.

Avalanche size distribution

The size of an avalanche is defined as the sum of the number of active neurons at each time step, from the beginning of
the avalanche till its end. Let S m denote the size of an avalanche starting from m seeds and Gm (z) denote the corresponding
generating function
Gm (z) ≡ hzS m i =∞Xz j Prob(S m = j).

(S.12)

j=0

Since the activity of network is assumed to be sparse, the avalanche that starts fromn (i)mo seeds consists of m independent avalanchesP(i)starting from a single seed. Therefore we can write that S m = m
i=1 S 1 , where S 1 is a set of i.i.d. random variables denoting
sizes of single-seed generated avalanches. At the level of the generating functions, this assumption leads to a simple expression
Gm (z) = [G1 (z)]m .

(S.13)On the other hand, we know how a seed neuron propagates the activity through the network in a single step: it activates mneurons with probability
pm =

λm −λ
e ,
m!

(S.14)

where λ is given by (10). This means that one seed generate an avalanche of size 1 + S m with probability pm , where 1 is fromthe first step (i.e., the seed) and S m is from the subsequent steps. Hence, the single-seed generating function can be calculated
as follows
G1 (z) =∞X
m=0

pm hz1+S m i = z

∞
X
m=0

pmGm (z).

(S.15)

6
We combine (S.15) with (S.13) and (S.14) and arrive at an implicit expression for the one-seed avalanche size generating function


∞




 S0 Xm
= z exp(λG1 (z) − λ).
(S.16)
p0 e +
pm [G1 (z)] 
G1 (z) = z 




m=1

Note that we used S 0 = 0 (no avalanche without a seed). In order to inspect the tail of the distribution of S 1 we introduce an
auxiliary function g() = 1 − G1 (1 − ) and expand the RHS of (S.16) assuming that   1 and g()  1 (valid for λ ≤ 1):!
λ2 g()2
1 − g() = (1 − ) 1 − λg() ++ O(g()3 ) ,(S.17)
2
which to the lowest order can be rewritten as, for λ , 1 1−λ√g() =  2, for λ = 1.

(S.18)The small  behavior of g() ∼  1/2 in the λ = 1 case translates into the tail behavior of the avalanche size density as
Prob(S 1 = s) ∼ s−3/2(S.19)

for large s.
B.

Avalanche lifetime distribution

Let T m be the lifetime (number of steps with nonzero activity) of an avalanche that starts from m seeds. By definition T 0 = 0
and T m ≥ 1 for m > 0. As before, we treat an avalanche from different seeds as independent, and thus the following identity
linking the survival probabilities holds
Qm (t) = 1 − Prob(T m ≤ t) = 1 − [Prob(T 1 ≤ t)]m = 1 − [1 − Q1 (t)]m .

(S.20)As in the case of the size distribution, we can unwrap the first step of the dynamics starting from a single seed, which gives
Q1 (t + 1) =∞X

pm Qm (t) = 1 −∞X

pm [1 − Qm (t)] .

(S.21)

m=0

m=1

We plug (S.20) and (S.14) into (S.21) and arrive at the following recursive relation
Q1 (t + 1) = 1 − exp [−λQ1 (t)] .

(S.22)

If λ > 1 there exists a non-zero fixed point corresponding to the non-zero probability of survival at t → ∞. In contrast, theactivity eventually dies out almost surely for λ ≤ 1. Assuming λQ1 (t)  1 the recursive relation simplifies to
Q1 (t + 1) = λQ1 (t) −λ2
Q1 (t)2 ,2(S.23)

which predicts an exponential decay for λ < 1. At the critical point λ = 1 and the decay is a power law. In that case the recursion
can be solved with an ansatz Q1 (t) = C/tδ , leading to δ = 1, as expected.
VI.

EXISTENCE OF A CONTINUOUS TRANSITION IN CAUCHY NETWORKS WITH A POSITIVE THRESHOLD

Let φ(x) be an activation function such that φ(x) = 0 for x below a positive threshold θ, and φ(x) ≈ C for sufficiently large x
Rb(i.e. x > m1 ). Without much loss of generality, we additionally assume that φ(x) ≥ 0 and a φ(x)dx < ∞. The integral in themean-field equation can be then decomposed as follows:1m(t + 1) =
π

Z∞
0

m1

φ(m(t)gz) 1dz≈
π
1 + z2Zm(t)gθm(t)gφ(m(t)gz) Cdz+π
1 + z2

Z∞
m1
m(t)g

1
m(t)g
dz
=
π
1 + z2

Zm1
θ

!
φ(y)
C
m(t)g
dy 
+ arctan. (S.24)
m1
m(t)g 2 + y2 π

7We expand (S.24) around m(t) = 0:
 m1

Z

φ(y)
C 
m(t)g  dy 2 + + O(m(t)3 ),m(t + 1) =
π 
m1 
y

(S.25)

θand conclude that the transition between quiescent and active state occurs at the critical point described by the equation
 m1

Z

g C 
φ(y)(S.26)
 dy 2 + = 1.
π
m1 
y
θ

Hence, this guarantees the existence of positive and finite critical g∗ that solves the above equation. It is possible to extend
these results to activation functions φ(x) that are non-zero around x = 0 and are non-saturating: the transition exists if φ(x) is
sufficiently superlinear around x = 0 and sublinear for large |x|.

[1] L. Molgedey, J. Schuchhardt, and H. G. Schuster, Physical review letters 69, 3717 (1992).
[2] N. Brunel, Journal of computational neuroscience 8, 183 (2000).
[3] 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).
[4] B. Derrida, E. Gardner, and A. Zippelius, EPL (Europhysics Letters) 4, 167 (1987).
[5] E. Goles, Discrete Applied Mathematics 13, 97 (1986).
[6] M. Massar and S. Massar, Physical Review E 87, 042809 (2013).
[7] G. Wainrib and M. N. Galtier, Neural Networks 76, 39 (2016).
[8] B. Luque and R. V. Solé, Physica A: Statistical Mechanics and its Applications 284, 33 (2000).
[9] C. v. Vreeswijk and H. Sompolinsky, Neural computation 10, 1321 (1998).
[10] L. Büsing, B. Schrauwen, and R. Legenstein, Neural computation 22, 1272 (2010).
[11] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM review 59, 65 (2017).
[12] C. Linssen, M. E. Lepperød, J. Mitchell, J. Pronold, J. M. Eppler, C. Keup, A. Peyser, S. Kunkel, P. Weidel, Y. Nodem, D. Terhorst,
R. Deepu, M. Deger, J. Hahne, A. Sinha, A. Antonietti, M. Schmidt, L. Paz, J. Garrido, T. Ippen, L. Riquelme, A. Serenko, T. Kühn,
I. Kitayama, H. Mørk, S. Spreizer, J. Jordan, J. Krishnan, M. Senden, E. Hagen, A. Shusharin, S. B. Vennemo, D. Rodarie, A. Morrison,
S. Graber, J. Schuecker, S. Diaz, B. Zajzon, and H. E. Plesser, “Nest 2.16.0,” (2018).
[13] P. Alstrøm, Physical Review A 38, 4905 (1988).
[14] M. A. Munoz, R. Dickman, A. Vespignani, and S. Zapperi, Physical Review E 59, 6175 (1999).
[15] G. Ódor, Reviews of modern physics 76, 663 (2004).