Critical Kauffman networks under deterministic asynchronous update

We investigate the influence of a deterministic but non-synchronous update on Random Boolean Networks, with a focus on critical networks. Knowing that ``relevant components'' determine the number and length of attractors, we focus on such relevant components and calculate how the length and number of attractors on these components are modified by delays at one or more nodes. The main findings are that attractors decrease in number when there are more delays, and that periods may become very long when delays are not integer multiples of the basic update step.


I. INTRODUCTION
Random Boolean networks are widely used as models for complex systems that consist of interconnected units that influence each other and have two states ("on" and "off"). Kauffman [1,2] used them as a simple model for gene regulation, but they can also be applied in a social and economic context [3,4], for neural networks and protein networks [5]. Recently it was shown that the idealized representation of genes as Boolean units is sufficient to understand the essential dynamics of certain real gene regulation networks. In these cases, one need not include rate equations for the concentrations of the molecules involved in the processes in order to identify the sequence of steps taken by such a system [6,7,8].
A random Boolean network (RBN) is a directed graph consisting of N nodes and kN links between them. The nodes have values σ i ∈ {0, 1} and receive input from k other nodes, which are chosen at random when the network is constructed. Each node i has an update function f i , which assigns to each of the 2 k states of its k input nodes an output 1 or 0. The update function of each node is chosen at random among all 2 2 k possible update function. All nodes are usually updated in parallel according to the rule The assignment of connections and functions to each node remains fixed throughout the whole time evolution, the model is therefore referred to as "quenched" [9]. The dynamics follows a trajectory σ(t) ≡ {σ 1 (t), . . . , σ N (t)} in configuration space which eventually leads to a periodically repeating sequence of configurations, called cycle, as the state space is finite and the dynamics is discrete. Such an cycle is called attractor if there is a set of transient states leading to it; these constitute the basin of attraction.
The dynamics can be classified [9] according to the way information spreads through the network.
1. In the frozen phase, all nodes apart from a small number (that remains finite in the limit of infinite system size) assume a constant value after a transient time. If in the stationary state the value of one node is changed, this perturbation propagates during one time step on average to less than one other node.
2. In the chaotic phase initially similar configurations diverge exponentially. Attractors are usually long, and a non-vanishing proportion of all nodes keep changing their state even after long times.
Critical networks are at the boundary between the frozen and chaotic phase [9], and neighboring configurations diverge only algebraically with time. Whether a network is frozen, critical or chaotic depends on the value of k and on the probabilities assigned to the different types of update functions. When all update functions are chosen with the same probability, networks with k < 2 are frozen, networks with k = 2 are critical, and networks with k > 2 are chaotic [10]. The usual synchronous way of updating is not very realistic [11] as natural systems are rarely controlled by an external clock. It is known that properties of attractors for synchronous dynamics can differ from those for asynchronous dynamics. For instance, for cellular automata part of the self-organization is closely tied to the synchronous updating [12]. For RBNs, the dynamics changes considerably when other updating schemes are chosen [13,14]. While in critical Boolean networks with parallel update the number of attractors increases superpolynomially with the network size [15], it becomes a power law for asynchronous stochastic update [16].
In this paper we will consider deterministic updating schemes that are not fully synchronous. This means that some nodes are less frequently updated than others. Such node-based delays can be motivated biologically: The expression of genes is not an instantaneous process, the transcription of DNA and transport of enzymes may take from milliseconds up to a few seconds. To each node i we assign a delay time τ i . The value σ i of node i is updated in time intervals τ i , and each node must be assigned an initial "phase" ϕ i < τ i (i.e., the time until the first update). The model is referred to as a Deterministic Random Boolean Network (DRBN). The system is deterministic as the succession of network states σ(t) is entirely defined by the initial condition σ(0) and the initial phases {ϕ i }. The case of parallel update, the so-called Classical RBN (CRBN), is a special DRBN with all τ i ≡ 1, ϕ i ≡ 0. The size of the state space Ω changes from |Ω| CRBN = 2 N to |Ω| DRBN = i 2 τi , when all τ i are integers.
The outline of the rest of this paper is the following: First, we review the concept of relevant components (Sec. II). This shows that the most frequent relevant components are simple loops, and less frequent are collections of loops with additional links within and between them. In the subsequent sections, we therefore study simple loops with one delayed node, simple loops with several delayed nodes, two loops with a cross-link and one delayed node, and a loop with one additional link and a delayed node (Sec. III-Sec. VI). In the conclusion, we discuss the consequences of our findings for the entire network, which is composed of several relevant components.
1. The state of frozen nodes becomes constant after some time. Interestingly, the nodes that become frozen are the same nodes most of the time, and they constitute the frozen core. The frozen core is identified by starting from nodes with constant functions and by iteratively identifying nodes that become frozen due to frozen inputs. Networks without frozen functions can also develop a frozen core [23], however, the mechanism is different. In critical networks, the frozen core comprises all but a proportion ∼ N −1/3 of nodes.
2. Relevant nodes are non-frozen and have different dynamics on different attractors, and they influence at least one relevant node [24]. They determine the attractors, and the number of relevant nodes scales in critical networks as N 1/3 [21].
3. Finally, there are the irrelevant nodes which are not frozen but are slaved by the relevant nodes.
The non-frozen nodes of critical networks essentially form a k = 1 network. This means that typically all but one input of a nonfrozen node are frozen. In critical networks, relevant nodes are arranged in O(ln N ) components, most of which are simple loops. Typically, there is only one component that is not a simple loop but has µ nodes with two relevant inputs [25]. There exist two possible components with µ = 1, and then more complicated components with µ > 1. If µ = 1, there are either two loops with an cross-link or a loop with an additional link, see Fig. 2. Two loops with a cross-link occur twice as often as a loop with an additional link [25]. Since the relevant components determine the long-term dynamics, we will in the following study their properties.

