Kauffman networks with threshold functions

0.13 MB·English

by  Florian Greil

Kauffman networks with threshold functions
Florian Greil and Barbara Drossel
Institut fu¨r Festk¨orperphysik, TU Darmstadt, Hochschulstraße 6, 64289 Darmstadt, Germany
(Dated: February 6, 2008)
We investigate Threshold Random Boolean Networks with K = 2 inputs per node, which are
equivalent to Kauffman networks, with only part of the canalyzing functions as update functions.
7
According to the simplest consideration these networks should be critical but it turns out that
0
they show a rich variety of behaviors, including periodic and chaotic oscillations. The results are
0
supported byanalytical calculations and computer simulations.
2

n
a RandomBooleannetworks(RBN) were introduced by neural networks,but also in models for genetic networks
J
S. Kauffman in 1969 [1, 2] to model the dynamics of ge- [7, 9, 10]. These functions represent four of the 12 can-

9 netic and metabolic networks [3], but they are also used alyzing update functions of Kauffman networks. Cana-

in a social and economic context [4, 5], and for neural lyzing functions are those non-frozen functions where at
]
n networks. Although Boolean models represent a strong least one value of a given input can fix the output of a
n simplificationofthe farmorecomplexreality,there exist node, irrespective of the value of the second input. All
-
s several examples where the modelling of a genetic net- nodes are updated in parallel according to the rule
i work by Boolean variables captures correctly the essen-
d
t. tialdynamicsofthe system[6, 7, 8]. For this reason,the σi(t+1)=fi({σj(t)})≡fi(σi1(t),σi2(t)). (2)
a study of RBNs remains an important step on the way
m towards understanding real networks. Node i depends on the nodes j, namely on node i1 and
- A random Boolean network is a directed graph with i2.
d randomly chosen links between N binary nodes. We de- Theconfigurationofthesystem~σ σ1,...,σN per-
n ≡{ }
note the state of a node with σ = 1 (we call +1 “on” forms a trajectory in configuration space. As the state
o i ±
c and 1 “off”), and the number of inputs per node with space is finite and the dynamics is discrete, some states
[ K. E−ach node i is assigned at random an update func- will occur again. If a cycle in state space has a set of

tion f . In this paper, we focus on the case K = 2 and transient states leading to it, it is called an attractor.
1 i
v on threshold functions Kauffman classified the dynamics of RBNs according
6 towhetheritischaotic orfrozen orcritical (asdescribed
7 in the review [11]). In a chaotic network, a perturbation
1 fi =sign σjcij sign(si) (1) at one node propagates on an average to more than one
≡
1 Xj node, leading to long attractors. In a frozen network,
0  
suchaperturbationpropagatesonanaveragetolessthan
7 wherethesumistakenoverthetwoinputnodesfornode
0 one node, and in the thermodynamic limit N only
/ i, and cij = −1 for inhibitory connections and cij = 1 a finite number of nodes are not constant after→a∞certain
t for excitatory connections. (This version of RBN was
a transient time. Critical networks are at the boundary
m used for instance in [9].) A connectionis excitatorywith
between these two types of behavior, with a perturba-
- probability p+ and inhibitory with probability 1−p+, tion of one node propagating on an averageto one other
d leading to the following table (Tab. I).
node. Therefore the difference between two almost iden-
n
o ticalinitialstatesincreaseslikeapowerlawintime. The
c Input f7 f11 f13 f14 number of nodes that are not frozen on all attractorsin-
v: (↓,↓) ↑ ↑ ↑ ↓ creases in a critical network as a power law N2/3 of
i (↓,↑) ↑ ↑ ↓ ↑ the system size N, as was found numerically i∼n [12] and
X
(↑,↓) ↑ ↓ ↑ ↑ analytically in [13, 14].
ar (↑,↑) ↓ ↑ ↑ ↑ It is the aim of this paper to show that this classifica-
probability (1−p+)2 p+(1−p+) p+(1−p+) p2+ tionbreaksdownforthesimpleclassofRBNsconsidered
here. Wefindamuchricherdynamicalbehaviorwithnot
TABLE I: The four possible update functions for the model onlyafrozenandachaoticphase,butalsowithtwotypes
used in this paper. The input configuration is given in the
of oscillating phases and several critical points.
firstcolumn,with↑denotingσi =1and↓denotingσi =−1.
Let us first apply the criticality condition in its sim-
The last row gives the probability for each function, the top
rowgivesthenameofthefunctionaccordingtotheKauffman plestversion: Forallfourupdatefunctions,theprobabil-
model. ity that the output changes if one input spin is flipped,
is 1/2. Since each node is on an average the input to
Inagreementwithotherauthors,wedefinesign(0)=1. two other nodes, a perturbation at one node propagates
Threshold functions are used not only in the context of on an average to one other node, and we should expect
2
the model to be critical. This is in agreement with the therefore conclude that the model is in the frozen phase
finding that K =2-RBNs that contain only canalyzing for p+ >1/2,that it is criticalfor p+ =1/2,andchaotic
update functions (but all of them with the same prob- for p+ <1/2.
ability) are critical [15]. However, this simple argument Thesameresultisobtainedbycalculatingthestation-
is based on the assumption that all four possible input ary value of the Hamming distance between two identi-
configurationsoccur equallyoften, whichmay be true at cal network realizations. The Hamming distance is the
the beginning of a simulation run, but may be wrong al- fraction of nodes for two configurations ~σ,~σ˜ that have
ready after one timestep. For this reason, Moreira and different values: D = (4N)−1 N (σ σ˜ )2. If we de-
i=1 i − i
Amaral [16] argued that the calculation should be per- notewithπ2 theprobabilitythPatanodechangesitsstate
formed such that the input configurations are weighted when both inputs are flipped, the time evolution of D is
with their frequencies in the stationary state. given by
Let us therefore next apply the rule given by Moreira
and Amaral and let us determine for what values of p+ Dt+1 =2Dt(1−Dt)π1+Dt2π2. (6)
it predicts that the model is frozen, critical, or chaotic.
We will then see later that the result is still not correct π2 in the stationary state is obtained by summing all 8
combinations leading to s 2 . It can be written as
for all parameter values. i ∈{± }
Wedenotewithbttheproportionofnodesinstateσi = π2 =1 2b(1 2p+)2+2b2(1 2p+)2 2p++2p2+ (7)
+1 at time t. In the thermodynamic limit, it changes − − − −
deterministically according to and is 1/2 for p+ = 1/2 and 1 for p+ = 1. If Dt is very
bt+1 = 1− b2t(1−p+)2+(1−bt)2p2++ small, we have
2bt((cid:2)1−bt)p+(1−p+) . (3) Dt+1 ≃2Dtπ1, (8)
(cid:3)
The expression in the square brackets is the probability
which allows for the growth of a small perturbation if
that an input combination leads to s = 2. In the
i − 2π1 >1orp+ <1/2,inagreementwithourresultabove.
stationary state, we have bt+1 =bt =b with The transition from a stationary value D = 0 to a sta-
tionary value D >0 occurs at the same point.
4p2+ 2p+ 1 5 12p++8p2+
b(p+)= − −2(1±q2p+−)2 . (4) 1
−
The sign in the numerator has to be chosen such that
b [0,1], therefore only the positive branch remains, see
∈
Fig. 1. For p+ =1/2,the denominatorvanishes, andthe
stationary solution of Eq. (3) is b = 3/4. For p+ = 1, 0.5
we have b = 0 and b = 1, with the first solution being
obviously unstable as it is destroyed by one node in the
sotfaTttheheeσdimy=enaa+nm1ni.cusTm.hbFeeorsreocpfo+nno=dds0eos,luwttoeiowhnhaiivscehabast=pabe(rl−etu1fir+xbead√tip5o)no/ina2tt. 0 ppDb12((((pppp++++))))
0 0.5 p+ 1
one node propagates is in the stationary state given by
2π1,withπ1beingtheprobabilitythatanodechangesits FIG. 1: The functions b(p+),π1(p+),π2(p+) and the station-
state when one input is flipped. We obtain it by adding ary value D(p+) vs. p+.
up the probabilities for those input configurationswhich
allowatransitionbetweenanoutput+1and−1andvice Fig.2showsDt fordifferentvaluesofp+ andforgiven
versa. This is true for half of the input configurations initial conditions. One can see that D approaches 0
t
leadingto si =0(the first4 termsin the followingequa- for large times if p+ > 0.5. Furthermore, one can see
tion) and for all input configurations for which si = −2 that Dt oscillates with period 2 for the smaller values of
(the last 4 terms): p+. Thisoscillationisanindicationthatthedynamicsin
the “chaotic” phase has some structure, which shall be
π1 = (1 p+)(1 b)(1 p+)b+p+bp+(1 b)+
− − − − investigated in the following.
p+b(1 p+)b+(1 p+)(1 b)p+(1 b)+ Let us therefore have a closer look at the supposedly
− − − −
(1 p+)b(1 p+)b+p+b(1 p+)(1 b)+ chaotic phase p+ < 1/2. We will see that the dynamics
− − − −
(1 p+)(1 b)p+b+p+(1 b)p+(1 b) is not chaotic at all for sufficiently small p+. The con-
− − − − siderations that have lead to our simple phase diagram
π1 = b+p+ 2bp+ (5)
⇒ − are flawed. The reason is that we have assumed that b
Forp+ =1/2,weobtainπ1 =1/2,for p+ =1,we obtain becomes for large times stationary for all p+. In order
π1 = 0. For p+ = 0, we obtain π1 = b 0.618. We to see that this need not be the case, let us first look at
≃
3
if
0.6
p =0.05
D(t) pp+++==00..0079 p+ >(3−√7)/4≡pcb ≈0.0886. (10)
p =0.11
+
p =0.40 Only above this value does the system have a stationary
+
0.4 p+=0.50 state with constantproportionsofnodes being “on”and
p =0.60
+
“off”.
We finally investigate in more detail the region p+ <
p , where the proportion of “on” and “off” nodes oscil-
cb
0.2 lates with period 2. For p+ = 0, every node oscillates
with period 2, and we have a global attractor of period
2. This need not necessarily be the case if b oscillates
withperiod2. Theattractorcouldbemuchlarger,while
0 the proportion of off and on nodes oscillates still with
0 20 40 t 60 period two. In order to determine for which parameters
an attractor with period 2 is stable, we performed again
FIG. 2: The time evolution of the Hamming distance D for
alinear stability analysis,but nowfor twotime steps to-
differentvaluesofp+whenbisnotstationary. Thecurvesare gether. We assume that the system is on an attractor
calculated according to Eqs. (3), (6) starting from D0 = 0.2
of length 2. Let there be every even time step a propor-
and b0 =0.5.
tion x of “on”-nodes and every odd step a proportion y.
The time evolutionis givenby Eq.(3), but now we com-
1
b bine two consecutive time steps. We flip one node and
t+1
look how the Hamming distance grows in comparisonto
0.8
the undisturbed systemafter two time steps. The condi-
tion that information spreading is critical is in the cycle
b
0.6 with period 2
1
0.4 π1(x) π1(y)= (11)
· 4
0.2 Combining Eq. (11) with the time evolution of x and y
we obtain three equations
FIG. 3: The m 0a 0p bt v 0s..2bt+ 10.=4 1− 0.b62t for 0p.8+ =bt0 1. The fixed xy ==11−[[xy22((11−pp++))22++((11−yx))22pp2+2+++22yx((11−yx))pp++((11−pp++))]]
point b is unstable as depicted by a sample trajectory. − − − − −
1
=(x+p+ 2xp+)(y+p+ 2yp+)
4 − −
the situation where p+ =0: We then have bt+1 =1−b2t. This system can be solved numerically and we obtain a
Thisisaone-dimensionalmapshowninFig.3. Thefixed critical value pcn = 0.0657. For p+ below this value a
point is unstable! Instead of having a stationary point perturbation at one node will die out and all nodes will
with a constant proportion of nodes in the two states, again blink with period 2. Above this value, attractors
the system oscillates between a configuration where all must be longer than 2.
nodes are switched on and a state where all nodes are We checkedtheseanalyticalpredictionsbyperforming
switched off. This is not chaotic dynamics at all, but computer simulations. In order to identify the transi-
very stable dynamics. In order to determine the range tion at p , we measured the median attractor length.
cn
ofp+values,for whichthe fixed pointvalue ofb is unsta- As shown in Fig. 4, we find with increasing network size
ble, we performed a linear stability analysis. The ansatz an increasingly sharptransition. Below the transitionat
bt+1(b+δbt) = b+δbt+1 leads in linear order in δb and pcn, the proportion of attractors of length 2 converges
for p+ <1/2 to to some nonzero value with increasing system size, in-
dicating that cycles of length 2 are stable. Above the
2 transition, the median attractor length increases more
δbt+1 = 2 bt(1 2p+) +(1 2p+)p+ δbt
− (cid:16) − − (cid:17) and more rapidly with increasing system size, indicating
= 1 5 12p++8p2+ δbt =:M δbt (9) a diverging median. Another finding is that attractors
(cid:18) −q − (cid:19) · become again shorter as the critical point p+ = 1/2 is
approached.
InthelaststepweusedEq.(4). Thefixedpointisstable Wealsoevaluatedthefrequencyofphasejumpsinb(t)
if the real part of M is smaller than 1, which is the case on the attractors. Each phase jump is a deviation from
4
1000 sense defined by Kauffman. For p+ > 1/2, the network
is in the frozen phase.
100 The lesson to be learned from this study is that the
dynamical behavior of Boolean networks can be much
10
richer than expected from simple considerations. There
median attractor length vs. p
+
isindeednoreasonwhysomemodelshouldnotalsoshow
1
0 0.1 0.2 0.3 0.4 0.5 globaloscillationswithhigherperiodsorperioddoubling
N=802 cascades in the temporal behavior of b(t). Real genetic
0.15 N=602
N=402 networkscanbe expected to haveanevenricherdynam-
N=202
0.10 ical behavior. If the simple classification into “frozen”,
“critical”and“chaotic”networksfailsalreadyintheran-
0.05
dom model presented in this paper, it will be even less
median jump frequency vs. p
+ suitable for real genetic networks, which have attractors
0.00
0 0.1 0.2 with very specific properties related to the function of
thenetwork. Amoresophisticatedwayofdescribingand
FIG.4: Numericalverification forthetransition pcn andpcb, classifying the dynamical behavior of Boolean networks
both marked by vertical arrows. The upper panel shows the
is therefore required.
median attractor length, the lower panel the median jump
frequency on each attractor candidate in dependence of p+. Thisworkwassupportedbythe Deutsche Forschungsge-
Each data point corresponds to5000 sample networks ofsize meinschaft (DFG) under Contract No. Dr200/4-1.
N with fixed p+ and two initial conditions per realization.
Thetimeevolutionislimitedto5000computationalstepsfor
both the transient and the attractor length. The median is
thereforeamoresensibleobservablethanthemeanvalue: Ifa
attractor candidatecannot beverified becausethetimelimit [1] S. A. Kauffman,J. Theo. Bio. 22, 437 (1969).
was already reached we can still calculate themedian. [2] S. Kauffman, Nature224, 177 (1969).
[3] S.Kauffman,C.Peterson,B.Samuelsson,andC.Troein,
inProceedings ofthe National Academy of Sciences USA
an oscillation with period 2. The result is shown in the (2003), no. 25 in 100, pp.14796–14799.
lower part of Fig. 4. With increasing system size, there [4] J. M. Alexander, in Philosophy of Sci-
is an increasingly sharp transition at p between zero ence Assoc. 18th Biennial Mtg (2002), URL
cb
http://philsci-archive.pitt.edu/archive/00001049/.
phase jumps and a finite proportion of phase jumps.
[5] M. Paczuski, K. E. Bassler, and A. Corral, Phys. Rev.
Wesummarizethedifferenttypesofdynamicalbehav- Lett. 84, 3185 (2000).
ior in the following diagram, Fig. 5. [6] S. Bornholdt, Science 310, 449 (2005).
[7] F. Li, T. Long, Y. Lu, Q. Ouyang, and C. Tang, Proc.
nonfrozen frozen Nat. Acad. Sci. 101, 4781 (2004).
b has period 2 b is constant in time [8] R. Albert and H. G. Othmer, Journal of Theoretical Bi-
ology 223, 1 (2003).
all nodeshave nodes oscillate differently or not at all [9] T. Rohlf and S. Bornholdt, Physica A 310, 245 (2002).
period 2
[10] S. Bornholdt and K. Sneppen, Proc. R. Soc. London B
p+ pcn pcb 0.5 267, 2281 (2000).
[11] M. Aldana-Gonzalez, S. Coppersmith, and L. P.
FIG.5: Overviewofthedynamicbehaviorofthemodelwith Kadanoff, Perspectives and Problems in Nonlinear Sci-
K =2 in dependenceof theparameter p+. ence pp.23–89 (2003).
[12] J. E. S. Socolar and S. A. Kauffman, Phys. Rev. Lett.
Toconclude, wehaveshownthatthe simplestThresh- 90, 068702 (2003).
oldRandomBooleanNetworkshowsthreedifferenttypes [13] B. Samuelsson and C. Troein, Phys. Rev. Lett. 90,
098701 (2003).
of phase transitions and not just the generally expected
[14] V. Kaufman, T. Mihaljev, and B. Drossel, Phys. Rev. E
transitionbetween a frozenand a chaotic phase. For pa- 72, 046124 (2005).
rameter values p+ < pcn, all nodes oscillate stably with [15] U. Paul, V. Kaufman, and B. Drossel, Phys. Rev. E 73,
period two. For pcn <p+ <pcb, the fractionof on-nodes 026118 (2006).
oscillate with period two, but attractors are longer. For [16] A.A.MoreiraandL.A.N.Amaral,Phys.Rev.Lett.94,
pcb <p+ <1/2,the dynamicalbehavioris chaotic inthe 218702 (2005).