the effect on gene i of the product of gene j. The matrix is generally not symmetric and diagonal elements, wii, represent self- regulation. The fraction c of nonzero entries in the W matrix is a parameter of the model and reflects the density of the network. The degree of a gene is represented by Ki and SK T~cN, 0ƒcƒ1 is called network connectivity [23,24]. Each network is a dynamical system, with state vector S(t)~ðs1(t),:::, sN (t)Þ repre- senting the expression levels of each gene at time t. The deterministic, discrete-time dynamics of S(t) are modeled by the set of nonlinear coupled difference equations si(tz1)~f " N # X wij sj (t) , j~1 1 ð Þ where f is a normalization function that prevents the system from diverging (Text S1). More specifically, f is a threshold function (either step or sigmoid; Figure S13), representing cooperative binding and saturation in gene expression. The network is updated synchronously (see [25–27] for asynchronous updates). We define Equation (1) as the development process (see [8,9] for an illustration of the model, as well as a discussion of the biological motivations and assumptions behind it). Since the state space of the model is finite and the dynamics deterministic, the system will eventually reach an attractor given an initial gene expression state. The attractor can either be a fixed point or a limit cycle. The simplicity of the model allows for evolutionary simulations, where a standard population genetic model of mutation, recombination and selection acts upon a population of gene networks, and the network’s state is taken as its phenotype. Despite their level of abstraction, Boolean networks have been highly successful, both at reproducing experimental results in different Most Networks in Wagner’s Model Are Cycling organisms [2,18,28,29], and allowing for theoretical predictions about the evolution of network properties such as robustness, evolvability, and many others [8,9,11,13–17,19–22,30–38]. The properties of this gene network model in the absence of evolution or any kind of selection have attracted considerably less attention in the genetic regulation literature [11,12,18,20,30,31]. On the other hand, the physics community has been studying some theoretical properties of RBNs for some time (see [6,7] for recent reviews). In RBNS, each node is assigned an update function that prescribes the state of the node in the next time step, given the state of its input nodes. This update function is chosen from the set of all possible update functions according to some probability distribution. Since each of the K inputs of a node can be on or off, there are M~2K possible input states. The update function has to specify the new state of a node for each of these input states. Consequently, there are 2M different update functions [7]. RTNs are boolean networks with threshold functions only. The update function is Equation (1) with f (x)~sgn(x) (Text S1). While some analytical results have been obtained for the general class of RBNs, they usually apply only under some restricted conditions, such as in the limit of very large networks, specific network connectivities, or combinations of boolean functions. It has been shown that some results derived under these assumptions break down when only a subset of boolean functions is considered [39–41]. This is the case for RTNs, and although interesting in their own right, theoretical work done with RTNs seems to have been limited [3,24,42,43]. Results As we can see in Figure 1, cycles seem to dominate the dynamics, independently of network size N and degree K . Figure 1. Cycles dominate the dynamics. Average stability (Equation (2)) for different network sizes N, degree K and topology. K~N means c~1. Equation (1) is solved up to n~108 times with wij*N (0,1), si*f{1,1g, f (x)~sgn(x) and sgn(0)~1. Noisy tails for K~2 result from insufficient samples (Figure S4). Our measure is binary (the outcome is either 0 or 1 for each trial) and, for that reason, we do not find it helpful to present a variance measure. Instead, we present full stability distributions for similar experiments in Figure 2. Boxed region represents the size of the genome- wide regulatory networks of E. coli [45] and yeast [47]. doi:10.1371/journal.pone.0034285.g001 PLoS ONE | www.plosone.org 2 April 2012 | Volume 7 | Issue 4 | e34285