III. LOOPS WITH ONE DELAYED NODE
The simplest relevant component is a loop consisting of nodes with k = 1 incoming edges and Boolean functions which either copy (⊕) or invert (⊖) the previous node's value. A constant function within a loop will freeze the whole loop, and can therefore not occur in relevant loops. A loop with n inversions can be mapped bijectively onto one with (n − 2) inversions by replacing two ⊖ with two ⊕ and by inverting the values of all nodes between these two couplings. It is therefore sufficient to distinguish loops with an even or an odd number of inversions, and we call them "even" and "odd" loops respectively. When discussing even loops, we consider loops with only ⊕ functions. To odd loops we assign one ⊖ function, the connection to a node with this function is called "twisted". If no node is delayed, an even loop with a prime number N ∈ P of nodes returns to its initial state after N time steps. If N is not prime, shorter cycle lengths exist. Irrespective of the updating scheme, there are two fixed points for an even loop, namely σ ∈ { 0, 1}. An odd loop with N ∈ P nodes returns to its initial state after 2N time steps. It has no fixed points, For synchronous update, the shortest attractor of an odd loop has period 2, with alternating 1's and 0's.
Let us now introduce a delay in the loop. Let node 1 be the delayed node, i.e., we choose a value 1 < τ 1 ≡ τ * ∈ N, while τ i = 1 for i > 1. This means that one gene needs much longer to be expressed than all the others. Since node 1 remains at the same value for τ * time steps, node 2 receives τ * times the same input, leading to blocks of size τ * travelling around the loop. When the head of a block arrives at node N , node N will have the value of this block for τ * time steps, and during one of these steps node 1 will be updated. In the following, we will consider even and odd loops separately.
Even loops. An even loop with one delayed node has two fixed points, just as the loop with no delayed node. The other attractors are characterized by blocks of length τ * traveling around the loop. Assume node 1 is updated at time 0. A block with the value of node 1 will start travelling around the loop, and the head of the block will arrive at node N at time N − 1. The next update of node 1 will be at time T = ⌊(N − 1 + τ * )/τ * ⌋ · τ * , where the Gaussian brackets ⌊x⌋ denote the largest integer less or equal x. The value of node 1 becomes the same as at time 0, and the same block travels around the loop again. The same consideration can be made for all starting times that are multiples of τ * , leading to the result that the state of the loop is repeated every T time steps.
After the transient time N −1, the loop has reached an attractor that contains ξ = T /τ * blocks, each of which either has only values 0 or values 1, and with the blocks that contains node 1 and N being shorter than the other ones, if N is not a multiple of τ * . Every attractor corresponds to a pattern of blocks travelling around the loop. If the number of blocks is a prime number, ξ ∈ P, the length A ⊕ of the attractors is identical to T , For τ * ≥ N there is only ξ = 1 block containing all nodes, and the fixed points are the only attractors. For ξ / ∈ P, the pattern of blocks can have a period that is a divisor of ξ, in which case the attractor length is shorter.
The number of different attractors, ν ⊕ , can be calculated from the number of different patterns, 2 ξ . Including the 2 fixed points this leads to If the number of blocks is not prime, ξ / ∈ P, the number of attractors increases as the length of some attractors is shorter.
Odd loops. Without loss of generality we assign the twisted edge of the odd loop to be in front of the delayed node. As in the synchronous case, there are no fixed point attractors for odd loops. Let again node 1 be delayed and updated at time 0. At time T , node 1 will have the opposite state as the original one. After time 2T , node 1 returns to its original state, which implies that the loop returns to its original state after 2T time steps, if it is on an attractor. If all nodes are identical initially, a single domain wall travels around the loop, and after T time steps all nodes are again identical, but with the opposite state. The shortest attractor has a period 2τ * , and it has alternating blocks. If the number of blocks ξ = T /τ * is a prime number, all other attractors have the period 2T , and the number of different attractors is If ξ / ∈ P, the number of attractors increases as the length of some attractors is shorter. If τ * ≥ N , there is only one attractor with period 2τ * .
Non-integer delays. Let us now consider the case that τ * is not an integer, but a rational or a real number. Real numbers can be approximated by a series of rational numbers, and we therefore consider the case τ * = r/s with two incommensurate integers r and s, with r > s. During r time steps, s blocks emerge from node 1, part of them of length r s , part of them of length r s +1. When the first block reaches node 1, the same sequence of blocks will emerge again only if node 1 is in the same phase at its next update as it was at its first update. Otherwise, the pattern of blocks will be changed at each circulation around the loop, until it starts repeating again after s (or 2s for an odd loop) circulations (or a divisor of it). It follows that for irrational values of τ * the dynamics never become exactly periodic but are quasiperiodic. Of course, for values τ * > N , the only attractor is a fixed point (for an even loop) or a state with only one domain wall (for an odd loop).

