Reverse Engineering of Genome-wide Gene Regulatory Networks from Gene Expression Data.

Transcriptional regulation plays vital roles in many fundamental biological processes. Reverse engineering of genome-wide regulatory networks from high-throughput transcriptomic data provides a promising way to characterize the global scenario of regulatory relationships between regulators and their targets. In this review, we summarize and categorize the main frameworks and methods currently available for inferring transcriptional regulatory networks from microarray gene expression profiling data. We overview each of strategies and introduce representative methods respectively. Their assumptions, advantages, shortcomings, and possible improvements and extensions are also clarified and commented.

INTRODUCTION

Transcriptional regulation plays crucial roles in protein synthesis and its dynamical responses to internal and external signals, such as development processes and environmental stimuli [1, 2]. The temporal and spatial levels of mRNA and ultimately protein abundance are actually controlled by transcriptional regulations in a cell [3]. A regulation system consisting of genes, RNAs, proteins, and other molecules constructs the complicated regulatory interactions during sequentially transcriptional, post-transcriptional, translational and post-translational processes, which structure into multiplex networks [4]. A transcriptional regulatory network generally refers to regulatory activities between regulators, e.g. transcription factors (TFs), and their targets, e.g. genes [1, 5]. A gene’s transcription will be initialized or terminated by the TF proteins binding to its promoter region generally at the 5’ upstream of the transcription start site. To some degree, the final expression abundance is mainly determined by the activation or repression of their regulatory relationships [2, 6, 7]. Without distinguishably considering the physical regulations, a gene regulatory network refers to a collection of gene-gene interactions corresponding to such regulatory relationships through their products, and the interactions in gene regulatory network denote this kind of regulations. In contrast, a transcriptional regulatory network represents the physical bindings and direct regulatory interactions between regulators and their targets [8]. It contains more concrete and specific regulatory information between TFs and genes. From a systematic perspective, genome-wide transcriptional regulatory networks in cells control gene expression dynamically and precisely in response to biological context specificities [9].

Identifying transcriptional regulatory networks is of paramount importance from deciphering transcriptional mechanisms to uncovering potential drug targets [10, 11]. Various network reconstruction methods have been proposed and they can be generally categorized as ‘bottom-up’ and ‘top-down’ methods. The traditional gene knockout experiments can be categorized as bottom-up methods, which firstly identify the detailed regulations between TFs and targets individually, and then summarize all these regulations to form a regulatory network. The genetic relationships between genes can be detected from the effected genes after knocking out some gene [12-14]. And a global gene regulatory network can be built up after collecting these experimentally identified genetic interactions. Alternatively, top-down methods refer to the emerging systems biology approaches of identifying the global regulatory interactions systematically and in parallel. They firstly acquire many potentially regulatory interactions and then validate each of them by additional experiments. For instance, ChIP-Seq technology makes the genome-wide identification of protein-DNA interactions possible [15, 16]. The regulatory elements of DNA-binding proteins such as TFs are identified from massively parallel sequencing [17]. A genome-wide regulatory network is then drafted from these identifications. The details of TF-target binding event in specific conditions are often checked by further experiments [18]. Microarrays are another type of systematic expression monitoring technologies, which measures the amount of mRNA produced during transcription by hybridization [19, 20]. The reconstruction or inference of regulatory network from microarray gene expression data is often called a reverse engineering process, which backwardly reasons the regulatory system from its observational behavior [21]. Recently, the reverse engineering of transcriptional control network from microarray data becomes very popular to revealing genome-wide regulations [21-25]. Numerous computational strategies have been proposed to reconstruct large-scale gene regulatory relationships from expression profiles [26-29]. Several papers [30-33] have summarized and compared the available strategies from different perspectives. For instance in [30], Emmert-Streib and colleagues presented a systematic overview and comparison study of the network inference methods. They conceptually categorized the existing methods from statistical learning perspective. In this review, we focus on these available computational methods by highlighting their assumptions, advantages, weaknesses, possible improvements and future research directions individually.

Computational methods of inferring transcriptional regulatory networks from expression data are highly motivated by the availability of genome-wide expression profiling data [34-37]. The activities of gene regulation are closely related to gene expression levels [6, 38]. Gene expression profiles of time series or perturbations indicate the dynamics and differences of genes and then imply the causal regulatory possibilities between them. Moreover, the individual gene pairs between regulators and target genes should also be considered with cooperative and systematic perspectives, such as co-regulations, competitive regulations of activators and repressors, and indirect genetic regulations [9, 37, 39]. A global transcriptional regulatory network is embedded with high interacting affinities between regulators and targets, which can be learned from transcriptomic data. And the details of individual regulatory events are hypothesized to be validated by further experiments [13, 40]. The top-down method generates a global view of regulatory relationships in form of network illustrating the context-dependent scenario of regulations. Existing computational methods of inferring regulatory networks are all to formulate the regulations into certain models with these measured expression values [23, 26, 27].

In this review, we firstly formulate the reverse engineering of transcriptional regulatory networks from transcriptomic profiles into a general framework, and then review the major available strategies developed to address this problem, e.g., correlation-based methods, Boolean network methods, Bayesian network methods, differential equation methods, and knowledge-based methods of integrating and evaluating prior regulations. We focus on introducing the assumptions and main ideas behind these strategies and their approximations in the modeling of regulatory systems. Then the current research directions and alternatives of deciphering regulatory network from expression data are discussed. A brief vision of reconstructing transcriptional regulatory networks from high-throughput expression profiling dataset is then concluded.

FRAMEWORK OF REVERSE ENGINEERING

The surge of microarray technologies provides unprecedented opportunities to measure genome-wide gene expression simultaneously [19]. Various strategies have been developed to infer the regulatory architectures from their corresponding gene expression profiles for transforming experimental data into regulatory knowledge [22]. The inferred networking linkages represent the regulatory relationships among these measured genes.

(Fig. ) illustrates the general framework of the reverse engineering of transcriptional regulatory networks from gene expression data. Essentially, transcriptional regulatory network reconstruction is to identify physical and genetic regulatory relationships between TFs and target genes from their expression profiles. Without distinguishing the difference between TF and its own gene, gene regulatory network is often used as an approximation to the transcriptional regulatory system. Since the abundance of TF protein is often not available, it is approximated by its gene s expression. Specifically, a transcriptional regulatory system is represented by a network, whose nodes refer to regulators and target genes and whose edges indicate their regulatory interactions. As shown in (Fig. ), from microarray gene expression data, such as profiles of time-series physiological processes or perturbation experiments of gene knockout or RNA interface, we reversely engineer the network structures and parameters, e.g., regulatory logic, causality and strength, from the measured gene expressions by developing models and algorithms. The measured genes are those nodes in the regulatory network, and the linkages and related parameters can be identified from the patterns underlying the gene measurements. The regulatory network and expression data are often represented by regulatory matrix and expression matrix , respectively, i.e.,

