Inference of Biological Pathway from Gene Expression Profiles by Time Delay Boolean Networks

One great challenge of genomic research is to efficiently and accurately identify complex gene regulatory networks. The development of high-throughput technologies provides numerous experimental data such as DNA sequences, protein sequence, and RNA expression profiles makes it possible to study interactions and regulations among genes or other substance in an organism. However, it is crucial to make inference of genetic regulatory networks from gene expression profiles and protein interaction data for systems biology. This study will develop a new approach to reconstruct time delay Boolean networks as a tool for exploring biological pathways. In the inference strategy, we will compare all pairs of input genes in those basic relationships by their corresponding -scores for every output gene. Then, we will combine those consistent relationships to reveal the most probable relationship and reconstruct the genetic network. Specifically, we will prove that state transition pairs are sufficient and necessary to reconstruct the time delay Boolean network of nodes with high accuracy if the number of input genes to each gene is bounded. We also have implemented this method on simulated and empirical yeast gene expression data sets. The test results show that this proposed method is extensible for realistic networks.


Introduction
In order to understand complex biological networks and pathways, we need to investigate global structures instead of individual behaviors since there are interactions and associations between genes. Due to the invention of high throughput technology, genome-wide expression profiles can be measured simultaneously [1]. However, it is still a great challenge to identify complex biological networks from genome-wide data because the number of gene interactions is huge [2]. In recent years, there has been a significant progress in research concerning genetic network models and network reconstruction problems.
Clustering and dimension reduction are important methods for grouping genes that have similar expression profiles [3,4]. In the framework of clustering, it is important to define the degree of similarity between genes. By the method of clustering, we can group genes that have similar expressions. However, we still cannot find the causal relationship between genes. Hence, apart from the relationship of similarity, we will also have to consider another causal relationship between genes.
There have been many methods proposed in the literature to tackle the problem of genetic regulatory network reconstruction. For instance, the steady state approach have been used to model gene regulatory networks [5]. In addition, the Bayesian network model is an important technique that has been studied extensively in the past two decades [6][7][8][9][10][11]. A Bayesian network is a directed acyclic graph (DAG) comprised of two components. The first component is comprised of nodes that correspond to a set of variables and a set of directed edges between variables with Markov properties. The second component describes a conditional distribution for each variable given its parents in the graph. model, nodes represent the gene expression states. The status of a gene is quantized to one of the two states: on or off, representing a gene as active or inactive respectively. The wiring of rules between nodes in the graph represents a functional link between genes and determines the expressions of target genes after giving a series of input genes. Under the structure of Boolean networks, the target gene is determined by a set of genes with specific rules. For each gene, if the indegree (i.e., the number of input genes to each gene) is bounded by a constant K, only O( log n) pairs of state transition are necessary and sufficient to reconstruct the original network with n nodes [27,28]. However, Boolean networks have been criticized for their deterministic nature. The assumption that every affected gene would be expressed immediately at the next time step may be unsound.
Another point of view of constructing genetic network is to focus on the indication the pairwise relationships between genes. Most of the works is to find the gene-pairs with similarity relationship [29][30][31][32][33]. The similarity of a gene-pair represents the two genes with the same expression or opposite expression. In 2005, Li and Lu proposed directed acyclic Boolean network and the statistical reconstruction method of SPAN to infer the pair wise relations of every element [34]. The proposed method can reconstruct Boolean networks from noisy array data by assigning an s-p-score for every pair of genes. In the study, they proposed another relationship between two genes: relationship of prerequisite under the Boolean network model. If gene A is a prerequisite for gene B, then the ''on'' status of gene A is necessary for the ''on'' status of gene B. Boolean implication network, with the similar aspect, investigated all Boolean implication between pairs of gene for large scale genome microarray datasets [35]. Following the model, Wang et al. [36] proposed a two step counting approach for constructing biological pathways with Boolean network. However, most of these methods only consider pair wise relationship in order to decrease the time complexity. Therefore, if the structure of network is a combination of a set of genes to affect another gene, the algorithms will lose some information and rules in the genetic network reconstruction.
In this study, we will consider a much more generalized model by combining the structure of the above two models. If a Boolean function with one or several genes is a prerequisite for a target gene, then the induction of the Boolean function with input genes is necessary for the expression of the target gene. Hence, the target will be influenced by the Boolean function with several input genes. However, the induction of the Boolean function may not activate the target gene immediately, but at a future time. Therefore, the target gene may not have been influenced right now and we will treat these relationships as time delay affection. In this study, we will infuse these additional relationships for more generalized systems.

