PHYSICAL REVIEW E 91, 032802 (2015)

Structure of attractors in randomly connected networks
Taro Toyoizumi*
RIKEN Brain Science Institute, Wako-shi, Saitama 351-0198, Japanand Department of Computational Intelligence and Systems Science, Tokyo Institute of Technology, Yokohama 226-8502, Japan

Haiping Huang†
RIKEN Brain Science Institute, Wako-shi, Saitama 351-0198, Japan(Received 14 November 2014; revised manuscript received 3 February 2015; published 6 March 2015)
The deterministic dynamics of randomly connected neural networks are studied, where a state of binary neurons
evolves according to a discrete-time synchronous update rule. We give theoretical support that the overlap ofsystems’ states between the current and a previous time develops in time according to a Markovian stochastic
process in large networks. This Markovian process predicts how often a network revisits one of the previously
visited states, depending on the system size. The state concentration probability, i.e., the probability that two
distinct states coevolve to the same state, is utilized to analytically derive various characteristics that quantify
attractors’ structure. The analytical predictions about the total number of attractors, the typical cycle length, andthe number of states belonging to all attractive cycles match well with numerical simulations for relatively large
system sizes.
DOI: 10.1103/PhysRevE.91.032802

PACS number(s): 84.35.+i, 02.50.Ey, 05.45.−a

I. INTRODUCTION

Neurons in the brain interact with each other in a heterogeneous and asymmetric way [1], producing complex neuronal
dynamics for information processing. In past decades, there
was a surge of research interests in randomly connectedneural networks [2–7]. Although their behavior is describedby simple deterministic equations, the resulting dynamics arerich, exhibiting fixed-point behavior, limit cycles, or highdimensional chaos. These networks are capable of generating
useful dynamic activity patterns after appropriate learning[8,9].
Simple models of neural networks [10–18] have been
explored to elucidate characteristics of their complex dynamics. In these networks, connections between binary neurons
are independently drawn from an identical distribution, andthe state of a network is updated simultaneously in discrete time steps without thermal noise. Thus, every initialconfiguration must evolve into an attractor, which is either
a fixed point or a limit cycle. Because a fixed point is alimit cycle of length l = 1, the whole state space is divided
into separated basins of attractions with heterogeneous cycle
lengths. Extensive numerical simulations were carried out toanalyze the typical cycle length and the number of cycles[12,19]. The typical length of the cycles was observed to
grow exponentially with the number of neurons n (such kindsof cycles are called chaotic attractors), and the total numberof attractors increases linearly with n. These quantities werealso analytically evaluated based on an empirical assumptionthat the dynamics loses memory of its nonimmediate past
[16].
In this work, we develop a dynamic mean-field theory tocharacterize the attractors of the asymmetric neural networkby extending the state concentration concept [17], recently*†

taro.toyoizumi@brain.riken.jp
physhuang@gmail.com

1539-3755/2015/91(3)/032802(10)introduced to characterize the robustness and quickness of a
network’s transient dynamics. Our analysis estimates the (cumulative) distribution for the cycle length of attractors, the total
number of attractors, and the volume of attractors in the state
space.
We remark that our work has threefold contributions
for understanding the statistical properties of the dynamics
of randomly connected neural networks. First, theoreticalsupport for the Markovian property of state concentration
dynamics (termed the annealed approximation in Ref. [16]) is
provided by computing the finite-size effect of the mean-field
theory by explicitly evaluating the quenched randomnessof network connections. Second, we provide a detailed
picture about how state concentration happens in randomly
connected neural networks. In particular, we quantify whatis the characteristic distance that typically leads to state
concentration and evaluate characteristic time scales underlying the state concentration dynamics. Finally, our theory
gives good consistency with numerical simulations on the
distribution of the cycle length, the typical cycle length, thenumber of cycles, and the total number of states belonging
to all attractive cycles. These three contributions complement
the previous studies [12,14,16,17,19] and provide deep insights towards the dynamics of randomly connected neuralnetworks.
The paper is organized as follows. In Sec. II, we define theneural network model and its dynamics. Mean-field analysis
is presented in detail in Sec. III. Results on the state concentration and statistical properties of attractors are discussed inSecs. IV and V, respectively. We summarize our results in
Sec. VI.
II. MODEL DEFINITION

We consider randomly connected neural networks consisting of n neurons (units). Each unit interacts with all
the other units with an asymmetric coupling. We use Jij torepresent the coupling strength from unit j to i, and Jij is032802-1©2015 American Physical Society

TARO TOYOIZUMI AND HAIPING HUANGPHYSICAL REVIEW E 91, 032802 (2015)

independent of Jj i (and others), and they follow the same
Gaussian distribution with zero mean and variance 1/n. The
state of neuron i (i = 1, . . . ,n) at time t + 1 (t = 0,1, . . . ) is
set according to the parallel deterministic dynamics in discrete
time steps by its input hi (t) as+1 (active state)(1)
σi (t + 1) = sgn [hi (t)] =−1 (silent state),

across different realizations of {Jij } contains information
about the structure of attractors. Let us define the overlap of two states, σ (t) = {σi (t)|i = 1, . . . ,n} and σ (s) ={σi (s)|i = 1, . . . ,n} along the same trajectory at different times
t > s, by1
qts ≡σi (t)σi (s).
n i=1nwhere the input is defined byhi (t) =nJij σj (t).

(2)

j =1