, ,

where entry is the regulatory interactions between the -th gene and the -th gene ( ), and entry represents the gene expression value of the -th gene ( ) at the -th experiment ( ). It is noted that refers to a sample or a time point with specific phenotype meaning. The process of reverse engineering is to determine the unknown elements of matrix from the known , which is a reverse strategy for reconstructing the underlying regulatory relationships of the system.

As illustrated in (Fig. ), there are four levels of clarity for the elements of , which answer different questions about the regulatory parameters respectively. Suppose there are two genes, and . From the available gene expression data , Level I inference is to determine whether there is a regulatory connection between and from data . Let and represent the regulatory interactions from to and that from to , respectively. Level I is to determine whether . The binary decision makings build the fundamental architecture of these regulations from gene expression data. Then, when we identify the causal influence from the regulator of TF to its target gene , Level II inference determines the edge direction and causality in the regulatory network, i.e., . In certain conditions or states, TF might activate or repress the transcription of a target gene, and the concentration of the target is then increased or decreased accordingly. The edge orientation underlying the Level Ⅲ regulatory relationship contains the type information of activation and repression, i.e., when activates , and when represses . More specifically, when we identify the regulation strength from to in the Level IV inference, such as , it provides concrete regulatory weight of its transcriptional dynamics. Level I inference is to reconstruct gene regulatory interactions, while the other inference levels contain more detailed information about transcriptional regulatory interactions, such as regulator and target, activation and repression, and concrete regulatory strength. The strong or weak regulation can then be relatively assessed when all the real numbers of regulatory strengths are determined. (Fig. ) shows the direct modeling of the regulation in an isolated gene pair and in a simple regulatory system respectively. The left graph refers to the regulation between and , while the right one shows the direct causality from to and the indirect influence transferring from . When the system contains a large number of genes, it is apparent that they are needed to be modeled in a systematic manner.

The intrinsic difficulties of transcriptional regulatory network reverse engineering come from several sources. Mathematically, one difficulty is the so-called curse of dimensionality, i.e., in the formation of expression matrix . For intensive cost, there are often a few samples ( ) of microarray that have been experimented, while thousands of genes ( ) have been tested simultaneously in each experiment [41]. From the statistical learning perspective, it is hard to infer a reliable solution of gene regulations from expression data [27]. Moreover, genome-wide regulatory networks tend to be sparse [34-36, 42], all of which result in the high likelihoods to achieve false positive regulations or low likelihoods to achieve false negative regulations [34-36]. Biologically, gene regulation is a complicated physiological process that contains some important steps, such as TF selectively binds to the upstream of the transcription start sites of certain genes to initialize the transcription. Thus, we often model the regulatory system by simplifying some mechanisms, such as cooperation or competition of the TF regulators [43]. Furthermore, the real environment of gene regulation is very dynamic with respect to temporal and spatial features. For example, the up-regulation of one gene encoding a TF can sequentially affect its downstream targets and some regulations can only take place in particular cell types [44, 45].

The reconstructed regulatory network is a graphical representation of transcriptional topology of both trans- and cis-regulations [46]. The static network structure is usually not efficient to describe the three-dimensional regulatory contexts in cells [47]. Moreover, the epigenetic regulations, such as DNA methylation [48], histone modification and nucleosome positioning [49], strongly influence transcriptional concentrations [50]. miRNAs are also regarded as crucial regulators in the post-transcriptional regulations [51]. The multiplex, hierarchical, heterogeneous regulatory processes are intensely cooperative to generate gene expression levels of mRNA abundance detected by microarray. At the same time, the microarray technique of measuring gene expression is still in its maturing period. The sample preparation, such as cell numbers [52], as well as data preprocessing alternatives including probeset design, background correction and normalization [53, 54], highly affect the quantitatively measured values. Furthermore, the cognate mRNA level is used to represent TF activity in the reverse engineering. The abundance mismatch between mRNA and protein also interfere with the inference of the regulation system [55]. These obstacles challenge the perfect reconstruction of regulatory relationship from expression data.

To address these difficulties of reverse engineering regulatory networks, numerous efforts have been devoted and many substantial regulations have been discovered by in silico methods and validated by traditional experiments [30-33, 56]. An international competition named DREAM (Dialogue for Reverse Engineering Assessments and Methods) has been initialized to catalyze the quantitative modeling of transcriptional network inferences [57, 58]. For evaluating the reconstruction performances, several types of measures have often been utilized, e.g., general statistical measures, functional consistency measures and network-based measures [30]. For widely-used statistical measures, the evaluations are often implemented by opening the expression profiling dataset and blinding the benchmarked network structure. After the transcriptional regulatory interactions are inferred from the data by some proposed method, the assessments are performed by comparing the identification results with the benchmarked network [23]. Compared to true regulations, these measures are employed to evaluate the predictions, e.g., sensitivity, specificity, accuracy, F-measure, and Matthews correlation coefficient [30, 59]. The tradeoffs between sensitivity and specificity are often presented by the receiver operating characteristic (ROC) curve. The area under the ROC curve (AUC) is often calculated for assessment [60]. Currently, many methods for reverse engineering regulatory networks have been available [22, 27, 30]. Instead of introducing them individually, we categorize them into several main streams of strategies and introduce their main ideas and philosophies.

EXISTING METHODS

Due to the difficulties mentioned above, the transcriptional regulatory network inferences are far from accurate and perfect [61], and almost all available methods have their own advantages and drawbacks [27, 61]. We summarize them into the following five categories, namely correlation-based methods, Boolean network methods, Bayesian network methods, differential equation methods, and integrative prior knowledge-based methods.

Correlation-based Methods

The first endeavor to identify the regulatory relationships in thousands genes measured in microarray is to investigate their pairwise correlations. If gene highly coexpresses with gene , that is to say, when gene s expression grows up, gene s expression grows up or down simultaneously, then the association between the two genes can be detected and modeled by some methods. The regulation can be inferred according to their transcriptional dependence. For multiple genes, clustering is often employed to identify the coexpressed genes [62, 63]. The genes in the same clusters or groups characterize similar expression patterns during physiological processes. They are often assumed to be regulated by the same or related TFs. Two correlation measures are widely used to detect the associated gene pairs, i.e., correlation coefficient [64] and mutual information [65].

The most popular linear correlation between two variables is Pearson s correlation coefficient (PCC). Suppose gene and gene have a series of measurements and , where , then the PCC between and is estimated by the sample correlation coefficient, i.e.,