IV. LOOPS WITH MULTIPLE DELAYED NODES
Next, we consider loops with multiple delayed nodes and integer delay times τ i ∈ N >0 . Loops with rational values τ i = r i /s i can be mapped on those with integer values of τ i by measuring time in units of the inverse of the least common multiple of all s i . In the following we will first look at two special cases before we focus on the general case where updates may occur in any order.
Sequential update. We choose τ i ≡ N and update the nodes in the order in which they occur on the loop, i.e., ϕ i = i − 1 (connection-wise (cw) update) or ϕ i = N − i (counter-connectionwise (cc) update).
For cw-update, all nodes of an even loop have the same value, after every node has been updated once. Thus, we have two fixed point attractors consisting of the two homogeneous configurations, σ ∈ { 0, 1}. For an odd loop, the attractor has a single domain wall that travels around the loop, and the period of the attractor is 2N .
For cc-update, node j is updated before node j − 1. Therefore, N update steps give the same result as 1 update step in the case of parallel update, and the results of Sec. III can be taken over.
Same delay, different phases. We now choose again τ i ≡ N , but we update the nodes in any order, i.e., the values of ϕ i are some permutation of the numbers 0 to N − 1. There are two classes of nodes, according to whether a node is updated before or after its predecessor. Nodes that are updated after their predecessor have after N time steps the same state as their predecessor. Such nodes and their predecessor are therefore part of the same effective node. The number N * of effective nodes is identical to the number of nodes that are updated after their predecessor. Let us give an example, for instance N = 6 and the updating order If we specify only whether a node is updated before (b) or after (a) its predecessor, this can be written as {a, b, a, b, a, a} leading to N * = 2. Once we have identified the effective nodes, we can map N time steps on such a loop with sequential update on one time step on a loop of size N * with parallel update. All results concerning attractor numbers and lengths obtained for loops with parallel update can then be transferred to loops with sequential update.
Different delay times. We now consider the general case where the delays τ i and the phases ϕ i can take any integer value. In order to determine whether the initial state of a given node influences the attractor, we proceed in the following way: We fix the state of this node, let us say, to 1, and we evaluate to which nodes this 1 propagates with time. In order to make sure that later on all 1s on the loop will be due to this initial 1, we set all other nodes to 0 and choose an even loop. When the chosen node is updated before its successor, the 1 is lost, and the initial state of this node does not affect the attractor. If the node is updated after its successor, the 1 has moved to the successor and is not yet lost. Next, we check whether the successor is updated before the 1 that is now there can propagate further. During the course of time, the 1 may spread to become a block of larger size, which continues to change its size with time. Now we consider the loop at times which are a multiple of τ * = lcm(τ i ). At these times, the phases of all nodes are the same as at the beginning. We wait until either the original node has again state 1 or until all 1s are gone. In the first case, the 1 will survive forever, and the initial state of the chosen node will consequently affect the attractors. In the second case, the chosen node does not affect the attractors.
If we repeat this procedure for every node, we will know the initial state of how many nodes N * will affect the attractors. We only consider these "relevant" nodes from now on, and we consider them only at times that are multiples of τ * . Let m be the number of relevant nodes through which each block moves during τ * time steps. If m and N * have no common divisor, we order the N * nodes in the sequence in which they are visited if the system is only considered at times that are multiples of τ * . Then we have mapped the task of finding the attractor number and length on a loop of size N with any rational delay times of the task of finding the attractor number and length in a loop of size N * . If m and N * have a common divisor l, we can map our system on l loops of length N * /l with parallel update and thus find the number and lengths of attractors.
For irrational delay times, the mapping on integer delays cannot be performed. Nevertheless, one can determine whether a 1 at a given node will eventually become a block that is so large that it will never vanish. If all delays have irrational ratios, there will eventually come a moment where all nodes are updated connection-wise,  and from then on there is only one block left.

