Compiling the parts list is the first step in developing any model of some complexity and is not always a trivial exercise. Simple parts lists of genes, transcription factors, promoters, binding sites and other molecular entities are useful means for assessing the network complexity and for comparing different organisms. Such parts lists can be the result of a genome-sequencing and annotation project, where the complete DNA sequence of an organism is determined and all (or at least many) genes and proteins are identified. A parts list could also be represented as a database of regulatory elements or it could be ontology terms of transcription regulation processes assigned to a set of genes.

Comparing such lists from different organisms can provide an indication of the complexity of the transcriptional machinery or be used to predict the presence or absence of particular metabolic pathways [4-6]. The number of known and predicted transcriptional regulators in eukaryotic organisms varies from about 300 in yeast to about 1000 in humans (Table 1).

Many publications address the computational identification of transcription factor binding sites for instance by analysing promoter sequences of coexpressed genes [7]. One approach is to search for short sequences that are overrepresented in the promoters of a particular group of genes (e.g. clusters of coexpressed genes or sets of longer sequences known to bind a particular transcription factor) in comparison to the promoter sequences of all other genes. This approach obviously depends on the availability of the sequences for many genes and their upstream regions. Such an approach was applied in a cell cycle study in Schizosaccharomyces pombe, where Rustici et al. showed that the presence or absence of consensus binding sites in the promoter regions corresponds to the cyclic expression pattern of the genes [8]. Genes with a peak expression at similar cell cycle stage often share similar sets of consensus binding sites.

However, the exact promoter regions are usually unknown and even the transcription start sites are only known for a few genes. Baker's yeast (Saccharomoyces cerevisiae) has a relatively small genome with short intergenic regions, considering about 600–1000 bp upstream of the translation start site (ATG) appears to be a good approximation for the promoter regions. In higher organisms like vertebrates the intergenic regions and thus the putative promoter regions are much larger than in yeast, therefore the identification of regulatory elements in the DNA sequence by computational means has turned out to be rather elusive. Some studies have focused on the computational analysis of higher-level organisation of transcription factor binding sites in promoters, such as frequently occurring combinations of known binding sites [9,10], or restricted the search for regulatory elements to conserved sequence regions, identified by genome comparisons, a method often referred to as phylogenetic footprinting [11]. However, phylogenetic footprinting does not always work, because the localisation and the binding sites themselves are not always conserved [12,13].

Transcription factor localisations can be identified experimentally, too. For example, individual binding sites can be detected using the "DNAse I footprinting assay"; proteins bound to the DNA protect it from degradation by DNAse I, therefore these regions can be analysed further [14]. Another common experimental method is the "electrophoretic mobility shift assay" (EMSA) sometimes called "band shift assay" or "gel retardation assay" – DNA fragments that are bound by protein move slower in an electrophoretic gel than unbound fragments [15,16]. These methods allow fine mapping of individual binding sites, but are very labour intensive. High-throughput methods such as the Chip-on-chip method (see additional file 2, which describes this methods in more detail) allow the genome-wide detection of binding sites for a transcription factor, but the spatial resolution and signal quality is limited. Furthermore, assigning transcription factors to their target genes based on the genomic localisation is difficult due to the size of intragenic and intronic regions and long range effects of some transcription factors.

Nevertheless parts lists provide the first impression of gene networks in different organisms and they are necessary to have before we continue by having a look at the network topology.

Once we know the transcription factors and their binding sites, we can describe the gene transcription regulatory networks by graphs with nodes corresponding to genes and edges to regulatory interactions [17]. (Note the difference between these discrete graphs and plots of mathematical functions also often referred to as graphs.) A short introduction to discrete graphs is given in additional file 1, for more details please consult for example Cormen et al. [18]. One important concept that we will use below is a representation of a graph by a so-called adjacency matrix, where the element a_ij in a row i and column j equals 1 (i.e., a_ij = 1), if node i is connected to node j, otherwise a_ij = 0. Graph representations have been used for various biological data sets ranging from protein-protein interactions networks to coexpression networks, they have been long used in mathematics, physics and computer science, and many aspects of graphs and their applications have been studied (e.g., [19-21]).