Boolean Network
Boolean networks were introduced by Kauffman (1969) forty years ago to represent genetic regulatory networks. First, we will review the definition of a Boolean network. A Boolean network G(V ,F ) is a directed graph consisting of two components: a set of nodes V~fv 1 ,v 2 , . . . ,v n g that corresponds to genes, and a list of Boolean functions F~ff 1 ,f 2 , . . . ,f n g that corresponds to the rule of interaction and combination of several genes. For every node v i [V , its expression is simplified to two levels: on and off, representing a gene as active or inactive. For every Boolean . . ,v i k are assigned to the node v i in the graph and represent the rules of regulatory mechanisms between genes. The expression of a gene is determined by the expression of the gene directly affecting it with a Boolean function. Therefore, the state of each node For each node v i , the gene expression state at time t is assumed to take either 0 (not-expressed) or 1 (expressed) and is expressed as y t (v i ). In a Boolean network, every gene expression profile at time tz1 is completely determined by the expression profile of a set of genes v i1 ,v i2 , . . . ,v i k at time t and the corresponding Boolean  For convenience, we converted the Boolean network G(V ,F ) to the wiring diagram G'(V ',F ') (See Figure 1) [37]. For each node v i [V , suppose v i1 ,v i2 , . . . ,v i k are the input nodes assigned to v i . Then we construct an additional node v ' i and connected the edge from v ij to v i' for each 1ƒjƒk. That is, the set of fv 1 , . . . ,v n g represents the gene expression profile at time t and the set of fv ' 1 , . . . ,v ' n g corresponds to the gene expression profile at time tz1. Hence we can treat the set of fv 1 , . . . ,v n g as the input values and the set of fv ' 1 , . . . ,v ' n g as the corresponding output values. Therefore, the output values of fv '

The Structure of Time Delay Boolean Network
In the previous subsection, we found that given the values of the node (V) at time t, the expressions at time tz1 will be updated immediately by specific Boolean function (F ). That is, for every gene v i [V, if the expression profile of a set of genes fv i1 ,v i2 , . . . ,v i k g at time t and the corresponding Boolean function f i is observed, the gene expression of v i at time tz1 is determined by ). However, in real genetic regulatory situations, the deterministic system has been criticized due to the existence of misclassification error and noise. In addition, some of the gene expression may result in time delay when the gene is influenced by one or several input genes. That is, the induction of Boolean function may not activate the target gene immediately, but in the future. Hence, it would have been much more flexible to use a non-deterministic network system. In this subsection, we will consider two relationships between the Boolean function and the target gene instead of the deterministic relation.
First, we will introduce the structure of time delay Boolean networks. Suppose there are n elements, v 1 ,v 2 , . . . ,v n in a Boolean network. For any elements v i with specific Boolean function f i , we have two kinds of pair wise relationship: prerequisite and similarity. We say that a Boolean function f i with specific k input genes v i1 ,v i2 , . . . ,v ik at time t is the prerequisite for the target gene v i at time tz1, if the on-status of Boolean function at time t is necessary for the on-status of gene v i at time tz1. This relationship is denoted In other words, if the Boolean function f i is not active at time t, gene v i will be inactive at time tz1. If it does not cause confusion, we will omit the notation of y and input genes as denoted by f i [v i . Moreover, for every gene v i , we use v v i as its dual (from 0 to 1, or from 1 to 0) in this paper. Therefore, for any Boolean function and target gene with a prerequisite relationship there are a total of two possible relationships: In this model, we do not consider the situation of a dual of Boolean function prerequisite to the target gene, that is Since for any Boolean function whose dual is a prerequisite to the target gene, there must exist another Boolean function that is a prerequisite to the target gene. For . Therefore, for the prerequisite relationship, we only consider the Boolean function that is a prerequisite to target gene and the dual of target gene.
The other type of relationship between Boolean function and target gene is similarity. We say that the Boolean function f i and  Table 2. Count profiles for the basic eight relationships without misclassification error.
target gene v i are similar if the status of the Boolean function and the status of the target gene are in the same expression, and we denoted this by f i ,v i . In the same way, we do not consider the situation of Boolean function similar to the dual of target gene such as f i , v v i in this study. Since if there is one Boolean function that is similar to the dual of target gene, there must exist another Boolean function that is similar to the target gene.
In the diagram, if a Boolean function f i is a prerequisite to v i , we draw a directed arrow from the vertex f i to v i and if f i is similar to v i , we use an undirected line to connect f i and v i .
In the model of time delay Boolean network we proposed, the output of the gene expression is not completely determined by the input state and Boolean function. The output expression may have more than one possible result in the time delay Boolean network. We illustrate the above construction by an example in Figure 2. It has three elements, one similarity and two prerequisite relationships. The possible outputs for every input state are listed in the right part of the graph. If we knew the network structure, some of the inputs would have more than one possible output expression in the time delay Boolean network.

Identification Algorithm
First, we only consider Boolean networks in which the maximum number of input genes is bounded by a constant K for every target gene, because it has been proven that the number of profiles required grows exponentially if K is not bounded [38]. For simplicity, we only show algorithms for the case of K~2. However, the algorithm can be intuitively generalized to any K in a straightforward way. For the inference of genetic network, we need to clarify the following questions for each target gene. N Which input genes will affect the target gene? N What kind of Boolean functions will be used for combining those input genes?
N What kind of relationship exists between the Boolean function and the target gene?
In this subsection, we propose an algorithm to clarify the above questions. The algorithm below is conceptually very simple since it simply uses output Boolean functions with input genes and relationships with target genes that are consistent with the data. First, for each output gene expression at time tz1 such as v ' i , we consider all the pairs of elements in V at time t, for instance v j and v h . Then we count the eight incidents of (v j ,v h ,v ' i ) being (0,0,0), (0,0,1), . . ., (1,1,1) from the sample and arrange them in a 2|4 table; see the left part of Table 1. We mark a cell ''+'' if the count is positive and mark it ''0'' otherwise.
For detecting whether there exists a Boolean function which is a prerequisite to the target gene, we will compare the 2|4 output table with the left four basic relationships in Table 2. We consider the basic relationships to be consistent with the output table if the position of 0 cell in the basic relationships is also 0 in the output table. By comparing the output table with the four basic relationships, we can find relationships that are consistent with the output tables. If there is more than one relationship that is consistent with the output tables, we would use the Boolean logic gate ''and'' to combine the Boolean function and transfer the result to another Boolean function. Hence, the final Boolean function is the prerequisite to the target gene. Similarly, by comparing the 2|4 output table with the right four basic relations in Table 2, we could get another Boolean function which is the prerequisite to the dual of target gene.
Moreover, if only one Boolean function occurred in above relationship, that is, if there is only one Boolean function that is the prerequisite to the target gene or the dual of target gene, we will treat that relationship as our final relationship between the Boolean function and the target gene. However, if both of the two prerequisite relationships happened (i.e. Af i and f ' i s:t: i ), we need to check whether these two relationships are in conflict. If the dual of f i is equivalent to f ' i , our conclusion for inference will be that f i is similar to the target gene (that is, f i *v ' i ); otherwise, we will treat it as if there is no relationship between the input genes and the target gene because we did not gather enough information to judge true relationships between v ' i and (v j ,v h ) at this moment. By the above identification procedure, we can find the corresponding input genes, Boolean function and their relationship for every target gene.

Identification Algorithm with Noisy Array
In previous subsection, we discussed an identification method for data without noise. In this section we will consider the situation of noisy array data. We assume that every element in the entry of (I ij , O ij ), j~1,2, . . . ,m switches to its reverse status with a misclassification probability p independently; that is O Ã ij~O ij with probability 1{p ; 1{O ij with probability p : Thus, the observed array (I Ã ij , O Ã ij ) contains misclassification error. Our goal is to reconstruct time delay Boolean network from noisy array of binary data (I Ã ij ,O Ã ij ). Similar to section 2, we assume that the maximum number of input genes is bounded by 2 for every target gene. We treat the data in the 2|4 table as a multinomial distribution with eight cells whose probabilities are q 000 ,q 001 , . . . ,q 111 as shown in the right part of Table 1, where q 000 zq 001 z . . . zq 111~1 . Similarly, we extract the data with misclassification error for every output gene and each pair of input genes as the 2|4 table. Now the observed data n 000 ,n 001 , . . . ,n 111 are not generated from the multinomial q 000 ,q 001 , . . . ,q 111 , but from another multinomial r 000 ,r 001 , . . . ,r 111 as shown in Table 3, where r 000 zr 001 z . . . zr 111~1 .
Because of the misclassification error, a portion of the samples of m 000 may change to the other seven cells. We use the notations of m 000,000 ,m 000,001 , . . . ,m 000,111 to represent the counts of eight cells changed from m 000 . Analogous notations are defined for m 001 ,m 010 , . . . ,m 111 . The splitting is shown in Table 4. Consequently, the generated probabilities (q 000 ,q 001 , . . . ,q 111 ) are calculated as follows: q i1i2i3,j1j2j3~p I(i,j) (1{p) 3{I(i,j) q i1i2i3 , where Here, we adopt the notation q i1i2i3,j1j2j3 analogous to m i1i2i3,j1j2j3 . The above parameters and splits are shown in Table 4. In the table, it is easy to find that the correspondence between two sets of counts and probabilities is the following: and ð3Þ For the complete data fm i1i2i3,j1j2j3 g, the log-likelihood is given by where q i 1 i 2 i 3 ,j 1 j 2 j 3 are those splitting probabilities. Since the complete data fm i 1 i 2 i 3 ,j 1 j 2 j 3 g are not observable, we use the EM algorithm to maximize the log-likelihood. In the E-step, the splitting counts of complete data fm i 1 i 2 i 3 ,j 1 j 2 j 3 g are evaluated by the conditional expectations using the current values of the parameters by the following formula E p,q 000 ,q 001 ,...,q 111 (m i 1 i 2 i 3 ,j 1 j 2 j 3 Dn j 1 j 2 j 3 )~n where i 1 ,i 2 ,i 3 ,j 1 ,j 2 ,j 3~0 ,1. One probabilities of q 000 ,q 001 , . . . ,q 111 are zero in those different hypotheses specified in Table 5. In the M-step, we maximize the conditional expectation of the loglikelihood for the complete data to obtain the maximum likelihood estimates (MLEs) of the parameters. According to the MLEs, we can compute the p-score for every pair of input genes and target gene, which are obtained by estimating for the misclassification probability under every prerequisite relationship. For the first step, we would like to determine the most probable relationships between every pair of input genes and one output gene. Next, we find the most probable Boolean function with pair input genes for every output gene and select candidate pairs of input genes and output gene for the watch list. Finally, we reconstruct a time delay Boolean network by integrating the relationship of those genes selected.
For one output gene v 0 i and a pair of input genes v j and v k , we define the p-scores p (vj i are, respectively, the maximum likelihood estimates of p under the triangular model: q 000~0 , q 010~0 , q 100~0 , q 110~0 , q 001~0 , q 011~0 , q 101~0 , q 111~0 .
According to the EM algorithm described above, we can evaluate the p-score for every output gene. We use the MLEp p to measure how well each hypothesis fits: the smaller the score is, the more likely that the corresponding hypothesis could be true.
If the samples are generated from a time delay Boolean network, p-score are quite useful for the discovery of true Table 5. The eight basic relationships and their probabilistic hypotheses and p-scores.

Relation
Hypothesis Scores relationships. Here we can consider the maximum compatibility criterion: to choose the maximum threshold value so that the selected relationships contain no conflicts [34]. We collect those relationships whose p-scores are smaller than a threshold. Known biological results are helpful for the determination of a threshold. For example, if we know the relationship (v 1 or v 2 )[v 3 is true, then the p-scores smaller than p (v1 or v2)[v3 should be in our watch list. As more relationships are included in the watch list, the more likely we are to observe incompatible ones. In general, we can choose the threshold that allows the maximum number of relationships with no conflicting relationships. Next we will demonstrate the method by illustration examples.

Theoretical Results
First, we will analyze the number of input/output pairs required for the network reconstruction of time delay Boolean network to be unique. The theoretical results of classical Boolean networks only consider the similar relationship [27,38,39]. The following results prove the theoretical results time delay Boolean networks that has a more flexible structure and consider both similar and prerequisite relationship.
Proposition 1. For all subsets of V with 2K genes, if all assignments (i.e., 2 2K assignments) of Boolean values appear in input expression patterns and all of its possible output expression patterns of the target gene are present, the identification of genetic network is determined to be unique, if it exists.
(Proof) Let z be any gene in V and suppose z is controlled by a Boolean function g(x i1 , . . . ,x i k ) with similarity or prerequisite relationship (i.e., g*z or g[z). If the Boolean function g is similar to z, the case is proved for the classical Boolean networks in Akutsu et al. (1998). Next, we consider the case of Boolean function g as a prerequisite to z. In this case, there must exist a specific input value fa 1 , . . . ,a k g for fx i1 , . . . ,x i k g such that z have two possible values 0 and 1. Hence, any other genes would not control z because all assignments of Boolean values are appearance. Let us illustrate the above statement by the example for the case of K~1 and K~2. If K~1 and x[z, when the input of x is 1, the outcome of z being both 0 and 1 will appearance. Therefore, given the input of x~1, the outcome of z is not deterministic no matter the value of any other gene y is 1 or 0. Hence, any other gene y would not affect gene z. If K~2 and g(x 1 ,x 2 )[z for some Boolean function g, there must exist an input (a 1 ,a 2 ) such that g(a 1 ,a 2 )~1. Then, for any other pair of gene fy 1 ,y 2 g where fy 1 ,y 2 g\fx 1 ,x 2 g~w, the outcome of z is not deterministic for any input of fy 1 ,y 2 g, if the input of fx 1 ,x 2 g is fa 1 ,a 2 g. In a case of fy 1 ,y 2 g\fx 1 ,x 2 g=w, we can prove that gene y i which does not belong to fx 1 ,x 2 g would not affect the gene z in a similar way.
Proposition 2. The probability that one sub-assignment with all of its possible results in the target gene does not appear among m random input expression pattern is at most 2(1{ 1 2 2Kz1 ) m . (Proof) For any fixed set of nodes fv i1 ,v i2 , . . . ,v i2K g, the probability that a sub-assignment v i1~vi2~. . .~v i2K~1 does not appear in one random input expression pattern is 1{ 1 2 2K . Thus, among the m random input expressions, the probability that v i1~vi2~. . .~v i2K~1 appears is t times is equal to In addition, the probability that all of the possible results in the target gene does not appear among t times input is smaller than ( 1 2 ) t{1 for 1ƒtƒm and equal to 1 for t~0. Hence the probability that one sub-assignment and all of its possible results does not appear among m random input expression is smaller than ) t{1 , and this can be bounded by 2(1{ 1 2 2Kz1 ) m by an algebra calculation.  Figure 1, we generate 100 samples with p = 0.05.

Samples Hypotheses Relation
Input Output q 000 = 0 q 010 = 0 q 100 = 0 q 110 = 0 q 001 = 0 q 011 = 0 q 101 = 0 q 111 = 0  Next we prove the main theorem. Theorem 1. For the identification of one time delay Boolean network of n nodes with maximum indegree ƒK, O(2 2Kz1: (2Kza) : log n) uniformly and randomly sampled input patterns are sufficient for exact inference with probability at least 1{ 1 n a . (Proof) We consider the probability that the condition of Proposition 1 is not satisfied under m random input expression patterns.
By Proposition 2, the probability that v i1~vi2~. . .~v i2K~1 with all of its possible results in the target gene does not appear among the m random input expression patterns is bounded by 2(1{ 1 2 2Kz1 ) m for any fixed set of nodes fv i1 ,v i2 , . . . ,v i2K g. Since the number of combinations of 2K nodes from a set of n possibilities is bounded by 2 2K : n 2K , the probability that the condition of Proposition 1 is not satisfied is at most 2 2K : n 2K : 2(1{ 1 2 2Kz1 ) m . It is not difficult to see that 2 2K : n 2K : 2(1{ 1 2 2Kz1 ) m vp holds for mw ln 2 : 2 2Kz1: (2Kz1z 2K log nz log 1 p ). Hence, we obtain the theorem by letting the non-identification probability p~1 n a . Next we develop an information theoretic lower bound on the number of input/output pairs needed for the identification of a time delay Boolean network.
Theorem 2. If the maximum indegree ƒK, at least V(2 K zK log n) input/output pairs are required for the identification of a time delay Boolean network in the worst case.
(Proof) The number of time delay Boolean networks is given by all the possible combination of Boolean function with k nodes from a set of n possibilities with all possible relations between Boolean functions with target node. Since there are V(n K ) possible combinations of input nodes, 2 2 K possible Boolean functions and 3 possible relations between Boolean function with each node, there are V((2 2 K : n K : 3) n ) Boolean networks whose maximum indegree is at most K. On the other hand, there are at most 2 n possible output patterns with one input expression pattern. Therefore, V(log 2 n ((2 2 K : n K : 3) n )) which is the same as V(2 K zK log n) input/output pairs are required in the worst case.

Example with Simulation and Real Data
We will illustrate our method by the example described in Figure 2. For the pair of samples consist of three elements list in the right part of Figure 2, we uniformly generated 100 input samples and their corresponding possible output samples with misclassification probability p~0:05. For the prerequisite relationship, if the status of Boolean function with input genes is on, then we allow the output value to have equal probability of on or off. The data can be arranged as input/output sample similar to that obtained from the microarray data with time. Namely, the input of each sample can represent the gene expression at time t and the output can represent the gene expression at time tz1. For each pair of input and output genes, we compute the 8 basic pscores that represent the 8 basic hypotheses in Table 5 for all of pair input genes and output genes. After the calculation, the simulation results of every p-score are listed in Table 6.
Beside the example with 3 elements, in order to shows the superiority of the proposed method can be applied to a larger network, a more comprehensive example with a larger network is given in Figure S1.
Next, we have to decide the threshold for choosing the relations. When we increase the threshold of the p-score, the relations whose p-score are smaller than the threshold will be chosen. Moreover, when the number is 0.138, the conflict occurs, since we have However, in our model, there are at most two genes that would affect an output gene. Therefore, we can choose 0.138 as our threshold and include relations whose p-score is smaller than the threshold. By these procedures, we can reconstruct the time delay Boolean network identical to Figure 2.
In the area of gene regulatory network study, Schuller has summarized regulatory cis-acting elements of structural genes of the nonfermentative metabolism and described the molecular interactions among general regulators and pathway-specific factors [40]. In the gene regulation of gluconeogenesis by Sip4 and Cat8 pathway, the carbon source control could be identified for the regulator Cat8; see ( Figure 6) in Schuller [40]. In this study, we apply our proposed approach to explore the expression profiles and show some exploratory result on the Cat8 pathway.
In order to demonstrate the effectiveness of reconstruction, we use the microarray expression dataset of yeast Saccharomyces cerevisiae produced by DeRisi et al. [1] and Spellman et al. [41]. In total, the data is comprised of 41 experiments after filtering out experiments with missing values. By these experimental microarray data sets, we can use our proposed method to reconstruct the biological pathway and the genetic regulation network result is shown in Figure 3. The result is consistent with the genetic network in the literature. That is, the restraint of Mig1 or activation of Snf1 is a prerequisite for the decreasing of Cat8. Moreover, the restraint of Snf1 or Cat8 is a prerequisite for the decreasing of Mls1. However, the negative similarity between Snf1 and Mig1 is undetectable in our current model.

Conclusions
In this paper, we have introduced the model of time delay Boolean network that generalizes the Boolean network model in order to cope with dependencies that have two kinds of relationships: similarity and prerequisite. The approach for reconstruction of genetic network inference from gene expression data relies on the assumption that the expression of a gene is likely to be controlled by a relatively small number (say k) of genes. Also, some bounds on the size of data required for the identification of the time delay Boolean networks under constant of indegree are stated and discussed. Moreover, the algorithm of the network reconstruction from noisy array data is developed.
One characteristic of a Boolean network is that all the variables in the graph are binary. If the data we observed is continuous or quantized to have more than two levels, we need to discretize them. For microarray data, the ratios of expression level would be one possible approach of discretization. That is, we can treat the gene as on (active) if the log-ratio of its expression is larger than zero. We treat it as off (inactive) otherwise. In general, biological background knowledge will be helpful for setting thresholds for discretizaion. On the other hand, if the samples are obtained from a time course, then we can consider the gene as on or off by detecting whether the gene is either increasing or decreasing with time.
The work in progress is aimed at evaluating the effectiveness of the described approach for inferring genetic networks from biological gene expression time series data. Besides that, implementation on some other real biological data is also an important task.
For the implement of the network reconstruction algorithm, the greatest complexity is the computation of p-score for each of the n! k!(n{k)! input elements and n output elements, where n is the number of elements and k is the number of indegree. It is an iterative algorithm to compute the MLE for the p-scores by EM procedure while the common practice is to set an upper bound for iterations in numerical implementation. Consequently, this keeps the O(n kz1 ) complexity for the computation of MLE. In addition, the sorting algorithm for the n! k!(n{k)! n data cost O(n kz1 log n) in terms of time. Hence, the overall time complexity for the network reconstruction is O(n kz1 log n) for this algorithm. Figure S1 An example of genetic network with 8 nodes.

Supporting Information
(PDF)