A Relative Variation-Based Method to Unraveling Gene Regulatory Networks

Gene regulatory network (GRN) reconstruction is essential in understanding the functioning and pathology of a biological system. Extensive models and algorithms have been developed to unravel a GRN. The DREAM project aims to clarify both advantages and disadvantages of these methods from an application viewpoint. An interesting yet surprising observation is that compared with complicated methods like those based on nonlinear differential equations, etc., methods based on a simple statistics, such as the so-called -score, usually perform better. A fundamental problem with the -score, however, is that direct and indirect regulations can not be easily distinguished. To overcome this drawback, a relative expression level variation (RELV) based GRN inference algorithm is suggested in this paper, which consists of three major steps. Firstly, on the basis of wild type and single gene knockout/knockdown experimental data, the magnitude of RELV of a gene is estimated. Secondly, probability for the existence of a direct regulation from a perturbed gene to a measured gene is estimated, which is further utilized to estimate whether a gene can be regulated by other genes. Finally, the normalized RELVs are modified to make genes with an estimated zero in-degree have smaller RELVs in magnitude than the other genes, which is used afterwards in queuing possibilities of the existence of direct regulations among genes and therefore leads to an estimate on the GRN topology. This method can in principle avoid the so-called cascade errors under certain situations. Computational results with the Size 100 sub-challenges of DREAM3 and DREAM4 show that, compared with the -score based method, prediction performances can be substantially improved, especially the AUPR specification. Moreover, it can even outperform the best team of both DREAM3 and DREAM4. Furthermore, the high precision of the obtained most reliable predictions shows that the suggested algorithm may be very helpful in guiding biological experiment designs.