In a directed graph (i.e., a graph where connections between nodes have a definite direction) we call genes (nodes) with outgoing edges (arcs) source genes. For a given source gene, we call the set of all genes with incoming arcs from that source gene its target genes. Regulatory relationships can be of various natures. For a specific model we need to define the precise meaning that we assign to the edges (Figure 2). For instance, an arc from a gene A to B may mean that source gene A is a transcription factor, which is known to bind to the promoter of target gene B. A rather different network will be obtained if an arc from A to B denotes the observation that the disruption (e.g., mutation) of source gene A changes the expression of target gene B. We will present examples for these types of networks in the next section. A widely studied type of molecular network of a different type is the protein-protein interaction network, where the nodes represent proteins and two proteins are connected by an undirected edge, if they bind to each other (Figure 2) [22]. A different network will be established by connecting genes based on their sequence similarity. Networks can also be built based on the co-occurrence of gene names in journal abstracts. If two gene names frequently occur in the same abstracts it is likely that they share some kind of functional relationship [23].

To illustrate knowledge that we can obtain from studying network topology, we will compare and combine information from two high throughput data sets for the yeast Saccharomyces cerevisiae. The first is obtained from chromatin immunoprecipitation experiments for transcription factors (ChIP network), while the second is obtained from microarray experiments on single gene deletion mutants (mutant network), for more detail on the experimental methods please refer to additional file 2. Microarray experiment measurements can be presented in a data matrix, where rows represent genes and columns particular experiments (hybridisations). For instance, in the ChIP experiment each column will correspond to a particular transcription factor (studied in the particular experiment), while in the mutant experiment each column will correspond to a particular mutant. In this way, not only rows, but also columns will correspond to genes. The measurement values are typically real numbers, such as intensity levels, expression levels, or p-values, depending on the data processing steps applied. By applying an additional data processing step, often called thresholding, we can transform these continuous values into discrete values (e.g., if we chose a threshold t, then we can replace any x_ij by 0, if x_ij <t, and by 1, otherwise). In this way we will transform the original measurement matrix for these experiments into an adjacency matrix defining a graph: two genes in this graph are connected, if the measurement value is higher than the chosen threshold.

The ChIP network is based on experimental data published by Harbison et al. [24]. As described in additional file 2 they used genomic tiling arrays to identify the genomic regions bound by transcription factors. The authors assigned each genomic region to one or two target genes based on proximity in the genome. Relative intensities of spots are the basis for an error model that assigns a probability score (p-values) to binding interactions, which we use for discretisation.

The starting point for the mutant network is the gene expression data matrix published by Hughes et al. [25]. Each experimental condition, which in our case is a particular gene deletion mutant, corresponds to a column and each gene corresponds to a row (for more detail see additional file 2). We discretise the data matrix using a gene-specific standard deviation estimate γ obtained from the error model proposed by Hughes et al. [25].

The size of the networks depends on the discretisation thresholds chosen (Table 2). The criteria used to choose these thresholds are rather subjective, for the following comparisons we focused on these thresholds: ChIP network p_t = 0.001 and mutant network γ = 2.5.

What do these two networks mean and how do they compare? In the ChIP network, an arc A→B means that the gene A codes for a transcription factor that binds to the promoter of gene B, while in the mutant network it means, that the mutation of A will change the expression level of B [23]. The ChIP network describes physical interactions, but it does not tell us anything about the effects of these interactions. The mutant network is similar to the one used in gene networks built by classical genetics means – we know that a mutation (perturbation) of the first gene has an effect on the second one, but it does not necessarily mean a direct physical interaction – there may be a long transcriptional or signalling cascade leading from the first gene to the second. In this way the first network is likely to contain direct interactions, while the second may include indirect interactions as well. However, it is possible that some of the 'direct' interactions of the Chip network are not biologically functional and thus may not be supported by mutant network. Most importantly it should be noted that the experimental conditions in the two experiments are not identical, which may result in considerable discrepancies between the two networks.