Therefore, by combining Eqs. (1)and (2), the dynamics are
summarized by σi (t + 1) = sgn[ j Jij σj (t)] in terms of theactivity, or equivalently by hi (t + 1) = j Jij sgn[hj (t)] interms of the input.
We later compare the dynamics of randomly connectedneural networks to that of random Boolean networks [20,21],
where each one of 2n states σ (t) = {σi (t)|i = 1, . . . ,n} israndomly mapped to another.
III. MEAN-FIELD ANALYSIS

We study the dynamical evolution of the overlap between
two states along a trajectory, expecting that its distribution

This overlap takes +1 if two states are the same and −1 if
one is the sign flip of the other. The overlap takes a discretevalue for a finite-size network, but can be approximatedas a continuous quantity in the large network size limit.
The mean-field theory provides the dynamics of this overlap
parameter and its fluctuation defined over the ensemble ofrandom {Jij } (see Appendix A). The stochastic dynamics ofthe overlap is well approximated for large n by a Markovian
process
Pt+1,s+1 (q) ≈W (q|q  )Pts (q  )dq  ,(4)where Pts (q) ≡ Prob(qts = q) is the probability of qts = q.
The transition probability is approximated for large but finite
n by a simple binomial distribution

n
1 + ϕ(q  ) n(1+q)/2 1 − ϕ(q  ) n(1−q)/2
W (q|q ) =22n(1 + q)/2
 
1 + q 1 + ϕ(q  ) 1 − q 1 − ϕ(q  )
ln
+
ln
,
≈ exp n H (q) +
2
2
2
2


(3)



where ϕ(q) ≡ (2/π ) arcsin q and H (q) ≡ − 1+q
ln 1+q
−
2
2
1−q
1−q
ln
.
Note
that
Eq.
(5)
summarizes
theprobability
22that n(1 + q)/2 out of n neurons take the same sign instates σ (t + 1) and σ (s + 1), given that n(1 + q  )/2 out ofn neurons take the same sign in the previous steps. The
binomial distribution in Eq. (5) suggests that the state overlap
for each neuron is approximately independent, occurring withprobability [1 + ϕ(q  )]/2 (see Appendix A for details).
A similar expression is obtained for random Boolean networks by replacing ϕ(q) with ϕBN (q) ≡ δq,1 , simply reflecting
the completely random nature of state transitions.
It is worth noting that, the dynamics of the overlap
becomes deterministic in the limit of large n according to
the central limit theorem, which is the so-called distancelaw [22–26], qt+1,s+1 = ϕ(qts ). In this equation, the equality
holds only at q = 0 and q = ±1, and otherwise |ϕ(q)| < |q|.
Hence, in the limit of large n, the overlap monotonicallyconverges to the stable solution of q = 0, implying that two
distinct states would never converge. On the other hand,√for finite n, the overlap fluctuates with amplitude ∼1/ n
about the deterministic solution (see detailed explanationsin Appendix A). Thus, the overlap can evolve from q < 1
to q = 1 in time, indicating that system’s state eventuallycomes back to one of the previously visited states in a finite
network.

(5)

The Markovian process of Eq. (4) sequentially provides
Pt+l,t (q) for t = 1,2, . . . for some positive time differencel = t − s given an initial distribution at t = 0. The initial
distribution is denoted by Pl,0 (q) ≡ Prob(ql,0 = q). Since theinitial state, σ (0), is selected randomly and independently from
{Jij }, we can set, without losing generality, the initial state tobe σi (0) = 1 for all i (see Appendix B). In this case, the initial
overlap of interest is expressed by1
σi (l)σi (0)n i=1nql,0 =1sgn [hi (l − 1)] .
n i=1
n

=

(6)

If l is small, Pl,0 (q) reflects the memory of the initial state
σ (0) and is hard to evaluate exactly. However, if l is large, themean-field result in Appendix A indicates that {hi (l − 1)|i =
1,2, . . . ,n} follows approximately a zero-centered independent Gaussian distribution with unit variance in the large
network-size limit. This means that the state overlap of Eq. (6)approaches a distribution centered around zero with variance∼1/n. In particular, Pl,0 (q) tends for large l to a binomial
distribution ( n(1+qnl,0 )/2 )2−n , where the probability of ql,0 = ±1032802-2

STRUCTURE OF ATTRACTORS IN RANDOMLY CONNECTED . . .

is approximately 2−n in the large network-size limit. We
confirm this property later with numerical simulations.

PHYSICAL REVIEW E 91, 032802 (2015)

In the remaining part of this section, we characterize inmore detail the quasistationary state in large n limit, fromwhich we extract the structure of attractors.
We first introduce an auxiliary notation

IV. STATE CONCENTRATIONIn this section, we consider how different states concentrate
in time. The Markovian dynamics of Eq. (4) are completely
characterized by the eigenvalues and eigenvectors of thetransition probability W [27]. Let fa (q) and λa ( 1) be theath eigenvector and eigenvalue of W , respectively. We rank
eigenvalues in a descending order, i.e., λ1  λ2  · · ·  λn+1(the number of possible values for q is n + 1). The distributionof the overlap is expressed by a weighted sum of theeigenvectors as
Pt+l,t (q) =n+1
(λa )t Aa fa (q),αt+l,t (q) ≡

+

a=1

1
1
0.1

1

3

1

4

0.01
0.001
5

10

15

20

n

FIG. 1. (Color online) The nontrivial eigenvalues λ3 and λ4 ofthe transition matrix W .

(8)where exp[nαt+l,t (q)]dq = 1 according to the normalization
constraint. With this notation, we can express the dynamics ofEq. (4) by
1
αt+l+1,t+1 (q) = ln W (q|q  )Pt+l,t (q  )dq n
1 + q 1 + ϕ(q  )ln≈ H (q) + max