V. TWO LOOPS WITH A CROSS-LINK
Now we consider a complex component consisting of a loop with N 1 nodes connected to a loop with N 2 nodes, see Fig. 2 (a). The node Σ is the one with two inputs, and its input nodes are labelled G 1 and G 2 . The first loop is either odd or even, the second loop can without loss of generality be chosen such that it has only ⊕-couplings, except at Σ. We insert no nodes between G 1 and Σ, as a system with m nodes on the cross-link can be mapped on a system with a direct link by connecting node number m (counted clockwise from G 1 ) directly with node Σ. We consider only the nontrivial cases where the coupling function of node Σ is a function that responds to both of its inputs. If we take into account that certain of these functions can be mapped onto each other by inverting the states of all nodes or by inverting the states on the first loop, we end up with three functions that are truly different. They are shown in Tab. I.
In the following, the results by Kaufman and Drossel [25] for components under synchronous update will be generalized to components with one delayed node.
Delayed node on first loop. If the first loop has a delay τ * , the value of G 1 can change only at times that are multiples of τ * , and the pattern of change is repeated after the attractor period of the first loop. If the first loop is on a fixed point, the second loop can be considered as an independent loop with a function at Σ, which depends on the fixed point value of the first loop. We therefore consider here the case that the first loop is not on a fixed point but provides a periodic input of period p 1 to Σ, with blocks of size τ * of identical bits. The second loop then behaves like an independent loop where the Boolean function at Σ changes after τ * steps, and where the pattern of changes is repeated periodically with period p 1 .
If f Σ = f rev , the second loop switches between truth and negation, for the (in)homogeneous canalizing function f ch (f ci ) the second loop changes between truth (negation) and constant value σ Σ = 1. The attractor period is at most 2p 1 N 2 . More detailed results for possible attractors of such a system under synchronous updating can be found in [25], the only difference being that p 1 is now related in a different way with N 1 . An interesting finding is that for the homogeneous canalizing function, the second loop becomes frozen on the value 1 if p 1 and N 2 have no common divisor. Furthermore, for the inhomogeneous canalyzing function, the first loop enslaves the second loop and imposes its period on it, if p 1 and N 2 have no common divisor.
Delayed node on second loop. We now proceed to the case where a single delayed node is on the second loop. The first loop behaves in the same way as for simple (synchronously updated) CRBN-loops. We focus on the sequence of values of the delayed node. We assign the delay to the node Σ: A system with the delayed node m nodes after Σ can be transformed into a system with the delay at Σ by rotating the first loop m nodes counterclockwise.
Node Σ responds to the input from G 1 only every τ * time steps. Let us denote with p 1 the period of the sequence of values of G 1 every τ * time steps, which is the period of the input sequence to Σ generated by the first loop at those times where Σ is updated. Let ξ = ⌊(N 2 − 1 + τ * )/τ * ⌋ denote the number of blocks of the second loop.
All results for the attractors on two loops with a crosslink and no delay can now taken over by replacing N 2 with ξ, by replacing nodes with blocks, and by taking τ * as the time unit. In particular, for a reversible function f Σ , the largest period is p 1 ξτ * . A homogeneous canalyzing function f Σ freezes the second loop on the value 1 if p 1 and ξ have no common divisor. Furthermore, for an inhomogeneous canalyzing function, the first loop enslaves the second loop and imposes its period (times τ * ) on it if p 1 and ξ have no common divisor.
General case. Now we consider the case of multiple (integer) delays on both loops. Let the first loop have a period p 1 and the second loop (if even and decoupled from the first loop) a period p 2 . The general system of two interconnected loops without any delay has been studied in [25]. There, it was shown that the attractor length lies between p 1 and 2p 1 N 2 /g, where g = lcm(p 1 , N 2 ). We can conclude that now the attractor length lies between p 1 and 2p 1 p 2 /g, where g is the greatest common divisor of p 1 and p 2 .
For a homogeneous canalyzing function f Σ = f ch , the second loop is frozen if p 1 and p 2 are incommensurable. The longest attractors occur for reversible functions, f Σ = f rev .