,

where and are the sample means of and , and and are the sample standard deviations of and . WGCNA (Weighted Gene Coexpression Network Analysis) is a representative method of building the gene coexpression regulatory network by employing PCC [66]. (Fig. ) shows its general framework [67]. Firstly, a clustering method such as hierarchical clustering is implemented to group thousands of genes into some clusters. In each cluster, the highly coexpressed genes are linked by correlation values. For example, when exceeds a defined threshold such as , a functional linkage between and is created in the resulting coexpression network. After the pairwise functional implications between any two genes are identified, a genome-wide network is built up. The simplicity underlies the method that makes it popular to analyze gene expression data, especially to build gene coexpression relationships [68]. Beyond the linear correlation metric of PCC, some rank-based correlations such as Spearman s correlation are also employed to detect the relationship between genes [69]. These correlations replace gene expression values to their relative ranks and then calculate the correlation coefficient between the two ranking lists.

Mutual information (MI) is often employed to measure the non-linear gene expression associations between pairs of genes [65, 70]. Generally, MI is an information-theoretic measure of the mutual dependence between two random variables. For two genes and , it is defined as

,

where two gene expression values construct two vectors, in which the elements denote their expression values in different samples respectively. and are the marginal probabilities of each discrete value in and in , respectively. is the joint probability of and . High MI value indicates that there may be a close relationship between the two genes, while low MI value implies their independence [60].

MI has been widely used to identify transcriptional regulatory relationships from gene expression data [71]. The quick and accurate estimation of MI is a crucial step in the reverse engineering because computing pairwise MI is nontrivial and quite time-consuming [72]. Similar to the PCC-based framework shown in (Fig. ), the available approaches compute the pairwise MI between all gene pairs and construct an association matrix. RN (Relevant Network) chooses the gene pairs when its MI value exceeds a given threshold of significant value [65, 70]. ARACNE (Algorithm for the Reconstruction of Accurate Cellular Networks) implements the data processing inequality on each connected gene triplet to remove the least significant edge in the MI relevant networks [73]. CLR (Context Likelihood of Relatedness) transforms the MI values into z-scores and connects the genes by employing a background sensitive estimator [74]. MRNET (Maximum Relevance Network) is built on the MI-based mRMR (minimum redundancy maximum relevance) feature selection method [75]. MINET presents a software package of MI estimators for inferring large-scale transcriptional regulatory networks [76]. By implementing these MI-based methods, some important transcriptional regulations have been revealed and validated [77, 78].

Unlike PCC and MI, maximum information correlation (MIC) is proposed to detect the strength of any type of linear or nonlinear correlations between genes [79]. MIC adopts binning as a scheme to apply MI to calculate the association between gene variables. It is defined as

,

where is the MI of and . , are the numbers of bins and bins divided, and the total number of bins is constrained to be less than some number . MIC defaults and is the sample size [79]. Although the effectiveness of MIC is controversial [80], it devotes an effort to identifying diverse types of gene relationships and indicates the importance of an association metric to identify genetic relationships [81].

The correlation or coexpression is a fundamental strategy to identify the regulatory relationships at the former Level I and Level IV inferences (Fig. ) and should be improved to be more reasonable in the reverse engineering [82]. Although it is found that the genes in the same grouped clusters tend to have similar functions, these genes might have no direct interactions with each other, and there is no any information to distinguish causal regulators and responsive targets. The built network is not directed (Level II) and without the causality of functional linkages (Level Ⅲ) [83], though it can be determined by additional information, such as annotated TFs [77]. Moreover, the clustering methods such as hierarchical clustering are highly dependent on the threshold chosen to cut the hierarchical tree (dendrogram). The number of clusters and chosen distance metrics also highly affect the resulting networks [81]. It is often assumed that there is modularity property in coexpression regulatory networks, which means dense connections between the genes within the same modules but sparse connections between genes in different modules [68]. The clusters form the building blocks of genome-wide regulatory networks. The linkages between these modules are often omitted in these available methods [84]. These functional linkages indicate the crosstalk and functional cooperation between these modules upon certain conditions [67, 85-87].

Another important issue of this type of methods is the isolated modeling of individual gene pairs as shown in (Fig. ). The regulatory effect from to can also be transferred from . The indirect regulations highly bias the inferred results [88, 89]. We should consider the degree of association with the removal of the effects from indirect regulations by controlling one or several other genes. Partial correlation coefficient can be employed to quantify the association between two genes when conditioning on other gene or genes [88]. For instance, conditioning on a gene or gene set , partial correlation between gene and gene is to measure the exact correlation between the parts of and that have no relationship with . The order of partial correlation coefficient is determined by the number of conditioned genes. Obviously, the mentioned PCC is the zeroth-order partial correlation coefficient. Theoretically, it can be raised to any arbitrary order. The first-order and second-order partial correlation is defined as and ,

respectively. In practice, it is difficult to calculate high-order partial correlation coefficient because of the curse of dimensionality. It is often estimated by developing some specific computation techniques in the reverse engineering of regulatory networks [89].

Similarly, conditional mutual information (CMI) measures the conditional dependency between two genes given other gene or gene set. The CMI of genes and given is defined as

CMI has been applied to reconstruct genome-wide regulatory networks [90-92]. The recently proposed MIC is also expected to be extended to calculate the conditional and partial versions for detecting more delicate and meaningful associations between genes [93].

Based on CMI, we proposed a reverse engineering method [60] by utilizing path consistency algorithm [94] to remove the edges with conditional independent correlation from the network. (Fig. ) shows the general framework of our PCA-CMI method. The main idea of PCA-CMI is to eliminate the edges with independent correlations recursively, i.e., from low to high order independent correlation until there is no edge that can be removed. Firstly, we began with a complete graph, in which all the possible regulations among these genes are contained. Secondly, for adjacent gene pair and , we calculated MI , i.e., zeroth-order CMI. We removed the edges between genes and if they have low or zero MI values. Thirdly, for adjacent gene pair and , we computed the first-order CMI conditioned on their adjacent gene . We removed the edge between them if they have low or zero CMI. The next step is to identify higher order CMI until there are no more adjacent edges to be eliminated [60]. Since it is also time-consuming to calculate CMI [60, 90], in our proposed algorithm, with the assumption of Gaussian distribution, CMI is estimated with Gaussian kernel probability density estimator [56].