Introduction
In the post-genomic era, one of the fundamental tasks is reconstructing gene regulatory networks (GRN) from experimental data and other a priori information. It is hoped that this reconstruction is helpful in both understanding cell functions and gaining additional insights about the processes of some complicated diseases that might lead to new target gene discovery. Recently, with the development of high-throughput technologies, such as DNA microarrays and mass spectroscopy, etc., it becomes possible to simultaneously collect thousands of gene expression data [1,2]. Stimulated by these technology advancements, a variety of different models and methods have been proposed for GRN reconstruction, such as Boolean networks [3,4], Bayesian networks [5,6], information theory based algorithms [7][8][9][10], ordinary differential equation (ODE) based methods [11][12][13], etc. In addition, some software packages, such as GeneNet, minet, etc., have been developed [8,9].
A challenge common to all these reverse-engineering methods is that in comparison with the dimension and complexity of a GRN, the collected experimental data is generally with a low SNR (signal-to-noise ratio) and the number of observations is not very large in every experiment. Another challenge is to evaluate the appropriateness of the assumptions adopted by these methods. To settle these problems, the Dialogue for Reverse Engineering Assessments and Methods (DREAM) project recently provided a set of benchmark networks that can be used to compare both advantages and disadvantages of different GRN topology inference methods [14][15][16][17]. Compared with other benchmark networks, one of the most attractive characteristics of the networks provided by the DREAM project is that they are extracted from actual biological networks and are able to represent some most important and typical biological modules. By far, it has become one of the most widely used benchmarks for GRN topology inference.
Several methods have been shown to be effective in inferring the structure of a GRN through participating the DREAM project. For example, the best performer of the DREAM3 subchallenges took an approach that firstly learns some Gaussian noise models from knockout experimental data, and then combines these results with those obtained through fitting time series experimental data to an ODE model [18]. The second place team of the DREAM3 Size 100 subchallenges utilized a mutual information (MI) based method and a so-called Inferelator 1.0 method, which takes sparsity of a GRN into account through penalizing the l 1 {norm of the kinetic parameter vector of the ODE model [19]. On the other hand, the DREAM project organizers applied the so-called Z-score to measure possibilities of the existence of a direct regulation from one gene to another gene [16]. Surprisingly but also interestingly, this simple method was proved to be placed at respectively the first (tie) and the third for the Size 100 subchallenges of DREAM3 and DREAM4.
From a statistical point of view, the Z-score based method is actually a t-test [20]. More precisely, to determine whether gene j has a direct regulation on gene i, it utilizes the absolute expression level variation (AELV) of gene i from the wild type after a perturbation on gene j. The larger the magnitude of this AELV is, the more unlikely that the change is due to measurement noise, and thus the larger the probability that gene i is directly regulated by gene j. This AELV, however, is sometimes not very effective in distinguishing a direct regulation from an indirect one, as possibilities can hardly be excluded that an indirect regulation causes an AELV larger in magnitude than some direct regulations [21]. To reduce estimation errors caused by indirect regulations, which is often called cascade errors [21,22], the best performer of the Size 100 subchallenges of DREAM4 suggested to refine the results of the Z-score based method through down ranking some feedforward edges [22]. Basically, the idea behind this treatment is to remove every direct regulation between two different genes in a GRN estimate, provided that it does not belong to a cycle and there exists another direct or indirect regulation between these two genes. This procedure has significantly improved the adopted estimation specifications, and therefore shown its efficiency in GRN topology estimations [22].
The results of [22] are encouraging. It seems, however, that further efforts are still required to make the estimation procedure applicable to practical problems, noting that as reported in [22], its prediction accuracy for some networks is still not very high, and thresholds exist that are significantly different from the recommended one but are capable of leading to a much better network structure estimate. In addition to this, our computational experiences with this method show that its precision-recall (PR) curve is still not very satisfactory for some networks. A detailed discussion on this issue is provided in the subsection of Further Discussions of the section of Results and Discussions.
To achieve a better GRN structure estimation, an innovative technique is proposed in this paper for GRN topology inference. The ideas behind the developed algorithm are relatively simple. That is, rather than absolute change, relative variation of gene expression level is adopted in measuring possibilities of the existence of a direct regulation between two different genes of a GRN. This algorithm consists of three major steps. That is, magnitude estimation and normalization of the relative variations, estimation of genes not regulated by other genes, modification of the normalized estimate for the magnitude of the relative variations and GRN topology identification. In the first step, relative expression level variation (RELV) of a gene is estimated using experimental data before and after another gene of the same GRN has been perturbed. This estimate is further normalized to reflect effort differences of regulating distinct genes. In the second step, on the basis of the estimated probability that the magnitude of the RELV of a gene is greater than a prescribed value, genes with a zero in-degree are estimated. Finally, in the third step, every normalized magnitude of the AELV of a gene with an estimated nonzero in-degree is adjusted to be greater than those with an estimated zero in-degree. Computational experiences with the Size 100 network inference subchallenges of both DREAM3 and DREAM4, as well as some other simulated large scale GRNs, show that this method can significantly outperform not only the Zscore based method, but also the best teams who utilized an integration of several widely adopted methods. The suggested method has also been integrated with the so-called down ranking method, which is recommended by the best network inference team of DREAM4. Once again, it has been confirmed through actual computations that this method is helpful in reducing cascade errors. The corresponding improvement, however, is not as significant as that to the Z-score based method. This means that some cascade errors have been reduced by the suggested method, which confirms from another aspect that the suggested method is really effective in distinguishing direct and indirect regulations of a GRN.
The outline of this paper is as follows. At first, the relative variation based estimation algorithm is illustrated. A technique is also provided that can integrate GRN topology prediction results using respectively steady state knockdown and knockout experimental data, as well as a procedure that integrates the method suggested in this paper with the so-called down ranking method. Afterwards, the proposed estimation method is assessed using the data sets of the Size 100 subchallenges of both DREAM3 and DREAM4. Variations of estimation performances with respect to parameters of the suggested method have also been reported, as well as estimation results using both steady state knockdown and knockout experimental data. In addition, estimation results are also given in which the suggested method is integrated with the socalled down ranking method. Finally, some concluding remarks are given about characteristics of the suggested method, as well as some future works worthy of further efforts.

Materials and Methods
Concerning a GRN with n genes, assume that measurement errors affect experimental data in an additive way, as well as that measurement errors with the expression level of gene i have an independent and identical normal distribution N(0,s 2 i ). Let x ji and x ½0 ji represent respectively the observed and the actual gene expression levels of gene i when gene j is knocked out or knocked down, and e ji the corresponding measurement error. Then, it is obvious from these representations that Moreover, denote by x ½wt i and x ½wt,0 i respectively the observed and the actual expression levels of gene i in the wild type, and g ji the steady expression level variation of gene i after the knockout/ knockdown of gene j. Then, from its definition, we have that g ji~x Define the RELV (relative expression level variation) of gene i resulted from a perturbation on gene j, denote it by d ji , as Note that in a GRN, every gene can usually be approximately assumed to be in one of the following two states, expressed and unexpressed states. Moreover, expression levels, that is, concentrations of the corresponding proteins or mRNAs, etc., of distinct genes usually take very different values, and sometimes these values may even have different orders [1,11,23]. These imply that when a gene is knocked out or knocked down, absolute variations of the expression levels of genes regulated by this externally perturbed gene may have very different magnitudes. These characteristics of g ji are not very attractive in GRN topology estimation, as they imply that the magnitude of g ji due to an indirect regulation may sometimes be significantly larger than that due to a direct regulation. On the other hand, when the RELV is utilized, the aforementioned problems can be partly overcome. More specifically, RELVs of every gene are roughly of the same order, which makes their comparisons more biologically significant than the AELV that is important in GRN structure estimations. In addition, if in a pathway of GRNs, every direct regulation, say, that from gene j to gene i, make a relative variation of the concentrations of the proteins or mRNAs, etc., related to the regulated gene i at most as large as that of the regulation gene j, then, it is obvious that in this pathway, the magnitude of every RELV due to an indirect regulation is certainly not greater than that due to a direct regulation. From these considerations, it appears that RELV is more attractive than AELV in GRN topology estimations.
But it is worthwhile to emphasize that self-activation usually exists in GRNs [23,24], which may make the RELVs of a pathway be amplified during cascade gene connections. This means that the aforementioned assumption may not be satisfied by every pathway of a GRN. To make things worse, for some genes of a GRN, there exist more than 1 directed pathways from one gene to another gene [14,24]. For example, when gene j directly regulates gene i (through its proteins), it is possible that gene k is also directly regulated by gene j and gene k directly regulates gene i further. Under such a situation, when gene j is externally perturbed, the RELV of gene i is due to both direct and indirect regulations. If one of these two regulations has an activation effect and the other has a suppression effect, then, their composite effects may significantly weaken that of the direct regulation and therefore result in an incorrect estimate using the above assumption, noting that generally, the RELV of a gene should be estimated from experimental data.
The above arguments show that although the RELV of a gene has some attractive properties in GRN topology estimations, it is still not very clear whether or not the adopted assumption is reasonable for most pathways of a GRN from a biological viewpoint. It appears, however, from our computational experiences, some of which are reported in the section of Results and Discussions, that this assumption may have some nice biological interpretations and is satisfied by many regulations existent in a GRN.
These arguments mean that relative concentration variation is, at least under many situations, able to differentiate direct and indirect regulations of a GRN, and is therefore more effective than the absolute one in GRN topology estimations, noting that indirect regulations are rich in a GRN. These arguments also imply that the larger the magnitude of d ji is, the more unlikely that the expression level variation of gene i after perturbing gene j is due to indirect regulations, and thus the larger the probability that gene i is directly regulated by gene j.

RELV Estimation
As x ½0 ji is generally not available from experiments, an estimate for d ji should be used in GRN topology inference. To obtain this estimate, the set to which d ji belongs most likely with a fixed probability is considered. Recalling that e ji is assumed to have a normal distribution N(0,s 2 i ), this is equivalent to compute the minimal interval that contains an estimate of d ji when the measurement noise e ji is assumed to be not greater in magnitude than l 1 s i for a fixed non-negative l 1 . On the basis of this observation and Equation (2), as well as the fact that both x ½wt,0 i and s i are always positive, it can be directly shown that Therefore, Note that in GRN topology inference, the sequence of probability of the existence of a direct regulation from one gene to another gene plays the most essential role. Moreover, it has been argued that the larger the absolute value of d ji , the higher the probability that gene i is directly regulated by gene j. Based on these considerations, the element with the maximal magnitude of the d ji s satisfying Equation (5) is taken as its estimate. Denote this estimate byd d ji . Then, from Equation (6), we have that and its value can be calculated from experimental data, provided that both s i and x ½wt,0 i are known. In practical applications, however, these two parameters are generally not available and should also be estimated from experimental data. If the set of genes that do not affect gene i, denote it by ND i , is known, then, some widely adopted estimates for s i and x ½wt,0 i are respectively as in which #( : ) stands for the element number of a set [25,26]. However, the set ND i is usually unknown before GRN topology inference, which invalidates adoption of the above estimates. On the other hand, it is now widely recognized that a large scale GRN usually has a sparse topology [11,13,22], which means that for most genes of a large scale GRN, #(ND i ) is very close to n which stands for the number of its genes. This means that in estimating s i and x ½wt,0 i , essential differences will not arise for most genes even if ND i is taken to be the whole set of the genes of a GRN. Based on these considerations, the following estimates are adopted in this paper for s i and x ½wt,0 i , in which differences have not been taken into account between a measurement for the expression level of gene i in the wild-type and those when some other genes have been knocked out and/or knocked down.
Although these estimates may be crude, they are widely adopted in GRN topology estimation and are capable of leading to a good network estimate [16,18]

RELV Normalization
Define a n|n dimensional matrix D with its j-th row i-th column element being the estimate of Dd d ji D when i=j and its diagonal element being zero, and denote its i-th column vector by d d d d i . Then, the above derivations make it clear that this matrix contains information about the probability of the existence of a direct regulation between any two different genes in a GRN. However, to infer the structure of a GRN from the aforementioned matrix D, an important fact must be taken into account. That is, in a GRN, some genes may be easily regulated by other genes, while regulations on some other genes may need more efforts [23,24]. As a matter of fact, even when every direct regulation of a pathway in a GRN satisfies that RELV of the regulated gene is not greater than that of a regulation gene, it is still possible that direct regulations to different genes lead to different magnitudes of these variations of the regulated genes, although under such a situation, a direct regulation of this pathway certainly makes the RELV of the regulated gene have a magnitude not smaller than that of every indirect one. Therefore, in order to obtain a good estimate from the matrix D about the topology of a GRN, an appropriate normalization is still required for the estimatedd d ji s among different genes.
Although it is still not very clear how to make a biologically significant normalization among the RELVs of different genes, as a primary study, it is suggested in this paper to use the p-norm of the vector d d d d i and the geometric average of its non-zero elements to achieve this objective, which are widely adopted in many fields like system analysis and synthesis, signal processing, etc., and have shown their efficacy [13,25,27]. While their effectiveness in GRN topology estimation is still not very clear theoretically, our computational experiences, part of them are given in the section of Results and Discussions of this paper, show that they are able to lead to an estimate much better than that without normalizations. More specifically, when the p-norm and the geometric average are respectively used in this normalization, the j-th row i-th column element of the matrix D, that is, It is worthwhile to note that this normalization does not change the diagonal elements, which is important as self regulation can hardly be identified from either knockout or knockdown experimental data. For presentation conciseness, the normalized matrix D using the vector p-norm and the geometric average is denoted respectively by D ½p and D ½g in the rest of this paper.

Estimation of Genes with a Zero In-degree
While normalization is helpful in balancing RELVs among different genes, another problem arises in GRN topology estimations. That is, this normalization usually leads to wrong estimates in network inference for genes not regulated by other genes. This is because that although for these genes, the corresponding computed d d d d ji s are generally very small, some of their normalized values may have a comparable magnitude with those of a gene that is regulated by other genes. Note that in a GRN, nodes with a zero in-degree, that is, genes that can not be regulated by other genes, extensively exist [16,23]. Therefore, special cautions must be taken to deal with them in inferring the structure of a GRN.
To distinguish genes that can be and can not be regulated by other genes, once again RELVs of a gene are considered when another gene is knocked out and/or knocked down, but in a different way. It is worthwhile to point that in principle, it is also possible to use Dd d ji D in estimating genes that are with a zero indegree. However, actual computations show that the corresponding estimate is not as effective as the following estimate. More precisely, for a prescribed l 2 [(0, 1, it is reasonable to regard that there exists a direct regulation from gene j to gene i when Dx . Otherwise, gene i is considered not to be directly regulated by gene j. As x ½0 ji can hardly be estimated with an acceptable accuracy from experimental data in actual GRN structure identification, probability is taken as a measure for the existence of a direct regulation from gene j to gene i. Recall that measurement errors are assumed to affect experimental data additively. From Equation (2) and the definition of g ji , it is obvious that Dx On the other hand, note that x ½0 ji stands for the actual expression level of gene i when gene j is externally perturbed, and therefore can not take a negative value. It is straightforward from this fact and Equation (1) that the following inequality should always be satisfied by the measurement error e ji .
Summarizing Equations (13) and (14), it can be declared that to guarantee the existence of a direct regulation from gene j to gene i, it is necessary and sufficient that Noting that measurement errors are also assumed to have a normal distribution N(0,s 2 i ), the above equation makes the probability computable for the existence of a direct regulation from gene j to gene i, provided that both s i andx x ½wt,0 i are available. In practical inference, s i and x ½wt,0 i are usually replaced by their estimatesŝ s i andx x ½wt,0 i that are provided in Equation (9). Denote the corresponding calculated probability byP P ji . Then, the above arguments make it clear that Moreover, the larger the P ji is, the higher the confidence is that gene i is directly regulated by gene j. LetP P i represent the maximum value ofP P ji when j varies over all the integers between 1 and n except i, that is,P P i~m ax 1ƒjƒn, j=iP P ji . The above arguments imply that ifP P i takes a small value, then, it is very possible that gene i is not regulated by any other genes of the GRN. In other words, the in-degree of this gene is equal to zero with a high probability.
To estimate genes that has a zero in-degree, both absolute value and relative largeness ofP P i are considered. This can make all the genes with an estimated zero in-degree have an estimate of the probability of being regulated by other genes not only very small, but also significantly smaller than that of every gene with an estimated nonzero in degree. More precisely, rearrangeP P i in an increasing order, that is,P P i i1 ƒP P i i2 ƒ Á Á Á ƒP P i in in which i i k = i i j when k=j and i i k [f1,2, Á Á Á ,ng. Define r k as r k~P P i i kz1 Let m denote the first integer with which r m takes the greatest value under the condition thatP P i im belongs to ½P min , P max for some prescribed P min and P max , then, all genes numbered as i i 1 , i i 2 , Á Á Á, and i i m are regarded not to be regulated by any other genes of the GRN.

RELV Magnitude Modification
With the above estimate about genes of a GRN that have a zero in-degree, the normalized RELV matrices D ½p and D ½g are modified, in order to get a better estimate about its structure. As genes estimated to be of a zero in-degree generally has a low probability of being regulated by other genes, its corresponding normalized RELVs must be adjusted to have a lower rank than those of genes that might be regulated by other genes. To achieve this objective, define d ½? 0 as the maximal magnitude of the normalized RELVs of the genes estimated to be not regulated by any other genes. That is, in which ? can be either p or g. With this value, the normalized RELVs are modified as follows, This modification makes every RELV of a gene possibly regulated by other genes greater in magnitude than any RELV of a gene estimated to be of a zero in-degree. Denote byD D ½? the n|n dimensional matrix with its j-th row i i kth column element beingd d ½?
j i i k , ?~p or g. Elements of this matrix are directly used to infer the structure of a GRN, according to the principle that the bigger the j-th row i-th element is, the higher the probability is that gene i is directly regulated by gene j.

Estimation Algorithm
In summary, to estimate the structure of a GRN, it is assumed in this paper that measurement errors in the expression levels of every gene have an independent and identical normal distribution, and affect experimental data additively. On the basis of the concept of the RELV of a gene, an algorithm is suggested in this paper for identifying direct regulations of a GRN. This algorithm consists of the following three main steps.
1. Using a prescribed l 1 and the estimates for s i (the standard variance of measurement errors) and x ½wt,0 i (the wild type expression level of gene i), which are given in Equation (9), as well as Equations (10) and (11), calculate the matrices D ½p or D ½g consisting of the normalized magnitudes of the estimates of the RELVs for every gene in a GRN. This is equivalent to construct the matrix D ½p or the matrix D ½g respectively as D ½p~ d d  (9) and (16), as well as a prescribed l 2 , calculate the estimate for the probability that gene i is directly regulated by gene j, that is,P P ji . ComputeP P i aŝ P P i~m ax 1ƒjƒn, j=iP P ji . RearrangeP P i into a monotonically increasing sequenceP P i i1 ƒP P i i2 ƒ Á Á Á ƒP P i in . Using some prescribed thresholds for the minimum and maximum ofP P i , say, P min and P max , determine the gene i i m that has aP P i im belonging to ½P min , P max and makesP P i i mz1 P P i im reach its maximum in the first time. Designate the in-degree of genes i i k , 1ƒkƒm, to be zero.

On the basis of Equations
3. Modify the matrices D ½p or D ½g according to Equations (17) and (18). Using elements of these modified matrices, queue possibilities of the existence of a direct regulation from the gene with the same number of the row to the gene with the same number of the column. The bigger this element is, the higher the confidence is for the existence.

Integration of Knockout and Knockdown Data
Currently, biological experiments can provide both steady state data and time series data. In addition, an experiment can be performed through knocking out a single gene, knocking down a single gene or simultaneously perturbing several genes. In this subsection, a method is proposed to integrate estimation results obtained respectively from steady state knockout and knockdown experimental data. Rather than to develop an efficient integration method, the purposes of this investigation are mainly to clarify characteristics of these types of experimental data in GRN topology estimations with the method suggested in this paper.
The integration method suggested in this paper is similar to cross validations [25,26]. That is, if a consistent estimate is obtained from different experimental data, then, confidence is strengthened about the correctness of this estimate. More specifically, if from both the knockout and the knockdown experimental data, the estimated magnitude of the RELV corresponding to a possible direct regulation has a large value, then, confidence about the existence of this direct regulation is increased. As a large scale GRN usually has a sparse topology, it appears reasonable to only modify a few RELVs for every gene in this integration. Moreover, as knockout experimental data is widely believed to be more informative than knockdown experimental data in GRN topology identification, for example, observations from the reported results of the DREAM project show that estimation performances with knockdown experimental data are usually significantly worse than those with knockout experimental data [16,22], a higher confidence is given to the knockout experimental data based estimates.
In this paper, it is suggested to modify the first 5 biggest RELVs of a gene obtained from knockout experimental data, provided that this gene is estimated to have a nonzero in-degree. More precisely, letD D ½?,ko andD D ½?,kd denote the n|n dimensional matrices consisting respectively of the modified normalized RELVs obtained from knockout and knockdown experimental data,d d their j-th row i-th column elements. Here, ? can be either p or g, and i,j~1,2, Á Á Á ,n. If a gene, say, gene i, is estimated to have a nonzero in-degree using the knockout experimental data, which is equivalent to i= [ i i 1 , i i 2 , Á Á Á , i i m f g , then, its modified normalized RELVs will be further adjusted according to the following procedures.

Integration with the Down Ranking Method
In reducing false positive errors in GRN topology inference, the so-called down ranking method has been proved to be very effective [22]. While the objectives of this method are almost the same as those of the algorithm suggested in this paper, different approaches have been utilized. Briefly, in the down ranking algorithm, it is assumed that an a priori estimation about the topology of a GRN has been obtained by some methods, and if a direct path between two genes is not in a cycle and there is another direct or indirect path between these two genes in the estimated GRN structure, then, the former direct path should be deleted. This idea has been further extended to the so-called strongly connected components. A detailed description can be found in [22].
In this subsection, a procedure is suggested to integrate the algorithm suggested in this paper with this down ranking method. The major proposes are to see whether performances in GRN topology estimation can be further improved, as well as whether the cascade errors reduced by the down ranking method can also be reduced by the method suggested in this paper. Taking into account characteristics of these two different methods, they are integrated in the following way, in which ? can be either p or g.
1. Compute elements of the matrix D ½? using Equation (11) and knockout experimental data. 2. For a given threshold value, say, c, a matrix D ½?,dr is obtained from the previously obtained matrix D ½? , using the down ranking algorithm. 3. For every i,j~1,2, Á Á Á ,n with i=j, compute the estimateP P ji using Equation (16) that stands for the probability of the existence of a direct regulation on gene i from gene j. Similar to that of the 2nd step of the estimation algorithm, estimate genes with a zero in-degree using these probabilities and some prescribed P min and P max . Modify the matrix D ½?,dr by the same token as that of the matrix D ½? , on the basis of (17) and (18). Denote the modified matrix byD D ½?,dr . 4. Queue possibilities of the existence of a direct regulations in the GRN according to the elements of the matrixD D ½?,dr , in the same way as that without method integrations in which the matrixD D ½p orD D ½g is used.
Note that in the estimation algorithm suggested in this paper, the matrix D ½? , ?~p or g, has been normalized which makes every element of this matrix belong to ½0, 1. This implies that when the down ranking method is applied to the matrix D ½Ã , a meaningful threshold value c should also belong to this interval. This is different from the situation when the Z-score based method is integrated with, in which the computed Z-scores between two different genes may vary in a much larger interval, that leads to a much bigger set for searching the optimal threshold value c.

Data Sets and Assessment Metrics
To illustrate the effectiveness of the developed inference algorithm, tests are performed on the Size 100 Network subchallenges of both DREAM3 and DREAM4, using the data set provided by the organizers. These subchallenges are designed to assess performances of an identification method for the structure of a large scale GRN [16]. They respectively contain 5 different benchmark networks which were obtained through extracting some important and typical modules from actual biological networks. There are three types of experimental data for each subchallenge, which are respectively knockout experimental data, knockdown experimental data, and time series experimental data. Predictions are compared with the actual structure of the networks by the DREAM project organizers using the following two different metrics in topology prediction accuracy evaluations. Moreover, for every network of DREAM3 and DREAM4, the p-values of the AUPR and AUROC specifications, which indicate the probability that random predictions would have the same or better performances, are computed, and finally a score is calculated using these p-values. More specifically, the logarithm of the geometric mean is calculated respectively for both the 5 AUROC p-values and the 5 AUPR p-values, and the score is taken as the absolute value of the average of these two logarithms. More detailed explanations can be found in [15,16] or the web site of the DREAM project at http://wiki.c2b2.columbia.edu/dream/. A larger score indicates a better performance of the adopted inference algorithm.
Noting that the GRN inference algorithms developed in the previous section are applicable only to steady state experimental data, concentrations of this section are focused on knockout and knockdown experimental data. As the suggested estimation algorithm without either data integration or method integration consistently gives much better performances when the knockout experimental data are used for the Size 100 subchallenges of both DREAM3 and DREAM4, which is in a good agreement with other methods reported by the participants of the DREAM project [16,19,22], the corresponding results are at first reported.

Performances for the Knockout Data
Using the knockout experimental data provided by the DREAM project organizer, GRN topology inference is performed for the Size 100 subchallenges of both DREAM3 and DREAM4. To investigate influences of different normalization on the prediction accuracy of the estimation algorithm, p~2 is firstly adopted for the p-norm based normalization, which is widely utilized in fields like system analysis and synthesis, signal processing, etc. [13,25,27]. Moreover, p~3:5 is also utilized which is found to be close to the optimal one for most networks of DREAM3 and DREAM4. In addition, the optimal p is also searched for the p-norm based normalization over the interval ½1, 101 for the Net3 network of DREAM4, and ½1, 11 for all the other networks, through an equally spaced sampling with 100 samples. This is because that actual computations show that for the Net3 network, the AUPR specification does not take its optimal value when the parameter p is restricted to the interval ½1, 11. In fact, it still increases around p~11. In this optimization, the desirable p is selected to be the sample that maximizes the AUPR specification. This is mainly because that due to some precision problems of the score computation method provided by the organizers, the computed p-value of some networks become zero which makes it impossible to compute the score of the corresponding estimation algorithm. These problems can also be understood from the results reported in Table 1, in which several computed p-values are zero. On the other hand, according to our computational experiences, significant improvement on the AUROC specification appears much more difficult. The results are provided in Tables 1 and 2, in which RELV ½2 , RELV ½3:5 , RELV ½pÃ represent respectively the results for the algorithm using the p-norm based normalization with p~2, p~3:5 and the optimal p; while RELV ½g those for the algorithm using the geometric average based normalization. With a little abuse of terminology, in the rest of this paper, these representations are used to indicate the suggested estimation method with the corresponding normalization, in order to avoid awkward statements.
In all these estimations, l 1~0 :01, l 2~0 :5 and P max~0 :2 are utilized. In addition, P min~0 :0005 and P min~0 :05 are respectively adopted for the subchallenges of DREAM3 and DREAM4. To compare prediction performances with the Z-score based method and the best team, the corresponding specifications are also included in these tables. It is worthwhile to note that the estimation accuracy specifications of the best team included here are obtained directly from the web site of the DREAM project. Their digit lengthes are different from the other results that are obtained through actual computations. The best values of the AUROC and the AUPR specifications for each network are written in boldface. In addition, relative performance variation is also provided for each network, immediately below the p-values of the estimation specifications. The first line (RPV-Z) gives results compared with the Z-score based method, while the second line (RPV-B) those with the best team. The averaged relative performance variation (ARPV) is provided immediately after the comparisons for each network. Furthermore, the optimal p for each network is given in parentheses in the last line of the RELV ½pÃ row. In the last column of Table 2, the obtained scores are also given for each method in the same row of their AUPR values. It should be pointed out that in Table 1, due to some technical issues with the software provided by the DREAM project organizers, the score can not be calculated for the best team of DREAM3 and is designated to be ? to facilitate comparisons with other methods, which is resulted from the high value of the AUPR specification for the Yeast3 network. This expression way is also adopted in other tables of this paper. As the AUPR specification of both the Z-score based method and the method suggested in this paper is better than the best team for the Yeast3 network of DREAM3, score comparisons among them are currently impossible and therefore the scores are omitted.
From these computation results, it is clear that although there are some performance differences among RELV ½2 , RELV ½3:5 , RELV ½pÃ and RELV ½g , they all show some accuracy improvements in GRN topology inference over the Z-score based method for most of the subchallenges. Especially, significant performance improvements over the best team of both DREAM3 and DREAM4 can even be seen with the estimation method using the p-norm based normalization. Moreover, improvements on the AUPR specification are more significant than the AUROC specification. On the other hand, it can be seen that when p is fixed to be 3:5, the performance is very close to that of RELV ½pÃ which utilizes the optimal p. But it is worthwhile to emphasize that in actual applications, there are still no methods for estimating the optimal value of the parameter p, which means that estimation performance comparisons with RELV ½pÃ are of little practical values. The purpose of providing results corresponding to the optimal p in these tables are only to make it clear that deviations of this parameter from its optimal value generally does not result in significant estimation performance deteriorations.
Results of Tables 1 and 2 also reveal that normalization indeed plays an important role in improving prediction accuracy. As the optimal parameter p for the p-norm based normalization usually can not be known in actual applications, discussions are concentrated on the results of RELV ½2 , RELV ½3:5 and RELV ½g . The obtained results show that among these three methods, although the 2{norm based normalization is widely utilized in fields like system analysis and synthesis, signal processing, etc., it seems more appropriate to use the 3:5{norm based normalization in GRN topology estimations. With this normalization, the suggested algorithm outperforms the Z-score based method almost in each adopted specification and in every network inference. Improvements in the AUPR specification are particularly evident with the Net1 network and the Net5 network of DREAM4, which are respectively greater than 14% and 13%. On the other hand, although RELV ½2 does not perform as well as RELV ½3:5 , it still yields better results than the Z-score based method for all the DREAM3 and DREAM4 subchallenges. These facts show that compared with AELV, RELV is indeed more RPV-Z: relative performance variation with respect to the Z-score based method; RPV-B: relative performance variation with respect to the best team; ARPV: averaged relative performance variation of the 5 networks.
{ RELV ½pÃ , which stands for the method with the optimal normalization parameter p, generally can not be applied in actual estimations. The purposes to include its inference results here are only to make it clear that significant estimation performance degradation does not occur when the parameter p deviates from its optimal value. *Due to some precision issues of the method suggested by the DREAM project organizers, these p-values can not be distinguished from zero in actual computations, which makes it impossible to compare scores of the adopted GRN topology estimation methods. doi:10.1371/journal.pone.0031194.t001 RPV-Z: relative performance variation with respect to the Z-score based method; RPV-B: relative performance variation with respect to the best team; ARPV: averaged relative performance variation of the 5 networks.

{
The purposes to include the inference results of RELV ½pÃ are completely the same as those of Table 1. That is, to clarify that deviation of the parameter p from its optimal value usually does not lead to significant estimation performance degradations. doi:10.1371/journal.pone.0031194.t002 effective in distinguishing direct and indirect regulations, and therefore reducing the so-called cascade errors in GRN topology inference.
In addition, compared with the best team of DREAM3, although the AUROC specification has become slightly worse for some networks, both RELV ½2 and RELV ½3:5 show improvement in the AUPR specification for every network, and the biggest improvement is greater than 20%. In comparison with the best team of DREAM4, although these two methods occasionally show some great performance decrements, for example, the AUPR specification of RELV ½2 for the Net5 network is about 14% lower than that of the best team, they still yield better results on average in both AUROC and AUPR specifications. According to the report in the web site of the DREAM project for the Size 100 subchallenges of DREAM3, the p{value of the AUPR specification corresponding to the Yeast3 network of the best team is very close to 0 and its score can not be calculated due to some precision difficulties [16]. This has been confirmed by the results reported in Table 1, in which several computed p-values are zero that is impossible in practice. As the AUPR specifications of RELV ½2 , RELV ½3:5 , RELV ½pÃ and RELV ½g with that network are all higher than that of the best team, the corresponding scores of these estimators can not be computed, either. For the DREAM4 subchallenges, based on the evaluation scripts provided by the DREAM project organizers, the score of the suggested method is computed for every adopted normalization which is also included in Table 2. These results make it clear that both RELV ½2 and RELV ½3:5 could have ranked the first place in the Size 100 subchallenges of DREAM3 and DREAM4. But it should be emphasized that these comparisons are only of some reference values, noting that all the participants of the DREAM project were completely blind to both the structure and the dynamics of the networks.
Concerning the subchallenges of DREAM4, note that the best team integrated their down ranking method with the Z-score based method, and the score improvement does not exceed 1:3 point. On the other hand, the scores of the methods RELV ½2 and RELV ½3:5 are respectively greater than this best team approximately 1:8 points and 4:0 points. These performance improvements appear not to be a small one. When the average ratio is considered for the subchallenges of DREAM3 about the improvements on the AUROC and the AUPR specifications, similar conclusions can also be achieved.
When p-values are directly used in comparing performances of these estimation algorithms, consistent conclusions can be achieved. For example, when the p-value of the AUPR specification is taken into account for the Net1 network of DREAM4, the values of the Z-score based method and the best team are respectively about 2:3709|10 17 times and 5:8786|10 6 times of that of the method RELV ½3:5 .
It appears also worthwhile to note that the best team of DREAM4 utilized an estimation method different from that adopted by the best team of DREAM3. The results of Tables 1 and 2 may imply that the method suggested in this paper shares advantages of different approaches, and overcomes to some extent their disadvantages. But a theoretically solid justification for this declaration is still under investigation, and further efforts are required to clarify the actual reasons behind these phenomena.
Note that in the estimation algorithm suggested in this paper, the step of estimating genes with a zero in-degree plays an important role. To see the effectiveness of the proposed method in this estimation, the number of genes estimated to be of a zero indegree is given in Table 3 for each network of DREAM3 and DREAM4, together with its actual value. In this estimation, l 2 , P min and P max are selected as the same as those adopted in obtaining the estimation results reported in Tables 1 and 2. In this table, an estimation error has also been given which stands for the number of genes that can be regulated by other genes but are estimated to be with a zero in-degree, which is called in this paper, with a slight abuse of terminology, also as a FN (false negative) error. Table 3 shows that the suggested method is really effective in estimating genes that can not be regulated by other genes. More detailed analyzes on the estimation results show that if an FN error occurs, then, the genes that are wrongly estimated to be of a zero in-degree are usually regulated by less than 2 other genes. Moreover, if a gene with a zero in-degree is wrongly estimated to be regulated by other genes, then, in the corresponding probabilities, say, P ji s, the number of values that are significantly greater than 0 is usually less than 2. These types of mistakes appear reasonable in GRN topology estimation, noting that a large scale GRN usually has a sparse structure and measurement errors may happen to make the estimated value for every RELV of a gene with a small in-degree indistinguishable from 0. On the contrary, measurement errors are also able to make a few estimated RELVs of a gene with a zero in-degree significantly different from 0.

Robustness of the Suggested Method
Recall that in the suggested GRN topology estimation algorithm, parameters l 1 , l 2 , P min and P max should be selected. While these parameters have some biological interpretations, their selection has not been completely settled from a theoretical viewpoint. It is therefore interesting to investigate how sensitive the estimation accuracy is to the variation of these parameters. As knockout experimental data is used, it appears reasonable to select l 2 as l 2~0 :5. On the other hand, P min~0 :0005=0:05 and P max~0 :2 also seem to be an appropriate choice, as a big relative change with a small P min does not result in a significantly large P i imz1 , and a great P max may lead to a large amount of mistakes of wrongly estimating a gene regulated by other genes as a gene with a zero in-degree. These arguments imply that in GRN topology estimation, selection of the parameter l 1 is more essential. To investigate influences of the parameter l 1 on the prediction accuracy of GRN topology inference, 100 samples are taken for this parameter which is logarithmically equally spaced over the interval ½10 {2 , 10. For every sampled parameter l 1 , values of AUROC and AUPR for each network of DREAM3 and DREAM4 are calculated with the suggested estimation algorithm using respectively the 2-norm and the 3:5-norm based normalizations. The difference between the obtained AUROC specification and that with l 1~0 :01, as well as the difference between the obtained AUPR specification and that with l 1~0 :01, are shown in Figures 1 and 2. In these calculations, l 2 , P min and P max are respectively fixed to be the same values as those used before.
From Figures 1 and 2, it can be seen that performances of the proposed algorithm do vary with the parameter l 1 . But these performances keep almost the same values if l 1 [½0:01,1. Moreover, except a few networks, these performances begin to decrease from l 1~1 . Consistent observations have also been found for the suggested inference algorithm with other p-norm based and the geometric average based normalizations. These results imply that in practical applications, it may not be very difficult to select an appropriate l 1 . In this paper, this parameter is usually chosen as 0:01.
To understand influences of different normalizations on GRN topology estimation accuracy, variations of the AUROC and AUPR specifications with the parameter p have also been investigated. The results are given in Figure 3. Note that in this figure, the parameter p for the Net3 network of DREAM4 should be modified. Its variation interval for this network is ½1, 101. Once again, to make the variations clearer, some particular values have been extracted from the calculated AUROC and AUPR specifications, which are given in detail in the caption of the figure. In these calculations, the parameters l 1 , l 2 , P min and P max are chosen as the same as those adopted before.
From Figure 3, it is clear that the adopted estimation accuracy metrics indeed vary with the parameter p. The optimal p that maximizes the AUROC specification is different from that maximizes the AUPR specification, and different network has a different optimal p. On the other hand, it is also clear from this  figure that although the optimal p is different for each network and each specification, significant specification change does not arise when the parameter p varies over a relatively large interval. More specifically, for each network, the variation of the AUROC specification is not larger than 0.01 in magnitude, and when p §2, the variation of the AUPR specification is not larger than 0.03 in magnitude. For some particular networks, such as Ecoli2, Yeast3 and Net4, the variation magnitude is much smaller. These observations suggest that in actual applications, it is not very difficult to find a suboptimal value for the parameter p. Particularly, p~3:5 appears to be an appropriate selection for every network of DREAM3 and DREAM4. This can also be confirmed from the results of Tables 1 and 2, which show that, compared with the results with the optimal p, significant performance degradation generally does not arise with the method RELV ½3:5 . It is worthwhile to note that p~3:5 is different from those that are widely adopted in system analysis and synthesis, in which p~1,2, or ? is used more extensively [25,27].
On the other hand, to investigate the validity of the suggested technique for estimating genes with a zero in-degree, the obtained i i m is perturbed to be i i mzj with j~+1,+2,+3,+4. This may simulate the situation under which i i m is different from its actual value due to estimation errors inŝ s i andx x ½wt,0 i , as well as the imperfectness of the adopted assumptions and numerical integration errors, etc. Through the aforementioned perturbations, the estimated number of genes with a zero in-degree can be changed respectively by +1,+2,+3,+4 with respect to that of the unperturbed one. The obtained results for the methods RELV ½2 and RELV ½3:5 are respectively shown in Figures 4 and 5. When other normalizations are utilized, consistent observations have  From Figures 4 and 5, it can be seen that estimation performances with some networks can become slightly better when i i m deviates from the value adopted in the suggested estimation algorithm. For example, both the AUROC and the AUPR specifications of the Ecoli1 network and the Net2 network are better when the gene numbered i i mz1 is also regarded to be of a zero in-degree, and the AUPR specification of the Yeast2 network and the Net5 network is a little higher when the gene numbered i i mz1 is also considered as a gene not regulated by other genes. However, these performance improvements are not very significant, and when all the networks are taken into account, it is still better to use i i m in GRN topology estimations. In addition, if there are small variations in i i m , significant performance decrement usually does not arise.

Performances for Integration of Knockdown and Knockout Data
In this subsection, the suggested method for integrating knockout and knockdown experimental data is applied to the Size 100 subchallenges of both DREAM3 and DREAM 4. As mentioned before, rather than to develop a high performance integration method, the major purposes to include these results are to clarify effectiveness differences of knockout and knockdown experimental data in GRN topology estimations when the suggested method is adopted. In order to compare estimation performances, results using the Z-score based method are also integrated with completely the same procedure, that are respectively obtained from the knockout and knockdown experimental data.
The computational results of the Size 100 subchallenges of DREAM3 and DREAM4 are given respectively in Tables 4 and 5, in which KD, KO and MIX stand respectively for the estimation results obtained from knockdown experimental data only, knockout experimental data only and both of them using the above integration algorithm. Due to space considerations, the reported results are restricted to those with respectively the 2-norm and 3:5-norm based normalization. When other normalizations are utilized, consistent observations have been obtained and the conclusions are similar. For comparisons, the results are also included that are obtained using the Z-score based method.
From Tables 4 and 5, it is clear that when applied to the DREAM4 subchallenges, the suggested integration procedure is able to improve estimation performances for both the Z-score based method and the estimation algorithm suggested in this paper. As a matter of fact, compared with the results using only knockout experimental data, although there is one network with which the AUPR specification has been slightly degraded when the method RELV ½3:5 is used, the final score of RELV ½3:5 has been increased by about 3.8 points. Furthermore, the scores of the method RELV ½2 and the Z-score based method have been increased more significantly, which are respectively about 5.6 and 7.0 points. These improvements seem not small, noting that the best team of DREAM4 integrated the Z-score based method with their down ranking method, but the obtained merits are less than 1.3 points. In addition, this integration method appears more effective for the Z-score based method. More specifically, under such a situation, for each network, every adopted specification has been improved, and the final score has been increased almost 10%. These observations are significantly different from those of [18][19][20], which indicated that when knockout experimental data are available, knockdown experimental data is of little values in GRN topology inference.
Similar conclusions can be achieved if the p-values of the obtained estimation specifications are directly compared.
However, when utilized in the DREAM3 subchallenges, the aforementioned integration procedure does not work very well either with the Z-score based estimation algorithm or the algorithm suggested in this paper. Compared with the results using only knockout experimental data, this integration even worsen almost every specification of each network. The reasons are still not clear which are worthy of further efforts. But from these observations, it is clear that compared with those of DREAM3, information in the data sets of DREAM4 about the structure of a GRN are more consistent which are respectively contained in the knockout and knockdown experimental data.
Note that although in DREAM3, only measurement errors are added into the simulated experimental data, variances of the measurement errors are assumed to be of the same value for every gene under all situations, no matter it is in the wild type, or when some genes of the GRN have been knocked out or knocked down. On the other hand, in DREAM4, external disturbances are added to both the simulated mRNA concentrations and the simulated protein concentrations, but both background noises and the fact that gene expression values are typically measured on a logarithmic scale have been taken into account in simulating these external disturbances. Such a treatment makes a simulated measurement error have a standard variance approximately proportional to a simulated actual value of the expression level of a gene [16,17]. Note that the magnitude of a knockout perturbation is twice as that of a knockdown perturbation. It can therefore be declared that compared with the knockdown experimental data of DREAM4, those of DREAM3 are more noisy, and hence less informative. These can also be seen from the differences between the AUROC/AUPR specifications using respectively only the knockout experimental data and the knockdown experimental data. As a matter of fact, it is clear from Tables 4 and 5 that for all the adopted estimation methods, compared with their counterparts in network topology estimations of DREAM4, the above differences are consistently larger in those of DREAM3, especially when the AUPR specification is considered. As the simulated data of DREAM4 are believed to be closer to actual biological experimental data than those of DREAM3 [14,17,20], it is hoped that the suggested integration method is helpful in practical GRN topology estimations. In addition to these, it is also clear from these tables that when only knockdown experimental data is utilized, the Z-score based method outperforms about 0:8 point the method suggested in this paper with the DREAM4 subchallenges. But when the DREAM3 subchallenges are coped with, the conclusions are completely the opposite, in which the methods suggested in this paper, no matter the method RELV ½3:5 or the method RELV ½2 , can obtain a score higher than the Z-score based method approximately 3:0 points.

Performances for Integration with the Down Ranking Algorithm
In this subsection, estimations are performed using the integration procedure proposed for the suggested RELV based inference method and the so called down ranking method. The corresponding results are given in Tables 6 and 7 when the pnorm based normalization with p~2 and p~3:5 are respectively used in these method integrations. The corresponding results are given in the rows started by RELV ½2,dr and RELV ½3:5,dr respectively. To compare the effectiveness of method integration, results obtained through integrating the Z-score based method and the down ranking method are also included, which are denoted by Z{Score ½dr . In these tables, only results with some typical and optimal values for the threshold of the down ranking method are included. In searching the optimal threshold value, the AUPR specification is once again taken as the cost function. Similar to Tables 1 and 2, estimation results with the optimal threshold, as well as those of Z{Score ½dr with c~2:5000, are included here only for some references. The major purposes for this inclusion are to clarify estimation performance degradations when the adopted threshold c deviates from its optimal value.
From Tables 6 and 7, it is obvious that the down ranking method is indeed helpful in improving estimation accuracy, no matter it is integrated with the Z-score based method or the method suggested in this paper. Moreover, compared with the AUROC specification, the AUPR specification has been improved more significantly. In addition, estimation performances for the DREAM4 subchallenges have been improved more greatly than those of the DREAM3 subchallenges.
An interesting observation from these tables is that when the optimal threshold value is adopted for the down ranking method, although the method RELV ½2 still outperforms the Z-score based method, performance differences among the methods RELV ½2 , RELV ½3:5 and the Z-score based method become smaller than those before the method integration. While this may mean that the down ranking method is more effective in improving the Z-score based method, it may also suggest that the method proposed in this paper is really effective in reducing the so-called cascade errors in GRN topology estimations.
The results of Tables 6 and 7 also indicate that although the algorithm suggested in this paper is able to reduce the so called cascade errors, there still exist some cascade errors that this algorithm fails to detect. This may possibly be due to the following three causes. One is the imperfectness of the experimental data in which several kinds of noises exist. One is the imperfectness of the adopted assumptions on measurement errors, which may have not appropriately reflected their actual distributions. The other one is that there may exist genes for which indirect regulations cause a RELV with a magnitude bigger than that caused by direct regulations.
On the other hand, it seems that the down ranking method is much more helpful in improving the prediction performance of the method RELV ½2 than that of the method RELV ½3:5 .
In applying the down ranking algorithm, a threshold value c should be chosen for extracting a primary estimation about the network structure from some computed weights or confidences about direct regulations between any two different genes of a GRN. There is, however, still no very solid theoretical guidance on how to suitably choose this threshold value [22]. As an example, it is reported in [22] that while c~2:0 is found through extensive numerical simulations to be the best selection for integrating with the Z-score based method, c~2:5 is more appropriate for the subchallenges of DREAM4. It is therefore interesting to investigate variations of estimation performances with this parameter. Due to space considerations and the fact that RELV ½2,dr outperforms RELV ½3:5,dr , discussions are restricted to the method with the 2norm based normalization. When the 3:5-norm based normalization is utilized, consistent observations have been obtained and the conclusions are similar. According to the results reported in [22], when the Z-score based method is to be integrated, the interval for the parameter c is selected to be ½0,5 in this paper. On the other hand, when the method suggested in this paper is to be integrated, this interval is chosen as ½0, 1. In these intervals, 100 equally spaced samples are used in searching the optimal c. Figure 6 shows these variations when the Z-score based method and the algorithm suggested in this paper with the 2-norm based normalization are respectively integrated with the down ranking method.
Variations of the AUROC and the AUPR specifications with the parameter c are shown in Figure 6. From this figure, it is clear that although the optimal value is different for each network and each specification, c~0:31 appears to be a good choice for the threshold value when the down ranking method is integrated with the estimation method suggested in this paper. The corresponding results for the Size 100 subchallenges of DREAM3 and DREAM4 are given respectively in Tables 6 and 7, together with those using the optimal c.
The results of Figure 6 are also consistent with the observations reported in [22]. That is, when the Z-score based method is integrated with the down ranking method, c~2:5 is more appropriate for the Size 100 subchallenges of DREAM4, although c~2:0 is generally believed to be the best selection.
From Figure 6, it is also clear that compared with the Z-score based method, estimation performances of the algorithm suggested in this paper is less sensitive to variations of the threshold around its optimal value, when they are respectively integrated with the down ranking method. This property is attractive in practical applications, recalling that it is still not very clear how to choose the optimal threshold value for a particular GRN and an experienced value usually deviates from the optimal one.

Further Discussions
As commented in [16], highly confident predictions in GRN topology estimations can become a good guidance to biological experiment designs. However, these predictions will be helpful only if their precisions are also sufficiently high. This requirement asks that a desirable estimation algorithm should have a PR (precision-recall) curve starting from the left upper corner, and decreasing monotonically and slowly with the increment of the recall rate. To see whether predictions made by the algorithm suggested in this paper share this property, the PR curve of the method RELV ½3:5 is shown in Figure 7 for each network of the Size 100 subchallenges of DREAM3 and DREAM4, which is based only on the knockout experimental data. To compare satisfaction degree about this requirement with the Z-score based method, its corresponding PR curve for each network is also included.
From this figure, it is obvious that for every network of DREAM3, when the recall rate is around 0, the prediction precision of the suggested estimation method is approximately equal to 1, and this prediction precision keeps large when the recall rate is less than some value. Moreover, this value is specially large for the Ecoli2 network. This suggests that for the DREAM3 Size 100 network inference subchallenges, predictions with a high confidence obtained by the suggested method are usually correct and are therefore helpful in the design of a follow-up experiment validation. This is different from the algorithm used by the best team and the second place team of DREAM3, which may not be very desirable in this aspect [16].
However, when applied to the DREAM4 subchallenges, the aforementioned properties do not hold for most of the networks. As a matter of fact, except the Net1 and Net5 networks, the PR { As noted in [22], c~2:5000 is obtained for Z{Score ½dr after a comparison with the actual network. On the other hand, the optimal c can hardly be obtained in actual estimations for each of Z{Score ½dr , RELV ½2,dr , RELV ½3:5,dr . The purposes to include their inference results here are only to clarify estimation performance degradations when an empirical parameter c is adopted. *Due to the some reasons as those of  The purposes to include the inference results of Z{Score ½dr with c~2:5000, Z{Score ½dr with the optimal c, RELV ½2,dr with the optimal c, RELV ½3:5,dr with the optimal c, are completely the same as those of Table 6. That is, to clarify estimation performance degradations when an empirical parameter c is adopted for these methods. doi:10.1371/journal.pone.0031194.t007 curve even does not start from the upper left corner. This means that there still exist some false positive errors among the estimated direct regulations whose existence is predicted with a high confidence by the suggested method. Furthermore, when the suggested method is integrated with the down ranking method, similar observations have been obtained. On the other hand, when the Z-score based method is utilized, consistent phenomena have been observed.
Nevertheless, a detailed analysis shows that concerning this requirement on GRN topology estimators, the Z-score based method does not outperform the method suggested in this paper, either. As an obvious example, in the DREAM4 subchallenges, when the Z-score based method is utilized, only the Net1 network has a PR curve starting from the upper left corner. More detailed comparisons are omitted, but it can be claimed from Figure 7 that the method suggested in this paper appears more helpful than the Z-score based method in guiding the design of a biological experiment to validate the actual existence of a predicted direct regulation.
When the Z-score based method is integrated with the down ranking method, which is adopted by the best team of DREAM4, the corresponding PR curves for these benchmark networks are very similar to those obtained from the Z-score based method. This implies that further enhancements are still required to make this integration applicable to practical GRN structure inferences.
Computations have been performed also on many other simulated large scale GRNs. The observed phenomena are consistent with what have been reported in this section.

Concluding Remarks
In this paper, an algorithm is developed for inferring GRN topology from steady state knockout/knockdown experimental data. Rather than the commonly used AELVs (absolute expression level variation), it utilizes RELVs (relative expression level variation) of a gene in gene knockout/knockdown experiments to measure possibilities of the existence of direct regulations among genes. Based on this variation, probability is estimated from experimental data for the existence of a regulation between two different genes of a GRN, which is further used to estimate whether or not a gene is regulated by any other genes. The estimated magnitude of the RELV of a gene is normalized and modified, on the basis of the estimation results about the existence Figure 6. Variations of the AUROC and AUPR specifications with the threshold value c. To make the variations clearer, the specifications shown are their deviations from those respectively with c~0:05 (for the Z-score based method) and with c~0:01 (for the algorithm suggested in this paper). doi:10.1371/journal.pone.0031194.g006 of direct regulations to it. These normalized and modified magnitudes are used in queuing the possibility of the existence of a corresponding direct regulation. A distinguished characteristic of this algorithm is that its computational complexity increases only quadratically with the number of genes in a GRN.
Computational results with the Size 100 subchallenges of both DREAM3 and DREAM4 show that this method can outperform not only the widely used Z-score based method, but also the best team of these subchallenges who used an integration of some well known methods. While these comparisons are only of some reference values, as all the DREAM project participants were completely blinded to both the structure and the dynamics of the simulated networks, it appears safe to claim that the suggested method is more efficient than the available methods in distinguishing direct and indirect regulations of a GRN. Integration with the so-called down ranking method show that the so-called cascade errors in GRN topology estimations can be further reduced. Precision analyzes show that highly confident predictions obtained by this method are usually more helpful in guiding designs of a biological validation experiment than those by the Z-score based method.
Further efforts along this line appear to test effectiveness of the suggested method with actual biological experimental data, to extend the suggested estimation method to biological experimental data in which several genes are simultaneously perturbed by external interferences, to give a more biologically significant normalization of the RELVs and selection of the parameters l 1 , l 2 , P min and P max , as well as to improve estimation accuracy of gene expression levels in the wild type and that of the variance of measurement errors. Challenges still remains there in reducing false positive errors among highly confident predictions, especially when the RELV of an indirect regulation is larger in magnitude than that of some direct regulations. It is also interesting to see whether some other structure information about a GRN, such as the power law, etc., can be helpful in making a more accurate prediction.
Dr.Soranzo and Dr. de la Fuente, for kindly providing the MATLAB codes of their user-friendly down ranking algorithm.
Availability The Matlab files for the method are available at http/sol;bioinfo.au. tsinghua.edu.cn/member/ylwang/Matlabfiles_REVL.zip.