First we can observe that both networks consist of one major component and almost all genes are part of it and are connected (this is true for a wide range of discretization thresholds). The degree distributions resemble roughly a power-law, i.e. most source genes have few target genes, while few source genes have many (Figure 3). Rung et al. discovered that the number of connections can indicate the functional class of the gene. Genes in the mutant network with many outgoing arcs (high outdegree) often encode proteins with regulatory functions, whereas genes with many incoming arcs (high indegree) are predominantly involved in metabolism [26]. Functionally related genes tend to be close in the networks, it is therefore possible to identify functionally related genes by comparing their neighbourhoods [23]. Manke et al. found directly interacting transcription factors and those, which are members of a protein complex, to occur together as putative DNA-binding modules more often than expected randomly [27].

When we compare and combine the mutant network and the ChIP network we immediately observe that their intersection is sparse. They share 102 edges connecting 13 from a total of 23 shared source genes to 93 distinct target genes (Table 3 and Figure 4). In the mutant network these 13 source genes are connected to 937 distinct target genes by 1157 edges (89 edges per source gene, 8.8% of the connections are in the intersection); whereas in the ChIP network they are connected to 631 distinct target genes via 924 edges (71 edges per source gene, 11.0% of the connections are found in the intersection). We used the hypergeometric distribution to compare target sets of all source genes [23] to identify source genes with target set intersections larger than expected by chance (Figure 5). We found that 9 of the 23 shared transcription factors Arg80p, Gcn4p, Hir2p, Mbp1p, Stb4p, Ste12p, Swi4p, Swi5p and Yap1p have significantly similar target sets in both networks (p < 0.05, 4000 genes in total) (Table 3). In these cases the transcription factor localisation might actually explain the changes in gene expression we see in the corresponding deletion mutants. One would assume that due to the nature of the experiments some effects in the mutant network are indirect effects that could be explained by a combination of direct connections in the ChIP network (Figure 6A). We find that indeed a number of connections in the mutant network can be explained by a combination of two edges in the ChIP network (Figure 6B).

Others have observed that when comparing different protein-protein interaction networks, their intersections are small, too. Only few interactions are reported by several experiments despite using just different methods to measure the same interactions [28]. Reassuringly, the shared interactions turned out to be more reliable than most of the data [28]. Here we work with experiments measuring different, if somehow related effects, but still the proportion of connections between functionally related genes is increased in the intersection of the two networks. In the intersection 40 of 102 (39.2%) connections connect genes that have the same cellular role in YPD, compared to less than 20% in the original networks (Table 2).

To summarise we can say that although the two networks are rather different, the part that is common to both is biologically more meaningful, and some of the indirect interactions of the mutant network can be explained by the direct interactions in the chip network.

Network topologies, particularly in yeast, have been widely studied and many interesting observations have been made. It has been proposed that the existence of highly connected genes (hubs) in a network might make networks more tolerant to random failure of network elements [29,30]. In protein-protein interaction networks it seems possible to classify hubs in combination with expression data: Han et al. [31] showed that hub proteins can be divided into two groups based on the level of coexpression between their neighbours in the network (the proteins directly connected to the hub proteins). Hubs with low coexpression seem to link functionally separate modules and removing these hubs leads to more rapid disintegration of the network [31]. However so far this has not been observed in transcription networks.

