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 [1]. 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[2].

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[3], the specification of the endomesoderm in sea urchin embryos[4], and dorsal-ventral patterning in the *Drosophila* embryo[5]. Although early studies of quantitative traits also revealed clues about such networks[6], 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[7]. 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[8], *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 [8] However, the Wagner model does not allow for differing concentrations of the transcription factors. Further studies [9]were stimulated by the potential of Wagner's original model. Siegal and Bergman[9] 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[10]. 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 [15]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[16]. 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[10] framework as Siegal and Bergman[9]. 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[20], 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 [21].

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[27]. 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[30] 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[38] to describe cellular differentiation and morphogenesis[39].

The discrete version of the Siegal and Bergman[9] 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 *S*
_{t+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[8] in order to allow detailed comparison between the continuous and discrete versions of the model.