q
2
2

(7)

where {Aa } is a set of initial coefficients that satisfies Pl,0 (q) =a Aa fa (q). Hence, as the time step increases, Pt+l,t (q)becomes progressively dominated by the components withlarge eigenvalues.
It is easy to see that W has two trivial eigenvectors
f1 (q) = δq,1 and f2 (q) = δq,−1 with degenerate eigenvalues
λ1 = λ2 = 1. Note that δq,q0 is the Kronecker delta function.
The third eigenvector f3 (q) is a nontrivial one and its
eigenvalue λ3 ≈ 1 − exp(−0.41n) exponentially approaches
1 with n (see Fig. 1 for the numerical result). The fourth
eigenvalue converges to λ4 ≈ 0.67 in the limit of large n.
The half-decay time of the ath component is describedby these eigenvalues and given by ta ≡ (ln 2)/(− ln λa ), or
equivalently, (λa )ta = 1/2. There is a clear gap between the
decay time of the third and fourth eigencomponents. This
result indicates that, for large n, the distribution of theoverlap must approach quickly a quasistationary state P∗ (q) ≡3
a=1 Aa fa (q) at around t4 ≈ 1.73 and stay unchanged until
t3 ≈ 0.69 exp(0.41n). In particular, the quasistationary state ischaracterized solely by f3 (q) except at q = ±1.
This analysis also suggests when the mean-field theory
breaks down—the theory is not applicable once the third
eigencomponent significantly decays at around t3 . That is,after an exponential time of exp(0.41n), the distribution of theoverlap becomes the linear combination of f1 (q) and f2 (q),i.e., every state becomes either the same or the sign flip of theothers. However, this never happens in a real system.

1
ln Pt+l,t (q),n1 − q 1 − ϕ(q  )ln+ αt+l,t (q  ) , (9)22where the Laplace method was applied in the second line
assuming large n. Note that, in the above expression, themaximizer q  of the second term is a function of q. In particular,
the well-defined asymptotic solution of Eq. (9), i.e.,
1 + q 1 + ϕ(q  )lnα(q) = H (q) + max
q
2
2
+

1 − q 1 − ϕ(q  )ln+ α(q  ) ,
2
2

(10)

with finite α(q) self-consistently provides the quasistationary
state. Note that Eq. (10) permits arbitrary discontinuity ofα(q) at q = ±1, reflecting that q = ±1 is the sink of theMarkovian process. However, in the following analysis, we
assume continuous α(q).
Next,wedefineindex
βt+l,t (qt+l,t ) ≡1ln
Prob(q|q=
1)
thatcharacterizes
the
t+l,t
t+l+1,t+1nprobability that two states σ (t + l) and σ (t) have overlap
qt+l,t before converging in the next step (qt+l+1,t+1 = 1). This
index is expressed, using the Bayes theorem, in terms of α by
βt+l,t (q  ) =1 W (1|q  )Pt+l,t (q  )lnn
Pt+l+1,t+1 (1)= ln

1 + ϕ(q  )+ αt+l,t (q  ) − αt+l+1,t+1 (1). (11)2This means that, for large n, most of the trajectories that leadto state concentration had an overlap specified by the peaklocation of βt+l,t , i.e., arg maxq  βt+l,t (q  ), in the previous step.
In the case of the randomly connected neural networks
studied here, βt+l,t (q) has two peaks. As shown in Fig. 2(b),one peak is located at q = 1 reflecting the monotonic increasein ln 1+ϕ(q)toward q = 1 and the other peak is located at2
q < 1 reflecting the peak of αt+l,t (q) at q = 0 in Eq. (11).
The q < 1 peak shifts to a larger positive value of q andits amplitude loses the dominance over the q = 1 peak as tincreases because αt+l,t (q) becomes blunt at large t [Fig. 2(a)].
The two peaks become comparable at around t4 . In finite-size
systems, the two peaks become indistinguishable once thedifference of the peak values becomes less than 1/n. The
result indicates that states concentrate mainly from q ≈ 0.5

032802-3

TARO TOYOIZUMI AND HAIPING HUANGPHYSICAL REVIEW E 91, 032802 (2015)

α

(a)
1.0

0.5

0.1

0.5

1.0

q

(b)

β
0.2

0.4

0.6

0.8

1.0

q

0.05

0.2

0.10

0.3
0.4

0.15

0.5

0.20

0.6
0.7

FIG. 2. (Color online) The Markovian dynamics of αt+l,t and βt+l,t in time. (a) The index αt+l,t characterizes the dynamics of distribution
P (qt+l,t ). The line color changes from the lowest curve (the orange curve at t = 0, i.e., αl,0 or βl,0 ) to the yellow, and finally to the gray (thetop curve at t = 10, i.e., αl+10,10 or βl+10,10 ). (b) The index βt+l,l characterizes the dynamics of distribution P (qt+l,t |qt+l+1,t+1 = 1). The result
indicates that states concentrate mainly from q ≈ 0.5 at the beginning but concentrate equally from q ≈ 0.75 and q ≈ 1 at the quasistationary
state. We used αl,0 (q) = H (q) − ln 2 as the initial condition assuming no correlations at starting points. The results hold for any l  1.