Luscombe et al. compiled data from ChIP-on-chip experiments for yeast to construct a network of 142 transcription factors, 3420 target genes and 7074 regulatory interactions [32]. To study the dynamics of this network they traced the paths from the target genes back to initial transcription factors, starting from target genes that are differentially expressed under particular conditions as demonstrated in previously published microarray experiments. Depending on the conditions, different sets of genes are expressed, leading to different sets of target genes as start points for the backtracking and different transcription factors along the path. This is consistent with results from ChIP-on-chip experiments and demonstrates that the topology of the gene regulatory network is not independent of the experimental conditions [24].

In a different line of investigations, Lee et al. [31], and Milo et al. [33], identified re-occurring structural elements (motifs) in the networks. They examined topological networks derived from ChIP-on-chip data for structures consisting of 3, 4 or more edges that occur in the original network more often than in randomised networks. Network motifs they identified to be significantly more frequent than in randomised networks included feed-forward and feedback loops. These motifs may partly be the result of gene duplications during genome evolution [34].

These are just some examples of possible analyses that can be performed on topology level. However, arguably the main reason to study the network topology is to prepare the ground for the next step of building more detailed models for gene networks. Before any logic or dynamic network model can be constructed we need to know which gene products interact and which are mutually independent. Even if we take the view that in the real-world network every gene is connected to every other gene to some degree, not all these connections are equally strong and a discretization step can be used to keep only the strong connections in the model. A complete network where everything is connected to everything may not be a practical approach to network modelling on a genome scale.

Arguably the most important question is if we can find modules, i.e. subnetworks that are relatively isolated from the rest of the network. If such modules are found, they can help us to use the reductionist approach later on by allowing modelling the parts of the network independently on a more detailed level. For instance, if we can build a dynamic model of an independent module, then we can perform simulations independently from the rest of the network. The existence of modules in biological systems has often been taken as an axiom [35]. However, a precise definition for what constitutes a module is elusive, and therefore this term has been used in various contexts [36]. In a graph representation it is natural to define a module as a 'relatively' isolated component, and indeed such components were found in protein-protein interaction networks. In contrast isolated components have hitherto not been found in the wiring diagrams of eukaryotic transcription regulation networks [26]. Several methods have been proposed to identify modules as groups of genes coexpressed under specific conditions [37,38], but there remain controversial opinions regarding the existence and nature of modules in gene networks [39,40]. Biologically meaningful pathways are sometimes used to define modules and to implement a reductionist approach. However, this easily breaks down in conditions when the pathways interact.

In general, data sets used for topological models have important limitations. While hundreds of organisms have been fully sequenced and many genes are identified relatively reliably, the data sets underlying most topological models are much less complete. Only a fraction of all protein-protein interactions in yeast have been tested; most large-scale experiments show high noise levels; and whereas the genome sequence is independent of particular growth conditions and (sometimes) is even conserved in fossils, data like protein-protein interactions and transcription factor localisations are condition dependent. In this context it is particularly important to note that some experimental methods are performed under conditions considerably different from the natural conditions in the cell. For example, the yeast-2-hybrid technology was used to determine protein-protein interactions between human proteins, yet the yeast cell provides a very different environment from a human cell [41,42]. This can be a considerable source for variation and systematic errors. Some experimental techniques are performed in test tubes, thus providing the most "unnatural" conditions. Unfortunately various limitations are unavoidable and we have to work with incomplete data for a limited set of conditions.

We can conclude that genome scale topological representations have helped us to make many interesting observations about the network properties, however the main question of finding well defined modules of such networks on topology level is still open.

Once we know the network topology, the next step is to study the rules of interaction between the different elements in the network. For instance, if a promoter consists of only one binding site for a transcription factor, we may want to know whether it is an activator or a repressor. If several transcription factors bind to a promoter, we need to know what each factor does, but also how these factors interact (Figure 7). Biological studies demonstrate that some promoters show combinatorial behaviour that can be approximated by Boolean functions (AND, OR, NOT and combinations of these), but in other cases the interaction is more complicated [43]. Linear functions, Boolean functions, decision trees, and Bayesian probability distributions have all been used to describe the network logic. We can distinguish between discrete functions and continuous functions. Discrete functions are based on the assumption that a gene can be in a finite number of states. In the simplest case we use only two different states to describe the activity of genes (e.g., expressed and not expressed). We can thus use Boolean functions to describe the interactions between transcription factors, e.g. "gene j is active, if transcription factor A AND B are bound to the promoter". It has to be stressed that such 'states' are only approximations of reality, that in the real world the interactions are not so well defined and are often fuzzy.

