BACKGROUND
high throughput technologies such as microarrays and yeast-two-hybrid assays have resulted in an explosion of biological data that can be represented as networks.
for example, microarray datasets can be analyzed as a coexpression network, in which nodes  represent genes and links  represent coexpression, the similarity of the level of expression of two genes over the samples in the study.
similarly, protein interaction data, such as that generated by yeast-two-hybrid assays, can be summarized as a network, with nodes representing proteins and edges representing physical interaction between two proteins.
genes and their products give rise to biological function through their interaction with each other and with other components of the cell.
the analysis of the above biological networks is therefore the natural way to understand cellular function on a genome-wide level.
in particular, we need a thorough understanding of the statistical properties of biological networks if we aim to make inferences, such as inferring evolutionary relationships between various networks, or separating signal from noise in imperfect network data.
erdős and rényi  <cit>  were the first to study the statistical properties of random graph models.
in their models , any edge between two vertices occurred independently of other edges with a constant probability p. in these graphs, however, the degree of a vertex  is a random variable with an approximately poisson distribution with λ = p, which is grossly at odds with most biological network observed to date  <cit> .
in real biological data, node degrees usually have heavy tail distributions  <cit> .
accordingly, in most statistical studies of biological networks, the null model is a random graph from the set having a degree distribution identical to that of the data, or a distribution in which the expected degrees are identical to those observed in the data  <cit> .
these models are themselves limited, because in addition to their degree distributions, biological networks show highly clustered connections  <cit>  and transitivity  <cit> .
indeed, it is difficult to assess which properties of a network would represent sufficient statistics that are biologically meaningful.
mahadevan et al.
<cit>  attempted to solve this problem by devising an increasing series of random network models they referred as the dk-series.
the distributions of the random networks are defined as uniform over the set of graphs having the same distribution of d-sized subgraphs as the observed network data.
particular cases of the series reduce to familiar distributions: the 0k distribution p <dig> is identical to the corresponding er distribution, which describes the average number of links per node.
the 1k distribution model tells us the expected degree of each node and assumes that the nodes are randomly connected conditional on the expected node degrees.
the 2k distribution p <dig> describes the interconnectivity of nodes with given degrees, maintaining the number ) of links between nodes of degrees k and k'.
the 2k distribution therefore preserves degree-degree correlations between nodes .
including still more connectivity information, the 3k distribution considers degree correlations among any  <dig> nodes, which include the transitivity of the network.
moving beyond pairs of nodes, various topological structures are possible.
for example, there are  <dig> different kinds of isomorphic structures for the 3k distribution.
increasingly larger subgraphs can be enumerated for d =  <dig> ,..., capturing increasingly complex features of a particular graph.
the dk-series is therefore an objective way to progressively include more features into a random graph model, just as each term in a fourier or taylor series progressively captures more details of a given function, and thus largely avoids the arbitrary selection of statistics that may or may not be sufficient or relevant to a particular process.
using this series as our starting point, we sought to evaluate the use of ever more inclusive dk distributions in the study of biological networks.
for four different biological networks, we trained dk models for d =  <dig>   <dig>   <dig>   <dig>  we first explored the properties of the models by evaluating their ability to predict in the observed networks the presence or absence of individual edges, as well as general network statistics.
we showed that the 2k model outperforms other models in predicting the presence and absence of the edges for both protein interaction and gene coexpression networks.
we then evaluated whether statistical significance against one of the models for subnetworks corresponded to biological significance.
we modeled our approach based on the scoring scheme used by tanay et al.
<cit> : they devised a pseudo-likelihood score for edges in a bipartite graph of genes and samples, in which edges occur according to a null model that corresponds to the 1k distribution, or to an alternative model, representing biological significance, independently with a high constant probability p =  <dig> .
they showed that this score results in improved accuracy in predicting functional gene groups, when compared with network density alone .
we reproduced their score using as the null distribution one of dk models for d =  <dig>   <dig>   <dig>   <dig>  we aimed to test the hypothesis that more inclusive distributions would result in a score for a set of nodes that is more indicative of biological significance, just as 1k was in the case of bipartite graphs.
we were surprised, however, to find that accuracy was only slightly increased with each successive dk distribution in the case of yeast protein interaction networks, while the 0k distribution  had the best predictive power in coexpression networks.
RESULTS
in this section, we first give a brief discussion of the dk models and the pseudo-likelihood methods for estimating the parameters.
next we study the accuracy of the dk models in predicting edges in the molecular networks.
then several statistics related to the networks are studied to evaluate if the random networks can approximate the observed networks.
finally, we evaluate if the dk models can be used to identify functionally homogeneous modules.
model description and parameter estimation
for each network, we created random graph models matching the 0k, 1k, 2k, and 3k distributions of the observed network.
for the 0k and 1k models, we took the degree sequence and number of edges as fixed properties of the network and thus defined the matching models.
in the 0k model, each edge occurs independently with probability , |e| being the number of edges and |v| the number of nodes in the real network.
in the 1k model, each edge occurs independently, conditional on the degrees k <dig> and k <dig> of its incident nodes, with probability p = min ,  <dig> - ε), for ε small .
for the 2k model, we calculated the probability p that two proteins with degree pairs  interact.
one intuitive approach is to estimate p by the fraction of interacting protein pairs among all the protein pairs with degrees , k <dig> ≤ k <dig>  however, for many degree pairs , the number of such protein pairs is small.
thus, the estimated value of p using this intuitive approach is not reliable.
to overcome this problem, we modeled p as a function of  and fitted the function using matlab.
details of the estimation method are given in the "methods" section.
the 3k model describes how protein triplets with degrees  interact with each other.
there are a total of eight possible interaction patterns among the three proteins as shown in figure  <dig>  as in the 2k model, directly estimating the eight probabilities corresponding to the interacting patterns is difficult due to the small number of protein triplets for many degree triplets .
thus, we reparameterized the probabilities to fit a logistic regression model, which is necessary to improve probability estimates for degree triplets for which we only have one or few examples.
details are given in the "methods" section.
the performance of predicting protein interactions using the dk distribution models
we began by studying the ability of the dk distribution models in their capacity to predict the presence and absence of links in the observed biological networks.
for each dk  distribution model, we predicted the probability that any pair of proteins interact using the estimated parameters obtained above.
given a cut-off threshold, protein pairs with interaction probability above the threshold were predicted to interact.
the predicted interaction pairs were compared with the observed protein interactions to study the ability of the dk distribution model to predict protein interactions.
to overcome the problem of extreme large number of non-interacting protein pairs over that of the interacting pairs, we randomly chose the same number of non-interacting protein pairs as the number of interacting pairs for the comparisons.
three different evaluation methods were used: the accuracy, the receiver operation curve , and the precision-recall curve.
figure  <dig> and figure  <dig> show the performance of the dk models in predicting protein interactions using the mips  <cit>  protein interaction data and a gene co-expression network based on gds <dig>  <cit> , respectively.
the performance of the 2k distribution model always outperformed the 1k distribution model, which in turn outperformed the 0k distribution.
the 3k distribution model performed comparably to the 2k distribution model in a coexpression network , but significantly worse in the case of a protein interaction network .
similar results were obtained for the biogrid  <cit>  protein interaction data  and other coexpression networks .
comparing statistical features of random networks from the dk models with that of the observed networks
we next studied if the random networks based on the dk distribution models approximate the observed interaction networks.
to achieve this objective, we generated  <dig> random networks based on the dk distribution models and calculated several statistical features of the random networks .
we studied five network statistical features as in  <cit> :
• λ1: average of the smallest eigenvalue of the laplacian of the graph matrix;
• λn-1: average of the largest eigenvalue of the laplacian of the graph matrix;
• d: average shortest distance between the nodes;
• σd: standard deviation of shortest distance between the nodes;
• r: average assortativity coefficients.
the elements of the laplacian matrix of a network are defined by  if node i with degree ki and node j with degree kj are connected and lij =  <dig> otherwise for i ≠ j, and lij =  <dig> if i = j. several other important network statistics  <cit> , e.g.
network resilience and performance, are tightly controlled by the smallest non-zero  and the largest  eigenvalues of the laplacian matrix.
therefore, we studied whether the corresponding eigenvalues of the dk random networks are close to that of the true network.
in addition, the distribution of the shortest distances between any two nodes provides information on how the nodes cluster together in the network.
we used two quantities, the mean and standard deviation of the shortest distances, to characterize this distribution.
finally, the assortativity coefficient of a network provides information on how nodes of different degrees link to each other.
although these five network statistical measures cannot fully describe the network of interest, they capture important network properties.
if the dk distribution models can approximate the true network well, these quantities in the dk random networks should be close to the corresponding values of the true network.
tables  <dig> and  <dig> give the average and standard deviation of the corresponding feature values from the  <dig> random networks for each of the dk distribution models based on the mips protein interaction network and gds <dig> coexpression network, respectively.
for the mips protein interaction data, the average values of the five network features for the random networks converge to the corresponding values of the observed network, indicating that the dk model converges to the true network as d increases.
on the other hand, for the coexpression network, the average values of the five network features of the 2k model are closer to the corresponding features in the true network than the 1k and 3k models.
the poor performance of the 3k model maybe due to the fact that the estimated parameters are not accurate due to the relative small number of the nodes in the coexpression networks.
similar tables for the other networks are provided in additional file 2: tables s1-s <dig>
λ1: average of the smallest eigenvalue of the laplacian of the graph matrix; λn-1: average of the largest eigenvalue of the laplacian of the graph matrix; d: average shortest distance between the nodes; σd: standard deviation of shortest distance between the nodes; r: average assortativity coefficients.
the quantity in the brackets indicate the standard deviation of corresponding metric
notations are the same as in table  <dig>
the performance of the dk distribution models for the identification of functionally homogeneous modules
our primary motivation of this study is to see if the more complex models, which can generally more accurately describe the observed network, are helpful in the identification of biologically functionally homogeneous modules.
statistical deviations from a suitable model would indicate evolutionary pressure and thus functional significance.
therefore, we can compare the functional relevance of each model by how well statistical deviations from the model correlate with the functional homogeneity of the corresponding nodes.
we designed scores from our models based on an alternative hypothesis that edges are present in a functional module with constant probability p. generally p should be close to  <dig> as most functionally homogeneous modules are highly clustered.
as in  <cit> , we chose p =  <dig>  in the main text.
to see the validity of our results for different values of p, we also changed p to p =  <dig>  and p =  <dig> .
our approach is similar to that used by tanay et al  <cit> , which they used a single null model  in the context of a bipartite graph.
with this score framework, we used a simulated annealing algorithm to find groups of genes with high scores, retaining every group encountered during the run of the algorithm and their scores under each of the null models.
finally, we called a gene group functionally homogeneous if it was enriched in at least one functional category from the gene ontology  <cit> .
we defined module enrichment by the hypergeometric test p-value, with a threshold of p < 10- <dig>  these gene groups were taken to be true positives, and the remaining gene groups were taken to be true negatives.
again, we varied this threshold from  <dig>  to 10- <dig>  and no qualitative changes in the results were observed, showing that our approach is robust to the parameter for calling functional homogeneity .
we then evaluated the four models by comparing how well they can predict functional homogeneity in the mips and biogrid yeast protein interaction networks, or in two different yeast gene coexpression networks, gds <dig> and gds <dig>
results based on the mips protein interaction network
we used the mips  <cit>  yeast protein interaction network to compare the ability of the dk distributions to predict functional homogeneity in gene groups of size  <dig>  note that the size of the gene groups cannot be too small.
otherwise it is very hard to distinguish functionally homogeneous modules from random gene groups.
the score for each gene group depends on the null model for the network.
we calculated the score for gene groups of size  <dig> using dk distribution model as the null network model, d =  <dig>   <dig>   <dig>   <dig>  a gene set was predicted to be functional homogeneous if the score is above a cut-off threshold.
our objective is to see which score functions can more accurately predict functional homogeneity.
therefore, we compare the predicted "functionally homogeneous" gene groups with the positive groups.
we measure the performance using accuracy, the roc curve, and the precision-recall curve as above.
note that the subnetwork scores defined based on dk models in equation  for d =  <dig>   <dig>   <dig> and equation  for d =  <dig> are not on the same scale and thus the prediction accuracy in figure 4a is not comparable for the same cut-off value of the subnetwork scores.
however, the maximum prediction accuracy for the dk models can be compared.
figure  <dig> shows the results based on the mips interaction data with p =  <dig>  and gene group size  <dig>  the results based on other combinations of p =  <dig> ,  <dig> ,  <dig>  and gene set size n =  <dig>   <dig> are given in additional file 1: figures s5-s <dig>  the corresponding results based on the biogrid  <cit>  protein interaction data are also given as additional file 1: figures s8-s <dig>
we found almost no difference between the performance of the different scores, with even the 0k model  performing only slightly worse than the rest.
figure 4a shows that the highest prediction accuracy for the 0k model is slightly smaller than that of the 1k-3k models, and the highest prediction accuracy of the 1k-3k models are similar.
similarly, the roc curve  and the precision-recall curve  of the 0k model are slightly lower than the corresponding curves for the 1k-3k models.
this result shows that the added information of degree correlations and transitivity do not influence the model enough to have a significant effect on the global prediction performance of functional homogeneity.
because the performance of the different dk distribution models are so similar, we hypothesized that the scores for the gene groups based on different dk models maybe highly correlated.
we tested this hypothesis by studying the spearman's correlation between the scores for the gene groups using different models.
the results are given in table  <dig>  it can be seen that they are indeed highly correlated.
similar results were obtained based on the biogrid protein interaction network .
gs: group size
results based on gene coexpression networks
we repeated the above analysis, this time using a yeast coexpression network built from a yeast microarray dataset, gds <dig>  downloaded from the ncbi gene expression omnibus  <cit> .
we constructed the network by calculating the pearson correlation coefficient  between the expression levels of every pair of genes.
then, each gene was represented as a node in the network, and we drew a link between two nodes if the pcc exceeded a certain threshold.
to make the network as comparable as possible to the mips network, we selected the threshold such that the degree of the most highly connected node in each network was the same .
we were again surprised to find that the 2k performance was similar to 1k, while the 3k model's performance was actually worse than the 1k and 2k models .
this observation maybe explained by the fact that the 3k model does not approximate the observed network well as shown in table  <dig>  perhaps most strikingly, the 0k model displayed the best performance, showing that a simple measure of density is a very good predictor of function in coexpression networks.
we also changed the threshold for the pcc between expression profiles to build the network to  <dig> , the gene group size n =  <dig>  and the parameter p to p =  <dig> ,  <dig>  in defining the score function.
the results are presented in additional file 1: figures s12-s <dig>  same qualitative results were obtained.
we also studied the spearman's correlation between the gene group scores for different dk models and the results are given in table  <dig>  although the scores based on 1k, 2k, and 3k models are highly correlated, they do not strongly correlate with the scores based on the 0k model.
similar results were obtained for different pcc cut-off thresholds .
gs: group size
we also performed the same analysis for another gene expression dataset, gds <dig>  <cit> .
the performance results for this dataset are presented in additional file 1: figures s19-s <dig> and additional file 2: tables s7-s <dig>  it should be noted that none of the dk models performed well in identifying functionally homogeneous modules based on this gene expression data set.
one potential reason is that the number of sampling points is only  <dig>  which is much smaller than that of gds <dig>  which has  <dig> sampling points.
thus the network constructed based on gds <dig> may not be reliable.
despite the drawbacks of this dataset, the conclusions from this dataset is qualitatively identical to those found for gds <dig>  this demonstrates the generality of our conclusions with respect to gene coexpression networks.
CONCLUSIONS
we studied the ability of dk distribution models to predict individual edges and functionally homogeneous modules in protein interaction and gene coexpression networks.
a pseudo-likelihood logistic estimation method was proposed to estimate the parameters in the dk distribution models.
we found that the 2k distribution model performs the best in predicting individual edges in both protein interaction and gene coexpression networks.
a pseudo-likelihood ratio score function was then defined to evaluate potential functional homogeneity based on the dk distribution models.
for yeast protein interaction networks, 1k, 2k and 3k models perform similarly and are slightly better than the 0k model in predicting functionally homogeneous modules.
the dk scores were very highly correlated for different d. this means that, between two different subgraph topologies, the variation in the denominator, the dk distribution likelihood, was small relative to that in the numerator, the constant-p likelihood.
in this case, most of the variation in scores between modules would be accounted for by the numerator.
the different probabilities between 1k, 2k, and 3k may be similar overall in the networks we studied.
for gene coexpression networks, the 0k model performs significantly better than the other models in predicting functionally homogeneous modules.
we noted that 0k, or density, performed remarkably well as a prediction method even in the yeast protein interaction network, being able to find extremely functionally homogeneous groups of genes .
this may simply reflect that highly dense subnetworks in a protein interaction network represent protein complexes, which are of necessity functionally homogeneous.
one future avenue of research could be to remove this type of functionally homogeneous modules from the data, since they are relatively uninteresting examples of functional homogeneity.
it may be that the subtle differences between the various dk distributions are useful to pick out homogeneous modules of more specific functions.

 high throughput technologies such as microarrays and yeast-two-hybrid assays have resulted in an explosion of biological data that can be represented as networks. for example, microarray datasets can be analyzed as a coexpression network, in which nodes represent genes and links represent coexpression, the similarity of the level of expression of two genes over the samples in the study. similarly, protein interaction data, such as that generated by yeast-two-hybrid assays, can be summarized as a network, with nodes representing proteins and edges representing physical interaction between two proteins. genes and their products give rise to biological function through their interaction with each other and with other components of the cell. the analysis of the above biological networks is therefore the natural way to understand cellular function on a genome-wide level. in particular, we need a thorough understanding of the statistical properties of biological networks if we aim to make inferences, such as inferring evolutionary relationships between various networks, or separating signal from noise in imperfect network data. erdős and rényi were the first to study the statistical properties of random graph models. in their models, any edge between two vertices occurred independently of other edges with a constant probability p. in these graphs, however, the degree of a vertex is a random variable with an approximately poisson distribution with λ = p, which is grossly at odds with most biological network observed to date. in real biological data, node degrees usually have heavy tail distributions. accordingly, in most statistical studies of biological networks, the null model is a random graph from the set having a degree distribution identical to that of the data, or a distribution in which the expected degrees are identical to those observed in the data. these models are themselves limited, because in addition to their degree distributions, biological networks show highly clustered connections and transitivity. indeed, it is difficult to assess which properties of a network would represent sufficient statistics that are biologically meaningful. mahadevan et al. attempted to solve this problem by devising an increasing series of random network models they referred as the dk-series. the distributions of the random networks are defined as uniform over the set of graphs having the same distribution of d-sized subgraphs as the observed network data. particular cases of the series reduce to familiar distributions: the 0k distribution p is identical to the corresponding er distribution, which describes the average number of links per node. the 1k distribution model tells us the expected degree of each node and assumes that the nodes are randomly connected conditional on the expected node degrees. the 2k distribution p describes the interconnectivity of nodes with given degrees, maintaining the number ) of links between nodes of degrees k and k'. the 2k distribution therefore preserves degree-degree correlations between nodes. including still more connectivity information, the 3k distribution considers degree correlations among any nodes, which include the transitivity of the network. moving beyond pairs of nodes, various topological structures are possible. for example, there are different kinds of isomorphic structures for the 3k distribution. increasingly larger subgraphs can be enumerated for d =,..., capturing increasingly complex features of a particular graph. the dk-series is therefore an objective way to progressively include more features into a random graph model, just as each term in a fourier or taylor series progressively captures more details of a given function, and thus largely avoids the arbitrary selection of statistics that may or may not be sufficient or relevant to a particular process. using this series as our starting point, we sought to evaluate the use of ever more inclusive dk distributions in the study of biological networks. for four different biological networks, we trained dk models for d = we first explored the properties of the models by evaluating their ability to predict in the observed networks the presence or absence of individual edges, as well as general network statistics. we showed that the 2k model outperforms other models in predicting the presence and absence of the edges for both protein interaction and gene coexpression networks. we then evaluated whether statistical significance against one of the models for subnetworks corresponded to biological significance. we modeled our approach based on the scoring scheme used by tanay et al.: they devised a pseudo-likelihood score for edges in a bipartite graph of genes and samples, in which edges occur according to a null model that corresponds to the 1k distribution, or to an alternative model, representing biological significance, independently with a high constant probability p =. they showed that this score results in improved accuracy in predicting functional gene groups, when compared with network density alone. we reproduced their score using as the null distribution one of dk models for d = we aimed to test the hypothesis that more inclusive distributions would result in a score for a set of nodes that is more indicative of biological significance, just as 1k was in the case of bipartite graphs. we were surprised, however, to find that accuracy was only slightly increased with each successive dk distribution in the case of yeast protein interaction networks, while the 0k distribution had the best predictive power in coexpression networks.  in this section, we first give a brief discussion of the dk models and the pseudo-likelihood methods for estimating the parameters. next we study the accuracy of the dk models in predicting edges in the molecular networks. then several statistics related to the networks are studied to evaluate if the random networks can approximate the observed networks. finally, we evaluate if the dk models can be used to identify functionally homogeneous modules. model description and parameter estimation for each network, we created random graph models matching the 0k, 1k, 2k, and 3k distributions of the observed network. for the 0k and 1k models, we took the degree sequence and number of edges as fixed properties of the network and thus defined the matching models. in the 0k model, each edge occurs independently with probability, |e| being the number of edges and |v| the number of nodes in the real network. in the 1k model, each edge occurs independently, conditional on the degrees k and k of its incident nodes, with probability p = min, - ε), for ε small. for the 2k model, we calculated the probability p that two proteins with degree pairs interact. one intuitive approach is to estimate p by the fraction of interacting protein pairs among all the protein pairs with degrees, k ≤ k however, for many degree pairs, the number of such protein pairs is small. thus, the estimated value of p using this intuitive approach is not reliable. to overcome this problem, we modeled p as a function of and fitted the function using matlab. details of the estimation method are given in the "methods" section. the 3k model describes how protein triplets with degrees interact with each other. there are a total of eight possible interaction patterns among the three proteins as shown in figure as in the 2k model, directly estimating the eight probabilities corresponding to the interacting patterns is difficult due to the small number of protein triplets for many degree triplets. thus, we reparameterized the probabilities to fit a logistic regression model, which is necessary to improve probability estimates for degree triplets for which we only have one or few examples. details are given in the "methods" section. the performance of predicting protein interactions using the dk distribution models we began by studying the ability of the dk distribution models in their capacity to predict the presence and absence of links in the observed biological networks. for each dk distribution model, we predicted the probability that any pair of proteins interact using the estimated parameters obtained above. given a cut-off threshold, protein pairs with interaction probability above the threshold were predicted to interact. the predicted interaction pairs were compared with the observed protein interactions to study the ability of the dk distribution model to predict protein interactions. to overcome the problem of extreme large number of non-interacting protein pairs over that of the interacting pairs, we randomly chose the same number of non-interacting protein pairs as the number of interacting pairs for the comparisons. three different evaluation methods were used: the accuracy, the receiver operation curve, and the precision-recall curve. figure and figure show the performance of the dk models in predicting protein interactions using the mips protein interaction data and a gene co-expression network based on gds, respectively. the performance of the 2k distribution model always outperformed the 1k distribution model, which in turn outperformed the 0k distribution. the 3k distribution model performed comparably to the 2k distribution model in a coexpression network, but significantly worse in the case of a protein interaction network. similar results were obtained for the biogrid protein interaction data and other coexpression networks. comparing statistical features of random networks from the dk models with that of the observed networks we next studied if the random networks based on the dk distribution models approximate the observed interaction networks. to achieve this objective, we generated random networks based on the dk distribution models and calculated several statistical features of the random networks. we studied five network statistical features as in: • λ1: average of the smallest eigenvalue of the laplacian of the graph matrix; • λn-1: average of the largest eigenvalue of the laplacian of the graph matrix; • d: average shortest distance between the nodes; • σd: standard deviation of shortest distance between the nodes; • r: average assortativity coefficients. the elements of the laplacian matrix of a network are defined by if node i with degree ki and node j with degree kj are connected and lij = otherwise for i ≠ j, and lij = if i = j. several other important network statistics, e.g. network resilience and performance, are tightly controlled by the smallest non-zero and the largest eigenvalues of the laplacian matrix. therefore, we studied whether the corresponding eigenvalues of the dk random networks are close to that of the true network. in addition, the distribution of the shortest distances between any two nodes provides information on how the nodes cluster together in the network. we used two quantities, the mean and standard deviation of the shortest distances, to characterize this distribution. finally, the assortativity coefficient of a network provides information on how nodes of different degrees link to each other. although these five network statistical measures cannot fully describe the network of interest, they capture important network properties. if the dk distribution models can approximate the true network well, these quantities in the dk random networks should be close to the corresponding values of the true network. tables and give the average and standard deviation of the corresponding feature values from the random networks for each of the dk distribution models based on the mips protein interaction network and gds coexpression network, respectively. for the mips protein interaction data, the average values of the five network features for the random networks converge to the corresponding values of the observed network, indicating that the dk model converges to the true network as d increases. on the other hand, for the coexpression network, the average values of the five network features of the 2k model are closer to the corresponding features in the true network than the 1k and 3k models. the poor performance of the 3k model maybe due to the fact that the estimated parameters are not accurate due to the relative small number of the nodes in the coexpression networks. similar tables for the other networks are provided in additional file 2: tables s1-s λ1: average of the smallest eigenvalue of the laplacian of the graph matrix; λn-1: average of the largest eigenvalue of the laplacian of the graph matrix; d: average shortest distance between the nodes; σd: standard deviation of shortest distance between the nodes; r: average assortativity coefficients. the quantity in the brackets indicate the standard deviation of corresponding metric notations are the same as in table the performance of the dk distribution models for the identification of functionally homogeneous modules our primary motivation of this study is to see if the more complex models, which can generally more accurately describe the observed network, are helpful in the identification of biologically functionally homogeneous modules. statistical deviations from a suitable model would indicate evolutionary pressure and thus functional significance. therefore, we can compare the functional relevance of each model by how well statistical deviations from the model correlate with the functional homogeneity of the corresponding nodes. we designed scores from our models based on an alternative hypothesis that edges are present in a functional module with constant probability p. generally p should be close to as most functionally homogeneous modules are highly clustered. as in, we chose p = in the main text. to see the validity of our results for different values of p, we also changed p to p = and p =. our approach is similar to that used by tanay et al, which they used a single null model in the context of a bipartite graph. with this score framework, we used a simulated annealing algorithm to find groups of genes with high scores, retaining every group encountered during the run of the algorithm and their scores under each of the null models. finally, we called a gene group functionally homogeneous if it was enriched in at least one functional category from the gene ontology. we defined module enrichment by the hypergeometric test p-value, with a threshold of p < 10- these gene groups were taken to be true positives, and the remaining gene groups were taken to be true negatives. again, we varied this threshold from to 10- and no qualitative changes in the results were observed, showing that our approach is robust to the parameter for calling functional homogeneity. we then evaluated the four models by comparing how well they can predict functional homogeneity in the mips and biogrid yeast protein interaction networks, or in two different yeast gene coexpression networks, gds and gds results based on the mips protein interaction network we used the mips yeast protein interaction network to compare the ability of the dk distributions to predict functional homogeneity in gene groups of size note that the size of the gene groups cannot be too small. otherwise it is very hard to distinguish functionally homogeneous modules from random gene groups. the score for each gene group depends on the null model for the network. we calculated the score for gene groups of size using dk distribution model as the null network model, d = a gene set was predicted to be functional homogeneous if the score is above a cut-off threshold. our objective is to see which score functions can more accurately predict functional homogeneity. therefore, we compare the predicted "functionally homogeneous" gene groups with the positive groups. we measure the performance using accuracy, the roc curve, and the precision-recall curve as above. note that the subnetwork scores defined based on dk models in equation for d = and equation for d = are not on the same scale and thus the prediction accuracy in figure 4a is not comparable for the same cut-off value of the subnetwork scores. however, the maximum prediction accuracy for the dk models can be compared. figure shows the results based on the mips interaction data with p = and gene group size the results based on other combinations of p =,, and gene set size n = are given in additional file 1: figures s5-s the corresponding results based on the biogrid protein interaction data are also given as additional file 1: figures s8-s we found almost no difference between the performance of the different scores, with even the 0k model performing only slightly worse than the rest. figure 4a shows that the highest prediction accuracy for the 0k model is slightly smaller than that of the 1k-3k models, and the highest prediction accuracy of the 1k-3k models are similar. similarly, the roc curve and the precision-recall curve of the 0k model are slightly lower than the corresponding curves for the 1k-3k models. this result shows that the added information of degree correlations and transitivity do not influence the model enough to have a significant effect on the global prediction performance of functional homogeneity. because the performance of the different dk distribution models are so similar, we hypothesized that the scores for the gene groups based on different dk models maybe highly correlated. we tested this hypothesis by studying the spearman's correlation between the scores for the gene groups using different models. the results are given in table it can be seen that they are indeed highly correlated. similar results were obtained based on the biogrid protein interaction network. gs: group size results based on gene coexpression networks we repeated the above analysis, this time using a yeast coexpression network built from a yeast microarray dataset, gds downloaded from the ncbi gene expression omnibus. we constructed the network by calculating the pearson correlation coefficient between the expression levels of every pair of genes. then, each gene was represented as a node in the network, and we drew a link between two nodes if the pcc exceeded a certain threshold. to make the network as comparable as possible to the mips network, we selected the threshold such that the degree of the most highly connected node in each network was the same. we were again surprised to find that the 2k performance was similar to 1k, while the 3k model's performance was actually worse than the 1k and 2k models. this observation maybe explained by the fact that the 3k model does not approximate the observed network well as shown in table perhaps most strikingly, the 0k model displayed the best performance, showing that a simple measure of density is a very good predictor of function in coexpression networks. we also changed the threshold for the pcc between expression profiles to build the network to, the gene group size n = and the parameter p to p =, in defining the score function. the results are presented in additional file 1: figures s12-s same qualitative results were obtained. we also studied the spearman's correlation between the gene group scores for different dk models and the results are given in table although the scores based on 1k, 2k, and 3k models are highly correlated, they do not strongly correlate with the scores based on the 0k model. similar results were obtained for different pcc cut-off thresholds. gs: group size we also performed the same analysis for another gene expression dataset, gds. the performance results for this dataset are presented in additional file 1: figures s19-s and additional file 2: tables s7-s it should be noted that none of the dk models performed well in identifying functionally homogeneous modules based on this gene expression data set. one potential reason is that the number of sampling points is only which is much smaller than that of gds which has sampling points. thus the network constructed based on gds may not be reliable. despite the drawbacks of this dataset, the conclusions from this dataset is qualitatively identical to those found for gds this demonstrates the generality of our conclusions with respect to gene coexpression networks.  we studied the ability of dk distribution models to predict individual edges and functionally homogeneous modules in protein interaction and gene coexpression networks. a pseudo-likelihood logistic estimation method was proposed to estimate the parameters in the dk distribution models. we found that the 2k distribution model performs the best in predicting individual edges in both protein interaction and gene coexpression networks. a pseudo-likelihood ratio score function was then defined to evaluate potential functional homogeneity based on the dk distribution models. for yeast protein interaction networks, 1k, 2k and 3k models perform similarly and are slightly better than the 0k model in predicting functionally homogeneous modules. the dk scores were very highly correlated for different d. this means that, between two different subgraph topologies, the variation in the denominator, the dk distribution likelihood, was small relative to that in the numerator, the constant-p likelihood. in this case, most of the variation in scores between modules would be accounted for by the numerator. the different probabilities between 1k, 2k, and 3k may be similar overall in the networks we studied. for gene coexpression networks, the 0k model performs significantly better than the other models in predicting functionally homogeneous modules. we noted that 0k, or density, performed remarkably well as a prediction method even in the yeast protein interaction network, being able to find extremely functionally homogeneous groups of genes. this may simply reflect that highly dense subnetworks in a protein interaction network represent protein complexes, which are of necessity functionally homogeneous. one future avenue of research could be to remove this type of functionally homogeneous modules from the data, since they are relatively uninteresting examples of functional homogeneity. it may be that the subtle differences between the various dk distributions are useful to pick out homogeneous modules of more specific functions. 