From a regulatory system perspective, linear regression methods identify the associations among genes comprehensively [95, 96]. Compared to the former correlation or partial correlation based methods, the regression methods model each gene by multiple predictors. They associate the expression of one gene to all the genes in the whole system and then identify these predictors by variable selection. So the cooperative regulatory relationships among genes can be identified simultaneously. Let denote a gene and be the gene set potentially regulate gene . Their relationship is modeled by a linear function, i.e., . The ordinary least squares, partial least squares and maximum likelihood methods can then be used to estimate the parameters of the linear system [97, 98]. Under the parsimony assumption, a regulatory network tends to be sparse [34, 36, 42]. Some variable selection method such as LASSO [99] and elastic net [100] are often employed to recognize the crucial regulators by the regularization techniques [101]. Specifically, LASSO minimizes the residual sum of squares subject to a bound on the -norm of the coefficients, i.e.,

where and Obviously, some coefficients may be shrunken to zero and the global linkages (coefficients) between these genes can be then inferred. We can find that the causal relationships or directions between these genes are embedded in the regression model. Regression combined with variable selection formulates the regulations into a systems biology approach to reconstructing the underlying genetic interactions from expression profiles. Apparently, regression-based methods achieve a sparse regulatory network and perform the four levels of regulation inferences shown in (Fig. ). For time course expression data, the vector autoregressive model is also employed to specify the gene expression value by a linear regression of those of earlier time points [97]. Similarly, Granger causality is modeled to learn time-lagged regulatory networks from time-course gene expression data [102, 103].

Boolean Network Methods

One of the main-stream strategies to reverse engineering transcriptional regulatory networks is based on Boolean networks. Boolean models treat the genes in a regulation system as logical elements [104]. It assumes that a single gene can be represented by a Boolean variable denoting whether it is expressed or not. The wiring of an element to one another corresponds to functional linkages between genes, and the Boolean rules determine the result of a regulatory signaling transduction given a set of input values [105, 106]. Boolean network provides a simple decision-making model of describing the regulatory mechanisms in a transcriptional system [104, 107, 108].

Specifically, a Boolean network is a directed graph , where the set of nodes representing genes. (Fig. ) shows a simple example. For each node , , a Boolean function is associated with it individually. The inputs of are from the specified parent nodes in to each node . The variable is Boolean and its value is often denoted as 0 or 1 which corresponds to the logical value True or False respectively. The logic operators AND , OR , and NOT are employed to define the Boolean operations in these genes [107]. At any given time , an expression pattern of names a state of a Boolean network, i.e., . The state at time point is determined by Boolean functions from the state , i.e., . The states of all nodes are updated according to their respective Boolean functions and all states transitions together correspond to a state transition of the regulatory network.

For representing the state transition, it is convenient to build a corresponding wiring diagram of a Boolean network as shown in (Fig. ) [106, 109]. For each node , let be the parent nodes of in . By introducing an additional node , we link an edge from ( ) to . Then in the resulting network. Apparently, the expression pattern of the additional node set is determined by individually and corresponds to the regulatory network state at the next time point. If we regard the expression patterns of the set as the input of , the expression patterns of are the output as shown in (Fig. ).

The reverse engineering of a Boolean network is to infer the Boolean functions at these nodes from expression data. When is known, the underlying network topology of regulations can be built spontaneously. An exhaustive search is to try out all Boolean functions on all combinations of out of genes. It is known to be an NP-complete problem and takes exponential time in the inference [105, 106]. So it is often tractable by employing certain computational techniques to avoid exponentially searching a consistent network structure with the observational data. When multiple network structures are found to be consistent with the gene expression data, more scoring metrics and assumptions can be defined to select one suitable regulatory architecture [26, 110].

Boolean network is a fundamental model of genetic system which identifies the network structure from a systematic perspective. It fulfills the Levels I, II and Ⅲ inferences of gene regulatory networks. The dynamic property and the simplicity in understanding and analyzing make it an attractive model of regulatory network reverse engineering. However, the binary and synchronous (i.e., the state of all genes updates to the next one at the same time) assumptions are not consistent with the true biological system [111]. To address these limitations, the discretization strategies and Boolean models have been extended in various ways to make them more biologically realistic and computationally tractable [26]. With the availability of gene expression data with larger sample size and higher quality, there have been approaches to introducing stochasticity to these models, such as probabilistic Boolean networks [112-114] in which the state transition diagram is stochastic. The generalized Boolean network models also try to cope with the shortcomings by enabling more sophisticated forms of logical update which allows asynchronous transition of elements [115].

Bayesian Network Methods

Definition. Bayesian network is a directed acyclic graph (DAG) representing a set of random variables and their joint probability distribution together with the family of conditional probabilities induced by the graph [116, 117].

Bayesian network is a typical probabilistic graphical model of causal inference in statistics. The general idea of learning Bayesian network structure from data is to evaluate each network structure with respect to the given data by defining a scoring function and to identify the optimal one according to the score [118]. The structure represents the conditional independence of these variables that facilitate their joint distribution to be decomposed. The graph is often assumed to follow the Markov property that each gene is independent of its non-descendents, given its parents in . By applying the chain rules of probability and the properties of conditional independency, the joint distribution on genes can be uniquely represented by the product form ,

where is the set of parents of in . In this way, each Bayesian network specifies the joint probability distribution over all genes down to the conditional distributions of the genes given their parents. As shown in (Fig. ), gene D is dependent on gene A and gene E, and independent on the other gene or genes. The global network probability is determined by the dependence structure between multiple interacting components.

The graphical representation consists of two distinct parts in reverse engineering transcriptional regulatory networks. The first component is a DAG representing the causal relationships of regulations (i.e., edges of set ) among a set of genes (i.e., nodes of set ). An edge exists from gene A to gene B if and only if A is a direct regulator of B. The second component is a set of parameter , which describes a conditional probability distribution of each gene, given its parent regulators. Taken together, the two components specify a probability distribution over the set of genes in , i.e., the network structure of regulations. Often, Bayesian scoring metric is derived to evaluate the posterior probability of a graph given the gene expression data , i.e.,

where is a constant which can be ignored [119]. In a Bayesian network framework, the calculation of the log marginal likelihood involves the probability of the data over all possible parameters assigned to . It is an NP-hard problem to select the maximum scored network structure given the data [117, 118]. Thus, the most probable network structure is generally implemented by approximating the posterior probabilities of the regulatory combinations heuristically [37, 118]. Bayesian network model becomes appealing for modeling causal relationships between these genes by selecting the most likely causalities in form of a DAG [9, 29, 119]. Some techniques have been developed to narrow down the search space to a tractable size. As an assumption, the basic form of Bayesian network cannot handle cyclic regulations and the temporal dynamic regulatory relationships [117]. Other alternatives have been proposed to extend the applicability of Bayesian network modeling, such as dynamic Bayesian network [120-124], module network [84] and state-space model [121, 125].

