Model transcriptional networks with continuously varying expression levels
© Carneiro et al; licensee BioMed Central Ltd. 2011
Received: 27 October 2011
Accepted: 19 December 2011
Published: 19 December 2011
At a time when genomes are being sequenced by the hundreds, much attention has shifted from identifying genes and phenotypes to understanding the networks of interactions among genes. We developed a gene network developmental model expanding on previous models of transcription regulatory networks. In our model, each network is described by a matrix representing the interactions between transcription factors, and a vector of continuous values representing the transcription factor expression in an individual.
In this work we used the gene network model to look at the impact of mating as well as insertions and deletions of genes in the evolution of complexity of these networks. We found that the natural process of diploid mating increases the likelihood of maintaining complexity, especially in higher order networks (more than 10 genes). We also show that gene insertion is a very efficient way to add more genes to a network as it provides a much higher chance of developmental stability.
The continuous model affords a more complete view of the evolution of interacting genes. The notion of a continuous output vector also incorporates the reality of gene networks and graded concentrations of gene products.
In the approximately ten years since the completion of the draft sequence of the human genome, researchers have become increasingly attuned to the many layers of complexity that underlie the mechanisms of life . Many new genes have been identified for transcription factors whose role is to activate or inhibit the production of other genes. The interplay between mutually interacting transcription factors defines a regulatory network that dictates the levels of RNA transcripts, signaling proteins, enzymes and other gene products. Such networks have emergent properties that are essential in every living system.
Understanding the organization and evolution of these networks has been a challenge because of their complexity. Experimental studies have been able to identify important roles of interacting regulatory networks, such as the ability of yeast to respond to environmental changes, the specification of the endomesoderm in sea urchin embryos, and dorsal-ventral patterning in the Drosophila embryo. Although early studies of quantitative traits also revealed clues about such networks, it has not been generally feasible to address the more general questions how they originate and evolve.
A mathematical model of mutually interacting transcription factors was first developed by Wagner. In this model, the level of expression of n transcription factors in an individual at time t is given by the values of the elements in a vector S with n entries. The mutual interactions between transcription factors are represented as an n × n matrix W whose elements w ij are real numbers. W defines the gene network, and the distribution of the non-zero entries of W specifies how connected the network is. Each row in W corresponds to a single transcription factor, and W ij represents the effect of transcription factor j on the production of transcription factor i. The vector S t of expression levels in an individual changes in ontology according to a characteristic time constant τ according to S t+τ = f(W × S t ), where f is a realvalued function applied element by element to the entries of W × S t . In Wagner's model, f (s i ) = -1, 0, or +1 according to whether the i-th element s i of W × S t is < 0, = 0 or > 0. The model therefore allows gene expression to be either completely repressed (-1), expressed at a basal level (0), or completely derepressed (+1). An individual is characterized by the final contents of its gene expression vector S. An individual is considered viable if, and only if, the final gene expression level vector converges to a stable state, meaning, even if the developmental process were continued [S t+τ = f(W × S t )] S would remain unchanged. Occasionally we use the term viable network, by which we mean a W matrix capable of yielding viable individuals from a subset of initial state vectors.
The gene expression patterns generated by this model were able to mimic certain aspects of Drosophila gene expression data  However, the Wagner model does not allow for differing concentrations of the transcription factors. Further studies were stimulated by the potential of Wagner's original model. Siegal and Bergman adapted the model to simulate the evolution of such transcriptional networks, which revealed that stabilizing selection in a population is sufficient to evolve robustness as predicted by Waddington's canalization model of development. In their model, Siegal and Bergman included parameters to dictate the shape of a sigmoidal function, which would allow for a continuous model to be explored, however in all their simulations the parameters were chosen such that the output vector would behave exactly like Wagner's.
This work was followed by additional studies focusing on such issues as network robustness to mutation[11–14]. Bergman and Siegal showed that knockout mutations (replacing an entire row and column with zeroes) would significantly increase the sensitivity to initial conditions, but later studies showed that the relationship between a gene's connectivity and its fitness effect upon knockout depends on its equilibrium expression level. In all these studies, the model was used with the limitation of Wagner's original -1,0,1 output vector. In the model discussed in detail below, we propose an alternate definition of deleterious mutation and analyze its effect on viability in a continuous output scenario.
Sexual and asexual reproduction as well as the coevolution of reproductive method and genetic architecture were analyzed by Azevedo et al[17, 18]. They showed that sexual reproduction evolves greater robustness than asexual reproduction, using the same canalization framework as Siegal and Bergman. But both studies rely on a version of the model in which transcription factors are either on or off and all mating is treated as haploid. In this study we show that, under the continuous model, the effects of mating in haploid and diploid populations are quite different.
The network structure of the Wagner model has been studied in many ways[7, 8, 13, 19–22], and the results suggest that large networks evolve to be sparse[16, 20, 23, 24] and modular[2, 21, 25, 26], The general argument holds that the cost of maintaining unnecessarily complex interactions is too large to be maintained, and once modules of interactions (i.e. smaller networks that achieve a viable final state) are formed, it is easier to combine them together than to evolve the same mechanisms de novo in new networks .
Other types of network models have also been developed. A system of coupled ordinary differential equations derived from the principle of chemical kinetics was used to describe genes and the concentration of their products. Boolean networks were also used to illustrate how gene interaction networks could be modeled together with their products assuming just two states, on or off[24, 28, 29]. Stochastic models based on the Gillespie algorithm[31, 32, 31]were motivated by recent experimental results that have demonstrated that gene expression can have a component affected by stochastic noise[33–35]. These models are limited to very small networks due to the inherent complexity of the algorithm[36, 37].
A continuous network model (that is, one in which the elements w ij in W and s j in the vector S are continuous real variables) might be capable of describing properties of the regulatory networks not present in the discrete model. One approach to a model of a regulatory network with continuous output was based on using artificial neural networks to describe cellular differentiation and morphogenesis.
The discrete version of the Siegal and Bergman model is limited by the concept of viability. Because of the bounding of the gene expression levels to discrete values (-1, 1), the development process St+1 = f(W × S t ) will discard every network in which the matrix W does not have an eigenvector pointing towards one of the vertices of the hypercube centered at the origin. In this sense most viable individuals in a population created using the discrete model will be essentially similar. In a continuous model, networks that have eigenvectors of the matrix W pointing to an edge or face of the hypercube can also yield viable networks. This is due to the broader range of acceptable values of the output vector. As a consequence, the continuous model should yield a more diverse population of developmentally stable matrices (a formal analysis of this point is presented in the additional file 1).
In this work, we describe another approach to a network model with continuous output generalizing that in Wagner (1996). The continuous variation in the abundance of the gene products creates additional complexity that allows a more complete description of the evolution of these networks. The model is intuitively appealing because different concentrations of transcription factors should affect gene expression quantitatively, resulting in different levels of activation and repression. The output vector would therefore be expected to be composed of elements that are continuous rather than discrete. We build on Wagner's model in order to allow detailed comparison between the continuous and discrete versions of the model.
Our model is similar to that of Siegal and Bergman in positing a gene regulatory network of n genes represented by the n × n matrix W, where w ij measures the extent to which the abundance of the product of gene j affects the production of the product of gene i. Given a vector S 0 of initial expression levels in an individual for each of the n genes, the individual is tested for "viability" in a stepwise process of "development" [St+1 = f(W × S t )] that tracks the amounts of the transcription factors in the individual over time. A viable individual is an individual that develops a stable output vector S.
The concept of stability is based on the development of the individual represented at time t by the vector S t . Each developmental step is modeled as the result of multiplication between the matrix W and the vector S t , yielding a new vector S t+1 , which is multiplied again by the same matrix W until the variation in S t+1 is less than some sensitivity constant σ when compared to previous values.
The multiplication of each row i of the matrix W by the vector S represents the interaction between every gene in the network. Each value in the matrix describes the type and strength of the interaction. Positive values mean activation or enhancement of production, whereas negative values mean repression or inhibition. Because the values in the matrix are continuous in the interval [-1, 1], the absolute value describes the strength of the interaction, be it positive or negative.
where lies in the interval [0, 1] and is the average of the expression levels over all times from τ- t to t. The number of genes in the network is n, and X and Y are vectors with n elements each. The number of steps (L) to reach stability is taken as a measure of path length. L max = 100 is the cutoff value in the number of development steps, after which a matrix is considered unstable. As a stability criterion we chose σ = 10-4. This model yields viable individuals in which the levels of the transcription factors are continuous, and is a straightforward extension of Wagner's model of transcription regulatory networks[8, 14]. Both models discard networks with imaginary eigenvalues of the matrix W owing to the cyclic behavior of the output vector and the lack of convergence to stable value of S.
which is a sigmoidal function centered at x = 0. Its curvature is determined by the constant a. Although the bounding effect of this function is obvious, the extent to which it spreads the expression levels is subtler and depends on the curvature.
We analyzed different values of a, and found that values greater than 2 result in a very large variance, forcing all values of the final state vector close to either -1 or 1, essentially making it a discrete model. Between the values of 1 and 2, however, we found that the values of the final state vector would be well distributed along the -1, 1 interval allowing us to use the model in a continuous fashion without affecting the number of viable networks we can yield through our development process. Hence for all simulations in this work, we used a = 1.5. This choice of a is virtually identical to the sigmoidal tanh function used by Kaneko (2011)
The continuous model gives us a more complete view of the evolution of interacting genes. It allows the addition of more genes to the networks and is more efficient in maintaining stability. The notion of a continuous output vector also creates a closer relationship with the reality of gene networks and gene products, where it is not sufficient to ascertain merely whether a gene is on or off. Gene product concentrations play an important role in determining the viability of individuals, and aids in the evolution and maintenance of complexity.
The different mechanisms to generate network complexity tested had a strong impact in the probability of yielding a viable individual, even in networks with many genes and especially for diploid mating. By "diploid mating", we mean that each element of the matrix W of the progeny network equals the arithmetic average of the corresponding elements in the parental matrices. The probability that a random network of 15 genes yields a viable individual is far smaller than that obtained by diploid mating between two viable networks suggests that sexual reproduction may be a key component in the evolution of complexity. The phenomena of insertion and deletion probably also play an important role in the evolution of complexity, given the high probability of a viable individual to remain viable after undergoing an insertion or a deletion.
The continuous model also gives insight into the mechanisms that regulate the evolution of complexity in a general setting that represents the concentration of the products of gene networks. The inclusion of gene product concentration brings the model closer to actual transcriptional networks, and perhaps gives us a better idea of the difficulty of obtaining complex networks that yield viable individuals in the real world.
Results and Discussion
The ICC tests whether the final output vectors of the n viable individuals are clustered together in small regions of the space [-1, 1]. In the equation for ICC, is the sample mean of the elements in the n-th individual, is the mean of all output vectors in the population, k is the number of genes in the network, and s2 is the variance of the elements among the n individuals. The value of ICC can be positive or negative, but a value close to zero means not correlated, whereas a value close to 1 or -1 means high intra-class correlation.
Generating Viable Individuals
Individuals were generated at random by drawing networks and initial vectors of gene expression levels from a uniform distribution on [-1, +1]. One possible interpretation of the initial vector is that it is the level of gene products passed by maternal inheritance to the zygote, where development of the embryo would begin.
The model depends on a set of initial conditions to start the developmental stage of the simulation using a randomly generated initial state vector. It is therefore unclear whether the viability of these individuals is determined by the choice of the initial output vector or by the wiring of the gene network.
To test the impact of the choice of an initial output vector we selected viable individuals and replaced their networks with randomly generated ones, however retaining the original initial state vector. We repeated this process 1000 times for each individual, generating a different network each time, while tallying the number of random W matrices that supported development into a viable individual. Analogously, we performed a similar test by keeping the network while randomly changing the initial state vector instead. These tests enabled an analysis of the relation between stability due to the input vector and stability due to the network.
Evolution of Complexity
Given the very low likelihood that a random W matrix with a high number of genes will be viable (that is, support development of random vectors to a stable state as defined by Equation ), we tested how easily complexity might evolve from combining ("mating") viable networks and producing a new network that may be interpreted as the "offspring" of two viable networks.
Mating was performed by defining a population of 1000 viable networks and mating two randomly drawn networks at a time. The "offspring" network was then tested for stability by iterating random initial vectors according to Equation 1. This process was repeated 1000 times to generate 1000 new progeny networks.
In the process of "haploid mating", a given gene is inherited at random from the network of either parent with equal probability. Accordingly, in the haploid mating process, we randomly selected individual rows from within the paternal or maternal network and copied them to create an offspring network. This process passes on parental genes without modification from one generation to the next. Repeating the selection process for each row yields a new offspring network with a random set of both parents' genes.
The initial state vector of the new offspring is chosen at random to equal a stable state of one of the parents. This procedure reflects the assumption that one of the parents would be passing on the general stable gene-product concentrations to its offspring, analogous to the interaction between an oocyte and its mother during the earliest stages of development.
It is interesting to draw a parallel between haploid mating in the discrete model and the continuous one. Haploid mating displays the same behavior in both models, with high efficiency in generating viable networks with a small number of interacting genes, but then efficiency falls off sharply as the number of genes increase. In the case of the discrete model, the efficiency drops to almost zero with 8 or more genes. In contrast, the continuous model maintains a more consistently slower drop with increasing number of genes, without ever reaching 0 even for networks of size 15.
Diploid individuals benefit from heterozygosity to modulate the effects of damage or deleterious mutations as well as from increasing diversity through the recombination events between the parental chromosomes. In the process of "diploid mating," each row in the W matrix of the progeny is calculated as the arithmetic mean of the corresponding rows in the W matrices of the parents. Biologically, this means that the effects on gene expression are additive, and effects due to dominance, overdominance, underdominance, epigenetics, parent of origin, and so forth are ignored. Taking the impact of each gene as the average of the impacts of this same gene in each parent tends to mitigate large negative or positive effects of the parental genes.
When applied to a set of 1000 viable networks, the diploid mating model generated viable progeny networks of up to 10 interacting genes in 19-32% of the iterations (Figure 4). This percentage is not as high as that in the haploid model, but diploid mating performs better as the number of genes increases. For networks with more than 10 genes, the number of viable offspring networks lies between 43-48%. The positive slope of the curve shows that the diploid mechanism with additive gene effects is very efficient in maintaining stability in complex networks.
A randomly generated network with 15 interacting genes has an 8.9% chance of being viable. When two viable individuals mate following the haploid-mating model, the likelihood of generating a viable network jumps to 22%, however diploid mating increases the likelihood to 47%. This increase may be due to the fact that these original two networks were already selected from a small pool of viable networks with 15 genes, and diploid mating maintains network stability better than haploid mating. We conclude that, while for any level of complexity (number of genes in the network) it is difficult to generate viable complex individuals at random, mating is relatively efficient in producing viable networks of the same level of complexity as those in the parents.
The difficulty in finding a viable network with more than 10 interacting genes prompted the question of whether increasing the number of genes of a viable network is more successful than generating a viable network at random. To answer this question we randomly inserted a gene into a viable network and developed stable state vectors to test whether stability was retained.
A gene insertion represents the phenomenon of a new gene being fully incorporated by the genome and interacting with the other genes in the network. In the inserted gene all interaction values are chosen at random from the uniform distribution [-1,1], and all pre-existing genes receive new randomly generated values for interaction with the newly inserted gene. The stable vector also receives a new randomly generated value, representing the initial concentration of the product of the new gene. The result is a new individual with an extra transcription factor that may or may not be viable when developed with the augmented network.
Similarly to the test with random insertions, the likelihood of obtaining a viable network after removing a gene was tested by deleting one gene at random from a viable network and developing viable individual state vectors to asses if it would remain viable. We performed 100 random deletions in each of the 1000 previously generated viable networks. A gene deletion comprises a row and column deletion in the network, plus an entry deletion for the corresponding gene product in the initial output vector.
For networks with few interacting genes, loss of a gene is critical, with very few networks remaining viable after a deletion. This result is compatible with the difficulty in finding viable networks when there are few interacting genes. With more complex networks the numbers are still high, for example, 67.9% for networks with originally 10 interacting genes, which is significantly greater than the 14% rate for randomly generated networks with 9 interacting genes. Deletion maintains 66.6% of the viable networks with 15 interacting genes.
We presented an alternative model to describe the development and evolution of gene transcription factors that allows for a continuous distribution of expression levels. This version of the model allows the study of more complex network (both in number of genes and degree of connectivity) given the additional classes of networks that yield viable individuals. The continuous model, however, makes it more difficult to define the concept of "neighboring networks" but this may be addressed by defining a threshold below which differences between networks define them as neighbors. Another limitation to our model is computing time, as the matrix multiplication in development and the tests for viability are more time consuming than in the discrete model.
We are grateful to Trevor Bedford and Bernardo Lemos for their helpful comments on the continuous output model. This work was supported by NIH grant GM07953 to DLH.
- Check H: Human genome at ten: Life is complicated. Nature. 2010, 464: 664-667. 10.1038/464664a.View Article
- Davidson EH, Levine MS: Properties of developmental gene regulatory networks. Proc Natl Acad Sci USA. 2008, 105: 20063-20066. 10.1073/pnas.0806007105.View ArticlePubMedPubMed Central
- Lee T, Rinaldi N, Robert F, Odom D, Bar-Joseph Z: Transcriptional regulatory networks in Saccharomyces cerevisiae. Science. 2002
- Levine M, Davidson EH: Gene regulatory networks for development. 2005, 102: 4936-4942.
- Stathopoulos A, Levine M: Dorsal gradient networks in the Drosophila embryo. Dev Biol. 2002, 246: 57-67. 10.1006/dbio.2002.0652.View ArticlePubMed
- Milkman RD: The Genetic Basis of Natural Variation. I. Crossveins in Drosophila Melanogaster. Genetics. 1960, 45: 35-48.PubMedPubMed Central
- Wagner A: Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization. Proc Natl Acad Sci USA. 1994, 91: 4387-4391. 10.1073/pnas.91.10.4387.View ArticlePubMedPubMed Central
- Wagner A: Does Evolutionary Plasticity Evolve?. Evolution. 1996, 50: 1008-1023. 10.2307/2410642.View Article
- Siegal ML, Bergman A: Waddington's canalization revisited: developmental stability and evolution. Proc Natl Acad Sci USA. 2002, 99: 10528-10532. 10.1073/pnas.102303999.View ArticlePubMedPubMed Central
- Waddington C, Kacser H: The strategy of the genes: a discussion of some aspects of theoretical biology. orton.catie.ac.cr. 1957
- Ciliberti S, Martin OC, Wagner A: Innovation and robustness in complex regulatory gene networks. Proc Natl Acad Sci USA. 2007, 104: 13591-13596. 10.1073/pnas.0705396104.View ArticlePubMedPubMed Central
- Ciliberti S, Martin OC, Wagner A: Robustness Can Evolve Gradually in Complex Regulatory Gene Networks with Varying Topology. PLoS Comput Biol. 2007, 3: e15-10.1371/journal.pcbi.0030015.View ArticlePubMedPubMed Central
- Huertasanchez E, Durrett R: Wagner's canalization model. Theoretical Population Biology. 2007, 71: 121-130. 10.1016/j.tpb.2006.10.006.View Article
- Wagner A, Wright J: Alternative routes and mutational robustness in complex regulatory networks. BioSystems. 2007, 88: 163-172. 10.1016/j.biosystems.2006.06.002.View ArticlePubMed
- Bergman A, Siegal ML: Evolutionary capacitance as a general feature of complex gene networks. Nature. 2003, 424: 549-552. 10.1038/nature01765.View ArticlePubMed
- Siegal ML, Promislow DEL, Bergman A: Functional and evolutionary inference in gene networks: does topology matter?. Genetica. 2007, 129: 83-103.View ArticlePubMed
- Lohaus R, Burch CL, Azevedo RBR: Genetic architecture and the evolution of sex. J Hered. 2010, 101 (Suppl 1): S142-57.View ArticlePubMed
- Azevedo RBR, Lohaus R, Srinivasan S, Dang KK, Burch CL: Sexual reproduction selects for robustness and negative epistasis in artificial gene networks. Nature. 2006, 440: 87-90. 10.1038/nature04488.View ArticlePubMed
- Conant G, Wagner A: Convergent evolution of gene circuits. Nat Genet. 2003
- Leclerc RD: Survival of the sparsest: robust gene networks are parsimonious. Mol Syst Biol. 2008, 4: 213-View ArticlePubMedPubMed Central
- Espinosa-Soto C, Wagner A: Specialization can drive the evolution of modularity. PLoS Comput Biol. 2010, 6: e1000719-10.1371/journal.pcbi.1000719.View ArticlePubMedPubMed Central
- Draghi J, Wagner GP: The evolutionary dynamics of evolvability in a gene network model. J Evol Biol. 2009, 22: 599-611. 10.1111/j.1420-9101.2008.01663.x.View ArticlePubMed
- Sackton T, Kulathinal R, Carneiro M: Population genomic inferences from sparse high-throughput sequencing of two. Genome Biol. 2009
- Andrecut M, Kauffman S: On the Sparse Reconstruction of Gene Networks. Journal of Computational Biology. 2008, 15: 21-30. 10.1089/cmb.2007.0185.View ArticlePubMed
- Misevic D, Ofria C, Lenski RE: Sexual reproduction reshapes the genetic architecture of digital organisms. Proc Biol Sci. 2006, 273: 457-464. 10.1098/rspb.2005.3338.View ArticlePubMedPubMed Central
- Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology. Nature. 1999, 402: C47-52. 10.1038/35011540.View ArticlePubMed
- Chu D, Zabet NR, Mitavskiy B: Models of transcription factor binding: sensitivity of activation functions to model assumptions. Journal of Theoretical Biology. 2009, 257: 419-429. 10.1016/j.jtbi.2008.11.026.View ArticlePubMed
- Kauffman S: The Origins of Order: Self-organization and Selection in Evolution (1993). emergence.org. 1993
- Kauffman S: Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology. 1969
- Ribeiro AS, Lloyd-Price J: SGN Sim, a stochastic genetic networks simulator. Bioinformatics. 2007, 23: 777-779. 10.1093/bioinformatics/btm004.View ArticlePubMed
- Gillespie D: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics. 1976, 22: 403-434. 10.1016/0021-9991(76)90041-3.View Article
- Gillespie D: Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry. 1977
- Blake WJ, KAErn M, Cantor CR, Collins JJ: Noise in eukaryotic gene expression. Nature. 2003, 422: 633-637. 10.1038/nature01546.View ArticlePubMed
- Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science. 2002, 297: 1183-1186. 10.1126/science.1070919.View ArticlePubMed
- Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda-infected Escherichia coli cells. Genetics. 1998, 149: 1633-1648.PubMedPubMed Central
- Elowitz MB, Leibler S: A synthetic oscillatory network of transcriptional regulators. Nature. 2000, 403: 335-338. 10.1038/35002125.View ArticlePubMed
- Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403: 339-342. 10.1038/35002131.View ArticlePubMed
- Vohradsky J: Neural model of the genetic network. J Biol Chem. 2001, 276: 36168-36173. 10.1074/jbc.M104391200.View ArticlePubMed
- Geard N, Wiles J: A gene network model for developing cell lineages. Artif Life. 2005, 11: 249-267. 10.1162/1064546054407202.View ArticlePubMed
- Kaneko K: Proportionality between variances in gene expression induced by noise and mutation: consequence of evolutionary robustness. BMC Evol Biol. 2011, 11: 27-10.1186/1471-2148-11-27.View ArticlePubMedPubMed Central
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.