The dynamics of critical Kauffman networks under asynchronous stochastic update

0.13 MB·English

by  Florian Greil

The dynamics of critical Kauffman networks under asynchronous stochastic update
Florian Greil and Barbara Drossel
Institut fu¨r Festk¨orperphysik, TU Darmstadt, Hochschulstraße 6, 64289 Darmstadt, Germany
(Dated: February 2, 2008)
We show that the mean number of attractors in a critical Boolean network under asynchronous
stochastic update grows like a power law and that the mean size of the attractors increases as a
5
stretchedexponentialwiththesystemsize. Thisisinstrongcontrasttothesynchronouscase,where
0
thenumberof attractors grows faster than any power law.
0
2
PACSnumbers: 89.75.Hc,05.65.+b,89.75.Hc

n
a
J RandomBooleannetworkswereintroducedin1969by figuration space, and their number was recently shown

5 Kauffman [1, 2] as a simple model for complex systems numerically[1,19,20,21, 22]andanalytically[23,24]to
consisting of units that interact via directed links. They growfasterthananypowerlawwiththe networksizeN.
]
n are used to model social and economic networks [3, 4], In contrast, we will show in the following that for asyn-
n neural networks, and gene or protein webs [5]. chronous stochastic update the mean number of attrac-
- ArandomBooleannetwork(RBN)is adirectedgraph tors grows as power law while their size increases like a
s
di with N nodes each of which takes a Boolean value σi ∈ stretched exponential with N. As the dynamics is no
0,1 . The number k of incoming edges is the same for longer deterministic, an appropriate definition of an at-
at. a{ll n}odes, and the starting points of the edges are cho- tractor must be given. An attractor is a subset of the
m sen at random. Usually these models are updated syn- configuration space such that for every pair of configu-
chronously, σ (t + 1) = c (σ (t),...,σ (t)), where rations on the attractor there exists a sequence of up-
- i i di,1 di,k
d the Boolean coupling function c of node i is chosen at dates that leads from one configuration to the other. In
i
n
random among a given set of functions and d denotes [7], such an attractor is called “loose attractor”. Start-
o i,j
the j-th input of node i. The configuration of the sys- ing from a random initial configuration, the system will
c
[ tem~σ σ1,...,σN thus performs a trajectory in con- eventually end up on an attractor.
≡{ }
figurationspace. Critical networksareofspecialinterest Let us first consider a set of nodes arranged in a
1
v [6]. Their dynamics is at the boundary between a frozen loop. Such loops occur as relevant components of criti-
1 phasewhereinitiallysimilarconfigurationsconverge,and calnetworks. Nontrivialdynamics occurs onlyif the two
8 a chaotic phase where initially similar configurations di- constant Boolean functions are omitted, the remaining
0 verge exponentially. Booleanfunctionsbeing“copy”( ),and“invert”( ). A
1 ⊕ ⊖
However, synchronous update is highly improbable in loopwith n inversions canbe mapped bijectively onto
0 ⊖
5 real networks, and it is used under the tacit assumption aloopwithn 2inversionsbyreplacingtwo withtwo
− ⊖
0 that going to asynchronous update will not modify the and by inverting the state of all nodes between these
⊕
/ essentialpropertiesofthesystem[7]. Buttherearegood two couplings. It is therefore sufficient to consider loops
t
a reasons to doubt the validity of this assumption. For in- withzeroinversions(“even”loops)andloopswithonein-
m
stance, for cellular automata it is well-known that some version(“odd”loops). The positionofthe -couplingin
⊖
- of the self-organizationis an artifact of the central clock the odd loop is called the twisted edge. For synchronous
d
[8]. For RBNs there is also recent evidence that devia- update, each configuration is on a cycle in configuration
n
o tions from synchronous update modify considerably the space and occurs again at most after N (2N) time steps
c attractors of the dynamics [9, 10]. for even (odd) loops. The number of cycles increases
:
v In this paper, we investigate a version of the model thereforeexponentiallywithN. Incontrast,mostconfig-
i where at each computational step one node is cho- urationsaretransientintheasynchronouscase,andonly
X
sen at random and is updated. A model with this two (one) attractors are left. The reason for this is that
r
a asynchronous stochastic updating scheme is often called a domain of neighboring nodes that have the same value
Asynchronous RBN (ARBN)[11, 12, 13], while the Clas- increasesordecreaseswithprobability1/N percomputa-
sical synchronous RBN is referred to as CRBN. ARBNs tionalstep. Thedomainsizethereforeperformsarandom
aremostlystudied numericallywiththe focus onvarious walk, and for an even loop no domain is left after of the
measures of stability [14, 15, 16, 17, 18]. ARBNs were order of N3 updates. The attractors are the two fixed
observed to be capable to generate an ordered behavior, points. For an odd loop, the nodes of a domain change
but the detailed properties of attractors have not been their state at the twisted edge, and the total number of
studied yet. domain walls is therefore odd. The attractor contains
We will show mostly analytically that the number of only one domain wall that moves around the loop, and
attractors changes completely when going from CRBNs the attractor comprises 2N configurations. The dynam-
to ARBNs. For CRNBs, attractors are cycles in con- ics on such a loop is closely related to the Glauber dy-
2
namics [25] of a one-dimensional Ising chain with cyclic der N−2/3, implying that of the orderN N−2/3 =N1/3
·
boundaryconditionattemperatureT =0,wherethedo- nodes have 2 nonfrozen inputs. Consequently, the pro-
mains also shrink and grow with a fixed rate and where portion N−1/3 of nonfrozen nodes (whether they are
the equal-time correlation function obeys a scaling form relevant or not) have 2 nonfrozen inputs. This result,
C(r,t)=f(r2t−1)[26]. Thedynamicsofanoddloopcan together with the other just mentioned features of the
be mapped ontothe dynamicsofanIsingchainwith one modelisconfirmedby(unpublished)studiesofourgroup.
negative coupling. It is a frustrated system in which not The other relevantnodes haveone relevantinput (as the
all bonds can be satisfied simultaneously. To conclude, second input comes from a frozen node). The remain-
by going from synchronous to asynchronous update, the ing non-frozen nodes (of the order of N2/3) are on trees
number of attractors of a loop is reduced from an expo- rootedin relevantnodes. Justas for thek =1 networks,
nentially large number to 1 or 2. This was also pointed thereareoftheorderofln(N)independentrelevantcom-
outin[27],whereadifferentasynchronousupdating rule ponents [21]. In contrast to the k = 1 networks, these
is used. componentsarenotalwayssimpleloops,butmaycontain
Letusnextconsidercriticalnetworkswithconnectivity severalnodeswithtworelevantinputs. Inordertoobtain
k = 1, where the Boolean coupling functions are again results for the number of attractors of the networks, we
“copy” and “invert”. Relevant nodes are those nodes have to investigate the attractors of such relevant com-
whose state is not constant and that controlat least one ponents.
relevant element [21]. They determine the attractors of Let us first consider relevant components that contain
thesystem. In[28],exactresultsforthetopologyofk =1 one node with two relevant inputs. These are two loops
networks are derived, and in [24] it is shown that the with one interconnection (# #-component), and a loop
−
numberofattractorsofacriticalk =1networkincreases with an extra link ( -component). The dynamics un-
⊘
fasterthananypowerlaw. Thenumberofrelevantnodes der synchronous update for such components is studied
scales as √N, and they are arranged in of the order of in [29], and it is found that the number of attractors in
ln(N) loops. The remaining nodes sit on trees rooted bothsystemsincreasesexponentiallywiththe numberof
in these loops. Under asynchronous update, each loop nodes. Withasynchronousupdate,thenumberofattrac-
hasatmosttwoattractors. Thenodesontreesrootedin tors becomes very small.
evenloopsarefrozenbecausetheloopisonafixedpoint. We discuss first two loops with one interconnection.
The nodes on trees rooted in odd loops can assume any The first loop is independent of the second loop, and its
combinationofstates,sinceonecanfindtoeachpossible attractor is either a fixed point (if the loop is even), or
state of a tree a sequence of updates that generates it. it has one domain wall moving around the loop. If the
Sinceeachloopisevenoroddwithequalprobability,the firstloopis onafixedpoint, itprovidesaconstantinput
mean number of attractors of networks with n relevant to the second loop, which therefore behaves like an even
loops is loop, or an odd loop, or a frozen loop. The system can
have at most three attractors. If the first loop is odd, if
1 n 3 n 3 ln(N/2) provides a changing input to the second loop, which can
2i = , (1)
2n (cid:18)i(cid:19) (cid:18)2(cid:19) ≃(cid:18)2(cid:19) therefore have an attractor that contains an arbitrary
Xi
and fluctuating number of domain walls. Consequently,
whichisapowerlawinN. Thenumberofnodesintrees a loop that has one external input can show one out of
is of the order of N. On average, half of the trees are four different types of behavior on an attractor:
rooted in odd loops. Consequently the mean attractor
1. The loop can be at a fixed point 0.
size increases exponentially with N.
Finally,weinvestigatethemostfrequentlystudiedcrit- 2. The loop can be at a fixed point 1.
ical networks with connectivity k = 2, where each of
the16possibleBooleancouplingfunctionsischosenwith 3. There is exactly one domain wall which moves
equalprobability. AllnodesapartfromtheorderofN2/3 around the loop.
nodes are frozen. This follows for instance from the fac-
4. Thenumberofdomainwallsinthe loopfluctuates.
tor ǫ3 in the last term of Eqn. (9) of [23], which implies
that only the fractionN−1/3 ofallnodes undergoa non- (Without loss of generality, we have assumed that all
frozen sequence of states in time. Numerical support for coupling functions for nodes with one input are “copy”.)
this result is presented in [20]. Nowweturntoaloopwithanextralink. Wecanagain
The number of relevant nodes scales as N1/3 [20], and assumethatallcouplingfunctions fornodeswithonein-
only a fraction N−1/3 of these relevant nodes have two put are “copy”. If this component has one fixed point
relevant inputs. This last statement follows from the (two fixed points), it is (they are) the only attractor(s).
evaluation of Eqns. (6) and (9) in [23] in saddle-point This isbecause onecanreachafixedpointfromanarbi-
approximation, where the width of the first term in the traryinitialstate by updating one node after anotherby
directionperpendiculartothelineofmaximaisoftheor- goingaroundthe loopinthe directionofthe links. After
3
at most two rounds the fixed point is reached. Only if twodomainwallsonthis section. Byrepeating the same
thecouplingfunctionforthenodewithtwoinputshasno sequence of updates, every even number of domain walls
fixed point, a more complicated attractor occurs. With- can be created in this section, and odd numbers can be
outlossofgenerality,wechoosethisfunctiontohavethe created by moving one domain wall out of the section.
output 1 if andonly if both inputs are 0. By considering If s is the number of sections, an upper bound for the
thepossibleupdatesequences,onefindsthatthecompo- numberofattractorsofthe componentis thereforegiven
nent can accumulate a large and fluctuating number of by 4s.
domain walls, as illustrated in Fig. 1. We checkedthis analyticalresult by computer simula-
tions. In order to make sure that we capture all attrac-
tors, we did a complete search of state space, which can
only be done for small networks. Starting from an ini-
tial state, we did N3 updates before assuming that the
a b c
systemis onanattractor,and we made sure that the re-
sults are not changed when the length of the initial time
periodis varied. Allstatesthatcanbe reachedfromthis
last state are on the same attractor as this state. All
other states that have been visited are marked as tran-
d e f sient states. Then we start with an unvisited state as
newinitialconditioninordertoidentifyfurthertransient
states and attractors. We constructed relevant compo-
nents by starting with one loop of a certain size, and by
iteratively inserting additional connections between two
randomly chosen nodes. The new connection contains
g h i a randomly chosen number of 1 to 4 nodes (such that
a section can contain two domain walls in its interior).
In these networks, the number of sections, s, is identi-
cal to the number of nodes with 2 inputs, µ. After each
insertion, we evaluated the number of attractors for dif-
FIG.1: Apossiblesequenceofconfigurationsofaloopwitha
ferent choices of coupling functions. This procedure was
cross-link showing how multiple domain walls are generated.
The coupling function of the node # with two inputs is such repeated more than 750000 times. The largest number
thatonlytwodarkinputsleadtolightoutput,allothercom- of attractors found in a system is shown in Table I as
binationgiveblackoutput. After(e)theproceduredescribed function of s=max(1,µ).
by (a) to (d) is repeated to obtain (f) and similarly for (g),
(h) and (i). µ νmax 4s realizations
0 2 1 227683
Equipped with these results, we now consider compo- 1 2 4 167370
nents with severalnodes with 2 inputs. We define a sec-
2 9 16 138541
tion to be a sequence of nodes starting at a node with
3 8 64 110263
twoinputs andending rightbefore anode with2 inputs.
4 23 256 73268
Suchasequencecanbranchandhaveseveralendpoints.
5 25 1024 40770
Clearly, the number of sections is the number of nodes
6 23 4096 15727
with 2 inputs; a simple loop is counted as one section.
A section is controlled by its first node, which is the one
TABLE I: Maximum number of attractors νmax as function
with 2 inputs. Just as for the loop with one external in- of the number of nodes with 2 inputs, µ, for networks with
put, a section can show on an attractor one out of the upto17nodes. Networkswithhigherµareprobedlessoften
four different types of behavior listed above. This is be- because if only short links are added there is no node left
causeallstatesthathavemorethanasingledomainwall which has not already two inputs.
inagivensectionmustbepartofthesameattractor. We
This leads us to the conclusion that a network con-
showthis bythe followingargument: Assumethatonan
sisting of the order of ln(N) relevant components, with
attractorthere occur twodomainwallsina section. The
componentnumberihavingµ nodeswith2inputs,can-
twodomainwallscanbedestroyedbyupdatingallnodes i
not have more than
betweenthe twowalls,suchthatthe domainenclosedby
the walls vanishes. A configuration with no wall on the ν =4max(1,µ1) 4max(1,µ2) ... 4max(1,µlnN) 4ln(N)+µ(2)
section (and with the state of all other sections unmod- · · · ≤
ified) is therefore also part of the attractor, and there attractors. ThisisapowerlawinN iftheprobabilitydis-
exists consequentlya waybackto the configurationwith tributionforthevalueofµbecomesindependentofN for
4
largeN. Indeed,aswehavementionedabove,eachofthe manmodelsincreasesasapowerlaw,andwethusregain
N1/3 relevant nodes has two (randomly chosen) relevant the original claim by Kauffman - albeit for models with
inputs with probability aN−1/3 (with some constant a). a different update rule than the original one.
Since this probability is independent for different nodes,
We thank V. Kaufman for useful discussions.
the value of µ is distributed for large N according to a
Poisson distribution with a mean a.
We thus have shown that in critical k = 2 networks
with asynchronous stochastic update, the number of at-
tractors grows as a power law in N, which is in strong [1] S. A. Kaufmann,J. Theo. Bio. 22, 437 (1969).
contrast to the synchronous case, where the number of [2] S. Kaufmann, Nature224, 177 (1969).
attractorsincreaseslikeastretchedexponentialfunction. [3] J. M. Alexander, in Philosophy of Science Assoc. 18th
We conclude with a discussionofthe size ofattractors Biennial Mtg (2002).
in these networks. There are of the order of N2/3 nodes [4] M. Paczuski, K. E. Bassler, and A. Corral, Phys. Rev.
on the trees rooted in the relevant components. These Let. 84, 3185 (2000).
[5] S.Kauffman,C.Peterson,B.Samuelsson,andC.Troein,
nodes in trees can adopt any configuration if the node
inProceedings ofthe National Academy of Sciences USA
they are rooted in can switch its state on an attractor.
(2003), vol. 100, pp.14796–14799.
Sinceanon-vanishingfractionofallrelevantnodesswitch [6] M. Aldana-Gonzalez, S. Coppersmith, and L. P.
theirstatesonanattractor,the sizeofthe attractorisof Kadanoff, Perspectives and Problems in Nonlinear Sci-
the order of ence pp.23–89 (2003).
[7] I.HarveyandT.Bossomaier,inFourthEuropeanConfer-
2N2/3 =exp N2/3ln2 . (3) enceonArtificialLife(ECAL97),editedbyP.Husbands
(cid:16) (cid:17) and I.Harvey (MIT Press, 1997), pp.67–75.
[8] T. E. Ingerson and R. L. Buvel, Phys.D 10, 59 (1984).
Thesizeoftheattractorsgrowslikeastretchedexponen-
[9] K. Klemm and S. Bornholdt (2004), cond-mat/0411102.
tial function and therefore faster than any power law.
[10] K. Klemm and S. Bornholdt (2003), q-bio/0309013.
Manyofourresultsholdalsoforotherkindsofstochas-
[11] C. Gershenson, in Artificial Life IX (MIT Press, 2004).
tic asynchronousupdate, for instanceif a certain(small) [12] C. Gershenson, J. Jan Broekaert, and D. Aerts, in Ad-
fraction of nodes is updated at each step, or if the time vancesinArtificialLife,7thEuropeanConferenceECAL
interval between two updates of a node is peaked at a 2003, edited by W. Banzhaf, T. Christaller, P. Dittrich,
J. T. Kim, and J. Ziegler (Springer, 2004), vol. 2801 of
value τ and Gaussian distributed around it. (The lat-
Lecture Notes in Artificial Intelligence, pp.615–624.
ter case describes our system for large N, when the net-
[13] C. Gershenson (2003), submitted to Phys. D.
work of relevant nodes is coarse-grainedsuch that of the [14] E. A.Di Paolo, BioSystems 59, 185 (2001).
order of N1/3 neighboring nodes are replaced by a sin- [15] P. Rohlfshagen and E. A. Di Paolo, BioSystems 73, 141
gle node that receives a delayed input from the previous (2004).
node.) Inthesemodifiedstochasticmodels,domainwalls [16] J. Hallinan and J. Wiles, in Ninth International Confer-
onanisolatedloopcanannihilate,butcannotbecreated ence on the Simulation and Synthesis of Living Systems
again, leading to the same attractors as with the com- (ALife9) (2004).
[17] B.MesotandC.Teuscher,inAdvances inArtificialLife.
pletelystochasticupdate. However,thestateofthetrees th
Proceedings ofthe 7 European Conference, ECAL2003,
rooted in the loops will be dominated by a few domain
edited by W. Banzhaf, T. Christaller, P. Dittrich, K. J.
wallswhenthedistributionofupdatetimesbecomesnar-
T., and J. Ziegler (Springer-Verlag, Berlin, Heidelberg,
row,withstateswithmoredomainwallsoccurringrarely. 2003),vol.2801ofLectureNotesinArtificialIntelligence,
Similarly, relevant loops that receive input from outside, pp. 367–376.
andrelevantcomponentswithnodeswithtwoinputswill [18] M.T.MatacheandJ.Heidel(2004), submittedtoPhys.
have attractors dominated by few domain walls, and the Rev. E.
[19] S. Bilke and F. Sjunnesson, Phys. Rev. E 65, 016129
actual size of attractors becomes in the thermodynamic
(2002).
limit N → ∞ smaller than the size obtained by consid- [20] J.E.S.SocolarandS.A.Kauffman,Phys.Rev.Let.90,
ering any possible sequence of updates.
068702 (2003).
The biological implications of these findings are clear [21] U. Bastolla and G. Parisi, Phys.D 115, 219 (1998).
and have been pointed out in [27]. Since biological net- [22] U. Bastolla and G. Parisi, Phys.D 115, 203 (1998).
worksdonothaveacompletely synchronousupdate, the [23] B.SamuelssonandC.Troein,Phys.Rev.Let.90,098701
numberofattractorsshouldbederivedfrommodelswith (2003).
[24] B. Drossel, T. Mihaljev, and F. Greil (2004), cond-
asynchronous update. In [9] it is found numerically that
mat/0410579.
thenumberofstableattractorsincreasessublinearly. At- [25] R. J. Glauber, J. Math. Phys. 4, 294 (1963).
tractors are called stable if they do not change when a [26] A. J. Bray, J. Phys. A 23, L67 (1990).
perturbation is added to a synchronous updating rule. [27] K. Klemm and S. Bornholdt (2004), q-bio/0409022.
The present paper extends this finding by showing that [28] H.FlyvbjergandN.J.Kjær,J.Phys.A21,1695(1988).
the number of attractors in critical asynchronous Kauff- [29] V. Kaufman and B. Drossel (2004), cond-mat/0410546.