Continuous functions use continuous values (real numbers) to represent the gene activity. Weights w_ij represent the interaction between genes i and j, which can be positive or negative. Thus the activity g_i of gene i can be calculated as the weighted sum of the activities of all n genes:

g_i = w_i1 g₁ + ... + w_in g_n

This approach assumes that the influence of one gene on another gene is linear. Note that the network topology will determine which of the weights w_ij are equal to 0 (i.e., if there is no arc from gene i to gene j in the network topology, then w_ij = 0). Like Boolean functions, linear functions are only approximations. For instance, it is not possible by linear functions to model a situation where the same transcription factor can play a role of an activator or repressor for the same gene, depending on the presence or absence of other transcription factors.

Although few promoters have been studied in great detail, there are excellent examples, such as the description of the promoter action logics of sea urchin developmental gene Endo16 [44]. The Endo16 promoter consists of almost 30 regulatory elements stretched over a region of 2.3 kb. Based on experimental data collected using modified promoter constructs Davidson and co-workers constructed a model expressed as an algorithm combining Boolean and linear functions. This algorithm takes as an input the occupancy information from 12 binding sites and outputs a value, that 'can be thought of as the factor by which, at any point of time, the endogenous transcription activity (...) is multiplied as a result of the interactions mediated by the cis-regulatory control system' [44]. Predictions of promoter manipulations based on this model have largely been confirmed in subsequent experiments. Extending their earlier work the group of Davidson compiled a regulatory network containing over 40 genes by the construction of a model that integrates extensive experimental evidence on early development of sea urchin embryos [45].

Recently Klamt et al. published an example for control logic networks [46], based on hypergraphs, which are an extension to the graphs described above. Several hyper-edges pointing to the same node represent OR relationships, but edges are allowed to combine to represent AND relationships. Weights on the edges distinguish positive and negative relationships. The authors provide a set of methods to analyse these networks, just to list a few examples: computation of all positive and negative signalling paths, computation of all positive and negative feedback loops and computation of minimal cut sets. These minimal cut sets report the smallest number of interventions necessary to force the network into a particular behaviour, for example, a minimal number of deletions necessary to block the activation of a particular downstream protein in a signalling cascade. These methods are implemented in the software tool CellNetAnalyzer and the example presented, a model of a signalling network for T-cell activation shows that these analyses are non-trivial for signalling networks of a typical size.

Soinov et al. used a supervised learning approach to build decision-tree-related classifiers. A decision tree is a predictive model (Figure 8). Soinov et al. built decision trees which allow us to predict the gene expression activity of a particular gene (leaf node) based on the expression data of other genes (interior nodes) [47]. Although the gene predictions are binary in this approach (the gene is predicted be "active" or "inactive"), this system utilizes continuous expression values, such as microarray data.

Bayesian networks provide a probabilistic framework for modelling gene regulatory networks [48-50]. Their graphical representation is a directed acyclic graph, where each node represents a variable and the edges represent dependencies. For a more detailed description of the application of Bayesian networks in gene expression analysis see the reviews by Pe'er [51] and Friedman [52]. Segal et al. applied a learning procedure based on probabilistic graphical models to networks consisting of groups of coregulated genes [53].