Based on the framework of Bayesian network, dynamic Bayesian network (DBN) introduces the time concept and models a stochastic temporal process of a set of random variables over time series [121-123]. It has been employed to describe the qualitative nature of the dependencies that exist between genes in a temporal process. The structure of a DBN is assumed to perform regulatory functions over discrete time points indexed by . Similar to the assumptions in Bayesian network, let be the gene expression vector of genes at time . For the time points , under the first-order Markovian assumption, i.e., is independent of for given , we thus have

in the time-course gene expression data [123]. As illustrated in (Fig. ), the underlying acyclic graph in Bayesian network can now be permitted to contain cycles. DBN model can explore the general network structure of gene regulations and overcome the shortcomings of the acyclic assumption and static network structure in Bayesian network learning models. A more complicated time-varying DBN model of describing the time-evolving network structures underlying the time series is also developed [126].

Differential Equation Methods

Differential equation formalisms including ordinary and partial differential equations have been widely used to describe and simulate dynamical systems in science and engineering. The powerful mathematical methods have been implemented to model the biochemical systems of metabolic processes and kinetic dynamics of genetic regulation processes [25, 26]. The regulatory interactions in form of network are revealed by the differential and functional relations between the time-dependent concentration variables [36, 127]. Here, we mainly introduce the ordinary differential equation (ODE) models in modeling transcriptional regulatory network. Partial differential equation (PDE) models contain the similar framework as ODE with more dynamic dimensions beyond the time in ODE [26]. ODE models directly consider the time differentiation and then the dynamics and causal relationships can be simultaneously identified in the four inference levels (Fig. ) of reverse engineering regulatory network.

In ODE models, the change rate of gene expression of a component in a regulatory system is modeled as a function of the concentrations of all the components. Mathematically, the general ODE model can be formulated as

,

where represents the gene expression values of genes at time point . The causal effects of gene expressions are embedded in the ODE system [128, 129]. Function describes the relationship between the first order derivative of X and the concentration of genes in the regulatory system. It is a linear or nonlinear function that describes the relationships between the change rate concentration of genes and their causal regulators. Specifically, a linear ODE model can be written as

,

where is the intercept and denotes the regulation effects of genes in the regulatory system on the rate of expression change of the -th gene.

The problem of network reconstruction from data is then transformed to identify the parameters in the ODE system. Traditionally, the least squares method and likelihood-based methods are implemented to find these parameters [27, 34]. Various techniques have also been employed to evaluate them [41, 130]. However, these methods are not effective for reverse engineering genome-wide regulatory networks. We and Lu et al. [128, 129] proposed an integrative pipeline to address the problem by introducing a two-step paradigm to identify these parameters effectively. The first step is to fit the mean curves of the gene expressions and then to estimate the derivative value respectively, i.e., , , , where is estimated continuously from the mean curve, is one of the set time points in range . Thus, the regulatory system becomes the following pseudo-regression model, i.e.,

,

where is the error term of estimation. Based on the parsimony assumption, the second step is to conduct the variable selection and estimation procedure by a regularization framework, such as LASSO [99] and SCAD [131], to shrink the variables as optimally as possible. The regulatory network is then reconstructed from the data when we identify the parameters of the formulated linear regression system. Original methods [128, 129] include a clustering procedure to divide these genes into groups with similar expression profiles, which helps to build a genome-wide network and simultaneously avoid the identifiability problem [132].

ODE is a directed network model and the dynamic feature of regulations is automatically and naturally quantified. In ODE models, gene regulations are modeled by derivative equations, which quantify the change rate of gene expression of one gene (dependent variable) in the system as a function of expressions of all related genes (independent variables) that refer to its regulators. In a transcriptional regulatory system, it is TFs that regulate the gene transcriptional processes. The abundance of TF proteins is the real independent variables. We usually have no such information and simply use the TF genes expression as approximation. Under such assumption, the reverse engineering of regulatory network becomes inferring the parameters of some specified functions such as the former linear function from gene expression data [128]. According to the differences between a mathematical modeling perspective and a statistical perspective lying in the network inference [133], ODE is to model the regulatory system but not to directly infer the regulatory network. The derivation equations are firstly assumed to describe the functional relationships among genes and their products. Then, the statistical techniques such as parameter estimation and variable selection are implemented to infer the regulatory architectures [128]. The resulting nonzero regulatory linkages construct a regulatory network. Time delay of the activation and self-degradation can also be flexibly integrated in the dynamical system by introducing certain terms in the differential equations, such as , where denotes a time delay and for the -th gene s self-regulation [26]. Compared to the former regression methods of modeling the mRNA concentrations of individual components in the system, ODE describes the derivatives of their concentrations. The strategies of parameter estimation are similar to each other.

Knowledge-based Methods

With the essential difficulties in the reverse engineering of regulatory networks, purely data-driven method is very difficult to identify genuine transcriptional regulations. It is hard to promise the effectiveness and efficiency of the reverse engineering only from gene expression profiles [22, 27, 134]. There are urgent requirements to develop novel methods that can utilize expression data in some alternative manners. At the same time, various prior knowledge of gene regulations from literature and genomic datasets can provide additional functional linkage information between genes, such as documented regulations [135, 136], TF binding sequence motifs in promoter region [45], ChIP-Seq data of protein-DNA binding [137] and protein-protein interactions [59]. These prior knowledge can be integrated together with gene expression data to identify transcriptional regulatory networks. Theoretically, the resolution space can be narrowed down to improving the identification significantly [138-140]. So it guides the inference in right direction and helps remove false positives in the predictions [141, 142]. Knowledge-based methods fall into two subcategories, the combination of prior knowledge and the evaluation of prior knowledge. We review them individually as follows.

The combination of prior knowledge is often implemented on the former reviewed reverse engineering methods. Bayesian network is one of the rational models to integrate prior knowledge in a principled manner to increase the inference reliability [140, 142]. According to the Markov assumption, the probability of a network structure can be decomposed as

,

where is the parents of in the DAG. The probability of a local regulatory structure is then calculated according to the structural knowledge priors,

.

The decomposition facilitates to incorporate the prior knowledge about regulatory structure into the network inference. Various techniques have been proposed to calculate these probabilities, i.e., and , as accurately and effectively as possible. Following a framework of statistical physics, [139] and [143] proposed an energy function to introduce the prior knowledge from multiple sources into the reverse engineering of regulatory network. Their main idea is to express the available prior knowledge in terms of network energy. Specifically, the prior knowledge about the regulatory relationship between gene and gene is represented by . Network energy of a network is then defined on the biological prior knowledge matrix. Then, a prior distribution over network structures is obtained by means of a Gibbs distribution [139]. The parameter of this distribution represents the weight associated with the prior knowledge relative to the gene expression profiles. In this way, the prior knowledge is integrated into a Bayesian network framework to learn the regulatory network structure. They achieved higher performance of inference in both simulated and real data [139, 143].

