This is the first time networks of the size of an organism have been simulated with this model. To the best of our knowledge, the largest network previously studied has N~400 genes [44]. We have extended this range to N~10,000 (Figure 1). This is comparable to E. coli (*1,800 genes identified in its regulatory network [45] or *4,300 annotated genes in the genome of the K- 12 strain [46]) and to yeast (3,420 network genes [47] or *5,8000 annotated genes in the Saccharomyces cerevisiae genome [48]). We have also seen that stability seems to decay much faster for degree K ~N than for smaller K . For networks of size N~65, for example, the probability of fixed points is already smaller than 10{4. More importantly, this number seems to approach 0 for larger N, unlike the fat tails of K~2. In other words, by choosing K ~N, one cannot study stable states for organism-size, genome- wide networks. Interestingly, dense networks are the most popular choice in the literature [16,22]. Our results have clear implications for the interpretation of previous studies [8,16]. By limiting their analysis to small dense networks with stable dynamics, most previous findings are limited in reach and may not apply to real biological networks. Figure 2 clearly suggests there are two types of networks in this model: stable and unstable. This bimodal nature of the stability distribution is not trivial. It is interesting that the stability of a random pair of matrix and initial state tells us how likely it is that the matrix is stable (or not) with any other random initial state. It is the network that determines stability, not the initial state. This implies that future studies using this model do not have to sample many different initial states to characterize network dynamics. Sampling pairs of networks and initial states, as we have done in most of our study, should suffice. These two classes of networks may have different topological properties. Since previous studies only use viable networks, most of the networks they allow to evolve are of the second type, i.e., stable. It is true that we do not expect biological networks to be random; the question is if we do want to start our simulations from a non-random set, are the chosen matrices biologically relevant? And, by selecting these and not others, what properties or biases are introduced into the simulations, prior to evolution? To the best of our knowledge, these questions have yet to be addressed. A great deal has been written about scale-free networks in biology [49,50]. In contrast, most networks we have studied here are regular, since this is the topology frequently used with gene regulatory networks [8,9], and thus the ones we are interested in characterizing. We have shown, however, that different topologies, including scale-free out-degree distributions, do not seem to change the overall results, in agreement with previous studies [22,34]. In Figure 3 we see that stability seems to decay faster with c for larger networks and in Figure S5 we have tried to show this explicitly finding that for small to intermediate size networks (N~10*20), the difference between the stability of dense and sparse networks is maximized. This is about the size of networks used in previous studies [11,33], where networks with different densities coexist in a population (i.e. adding or deleting connections is allowed). We believe this stability difference should be taken into account in analysis of such studies. Further, network stability does not seem to depend strongly on the nature of the matrix weights, either binary or real-valued (Figure 5). An important parameter of the model is how cells respond to their input signals. That is to say, is gene regulation a switch-like process, or is the response graded? While the former is implemented by a step function, the latter is modeled as a sigmoid curve (Figure S13). One could argue that a switch-like mechanism already introduces a lot of robustness to the model, in the sense that most changes do not produce a visible change in the Most Networks in Wagner’s Model Are Cycling phenotype (for xw0 or xv0 in the figure). Sensitivity sharply increases, however, at DxD?0, where the discontinuity occurs. Close to zero, very small changes in the regulatory inputs of a gene can quickly turn it on or off. In this DxD?0 region, the system is clearly not robust. The opposite can be said of a not very steep sigmoid function. By allowing continuous expression values, small changes in x lead to small changes in a gene’s state s~f (x). This choice, however, allows changes in positive or negative inputs around DxDw1 (or even DxDw5) to still produce visible changes in phenotype. In this sense, the system is less robust. We have shown that the two choices are only equivalent, at least in terms of stability, for large a, where the sigmoid behaves almost like a switch (Figure 6). With small a, for which the sigmoid is less steep and behaves more like a gradient, the choice between a step [8,14] or a sigmoid function [9,16] makes a difference (see Text S2 and Figure S12 on how to deal with f (0) in the discrete case). As already mentioned, a lot of work has been done on analytical properties of Random Boolean Networks [23,51] and it has been suggested that some properties of RTNs do not follow analytical results derived for the general class of RBNs [39–41]. The dynamics of RBNs with canalizing functions only, for example, seem to be dominated by fixed points [52]. We have shown here that this is clearly not the case for RTNs. It has also been suggested that the attractor length distribution of RBNs follows a power-law [7,40]. Again we have shown this is not true for RTNs, which instead seem to follow an exponential decay for a range of parameters (Figure S10 and Figure S11) perhaps slightly slower than initially suggested [43]. Methods The attractor reached by the dynamical system (1) is uniquely defined by the matrix W and the initial state S(0), and is either a fixed point or a cycle. Let n represent the number of pairs of W and S(0) for which Equation (1) is solved, and nf ƒn the number of times the attractor is a fixed point. To estimate the probability of reaching a fixed point steady state within the framework of this model, we generated up to n~108 random pairs of matrices and initial states, for different network sizes N and degrees K , and measured stability~nf =n, 2 ð Þ which takes values between 0 and 1. The estimated probability of cycles is given by 1{stability. As mentioned before, most of the evolutionary studies done with this model vary in the parameters used in Equation (1). We estimate the dependence of stability measured by Equation (2) on most of these parameters. We list in Table S1 all the variations we have studied along with the corresponding figures and references. Special relevance is given to the range of parameters used in previous studies. For this reason, we mostly study small [8,9] and dense [16,22] networks. Real gene networks, however, can be quite large [45,47] and appear to be sparsely connected, with an average of two transcriptional regulators per gene [32]. Stability estimates are also shown for these more realistic topologies [49,53]. More detailed methods and algorithms are included in Text S1. Supporting Information Text S1 Supporting Methods. More details on experimental procedures and algorithms for solving Equation (1) for both the discrete and continuous cases. (PDF) PLoS ONE | www.plosone.org 6 April 2012 | Volume 7 | Issue 4 | e34285