Table Of Content

Kauffman networks with threshold functions Florian Greil and Barbara Drossel Institut fu¨r Festkörperphysik, 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 stationary 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 condition 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, indicating 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).