Based on an ODE model, we proposed a method of linear programming (LP) to integrate prior knowledge in the reverse engineering of regulatory network [138]. The main idea is to build an LP model to minimize the association gap between gene expression data and network structure with constraints of the priori of regulatory relationships, and then to solve the LP to obtain the integrated regulatory network.

Specifically, given an experiment with genes and samples, the gene expression matrix is , where is the expression level of the -th gene in the -th sample. We employed an ODE model to quantify the rate of change of gene expression as a function of the expression of other genes [138]. Due to the unclear structures of regulatory system and data scarcity [41, 95, 138], we used the simplest linear additive models:

for , where the state variable is the mRNA concentrations of gene at time point , is the self-degradation coefficient, is the external stimuli, which is set to 0 when there is no external input, and represents the error and noise. describes the type and strength of the effect of gene on gene , whose positive, zero or negative values indicate the activation, naught or repression regulatory relationships between them respectively. For simplicity, we set . Hence, the equations can be described as:

where After we approximated by and neglected the error part, the linear additive model becomes

where . Instead of solving the equations by singular value decomposition (SVD) technique [41, 95, 127, 138], we derived a sparse regulation network [36, 56] based on an LP model. At the same time, more and more prior knowledge of gene regulatory network can be obtained from various sources. For example, if we know that gene and gene are interactive with the rule that activates , such priori should be guaranteed in the inference procedure and the inferred network should contain such information as activates .

In our LP model [138], the objective function is to minimize the number of gene connections to realize the sparseness of the inferring network, and the constraints are the linear additive equations and the prior knowledge of some local network structures. The model is described as

There are variables and constraints. It is equivalent to solve a canonical LP:

Clearly, there are variables and constraints. In the canonical form, the linear objective function can be defined as:

The sparseness and the prior knowledge for regulatory network are represented in the objective function and in the constraints of the LP model, respectively. When we let and , the objective function becomes

Hence, a sparse regulatory network is achieved from gene expression data by minimizing these regulatory strength coefficients with the constraints of the prior knowledge about the gene relationships. Generally, there are three kinds of the prior knowledge about the functional relationship between gene and gene ; activates/represses ( ), has no any relationship with ( ), and has some relationship with , but unclear of positive or negative regulation ( ). These prior knowledge are reflected in the constraints by the defined rules. If gene is an activator of gene ( ), we set as a constraint in our LP model. Conversely, if gene represses gene ( ), we set . If gene has no any relationship with gene ( ), we set as a constraint. If it is unclear which one is an activator or repressor ( ), we set the constraint as and . By solving the two LP models with the two constraints respectively, we selected the sparser solution as the inferred network [138].

Evaluating Prior Regulations

Due to the complexity of gene regulation and the difficulty of network inference from expression profiles, reverse engineering cannot easily identify genuine regulatory relationships [27, 134]. An amount of knowledge about gene regulations has been deciphered by decades of endeavors [41, 144]. Alternatively, we can evaluate the knowledge-based gene regulations documented in literature and databases and filter out the activated regulations in certain biological conditions and phenotypes. The screening evaluation procedure provides direct evidence for highlighting the condition-specific regulatory network in biological system [91, 134, 144]. Based on the available or predefined regulatory networks, the consistency between architecture and expression are measured, and the most rational network structure with the expression data can be revealed [145, 146]. In the evaluation strategy, each of the reference networks is assessed by measuring the correspondence between network structures and gene expression profiles. The comparison of matching significance in these knowledge-based regulatory networks can identify the responsive regulatory networks of certain conditions and phenotypes.

Network structure determines the regulatory functionality and robustness [147, 148]. The new forward-like engineering of matching network structure with gene expression data provides more alternatives to investigate the regulatory relationships. The original paper in this direction was published in [144]. The authors proposed a Gaussian graphical model to represent the causal relationships of regulatory network architecture and defined a graph consistency probability to measure the goodness of fitting between network and data. However the directed acyclic graph assumption limits its generality and applicability. Collaborating with the senior author of the original work, we introduced a DBN model to handle general regulatory networks [134]. Specifically, by recursive factorization, the joint probability distribution of a certain directed network architecture is represented as a product of the individual density functions conditioned on their parent variables [134, 144], i.e.

in graph . Let be the gene expression of genes at time point . Thus, for , under the first-order Markovian assumption that is independent of for given , we have in the time course data. Assume

,

where

where is the regulatory coefficient of is the error vector and with . According to linear assumption [125, 134], the log-likelihood function

.

Although the binary regulatory relationship between gene and gene is available in the priori, the details of activation ( ), repression ( ), no regulation ( ), as well as the Level IV of regulatory strength are unknown, especially in specific conditions. So we employed a quadratic programming (QP) to calculate the likelihood value by optimizing the coefficients , i.e.,

The constraints in the QP represent the regulatory strength between and . Based on the log-likelihood value, the significance of a network architecture was evaluated by a random sampling process [134, 145, 146]. For each regulatory network, we randomly generated (e.g. 2000) networks by rewiring the same number of regulations in the nodes of the evaluating network. An empirical p-value is calculated to evaluate its statistical significance, i.e.,

,

where is a random network, is the maximum log-likelihood value of the random network and the evaluating network . The evaluation provides a powerful alternative to identify responsive regulatory networks in certain dynamics of environment and condition [134].

Apparently, the knowledge-based regulatory relationships among these genes are not complete and the reference network library should be as complete as possible. To the ends, [149] and [150] have developed methods to integrate inference and evaluation in the same framework by completing the gene network with modifications so that the resultant network achieves more consistency with the gene expression data. The missing regulations can be identified from initial incomplete prior network. Due to the difficulties of pure data-driven inference of regulatory network, the alternatives of combining prior knowledge and evaluating prior gene regulations show promising research directions to investigate transcriptional regulatory network from gene expression data [134].

DISCUSSION AND CONCLUSION

In this review, we summarized the state-of-the-art methods of reverse engineering transcriptional regulatory networks from gene expression data and categorized them into several general frameworks, i.e., correlation-based methods, Boolean network methods, Bayesian network methods, differential equation methods and knowledge-based methods. (Table ) lists these strategies and their typical methods. Some methods implement hybrid models and employ several computational techniques to reversely engineer regulatory networks [41, 83]. These methods such as REVEAL [105], BC3NET [151] and GENIE3 [152] can be classified into multiple categories. For simplicity, we only categorized