at the beginning but concentrate equally from q ≈ 0.75 andq ≈ 1 at the quasistationary state.
These dynamics of the state overlap reflects the specificstructure of attractors as we shall show below. In contrast to theabove situation, for trivial dynamical systems that converge toa unique fixed point (e.g., hi (t + 1) = [1 + hi (t)]/2), βt+l,t (q)has a unique peak, which tends to approach q = 1 at large t,
indicating that most states concentrate from nearby locations.
On the other hand, in random Boolean networks, statesconcentrate randomly from any overlap values. Because moststates are orthogonal to each other for large n, states mainly
concentrate from q ≈ 0 (see Appendix C).

V. STATISTICAL PROPERTIES OF ATTRACTORSIn this section, we analytically describe the statistical properties of attractors for randomly connected neural networks
using the state concentration probability [17]. The state concentration probability pt+1,s+1 characterizes the conditional
probability of σ (t + 1) = σ (s + 1) given that no states up
to time t along the trajectory are the same or the sign
flip of the others. Because of the symmetry, pt+1,s+1 alsocharacterizes the probability of σ (t + 1) = −σ (s + 1) given
the same condition. Hence,
pt+1,s+1 ≡ Prob(qt+1,s+1 = ±1|{qt  ,s  = ±1|t   t,s  < t  }).
(12)This state concentration probability is further approximatedunder the Markovian approximation of Eq. (4) by
pt+1,s+1 ≈W (qt+1,s+1 = 1|qts )P (qts )dqts
qts =±1= exp[nαt+1,s+1 (1)],(13)

which directly follows from Eq. (9). Note that, based on theconsideration of the previous section, we used in the secondline that the result is not sensitive to the exclusion of q  = ±1
from the integral for large n. This is because the maxq  in
Eq. (9) is insensitive to its argument at q  = ±1 unless theinitial distribution Pl,0 (q) is sharply peaked at q = ±1, which
is not the case here [cf. Eq. (6)].
Hence, based on the Markovian property, the probability
that the dynamics starting from σ (0) comes back for the firsttime to σ (0) after l steps without visiting any sign flip ofpreviously visited states is described for large n byP̃ (l) ≡ Prob({q1,0 = ±1},{q2,s = ±1|s = 0,1}, . . . ,{ql−1,s = ±1|s = 0,1, . . . ,l − 2}, ql,0 = 1)1l−2

(1 − 2p2,s ) · · ·= (1 − 2p1,0 )s=0= pl,0 exp

l−1 
t−1


(1 − 2pl−1,s )pl,0s=0ln(1 − 2pt,s ) .

(14)

t=1 s=0
Note that, in the second line of Eq. (14), the factor t−1
s=0 (1 −2pt,s ) describes the probability that the state makes a transition
at time t to a state distinct from {±σ (s)|s = 0,1, . . . ,t − 1}.
The final factor, pl,0 , describes the probability of coming backto the initial state σ (0) after l steps.
Altogether, the probability that a certain state, σ (0), belongs
to a cycle of length l [revisiting σ (0) for the first time after lsteps] is described for large n by [16]P̃ (l)(odd l)P (l) =P̃ (l) + P̃ (l/2) (even l).

(15)

Notably, the probability takes different expressions for odd
and even l. If l is odd, Eq. (14) directly gives the probability.
If l is even, there are two separate kinds of contributions

032802-4

STRUCTURE OF ATTRACTORS IN RANDOMLY CONNECTED . . .

PHYSICAL REVIEW E 91, 032802 (2015)

steps from the initial value of pinit ≈ exp(−0.69n) to theasymptotic value p∞ ≈ exp(−0.46n).
Therefore, P̃ (l) of Eq. (14) can be further approximatedusing pinit and p∞ by
 l−1 t−1FIG. 3. There are two kinds of limit cycles if the cycle length lis even. (a) In the first kind of cycles, the cycle closes without evervisiting the sign flip of previously visited states. (b) In the second
kind of cycles, the state first makes a transition to the sign flip of theinitial state after l/2 steps, i.e., σ (l/2) = −σ (0). If this happens, thecycle must close after l steps.

depicted in Fig. 3. The first contribution is from cycles thatclose without ever visiting the sign flip of their history. The
second contribution is from cycles that involve a transition
at step l/2 to the sign flip of their initial state, which thenguarantees that the cycle closes in l steps.
The final step is to evaluate the state concentration probability pt,s . The initial state concentration probabilities are simply
given bypl,0 ≈ 2−n ≡ pinit

p∞ ≡ lim pt+l,t = exp [nα(1)]

(17)

t→∞

for any l  1, where α(1) = −0.46 from Eq. (10). That is, thestate concentration probability quickly converges in several

2
n=15
n=20
n=25

2 n pl,0

1.5
1
0.5

0

10

20

30

√
where τ ≡ −2/ ln(1 − 2p∞ ) is the characteristic cycle
length that grows exponentially with the system size, consistent with the numerical observations [19]. Note that, inthe first line of Eq. (18), we used the relationship that|p∞ − pts |  |p∞ − pinit | (for any t and s; see Fig. 2) to upperbound the deviation of pts from p∞ . To make the contribution∞ −pinitof the O(2tc l p1−2p) negligible, the approximation in the∞second line assumesp∞ − pinit−(1 − 2p∞ ) ln(1 − 2p∞ )1 − 2p∞.
and l2tc (p∞ − pinit )(16)

for large n and l as discussed in Sec. III. Although thisapproximation is inaccurate for l < 10, it becomes accurate
for large n over a wide range of l that includes the typicalcycle length (Fig. 4).
On the other hand, the state concentration probability pt+l,t
at t  1 is computed sequentially by Eqs. (9) and (13). In
particular, this probability quickly converges within several
steps [tc ≈ 5; see Fig. 2(a)] to the quasistationary value
of