VI. LOOP WITH ONE ADDITIONAL LINK
The other complex component with one node with two inputs is a loop with N = L + M + 2 nodes and one additional link. We call again the node with two inputs Σ, and its inputs G 1 and G 2 , see Fig. 2 (b). The link from G 2 to Σ can be treated as a direct link: A system with n < L nodes in the additional link can be mapped onto a system with a direct link by connecting node (M + 1 + n) to Σ (if we neglect delays). We consider five update functions at Σ, compare Tab. II (we now use the common decimal representation as identifiers for the functions). The other canalizing or reversible functions yield the same result, one only has to invert the output values of the truth table. Without loss of generality, all other functions in the loop are copy functions.
Let there be one delayed node in the component. We distinguish two cases according to the position of the node with update period τ * > 1: 1. The delayed node lies on the first M + 1 nodes (including G 1 ). Without loss of generality the delay can be shifted to node Σ.
2. The delayed node is in the chain of nodes between G 1 and Σ and can then be shifted to G 2 .
In the first case, the component can be reduced to a network of effective nodes by looking at the network only every τ * time steps. Each effective node corresponds to a block of τ * nodes which are at the same state for t mod τ * = 0. The results for the synchronous case (as studied in [25], Sec. 4) hold for the effective vari-ablesÑ,M , ⌈x⌉ denotes the smallest integer greater or equal to x: In the following we will study the second case, where G 2 is the delayed node. Canalyzing function at Σ. If there is a canalyzing function at Σ, there exist at least two types of attractors. The first type of attractors is obtained by requiring that G 2 does never have its canalyzing value at the moment when Σ is updated. (The canalyzing value is 0 for f 13 and f 2 and 1 for f 1 and f 14 .) In this case the loop consisting of the M + 2 nodes from Σ to G 1 is an even loop for f 13 and f 14 and an odd loop for f 1 and f 2 . Just before node Σ is updated, the state of node G 2 and the state of all nodes that will in nτ * time steps determine the state of G 2 must have a value such that Σ never has its canalyzing value (n is any positive integer). For functions f 13 and f 14 , this condition fixes the    entire component at the same value if M + 2 and τ * have no common divisor. If their greatest common divisor g is larger than 1, only the value of every g th node on the loop of length M + 2 is fixed by this condition. The number of attractors of the first type is therefore that of an even loop with (M + 2)(g − 1) nodes and with no delays. For functions f 1 and f 2 , the condition that G 2 does never have its canalyzing value can only be satisfied if τ * /g is even. The number of attractors of the first type is then that of an odd loop with (M + 2)(g − 1) nodes and with no delays. Compared to a component with no delays, this new type of attractors increases the mean attractor length if the canalyzing function is f 14 . Without delays, only very short attractors (apart from two fixed points) can occur [25] for f 14 , and the increase in attractor length for values τ * > 1 is clearly visible in Fig. 3.
On all other attractors, G 2 has at least sometimes its canalyzing value. The second type of attractors referred to above is obtained if M + 2 is a multiple of τ * . If the loop of M + 2 sites consists of blocks of size τ * , the dynamics can be mapped on that of an effective component withÑ = N τ * nodes and withM = M+2 τ * − 2 with no delay but time steps of length τ * , and the results of [25] can be taken over. In addition to attractors consisting only of homogeneous blocks, further attractors can be constructed by realizing that only one bit in each block is the one that triggers node G 2 . The value of bits that do not trigger node G 2 does only matter when the block reaches G 1 : If at this moment Σ is not canalyzed by G 2 , an inhomogeneous block will not be homogenized, but copied or inverted to node Σ. An inhomogeneous block can therefore survive forever if the blocks that are at G 2 at the moment where the inhomogeneous block is at G 1 do not have their canalyzing value. However, this implies that these noncanalyzing blocks are copied to Σ from G 1 , which is only possible for f 13 . Indeed, from [25] we know that a periodM + 2 of attractors on the effective component is only possible for this function, unlessÑ has special values.
Finally, let us look for nontrivial attractors that can occur even if M + 2 and τ * have no common divisor. Let us choose the function f 2 . An isolated block of size τ * of 1s in the initial state will survive forever since there will be a 0 at G 1 while this block is at G 2 . In fact, there exist a multitude of such attractors where there is a 0 at the right position at distance L behind a block of 1s. This explains why the number of attractors found numerically for τ * > 1 is larger for f 2 than for the other canalyzing functions (see Fig. 3). The length of these attractors can be larger than N + τ * − 1, as can also be seen in Fig. 3.
Function f 9 . f 9 is a reversible function, i.e., if one of the inputs changes its value the output changes, too. If there are no delays, the dynamics is reversible, and therefore all states are on cycles. There is 1 fixed point σ = 1. A striking feature of the synchronous case is that cycles of the order 2 N exist [25].
The exhaustive numerical attractor search suggests that the maximal attractor length can be approximated by τ * times the maximal attractor length in the nondelayed case for odd τ * .

VII. CONCLUSION
In this paper, we studied the influence of a deterministic but non-synchronous update on critical Kauffman network by introducing node-based delays. The dynamics of critical networks can be derived from the dynamics on its relevant components, most of which are simple loops, and some of which have a few nodes with two in-puts. For this reason, we have studied in this paper the three simplest types of relevant components. Not surprisingly, delays typically increase the attractor lengths and reduce the attractor numbers. New types of attractors emerge in the presence of delays. The basins of attraction are naturally larger and thus the path to the attractor becomes more robust. If all delays are randomly chosen real numbers, loops are most likely to be frozen or on a single attractor. Similarly, more complex components with real delays should have far less attractors than for parallel update.
It will be interesting to see how these results are affected when nonrandom networks, such as real gene regulation networks, are considered. Clearly, they are not updated in parallel. Some networks, such as in budding yeast [7] appear to be very robust with respect to the introduction of delays. This means that their choice of connections and functions is such that the update sequence does not matter much. It remains to be seen if this is a general feature of all those networks than can be described by using a Boolean idealization.