them into one of them. For instance, REVEAL also employs mutual information technique beyond Boolean network, so it can also belong to correlation-based methods. In (Table ), some methods such as TRNInfer [153] and Inferelator [154] reconstruct the four levels of transcriptional regulatory relationships, while others such as PCA-CMI [60] and WGCNA [68] generally identify gene regulations without direction information.

After the emergence of high-throughput microarray techniques, great efforts have been undertaken to infer transcriptional regulatory networks from gene expression profiles. Because of the complexity of gene regulations, it is still a challenging task to infer genome-wide regulatory networks from expression data by mathematical modeling [179]. Various computational methods have been proposed to interpret gene expression data and decipher the regulation mechanism of controlling gene expression. The reviewed methods are very useful for providing the quantitative models of harnessing the perturbation and time series of gene expression datasets and identifying the causal relationship of transcriptional regulations. In turn, the endeavors of coupling the regulatory interaction between genes imply the paramount importance of gene regulations in the study of genomics and genetics. It is difficult to assess these methods and select the best one that supersedes all the others by defining some benchmark standards [23]. The details and assumptions in the modeling of real regulation systems as well as the gene expressions in specific conditions and phenotypes determine the superiority of each method. The simple model as Boolean network can reveal critical implications in transcriptional regulation systems [107].

Besides the methods reviewed above, some other methods such as supervised learning [176, 177], feature selection [152, 178], neural network [174] and Petri net [175] methods have also been proposed to address the problem of learning transcriptional regulatory network from gene expression data. Most of these miscellaneous methods are heuristic for mining the relationship between genes from expression profiles. The availability next generation sequencing (NGS) technologies, e.g., RNA-Seq [180], can generate transcriptomic data of higher quality. Theoretically, these reviewed methods can be extended easily to reverse engineering transcriptional regulatory networks from RNA-Seq data. At the same time, identification of the causal regulatory mechanism of gene expression dynamics from gene expression data is constrained from the assumptions and approximations in the models. For instance, time delay between the activation of a TF and its downstream target genes widely exists in the regulatory relationships, which has not been well considered in the available methods [123, 181]. Also, the dynamics of regulation has not been modeled sufficiently, i.e., the regulation strength between TF and targets are always time-varying with temporal features [126]. The reviewed methods can be extended to integrate these important regulatory features into the models of causal regulatory relationships between genes. The reverse engineering methods will become more and more sophisticated for modeling transcriptional regulatory systems as comprehensively as possible.

Beyond utilizing gene expression data, an important research direction in building transcriptional regulatory network is to predict the interaction between TF protein and DNA by machine learning methods. Currently, one of the most important problems in the predictions is how to effectively formulate a biological sequence with a discrete model or a vector, yet still keep considerable sequence order information. All the existing operation engines, such as covariance discriminant (CD) [182, 183], neural network [184], support vector machine (SVM) [185, 186], random forest [187], conditional random field [188], nearest neighbor (NN) [189]; OET-KNN [190], Fuzzy K-nearest neighbor [191, 192], ML-KNN algorithm [193], and SLLE algorithm [183], can only handle vector but not sequence samples. However, a vector defined in a discrete model may completely lose all of the sequence-order information [194]. To avoid completely losing the sequence-order information for proteins, the pseudo amino acid composition or Chou’s PseAAC was proposed [194, 195]. Ever since the concept of PseAAC was proposed in 2001 [194], the approach of representing protein/peptide sequences has been widely used in all the areas of computational genomics [196]. Moreover, the concept of PseKNC (Pseudo K-tuple Nucleotide Composition) was introduced to deal with DNA/RNA sequences in computational genomics, such as for identifying nucleosome positioning [182] and predicting recombination spots [197]. The development of PseAAC for protein sequences and PseKNC for DNA sequences will highly facilitate the prediction of transcriptional regulatory interactions between TF protein and DNA only from sequence information [187, 198, 199].

The challenges of reverse engineering are not only from the information availability, but also from the complexity of regulation system [200]. The measured gene expression levels are not merely determined by the activity of its transcriptional regulators. Post-transcriptional regulations (microRNA silencing [51]) as well as epigenetic modifications (DNA methylation [48] and histone modification [201]) on the gene sequences also highly affect the levels of gene expression. The sequential and combinatorial regulations of gene expression among epigenetic factors, TFs, microRNAs should be considered systematically in reverse engineering regulatory systems when these genomic datasets are available [202]. In the future, the heterogeneous regulatory system with multiple genetic and epigenetic factors should be modeled to integrate transcriptional and post-transcriptional regulations. The integration of genomics, transcriptomics, proteomics datasets, such as ChIP-Seq, protein-binding motifs, gene expression, miRNA abundance, ratios of DNA methylation and chromatin modification, and prior knowledge of regulation, from multiple levels and various aspects of gene regulations provides a possible solution to reconstruct context-specific gene regulations [47]. The networks inferred from various levels can crossly validate each other for accurately identifying gene regulations underlying the whole system. Furthermore, the contradicted identifications in these inferences should be analyzed carefully. They might be caused by the noisy datasets, unrevealed regulatory mechanisms, and specific phenotype associations. Reverse engineering of transcriptional regulatory networks by integrating multiple datasets is a very important research direction [41, 153]. Consistent regulatory relationships at multiple levels shed a brilliant light on the gene expression dynamics in response to various internal signals and external stimuli.

In conclusion, a genome-wide inference of transcriptional regulatory networks from gene expression data provides a promising way to decipher the large-scale causal regulatory relationships among genes. Model-based computational methods of harnessing genomic data facilitate the discovery and revolutionize the research of gene regulation. We summarized the advantages and commented on the improvement possibilities of addressing the disadvantages of these methods individually. The assumptions of modeling the spatial and temporal gene regulations will become more and more reasonable with the accumulation of knowledge about gene regulations. The models will also become more and more close to the real complexity of gene regulation when we obtain better gene expression data with enough sample size and dedicated experiment design, multilevel biological processes, higher quality of expression signals, and systematic perspectives. Knowledge-based methods of integrating existing priori and gene expression seem to be powerful and flexible to decipher the genuine transcriptional control circuits in regulatory systems.

Table 1.

