Most Networks in Wagner’s Model Are Cycling Ricardo Pinho1,2*, Elhanan Borenstein3, Marcus W. Feldman1 1 Department of Biology, Stanford University, Stanford, California, United States of America, 2 PhD Program in Computational Biology, Instituto Gulbenkian de Cieˆncia, Oeiras, Portugal, 3 Department of Genome Sciences, University of Washington, Seattle, Washington, United States of America Abstract In this paper we study a model of gene networks introduced by Andreas Wagner in the 1990s that has been used extensively to study the evolution of mutational robustness. We investigate a range of model features and parameters and evaluate the extent to which they influence the probability that a random gene network will produce a fixed point steady state expression pattern. There are many different types of models used in the literature, (discrete/continuous, sparse/dense, small/large network) and we attempt to put some order into this diversity, motivated by the fact that many properties are qualitatively the same in all the models. Our main result is that random networks in all models give rise to cyclic behavior more often than fixed points. And although periodic orbits seem to dominate network dynamics, they are usually considered unstable and not allowed to survive in previous evolutionary studies. Defining stability as the probability of fixed points, we show that the stability distribution of these networks is highly robust to changes in its parameters. We also find sparser networks to be more stable, which may help to explain why they seem to be favored by evolution. We have unified several disconnected previous studies of this class of models under the framework of stability, in a way that had not been systematically explored before. Citation: Pinho R, Borenstein E, Feldman MW (2012) Most Networks in Wagner’s Model Are Cycling. PLoS ONE 7(4): e34285. doi:10.1371/journal.pone.0034285 Editor: Matjazˇ Perc, University of Maribor, Slovenia Received February 10, 2012; Accepted February 25, 2012; Published April 12, 2012 Copyright: 2012 Pinho et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: The Ph.D. Program in Computational Biology is sponsored by Fundac¸a˜o Calouste Gulbenkian, Siemens SA, and Fundac¸a˜o para a Cieˆncia e Tecnologia (fellowship SFRH/BD/33531/2008). The research was also supported in part by National Institutes of Health grant GM28016 and National Science Foundation award CNS-0619926. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The Ph.D. Program in Computational Biology is sponsored by, among others, Siemens SA. There are no competing interests or relationship between the authors and this commercial funder, besides that. This does not alter the authors’ adherence to all the PLoS ONE policies on sharing data and materials. * E-mail: rpinho@stanford.edu Introduction Gene regulatory networks have been studied intensively in recent years, both by physicists and biologists, who have provided different insights into this important field [1,2]. We present numerical simulations to investigate stability in large Random Threshold Networks (RTNs) [3]. This spin glass or neural network-type model [4] represents a subclass of Random Boolean Networks (RBNs) [5]. We consider all attractor types rather than only those networks that have fixed points. Thus our framework is that used by physicists [6,7] rather than that used by biologists [8,9]. Wagner [8,10] introduced a version of this gene network model to study the evolution of genetic robustness. The gene network is represented by a dynamical system whose state variables correspond to expression levels of the network’s genes. A network is said to be robust if it retains the same expression state after mutation. The transient from an initial state to an attractor represents a developmental process and the fixed point attractor represents the phenotype. For this reason, fixed points have been traditionally considered developmentally stable, and networks that have cycling dynamics are not allowed to survive. This requirement of developmental stability [9] can be viewed as viability selection [11]. Implementations of the model for evolutionary simulations have varied in parameters such as network size, connectivity, normal- ization function, and whether the components of both the state vector and the matrix are discrete or continuous [8,9,11–22]. In principle, all of these parameters may influence the dynamics of the model and, consequently, the results of evolutionary simulations. It is well known, for example, that prior to evolution, smaller networks are more robust to mutations than larger ones, and that this relationship reverses after selection [8]. Our goal is to systematically explore how changes in all of these parameters affect the probability of fixed points, in the hope of motivating discussion of the relevance of this model for evolutionary analysis. To this end,we focus our attention solely on the gene network model itself, on what has been called developmental dynamics, without evolution. We generate millions of random networks of size up to N~10,000 and measure the probability of fixed point dynamics for most of the different parameterizations reported in the literature, and show that cycles always dominate network dynamics. Fixed point steady states are the exception, not the rule in this gene network model. We also show that stability, defined as the probability of fixed point dynamics, decreases with network size and density. Stability distributions are bimodal: some matrices are always stable independently of the initial state, while others never reach a fixed point. Other measurements like period-size distributions show further deviation of the properties of this network model from those of the general class of RBNs. The layout of the paper is the following: we finish the Introduction by presenting a more detailed version of the Model; we next present our Results (more of which are detailed in Text S2 and Supporting Information Figures), followed by a Discussion; we conclude with a short Methods section where we present a summary of each model variant used to produce our Figures (Table S1; more detailed methods can be found in Text S1). Model The model consists of an interaction network of N genes, represented by an N|N matrix, W , whose elements, wij, denote PLoS ONE | www.plosone.org 1 April 2012 | Volume 7 | Issue 4 | e34285