0



p∞ − pinitP̃ (l) = pinit exp
ln(1 − 2p∞ ) + O 2tc l
1 − 2p∞t=1 s=0 2lln(1 − 2p∞ )≈ pinit exp
2
 2
l
(18)
= pinit exp − 2 ,
τ

40

50

l
FIG. 4. (Color online) Simulation results of the initial state concentration probability pl,0 as a function of l when n varies. The results
are obtained based on statistics collected from 2 × 109 networks. The
initial state is always set to σi (0) = 1 for all i (see Appendix B). The
sampling error increases with n because of exponential increase of
the state space. Note that p1,0 = 2−n is an exact result.

4tc

l
(19)
2

The first condition in Eq. (19) requires that − l2 ln(1 −∞ −pinit
2p∞ )2tc l p1−2p, while the second condition ensures that∞p∞ −pinit
1. The range of l specified by Eq. (19) is
2tc l 1−2p∞
roughly 10lexp(0.46n)/(2tc ) at n > 10. Hence, thecharacteristic cycle length τ ≈ exp(0.23n) is well within thisrange. Incidentally, τ is known to also characterize the typicaltransient time scale to enter a limit cycle [16].
The probability of observing a cycle of length l is given byn
P (l)/(lZ) with a normalization constant Z ≡ 2l=1 P (l)/ l.
In this expression, the probability, P (l), of a state belonging
to a cycle of length l should be divided by l to provide thecycle length probability since all states within a cycle sharethe same cycle length. Note that the normalization constantZ represents the probability of a state belonging to a cycle(attractor). Figure 5(a) shows the comparison of numericallyobtained cycle length probability with its theoretical estimate.
Numerical details to collect the statistics of the attractors
are given in Appendix D. The theory nicely captures thisprobability at around the characteristic cycle length, including
the difference in probability for odd and even cycle lengths, asn becomes large. However, the deviation is large for nontypicall in finite networks. The cumulative distribution of cyclelength is similarly obtained by F (l) ≡ (1/Z) ll  =1 P (l  )/ l  ≈ll/2[ 1 P̃ (l  )/ l  dl  + 1 P̃ (l  )/(2l  )dl  ]/Z. The comparison of
F (l) with the numerical results is shown in Fig. 5(b). The
discrepancy tends to become small for larger n [see the inset
of Fig. 5(b)].
The first moment (mean value) and the second momentof the distribution can be computed analytically as well.

032802-5

TARO TOYOIZUMI AND HAIPING HUANGPHYSICAL REVIEW E 91, 032802 (2015)

FIG. 5. (Color online) (a) Probability distribution of cycle lengths. (b) Cumulative distribution of cycle lengths. The numerical data is
obtained from 1000 samples for n = 10, 500 samples for n = 15, and 200 samples for n = 17. The inset shows an enlarged view at small cycle
length.

Their values are given by√
4 πτ [1 − erf(1/τ )]
l =
,
−t
∞
3 1/τ 2 e t dt
l2 =2τ 2 e−1/τ

(20)

2

(21)

,
∞ e−t
1/τ 2 t dt∞−twhere erf(x) = √2π 0 e−t dt and 1/τ 2 e t dt  −γE − α(1)nin the large n limit, where γE = 0.5772 is the Euler constant.
The theoretical predictions are compared with the numericalresults in Fig. 6. The exponential growth of the typicalcycle length is verified, which suggests that chaotic attractors
exist in the state space of a randomly connected neuralnetwork.
x

2

Following the same spirit, one can derive the number of
attractors as 2n Z  −3α(1)n/4 − 3γE /4, from which a linear
dependence [12,14,16] is confirmed [see also Fig. 7(a); thelinear fit gives the slope 0.360 ± 0.010, compatible with thetheoretical value 0.342). Another interesting quantity is thenumber of the attractive states Natt belonging to all cycles(e.g., a cycle of length l has l attractive states), which is
expected to grow exponentially with the network size n. Thisnquantity is evaluated by our theory as Natt = 2n 2l=1 P (l),and can be quantified as the growth rate (entropy density) s =limn→∞ n1 ln Natt . In the large n limit, we obtain s = −α(1)/2,
which is compared with the numerical results at finite n. As
shown in Fig. 7(b), as n increases, s decreases, approaching
the asymptotic limit 0.2277.
The deviation at small n (or at l far from the characteristic
length; see Fig. 5) comes from three approximations. One isFIG. 6. (Color online) The first (mean) and second moment of the cycle length distribution. Theoretical predictions and numericalsimulations are compared. The results are averaged over many random realizations of the networks (from 1000 samples for n = 10 to
100 samples for n = 18).
032802-6

STRUCTURE OF ATTRACTORS IN RANDOMLY CONNECTED . . .

PHYSICAL REVIEW E 91, 032802 (2015)

FIG. 7. (Color online) (a) Linear dependence of the number of cycles on network size n. (b) Entropy density of the attractive states definedby s = n1 ln Natt . As n increases, the numerical data approach the theoretical prediction.

Eq. (16), which becomes invalid at small l where pinit also
depends on l, but Eq. (16) becomes reasonable for large l (asoccurs in our case where the typical cycles are long). The
second one is Eq. (18). This approximation is valid in therange of l specified by Eq. (19). Note that this condition isconsistent with the numerical results shown in Fig. 5. The lastapproximation is Eq. (4), which breaks down for small n at
which two or more time-steps memory should be considered.
In the large network size limit, these approximations become
exact and the dynamics can be described by a Markovian
process in terms of the state overlap. Thus, as we focus on the
structure of attractors at a relatively large but finite n, theseeffects are not significant.