Some available strategies and their representative methods for inferring regulatory networks from gene expression profiles. Their supporting websites and original publications are also shown. Some R packages (http://cran.r-project.org) for Bayesian learning and differential equation parameter identification are also shown. In each category, the methods are ordered alphabetically.

Category	Method	Website	Reference
Correlation-based methods	ANOVA	http://www2.bio.ifi.lmu.de/˜kueffner/anova.tar.gz	[155]
	ARACNE	http://wiki.c2b2.columbia.edu/califanolab/ index.php/Software/ARACNE	[56, 73]
	CLR	http://cran.r-project.org/web/packages/parmigene	[74]
	C3NET	http://cran.r-project.org/web/packages/c3net/index.html	[71]
	GLMNET	http://cran.r-project.org/web/packages/glmnet/	[99, 100]
	grangerTlasso	http://www.biostat.washington.edu/~ashojaie/	[103]
	MINET	http://cran.r-project.org/web/packages/minet/	[76]
	MRNET	http://penglab.janelia.org/proj/mRMR/	[75]
	ParCorA	http://www.comp-sys-bio.org/software.html	[88]
	PCA-CMI	http://csb.shu.edu.cn/subweb/grn.htm	[60]
	Relevance Network	http://buttelab.stanford.edu/start	[65, 70]
	Schafer and Strimmer	http://strimmerlab.org/software.html	[89]
	Simone	http://cran.r-project.org/web/packages/simone/	[156]
	Stuart et al.	http://cmgm.stanford.edu/~kimlab/multispecies/	[64]
	WGCNA	http://labs.genetics.ucla.edu/horvath/ CoexpressionNetwork/Rpackages/ WGCNA	[66, 68]
Boolean network methods	Akutsu et al.	http://www.bic.kyoto-u.ac.jp/takutsu/members/takutsu/	[106, 109]
	Antelope	http://turing.iimas.unam.mx:8080/AntelopeWEB/content/about.jsp	[157]
	BoolNet	http://cran.r-project.org/web/packages/BoolNet	[158]
	BooleSim	https://github.com/matthiasbock/BooleSim	[159]
	Handorf and Klipp	http://code.google.com/p/libscopes/wiki/Paper2011	[160]
	Modent	http://acgt.cs.tau.ac.il/modent/	[161]
	REVEAL	Not available	[105]
	Shmulevich et al.	http://shmulevich.systemsbiology.net/	[112, 113, 115]
Bayesian network methods	ARTIVA	http://cran.r-project.org/web/packages/ARTIVA/index.html	[162]
	BC3NET	http://cran.r-project.org/web/packages/bc3net/index.html	[151]
	Beal et al.	http://www.cse.buffalo.edu/faculty/mbeal/	[125]
	BNFinder	http://bioputer.mimuw.edu.pl/software/bnf	[163]
	BNLEARN	http://cran.r-project.org/web/packages/bnlearn	[164]
	BNT	http://code.google.com/p/bnt/	[120]
	Frideman et al.	http://www.cs.huji.ac.il/labs/compbio/expression/	[117, 119]
	GeneNet	http://cran.r-project.org/web/packages/GeneNet	[83]
	G1DBN	http://cran.r-project.org/web/packages/G1DBN/index.html	[165]
	GlobalMIT	https://code.google.com/p/globalmit	[166]
	Module network	http://ai.stanford.edu/~erans/module_nets/	[84]
	TESLA	http://sailing.cs.cmu.edu/tesla/index.html	[126]
	SSM	http://www.chems.msu.edu/groups/chan/ssm.zip	[167]
Differential equation methods	Chen et al.	Not available	[127]
	deSolve	http://cran.r-project.org/web/packages/deSolve	[168]
	D'haeseleer et al.	Not available	[95]
	D-NetWeaver	https://cbim.urmc.rochester.edu/software/d-netweaver/	[128, 129]
	GRNInfer	http://doc.aporc.org/wiki/Software	[41]
	Inferelator	http://bonneaulab.bio.nyu.edu/software.html	[154]
	Tegner et al.	http://www.bu.edu/bme/people/primary/collins/	[34, 36]
	TRNInfer	http://www.sysbio.ac.cn/cb/chenlab/software.htm	[153]
	Wahde and Hertz	http://www.nbi.dk/~hertz/	[169]
Knowledge-based methods	Banjo	http://www.cs.duke.edu/~amink/software/banjo	[142]
	BNP	http://research.bioe.bilgi.edu.tr/bnp/	[170]
	Greenfield et al.	http://bonneaulab.bio.nyu.edu/software.html	[171]
	Hill et al.	http://mukherjeelab.nki.nl/DBN	[172]
	Linear programming	http://doc.aporc.org/wiki/Software	[138]
	Liu et al.	http://doc.aporc.org/wiki/Software	[134]
	Network energy	Not available	[139, 143]
	Network Screening	http://www.molprof.jp/~horimoto/	[144]
	PLASSO	http://nba.uth.tmc.edu/homepage/liu/pLasso	[173]
Miscellaneous methods	GENIE3	http://homepages.inf.ed.ac.uk/vhuynht/software.html	[152]
	Neural network	http://www.me.chalmers.se/~mwahde	[174]
	Petri net	http://dnagarden.hgc.jp/en/doku.php/software	[175]
	Supervised learning	http://cbio.ensmp.fr/sirene	[167, 177]
	TIGRESS	http://cbio.ensmp.fr/tigress	[178]
Differential equation methods	Chen et al.	Not available	[127]
	deSolve	http://cran.r-project.org/web/packages/deSolve	[168]
	D'haeseleer et al.	Not available	[95]
	D-NetWeaver	https://cbim.urmc.rochester.edu/software/d-netweaver/	[128, 129]
	GRNInfer	http://doc.aporc.org/wiki/Software	[41]
	Inferelator	http://bonneaulab.bio.nyu.edu/software.html	[154]
	Tegner et al.	http://www.bu.edu/bme/people/primary/collins/	[34, 36]
	TRNInfer	http://www.sysbio.ac.cn/cb/chenlab/software.htm	[153]
	Wahde and Hertz	http://www.nbi.dk/~hertz/	[169]
Knowledge-based methods	Banjo	http://www.cs.duke.edu/~amink/software/banjo	[142]
	BNP	http://research.bioe.bilgi.edu.tr/bnp/	[170]
	Greenfield et al.	http://bonneaulab.bio.nyu.edu/software.html	[171]
	Hill et al.	http://mukherjeelab.nki.nl/DBN	[172]
	Linear programming	http://doc.aporc.org/wiki/Software	[138]
	Liu et al.	http://doc.aporc.org/wiki/Software	[134]
	Network energy	Not available	[139, 143]
	Network Screening	http://www.molprof.jp/~horimoto/	[144]
	PLASSO	http://nba.uth.tmc.edu/homepage/liu/pLasso	[173]
Miscellaneous methods	GENIE3	http://homepages.inf.ed.ac.uk/vhuynht/software.html	[152]
	Neural network	http://www.me.chalmers.se/~mwahde	[174]
	Petri net	http://dnagarden.hgc.jp/en/doku.php/software	[175]
	Supervised learning	http://cbio.ensmp.fr/sirene	[176, 177]
	TIGRESS	http://cbio.ensmp.fr/tigress	[178]