There are situations where neither Boolean rules nor linear functions are powerful enough to express the control logics: transcription factors might bind competitively, if one factor is bound, the other one is excluded, as is the case for example in the phage λ switch between lysis and lysogeny [54]. In some cases, transcription factors have to form homodimers or heterodimers to be fully functional. The transcription factors might have to bind sequentially or might act synergistically. In these situations it might be necessary to use more complex functions (here this would be solved by Boolean circuits with memory or delay). It remains an open question what is the minimum repertoire of functions to describe regulatory logics.

The knowledge of the parts list of a network, its topology and the control logics are necessary requirements in order to expand the model to capture dynamic changes during time. Compared to the approaches above, the dynamic models can be described as 'classical' approaches to gene network modelling, as many of them have been developed and studied long before the current genome era. Typically, they are relatively small, involving only a few genes. They aim at describing and often simulating the dynamic changes in the state of the system and predicting the network's response to various environmental changes and stimuli.

Various dynamic models have been proposed. Greller and Somogyi subclassified them [55] as follows: "Dichotomies for framing our thinking on how to best represent a particular biological network problem include the following distinguishing attributes: quantitative versus qualitative measurements; logical versus ordinal variables (e.g. Boolean versus abundances); deterministic versus probabilistic state transitions (e.g. differential equations versus hidden Markov); deterministic versus statistical overall system description (e.g. vector field versus Bayesian belief network probability distributions); continuous versus discrete state (e.g. continuous intensities or concentrations versus low, medium and high); nonlinear versus linear elementary interactions and state update rules (e.g. multiplicatives, sigmoids or non-monitonics versus linear ramps); high-dimensional versus low-dimensional (e.g. >> 100s of variables versus << 100 variables); stochasticity present and profound versus absent or present as nuisance noise (e.g. probabilistic state transitions versus small amplitude errors); measurement error substantially corrupting and obfuscating versus negligible distortion." In the following sections we describe several approaches following the discrete to continuous model axis. The discrete model approaches we consider include Boolean network based models [56-58] and Petri nets [59-62], the dynamic systems are based on difference or differential equations [63-65]. We will then discuss hybrid models, which combine discrete and continuous elements [66-68].

Reverse engineering refers to an approach where one starts from data and tries to design a model that fits the data (semi-) automatically within the given model class, without additional prior hypothesis about the biological system. The model derived from the data is judged by the results of simulations compared to new experimental data. For example, one could use a gene expression data set to construct a particular gene network model that is consistent with the data. Inconsistencies between simulated data generated using this model and new data, that has not been used to construct the model, indicate shortcomings of the model. These inconsistencies can be used to choose between alternative models, or to improve the model. However, reverse engineering is possible only (1) if we have chosen an appropriate model class (in the sense that the desired properties of the real world network can be described in it), and (2) if we have enough quantitative data describing the behaviour of the system. Of course, even if the answers to these two questions are positive, reverse engineering is still a difficult problem, and few efficient algorithms are known. The methods chosen for reverse engineering depend crucially on the kind of modelling technique used. Quantitative models are normally more demanding than qualitative models. Dynamic models contain many parameters, and detailed experimental data are required to work out the parameters.

Miyano et al. have proposed algorithms to infer Boolean networks [67,72] and Friedman et al. developed methods to extract probabilistic graphical models, such as Bayesian networks from experimental data [49,52]. Tegner et al. proposed an approach for the reverse engineering of dynamic gene networks based on integrating genetic perturbations [97]. They identified " [...] the network topology by analysing the steady-state changes in gene expression resulting from the systematic perturbation of a particular node in the network." [97]. However, they only apply their approach to simulated data and to a comparatively small biological system consisting of only 5 genes.

A powerful approach to test our understanding of gene regulatory networks is to build new networks from scratch in an approach called synthetic biology. Predictions of small models have been successfully tested experimentally using specifically engineered control circuits, such as feed forward loops [98] and feedback loops [99-103]. In a sense this is reverse engineering of a real world network. For a more detailed description see the reviews by Kaern et al. and Ball [104,105].