VI. CONCLUSION

In this work, we studied the deterministic dynamics of arandomly connected neural network and proposed a simple
Markovian stochastic process to describe the evolution of theoverlap of two states along the dynamics trajectories. The
properties of the state concentration can be studied by a meanfield computation, and furthermore, the theoretical cumulative
distribution of cycle length is compared with the numericalsimulation results. The typical length of cycles is predicted
and observed to grow exponentially with the network size. The
number of attractive states on all cycles has also an exponential
growth with the network size, and its typical value can also be
predicted by our theory.
Our theory should have the potential to be generalized totreat more complex situations, e.g., couplings between neurons
are correlated, where one time-step memory is not enough todescribe the dynamics and strong memory effects induced byretarded self-interaction could be incorporated by introducing
a back-action field (two time-steps memory) [18]. The current
analysis is also restricted to the parallel type of dynamics,
whereas the sequential (asynchronous) dynamics seems to
be more natural, and our current method may apply to thistype of dynamics, although the computation will become more
complicated. However, the statistical properties of attractors

would not change qualitatively, as expected from numericalsimulations [12].
Our work is expected to provide insights towards understanding how the neural network processes information andstores temporal sequences [28], which will be left for future
study.

ACKNOWLEDGMENTS

We are grateful to Ugo Bastolla, Naoki Masuda, Hiroyasu
Ando, and Shun-ichi Amari for useful discussions. This work
was supported by RIKEN Brain Science Institute and the Brain
Mapping by Integrated Neurotechnologies for Disease Studies(Brain/MINDS) by the Ministry of Education, Culture, Sports,
Science and Technology of Japan (MEXT).

APPENDIX A: DYNAMIC FUNCTIONALINTEGRAL METHOD

We compute the dynamics of the state overlap using the
dynamic functional integral (mean-field) method (see [5] forsimilar calculations). In this section, we express time indicesas lowercase characters, e.g., hi (t) = hit , and follow theconvention that summations are neglected
if the same indices
appear twice in an expression, e.g., j Jij σj t = Jij σj t .
Let us first define the ensemble of state trajectories,h = {hit |i = 1, . . . ,n, t = 0,1, . . . }, averaged over different
networks:δ(hit − Jij σj t )P (h) ≡
i,t

,

(A1)

J

where δ(·) is the Dirac delta function and [·]J is the averageover the random couplings. In the following, we denote by ·
an average with respect to P (h).

032802-7

TARO TOYOIZUMI AND HAIPING HUANGPHYSICAL REVIEW E 91, 032802 (2015)

The joint distribution of the overlap q = {qts }t>s is1δ qts − σj t σj s
P (q) ≡
n
t>s


=

δ(hit − Jij σj t )

dhit
i,t


=

=

=

=

=

i,t

dhit d ĥit
2π
i,t

J

[exp(i ĥit hit − i ĥit Jij σj t )]J

dhit d ĥit
2π
i,t





1
1
δ qts − σj t σj s
exp i ĥit hit − qts ĥit ĥis2nt>snd q̂ts
2πt>s
nd q̂ts2πt>s

dhit d ĥit
2π
i,t



1
exp ni q̂ts qts − i q̂ts σj t σj s + i ĥit hit − qts ĥit ĥis2

exp [−nf (q,q̂)] ,(A2)dH exp L,1
L ≡ i ĥt ht − qts ĥt ĥs − i q̂ts σt σs ,2

(A3)

and dH ≡ t (dht d ĥt )/(2π ). In the derivation, we haveused the Fourier transformation of the delta function δ(x) =d x̂/(2π ) exp(i x̂x) for each delta function in Eq. (A2), andwe have taken the average over the independent Gaussianvariables {Jij } of mean 0 and variance 1/n. For large n, thedistribution of Eq. (A2) is well approximated by a Gaussian
distribution, where the peak is specified by the saddle-point
equations:

0=

∂f
1
= −i q̂ts + ĥt ĥs L ,
∂qts
2
∂f
= −qts + σt σs L ,∂i q̂ts(A4)

with average ·L ≡ ( ·eL dH )/( eL dH ). We can easily seethat q̂ = 0 is a solution of Eq. (A4) [5]. Hence, if q̂ = 0, theaverage ·L0 ≡ ·L |q̂=0 is an average over Gaussian h of mean
ht L0 = 0 and covariance ht hs L0 = qts , which simplifies thesaddle-point equation of Eq. (A4) in terms of a closed-form
expression of q by
qt+1,s+1 = σt+1 σs+1 L0




=
DxDy sgn(x)sgn qts x + 1 − qts2 y= ϕ(qts ),1
δ qts − σj t σj s
n
t>s





1
1
ĥit ĥis σj t σj sδ qts − σj t σj s
exp i ĥit hit −
2n
n
t>s

f (q,q̂) ≡ −i q̂ts qts − ln

0=



1
δ qts − σj t σj snt>sdhit d ĥit
2πi,twhere we have defined the action





(A5)

with ϕ(q) ≡ √(2/π ) arcsin(q) and a Gaussian measure Dx ≡exp(−x 2 /2)/ 2π dx. Note that |ϕ(q)|  |q| and the equality