recent studies have shown that the degree distribution of the nodes is not an adequate statistic in many molecular networks. studying the statistical properties of these networks in turn allows us to infer biological function. many aspects of biological functions can be modeled by biological networks, such as protein interaction networks, metabolic networks, and gene coexpression networks. our use of dk models showed that incorporation of degree correlations could increase predictive power in some contexts, albeit sometimes marginally, but, in all contexts, the use of third-order degree correlations decreased accuracy. the approach was used to analyze the mips and biogrid yeast protein interaction networks, and two yeast coexpression networks. we showed that 2nd order degree correlation information gave better predictions of gene interactions in both protein interaction and gene coexpression networks. however, in the biologically important task of predicting functionally homogeneous modules, degree correlation information performs marginally better in the case of the mips and biogrid protein interaction networks, but worse in the case of gene coexpression networks. we sought to extend this statistic with 2nd and 3rd order degree correlations and developed a pseudo-likelihood approach to estimate the parameters. however, it is possible that other parameter estimation methods, such as maximum likelihood, will show the usefulness of incorporating 2nd and 3rd degree correlations in predicting functionally homogeneous modules.
 many aspects of biological functions can be modeled by biological networks, such as protein interaction networks, metabolic networks, and gene coexpression networks. studying the statistical properties of these networks in turn allows us to infer biological function. complex statistical network models can potentially more accurately describe the networks, but it is not clear whether such complex models are better suited to find biologically meaningful subnetworks.  recent studies have shown that the degree distribution of the nodes is not an adequate statistic in many molecular networks. we sought to extend this statistic with 2nd and 3rd order degree correlations and developed a pseudo-likelihood approach to estimate the parameters. the approach was used to analyze the mips and biogrid yeast protein interaction networks, and two yeast coexpression networks. we showed that 2nd order degree correlation information gave better predictions of gene interactions in both protein interaction and gene coexpression networks. however, in the biologically important task of predicting functionally homogeneous modules, degree correlation information performs marginally better in the case of the mips and biogrid protein interaction networks, but worse in the case of gene coexpression networks.  our use of dk models showed that incorporation of degree correlations could increase predictive power in some contexts, albeit sometimes marginally, but, in all contexts, the use of third-order degree correlations decreased accuracy. however, it is possible that other parameter estimation methods, such as maximum likelihood, will show the usefulness of incorporating 2nd and 3rd degree correlations in predicting functionally homogeneous modules. 