holds only at q = 0 and q = ±1. Hence, unless q = 1 initially,
the overlap rapidly converges in a few steps to zero in the
n → ∞ limit.
The order parameters fluctuate around the saddle-point
solution of Eq. (A5) for finite n. This fluctuation of q andq̂ is characterized to the leading order by the Hessian matrixof f , i.e.,
Ats,t  s  ≡∂ 2f
1 ∂ ĥt  ĥs  L0=,∂qts ∂qt  s 
2∂qts

Bts,t  s  ≡∂ 2f∂ σt  σs  L0= −iδt,t  δs,s  + i,∂qts ∂ q̂t  s ∂qts

Cts,t  s  ≡∂ 2f= σt σs σt  σs  L0 − σt σs L0 σt  σs  L0 ,∂ q̂ts ∂ q̂t  s (A6)

for t > s and t  > s  , where the Hessian matrix is evaluated at
the saddle-point solution of the order parameters, i.e., q̂ = 0
and the solution of Eq. (A5).
In the current setup, the Hessian matrix is simply given byAts,t  s  = 0,
Bts,t  s  = −iδt,t  δs,s  + iϕ  (qts )δt  ,t+1 δs  ,s+1 ,(A7)

Cts,t  s  = δt,t  δs,s  + O(q 2 ),for t > s and t  > s  , where ϕ  (q) ≡ dϕ(q)/dq. Note that theO(q 2 ) contribution in Cts,t  s  can be more explicitly estimated,for example by applying Plackett’s approximation [29]. Here,
we would like to evaluate the (n multiplied) covariance of theoverlap parameter, Ãts,t  s  ≡ n Cov[qts ,qt  s  ]. By applying thematrix inversion lemma, we find that its inverse is032802-8

Ã−1 = A − BC −1 B T = (iB)(iB)T + O(q 2 ).

(A8)

STRUCTURE OF ATTRACTORS IN RANDOMLY CONNECTED . . .

This relation indicates that for small q the linear combination,√η ≡ n(iB)T δq,(A9)of the fluctuation of the overlap parameter, δq, is white
Gaussian random variables. To see this, one can apply thetransformation of variables and find that

 n

√
P (η) = δ[η − n(iB)T δq] exp − δq T Ã−1 δq dδq
2


1 T
(A10)
∼ exp − η η .
2Thus, Eq. (A9) indicates that the finite-size fluctuations of theorder parameter are described by1√ ηts = [δt  ,t δs  ,s − ϕ  (qt  s  )δt,t  +1 δs,s  +1 ]δqt  s n= δqts − ϕ  (qt−1,s−1 )δqt−1,s−1 .

(A11)

Altogether, summarizing that qts = ϕ(qt−1,s−1 ) in the n →∞ limit and that the finite-size correction is described by
Eq. (A11), we obtained, for finite n,1
qts = ϕ(qt−1,s−1 ) + √ ηts + O(q 2 ),
n

(A12)

which is a simple Markovian process that involves white
Gaussian noise of variance 1/n.
Recalling the definition of the overlap parameter, qt+1,s+1 =sgn(hit )sgn(his )/n, and that hi with different i tend to become
independent in the n → ∞ limit, we know that the overlap
parameter must be distributed approximately according toa binomial distribution. Extrapolating this observation, theresult of Eq. (A12) is consistent with the Markovian dynamicsof
P (qt+1,s+1 ) = W (qt+1,s+1 |qts )P (qts )dqts(A13)with the binomial transition probability 
2m − n qts
W qt+1,s+1 =
n 
 

 
1 + ϕ(qts ) m 1 − ϕ(qts ) n−m
n
,
=
2
2
m

PHYSICAL REVIEW E 91, 032802 (2015)

To see this point, we consider a simple transformation ofvariables,σ̃i (t) = σi (t)σi (0).

(B1)

The state overlap is also described in terms of these transformed variables by qts = (1/n) i σ̃i (t)σ̃i (s) and the initialstate is given by σ̃i (0) = 1 for all i.
These transformed variables follow the same update rule as
the original one,⎛⎞J˜ij σ̃j (t)⎠ ,(B2)σ̃i (t + 1) = sgn ⎝jexcept that the coupling matrix is given by J˜ij = σi (0)Jij σj (0)instead of Jij . Notably, the distribution of {J˜ij } is the same asthat of {Jij } as long as σ (0) is chosen independently of {Jij }.
Therefore, to study the dynamics of the state overlap, we canalternatively study the dynamics of these transformed variables
with the initial condition {σ̃i (0) = 1|i = 1,2, . . . ,n}.
APPENDIX C: THE DYNAMICS OF THE STATE OVERLAPIN RANDOM BOOLEAN NETWORKS

The dynamics of the state overlap in random Boolean
networks is described by Eq. (9) with ϕBN (q) = δq,1 , which
is simply
H (q) − ln 2 + αts∗(q = 1)
αt+1,s+1 (q) =(C1)max{αts (1), − ln 2 + αts∗ } (q = 1),
where αts∗ ≡ maxq  =1 αts (q  ). Let us assume that there is noperfect overlap of states initially, i.e., Prob(ql,0 = 1) = 0.
∗∗This means that αl,0 (1) < − ln 2 + αl,0and limn→∞ αl,0= 0,because the initial overlap distribution Pl,0 (q) = exp[nαl,0 (q)]must be normalized. Thus, the dynamics of Eq. (C1) converges
in one step to a stationary solution
α(q) = H (q) − ln 2.

(C2)

Moreover, we have from Eq. (11)(A14)

where m indicates the number of units taking the same state attime t + 1 and s + 1.
In summary, this result shows that the Markovian dynamics
of Eq. (A13) provides a good approximation of the dynamics
of the overlap parameter once O(q 2 ) terms become negligible
near the stationary state.
APPENDIX B: THE DYNAMICS OF THE STATEOVERLAP DOES NOT DEPEND ON THE
INITIAL STATE σ (0)A specific choice of the initial state σ (0) is not important tostudy the dynamics of the state overlap for a random ensembleof networks as long as σ (0) is selected independently of the
network connections {Jij }. Without losing generality, we canset σi (0) = 1 for all i.

1 + ϕBN (q  )+ H (q  )β(q  ) = ln2
H (q  ) − ln 2 (q  = 1)=
0(q  = 1).

(C3)This indicates that states mainly concentrate from q = 0 if
they do not already concentrate.
This analysis also provides important information about
the eigenvalues of the transition matrix W at the large
network size limit. The first eigenvalue is trivial, λ1 = 1, withthe eigenfunction f1 (q) = δq,1 , indicating that states neverseparate once they concentrate. The second eigenvalue, λ2 ≈1, is a nontrivial one that corresponds to the quasistationary
state with the eigenfunction f2 (q) = exp[nα(q)], where α(q)is given by Eq. (C2). The other eigenvalues λk for k  3 areall zero because the distribution of the overlap converges in asingle step to the quasistationary state. Furthermore, the fact

032802-9

TARO TOYOIZUMI AND HAIPING HUANGPHYSICAL REVIEW E 91, 032802 (2015)that state concentration happens with probability 2−n at eachtime step suggests that λ2 = 1 − 2−n .
APPENDIX D: SIMULATION DETAILS OF
THE DYNAMICS

The total number of states is 2n . They form a state set
called S. We also denote a path set P recording the states on a
dynamics trajectory. Only the state index is stored in both sets.
Step 1. Choose the first state σ 1 in S as a starting point
for the parallel dynamics, and remove this state from S at thesame time.

[1] H. Ko, S. B. Hofer, B. Pichler, K. A. Buchanan, P. J. Sjöström,and T. D. Mrsic-Flogel, Nature (London) 473, 87 (2011).
[2] K. Rajan, L. F. Abbott, and H. Sompolinsky, Phys. Rev. E 82,
011903 (2010).
[3] H. Sompolinsky, A. Crisanti, and H. J. Sommers, Phys. Rev.
Lett. 61, 259 (1988).
[4] C. van Vreeswijk and H. Sompolinsky, Science 274, 1724
(1996).
[5] T. Toyoizumi and L. F. Abbott, Phys. Rev. E 84, 051908 (2011).
[6] S. Ostojic, Nat. Neurosci. 17, 594 (2014).
[7] J. Aljadeff, M. Stern, and T. Sharpee, Phys. Rev. Lett. 114,088101 (2015).
[8] D. Sussillo and L. Abbott, Neuron 63, 544 (2009).
[9] R. Laje and D. V. Buonomano, Nat. Neurosci. 16, 925 (2013).
[10] G. Parisi, J. Phys. A: Math. Gen. 19, L675 (1986).
[11] E. Gardner, B. Derrida, and P. Mottishaw, J. Phys. 48, 741(1987).
[12] H. Gutfreund, J. D. Reger, and A. P. Young, J. Phys. A: Math.
Gen. 21, 2775 (1988).
[13] H. Rieger, M. Schreckenberg, and J. Zittartz, Z. Phys. B 74, 527(1989).
[14] M. Schreckenberg, Z. Phys. B 86, 453 (1992).

Step 2. σ 1 evolves to σ  by one step of the parallel
dynamics (all neurons’ states are updated for one time).
Step 2.1. If σ  ∈ S, remove it from S, put the index ofσ into P, and continue to perform the parallel dynamics,i.e., let σ 1 = σ  , then go to step 2.
Step 2.2. Otherwise, compare σ  with the one in set
P and if they coincide with each other, a new cycle isidentified and the length is recorded at the same time, then
go to step 3; otherwise, no new cycle is found and go tostep 3.
Step 3. Go to step 1 until the set S becomes empty.

[15] A. Crisanti, M. Falcioni, and A. Vulpiani, J. Phys. A: Math. Gen.
26, 3441 (1993).
[16] U. Bastolla and G. Parisi, J. Phys. A: Math. Gen. 30, 5613
(1997).
[17] S.-I. Amari, H. Ando, T. Toyoizumi, and N. Masuda, Phys. Rev.
E 87, 022814 (2013).
[18] H. Huang and Y. Kabashima, J. Stat. Mech.: Theory Exp. (2014)
P05020.
[19] K. Nützel, J. Phys. A: Math. Gen. 24, L151 (1991).
[20] S. A. Kauffman, J. Theor. Biol. 22, 437 (1969).
[21] S. A. Kauffman, Nature (London) 224, 177 (1969).
[22] B. Derrida and Y. Pomeau, Europhys. Lett 1, 45 (1986).
[23] B. Derrida and D. Stauffer, Europhys. Lett. 2, 739 (1986).
[24] B. Derrida, J. Phys. A 20, L721 (1987).
[25] S.-I. Amari, Kybernetik 14, 201 (1974).
[26] K. E. Kürten, Phys. Lett. A 129, 157 (1988).
[27] R. G. Gallager, Stochastic Processes: Theory for Applications (Cambridge University Press, Cambridge, UK,2014).
[28] T. P. Vogels, K. Rajan, and L. F. Abbott, Annu. Rev. Neurosci.
28, 357 (2005).
[29] R. H. Bacon, Ann. Math. Statist. 34, 191 (1963).

032802-10