BACKGROUND
a major challenge in systems biology is to infer mechanistic gene interactions from high throughput microarray data  <cit> .
underlying this challenge is the problem to find causal or regulatory relationships between genes.
robust solutions to this problem would provide us with a transcriptomic map of a genome that allows us to accurately predict the effect of gene perturbations.
previous efforts to detect mechanistic relationships from gene expression data can be broadly divided into linear correlation and probability based methods.
linear correlation based methods, such as clustering  <cit> , correlation networks  <cit>  and graphical gaussian models  <cit> , have a long and fruitful history in statistical modeling and bioinformatics.
these linear methods are computationally fast and relatively easy to interpret.
however, a key limitation with these methods is that they assume linear relationships between variables.
while some components of any transcriptional regulatory network are linear, nonlinear events such as or, and, and xor type transcriptional regulation are relatively commonplace  <cit> .
these nonlinear interactions would not be captured with a linear model, leading to spurious relationships between variables.
probability based methods have also been used to detect relationships between genes.
these probability methods include probabilistic boolean networks   <cit> , bayesian networks  <cit>  and mutual information networks  <cit> .
probability based methods can capture both linear and nonlinear regulatory relationships and are noise tolerant.
however, many of the current probability based tools used in systems biology suffer from the following three limitations:  data discretization  <cit> ,  pairwise testing  <cit> ,  emphasis on correlation over causality  <cit> .
to transform continuous data into a more easily computable form, most probabilistic methods require the data to first be discretized into a finite number of bins, such as high, medium, and low  <cit> : the number of bins used in discretization is difficult to choose, and is generally selected at some consistent yet arbitrary point.
unfortunately, different binning procedures can produce different analysis results  <cit> , suggesting that the act of binning alone introduces errors into the analysis.
methods that search for pairwise associations only focus on a single relationship between regulator and target at a time.
pairwise association networks have been created using classical mutual information  <cit> .
however, simple pairwise relationships are likely less common than multivariate relationships in real biological systems, as the expression of most genes is regulated not by a single gene but more likely by multiple genes.
methods that allow multivariate interactions such as bayesian networks or some fuzzy logic approaches  <cit>  are inherently superior in this respect.
a final challenge in creating mechanistically predictive transcriptional models is the ability to identify not just correlative but also causal models.
for the definition of causal relationship, we adopt the notion of probabilistic causation  <cit> .
informally, event a  causes event b if and only if a's occurrence alters  the probability of b. this sometimes reflects imperfect knowledge  of a deterministic system but more frequently suggests a stochastic nature of the causal system under study.
although difficult, causal relationships have been learned properly from non-sequential observational data  <cit> .
probabilistic graphical modeling methods like bayesian networks have been used to infer causal models from gene expression data  <cit> .
however, many probabilistic approaches are able to make correlative networks but not necessarily causal networks  <cit> .
their multivariate scoring metrics such as conditional probability and mutual information are still generalized two-way correlation between the target and the parent set.
similar to the classical two-way metrics, these generalized correlations alone cannot differentiate between a causal versus confounding  parent set.
true causal relationships like genetic regulation feature positive higher order interaction  <cit> , the non-additive effect above the sum of the lower order interactions  <cit> .
for instance, for regulation involved two regulators such as or, and, xor type relationships, two regulators together account for much more in the target than they individually can .
intuitively such non-additive effect can be described as coordination or synergy between parents .
on the other hand, confounding models commonly have no or negative higher order interaction .
we propose that with such high order interaction considered, we can better differentiate true causal model versus confounding models.
in this work, we demonstrate a novel algorithm that attempts to overcome all three limitations using a continuous high order mutual information based scoring metric we call mi <dig> .
note that continuous two-way mutual information has been described previously  <cit> .
high order interaction information  has been employed to model complex interactions  <cit> .
however, both two-way mutual information and high order interaction information are symmetric and as such unable to make causal statements.
mi <dig> combines 3rd order interaction information with the asymmetric mutual information between target and regulator set to account for the direction of regulation.
mi <dig> is novel as a combinatorial probabilistic metric and an integrated statistical learning method.
in this work, we compare mi <dig> to other probability based methods quantitatively and qualitatively using synthetic data where the true model is known.
next we apply mi <dig> and control methods to reconstruct regulatory networks centered at the transcription factor myc from a published high throughput microarray dataset  <cit> .
the learning results are then evaluated numerically and biologically.
learning myc centered transcriptional regulatory network represents an ideal test case for mi <dig> as myc is a well characterized transcriptional regulator that acts in tandem with a finite set of co-effectors and regulates the expression of a large group of genes  <cit> .
myc has been well investigated  <cit>  and online databases of myc targets  <cit>  are available for validation purpose.
despite these efforts, many cofactors and targets remain unidentified, and corresponding regulatory mechanisms unknown  <cit> .
as a result, an integrated understanding of myc dependent transcriptional regulation has remained out of reach  <cit> .
in this study, we use mi <dig> to derive an accurate transcriptomic map surrounding myc from the same gene expression dataset used to identify myc targets  <cit> .
the approaches used here are general and can be directly used for any transcriptional regulator given sufficient gene expression data.
RESULTS
mi <dig> validation with synthetic data
we validated mi <dig> against other commonly used methods listed in table  <dig>  including a discrete version of mi <dig> , two-way mutual information  and a log conditional probability score used in bayesian network  learning.
each control method represents one of the three major issues described in the introduction, with details given in table  <dig>  learning was carried out using data sampled from a synthetic regulatory network, described in additional file 1: supplementary figure  <dig> and supplementary table  <dig>  where the true network structure is known.
we learned the best two-parent regulatory model  for each dependent node  by exhaustively searching through each possible parent set and scoring with each metric.
note that each control method compares to and validates mi <dig> in one of the three major aspects described in introduction: data discretization ; pairwise testing ; emphasis on correlation over causality .
all scores are calculated based on continuous nonparametric probability density estimation, except dmi <dig> based on discretization using  <dig> bins of equal size.
† in this paper, log conditional probability and bn are used interchangeably.
# performance rank for real data experiment is based on qualitative comparison.
the resulting best scoring network from a representative experiment is shown in figure  <dig>  using the mi <dig> score, we recovered the true models for all dependent variables with exactly two parents, including u <dig>  u <dig> and u <dig>  for variables with fewer or more than two parents, i.e.
u <dig>  u <dig> and u <dig>  mi <dig> detected the best two-parent representative of the true models.
continuous mi <dig> outperformed dmi <dig> as dmi <dig> identified poor models for u <dig>  u <dig>  and u <dig>  the bn tended to select confounding yet correlative models with low or negative coordination  between the two parents.
for example, the bn score selected u2+u <dig> and x3+u <dig> over x1+x <dig> as the top  <dig> models for u <dig>  therefore, the coordinative component in mi <dig> is necessary to differentiate the true parent set from the confounding set.
compared to mi <dig>  mi <dig> as well as log conditional probability consistently gave more accurate models whenever there was a difference, demonstrating their advantage in capturing higher order relationships.
the existence of two way edges or edges with reversed direction showed that mi <dig> could not identify direction of causality between variables.
in addition, the two parents for nodes u <dig>  u <dig>  u <dig> and u <dig> picked by mi <dig> have highly negative coordination with each other.
these results demonstrate that, among the methods tested, mi <dig> most accurately identified the underlying regulatory network for both linear and nonlinear relationships between variables .
next we quantitatively compared the performance of mi <dig> to other commonly used methods in terms of both sensitivity  and precision   <cit> .
in figure  <dig>  only sensitivity curves are shown because the precision curves are essentially the same but shifted.
figure 3a provides the absolute performance, while 2b shows the relative performance.
the relative performance is a more meaningful comparison, given that the number of parents was fixed, although both results are quite similar.
the absolute sensitivity and precision mi <dig> algorithm achieved were  <dig>  and  <dig>  respectively , and the relative levels are both  <dig>  .
in this comparison, mi <dig> consistently outperformed dmi <dig> across all different sample sizes.
also mi <dig> was more robust than dmi <dig> in that the sensitivity and precision curves have smaller error bars .
in addition mi <dig> always outperformed the correlative bn.
mi <dig> consistently demonstrated the lowest performance by a large margin as long as the sample size was greater than  <dig>  all methods reached a plateau at ~ <dig> samples, indicating that the  <dig>  sample default used in this paper is appropriate for all  <dig> methods to learn two parent regulatory models .
finally, all four methods were ranked in terms of performance in table  <dig>  overall, mi <dig> always gave the highest true positive and the lowest false positive rate, and significantly outperformed all control methods , details of statistical tests shown given in additional file 1: supplementary table  <dig>
while experiment with the above-mentioned small synthetic network clearly proved the principles of mi <dig> method, we also scaled up to a synthetic network that has the size of large gene regulatory networks yet still allows exhaustive search of two-parent models.
this large synthetic network with  <dig> nodes and  <dig> edges was created by tiling  <dig> copies of the small network  plus  <dig> cross-tile edges  <cit> .
experiments with data sampled from this large synthetic network yielded similar results  to those from the small synthetic study above.
mi <dig> applied to high throughput microarray data
we used mi <dig> and control methods to infer regulatory network centered at myc transcription factor from a human b cell microarray dataset.
note that the same dataset had been generated and used for identifying myc target genes by another group  <cit> .
instead of doing an exhaustive search of co-regulator pairs for each target as in the synthetic data, we fixed one of the regulators to be myc and the target to be a known myc target, and searched for the second regulator.
this constraint imposed by our specific biological focus made the analysis more tractable and our results more testable, because we only need to select and test the second regulator .
notice that this simplified problem is a sub-case of the synthetic problem.
we are still using the same scoring metrics  and following the same procedure , except that one parent node is fixed by introducing extra literature data.
in this sense, all methods are still comparable.
experiments with synthetic data showed that such simplification does not change the final results as long as we are introducing a real parent of the target with enough marginal dependency, i.e.
i >  <dig> , for mi <dig>  dmi <dig> and bn.
for mi <dig>  fixing r <dig> = myc does change the results, but it makes sense when taken as prior knowledge introduction.
we pre-filtered myc targets, t, with i ≥  <dig>  to prevent bias upon fixing r <dig> = myc, and to speed up analysis similar to candidate parent set selection in the sparse candidate algorithm  <cit> .
the verified targets were retrieved from the myc target gene database  <cit>  available online  <cit> .
after pre-filtering using the constraint i ≥  <dig> ,  <dig> myc targets remained as shown in additional file 1: supplementary table  <dig>  for each filtered target of myc we selected top  <dig> cofactor  models using mi <dig> or control methods.
because for each target gene, there are usually multiple models which score almost the same and are equally interesting biologically.
for example, several coregulated myc cofactors are involved in regulation of a target gene, any one of them can be selected as the best r <dig>  or multiple genes in a pathway/complex represent the same regulatory action equally well, all of them are sensible coregulators for a myc target.
this is slightly different from the synthetic experiment, where only there is  <dig> true or best model for each target.
nonetheless, keeping only top  <dig> model led to almost the same lists of most frequently selected cofactor  as the list based on top  <dig> models , except that the number of targets mapped to individual cofactors was smaller.
all other comparisons between mi <dig> and control methods led to the same results when top  <dig> models were used .
top  <dig> highest scoring cofactors are counted for each target.
cofactors in bold font are involved in myc dependent or general transcriptional regulation, those in italics are in the list of  <dig> verified myc targets with i ≥  <dig> .
mi <dig> and dmi <dig> selected models with significant coordination i, whereas the bn and mi <dig> selected models with high two-way dependency or i  is constant because r <dig> is fixed to myc) shown by figure  <dig>   <dig>  models inferred by all methods showed distinct patterns when plotted in three dimensional space .
these patterns suggest that two parents together explain the target expression well.
the difference is that bn and mi <dig> models showed distinct two dimensional patterns as well , while the mi <dig> and dmi <dig> models did not.
what mi <dig> and dmi <dig> captured are 3-way interactions in that neither of the two parents alone can describe the target well enough.
in contrast, the relationships bn and mi <dig> captured are essentially two-way, and as such do not require both parents.
this outcome is not surprising in that the mi <dig> metric favors strong three way interactions, while the bn and mi <dig> methods have no such favor and as such would be expected to include confounding two-way models more frequently.
bn and mi <dig> models had low or negative 3-way coordination, and are likely confounding models.
the relationship r2~r <dig> is similar to t~r <dig> and t~r <dig> follows a nearly perfect linear pattern .
such high similarity between r <dig> and t is unlikely true regulation but rather coregulation  relationship when considering various other factors affecting the target gene expression that are not counted by transcription level of the regulator, such as mrna to protein translation, protein modification, and localization of the regulator.
in other words, real regulators do not correlate so well with their targets .
we expect that the r <dig> factors predicted by the bn and mi <dig> methods is most often another myc target tightly coregulated with t instead of a coregulator, and indeed many top r <dig> are myc targets .
next we collected the top  <dig> cofactors and ranked each cofactor according to its frequency of being selected.
table  <dig> lists the top  <dig> most frequently selected cofactors using the four methods.
transcriptional regulatory networks centered at myc were constructed based on the top  <dig> cofactors and corresponding targets, as shown in figure  <dig>  literature validation was focused on these top  <dig> cofactor lists .
top  <dig> cofactor based transcriptional regulatory networks constructed by mi <dig> and dmi <dig> were larger and covered more myc targets than the networks created by bn and mi <dig> .
out of  <dig> myc targets, mi <dig> places  <dig> % of these targets while dmi <dig> places  <dig> %, bn places  <dig> %, and mi <dig> places  <dig> % of the targets.
in other words, more myc target genes are regulated by the top  <dig> mechanisms inferred by mi <dig> or dmi <dig>  which is more consistent with the current mechanistic understanding of myc dependent transcription that myc regulates a large number of targets   <cit>  as a global transcriptional regulator yet only interacts with a small set of cofactors   <cit> .
biologically, top  <dig> myc cofactor list selected by mi <dig> was more consistent with the literature than the lists created by the control methods .
seven out of ten mi <dig> top myc cofactors are involved in myc dependent or general transcriptional regulation.
gcn5l <dig> , ash2l, mizf, cbx <dig>  are chromatin structure modifiers, which change chromatin structure around target genes through chemical modification hence activate or repress their transcription.
chromatin structure modification by gcn5l <dig> and similar enzymes is a well documented mechanism for myc dependent transcriptional regulation  <cit> .
znf <dig>  <cit>  and mizf  <cit>  are transcriptional factors.
csk phosphorylates and activates gsk-3beta directly  <cit>  and indirectly  <cit> , while gsk-3beta phosphorylates, deactivates myc and promotes its degradation  <cit> .
shoc <dig> complexes with ras and raf and enhances map kinase activation  <cit> , which in turn positively regulates myc stability/activity by phosphorylation  <cit> .
in contrast, only  <dig> ,  <dig>  and  <dig>  top  <dig> myc cofactors selected by dmi <dig>  bn and mi <dig> respectively are transcriptional regulators based on gene ontology and literature.
moreover,  <dig>  and  <dig>  top  <dig> myc cofactors selected by bn and mi <dig> respectively are actually from the  <dig> verified myc targets.
other top  <dig> cofactors selected by bn and mi <dig> are likely 'unverified' myc targets given that they either share function annotations or have similar expression profile with these questionable cofactors.
in other words, bn and mi <dig> frequently produced confounding models where target genes were mistaken as myc cofactors, while mi <dig> and dmi <dig> produced no confounding models.
in figure 6d, the two-way edges between red nodes suggest that mi <dig> not only confounded coregulators with targets, but also failed to tell the causal direction of the relationships.
combined with numerical comparison in figure  <dig>   <dig>  these biological results show that unlike bn and mi <dig> scores, mi <dig> score effectively differentiates true causal models from confounding models because it takes the interaction between regulators into account.
discussion
in this study, we have used mi <dig> to identify mechanistically plausible relationships from gene expression data.
for synthetic data, mi <dig> recovered all true two-parent models, or the best representatives of the true models, and showed superior performance over the commonly used probability based methods including bayesian networks and classical two-way mutual information.
for experimental data, myc cofactors identified by mi <dig> are either true or strongly supported by literature, while cofactors identified by control methods make little sense.
notably, the same microarray dataset has been used to identify myc targets based on two-way mutual information  <cit> .
mi <dig> uses three strategies to improve its predictions.
first, mi <dig> does not require data discretization, and as such retains more of the information in the data.
this continuous method enhanced the learning quality significantly, as shown by the synthetic example in figure  <dig>   <dig>  second, we extended classical two-way mutual information to three-way, which allows mi <dig> to capture more complex relationships between regulators and targets.
third, the mi <dig> score considers high order interaction or coordination and better differentiates causal relationships from confounding relationships as was shown by both the synthetic and myc problems .
myc cofactors predicted by mi <dig> details agree with the established literature.
notably, four of the top  <dig> cofactors selected by mi <dig> are chromatin structure modifier genes, suggesting that chromatin structure modification is the primary mechanism for myc dependent transcriptional regulation.
this inference is directly supported by the independent experimental results of knoepfler et al , which provides further evidence of the role of myc on chromatin structure modification via histone acetylation and methylation.
among the top myc cofactors identified by mi <dig>  gcn5l <dig>  <cit> , csk  <cit> , and shoc <dig>  <cit>  are known or presumed coregulators for myc transcriptional activity.
all other seven myc cofactors selected by mi <dig> are novel, although their connections to myc or transcription are well documented.
all these results demonstrate that mi <dig> is an accurate and powerful method to infer regulatory models from microarray data.
in contrast, top myc cofactors inferred using control methods make much less sense biologically.
fewer of them are known transcriptional regulators and none of them is directly connected to myc function.
the fact that multiple myc targets were mistaken as top myc coregulators suggests that bn and mi <dig> methods have difficulty inferring true causal relationships from high throughput gene expression data.
generally speaking, it is sensible that some myc targets can be its cofactors as seen in feedback loops.
however, it is not likely that these myc targets taken as co-regulators are real co-regulators because of feedback loops, since almost all of them are not functionally related to transcriptional regulation or myc regulation activity.
similar confounding regulators were selected by control methods in the synthetic example .
figure  <dig>   <dig> show why such confounding models occurred.
there are likely feedback loops in myc regulation, however these feedback relationships could only be identified with knockout data or time series data.
in this work we only consider the general case where non-sequential observational gene expression data are available.
learning from high throughput microarray data was different from learning from the small synthetic dataset.
differences between methods were larger for the microarray data , compared to the synthetic experiment .
for the microarray data, mi <dig> and dmi <dig> were closer, whereas for the synthetic data bn and mi <dig> were closer .
this change in ranking suggests that the coordinative component was more significant than the difference made by using continuous versus discrete metric  or 3-way versus 2-way metric  for microarray data, but not for synthetic data.
these differences between microarray data and synthetic data can be ascribed to the fact that large numbers of highly correlative confounding models exist for the microarray data due to the large number of variables , especially coexpressed genes, while the synthetic data contained relatively fewer possible confounding models.
the high order mutual information framework presented here is generally applicable, although we have only described and used three-way mutual information.
the same set of strategies can be used to model arbitrarily high order relationships.
to learn a regulatory model with d dimensions or nodes  by exhaustive searching through a system with v variables, we need ~10*5d data samples for nonparametric probability density estimation  <cit> , and computation time is o. although 10*5d is conservative compared to sufficient sample size indicated in the performance curve, ~ <dig> for d =  <dig> , undoubtedly, both the required dataset size and computational time exponentially increase with d. therefore, 4-way or 5-way relationships require more samples than currently available microarray chips.
through the use of mi <dig> we have demonstrated that tailored probability based metrics can outperform more standard methods used in systems biology for identifying mechanistic regulatory relationships.
we expect that future enhancements to these scoring metrics are possible to identify larger sets of regulators while making fewer assumptions during the analysis.
CONCLUSIONS
mi <dig> is a novel method for learning probabilistic graphical models and addresses three common issues in previous methods simultaneously:  handling of continuous variables,  detection of more complex three-way relationships and  better differentiation of causal versus confounding relationships.
simulation experiments show that mi <dig> consistently and significantly outperforms the frequently used control methods such as discrete probabilistic networks, two-way mutual information networks and bayesian networks.
when applied to a human b cell microarray dataset, mi <dig> recovered cofactors for myc transcription factor and revealed the major regulatory mechanisms involved in myc dependent transcriptional regulation, which are directly verified or strongly supported by literature.
overall, mi <dig> provides a powerful method for inferring mechanistic relationships underlying biological and other complex systems.

 a major challenge in systems biology is to infer mechanistic gene interactions from high throughput microarray data. underlying this challenge is the problem to find causal or regulatory relationships between genes. robust solutions to this problem would provide us with a transcriptomic map of a genome that allows us to accurately predict the effect of gene perturbations. previous efforts to detect mechanistic relationships from gene expression data can be broadly divided into linear correlation and probability based methods. linear correlation based methods, such as clustering, correlation networks and graphical gaussian models, have a long and fruitful history in statistical modeling and bioinformatics. these linear methods are computationally fast and relatively easy to interpret. however, a key limitation with these methods is that they assume linear relationships between variables. while some components of any transcriptional regulatory network are linear, nonlinear events such as or, and, and xor type transcriptional regulation are relatively commonplace. these nonlinear interactions would not be captured with a linear model, leading to spurious relationships between variables. probability based methods have also been used to detect relationships between genes. these probability methods include probabilistic boolean networks, bayesian networks and mutual information networks. probability based methods can capture both linear and nonlinear regulatory relationships and are noise tolerant. however, many of the current probability based tools used in systems biology suffer from the following three limitations: data discretization, pairwise testing, emphasis on correlation over causality. to transform continuous data into a more easily computable form, most probabilistic methods require the data to first be discretized into a finite number of bins, such as high, medium, and low: the number of bins used in discretization is difficult to choose, and is generally selected at some consistent yet arbitrary point. unfortunately, different binning procedures can produce different analysis results, suggesting that the act of binning alone introduces errors into the analysis. methods that search for pairwise associations only focus on a single relationship between regulator and target at a time. pairwise association networks have been created using classical mutual information. however, simple pairwise relationships are likely less common than multivariate relationships in real biological systems, as the expression of most genes is regulated not by a single gene but more likely by multiple genes. methods that allow multivariate interactions such as bayesian networks or some fuzzy logic approaches are inherently superior in this respect. a final challenge in creating mechanistically predictive transcriptional models is the ability to identify not just correlative but also causal models. for the definition of causal relationship, we adopt the notion of probabilistic causation. informally, event a causes event b if and only if a's occurrence alters the probability of b. this sometimes reflects imperfect knowledge of a deterministic system but more frequently suggests a stochastic nature of the causal system under study. although difficult, causal relationships have been learned properly from non-sequential observational data. probabilistic graphical modeling methods like bayesian networks have been used to infer causal models from gene expression data. however, many probabilistic approaches are able to make correlative networks but not necessarily causal networks. their multivariate scoring metrics such as conditional probability and mutual information are still generalized two-way correlation between the target and the parent set. similar to the classical two-way metrics, these generalized correlations alone cannot differentiate between a causal versus confounding parent set. true causal relationships like genetic regulation feature positive higher order interaction, the non-additive effect above the sum of the lower order interactions. for instance, for regulation involved two regulators such as or, and, xor type relationships, two regulators together account for much more in the target than they individually can. intuitively such non-additive effect can be described as coordination or synergy between parents. on the other hand, confounding models commonly have no or negative higher order interaction. we propose that with such high order interaction considered, we can better differentiate true causal model versus confounding models. in this work, we demonstrate a novel algorithm that attempts to overcome all three limitations using a continuous high order mutual information based scoring metric we call mi. note that continuous two-way mutual information has been described previously. high order interaction information has been employed to model complex interactions. however, both two-way mutual information and high order interaction information are symmetric and as such unable to make causal statements. mi combines 3rd order interaction information with the asymmetric mutual information between target and regulator set to account for the direction of regulation. mi is novel as a combinatorial probabilistic metric and an integrated statistical learning method. in this work, we compare mi to other probability based methods quantitatively and qualitatively using synthetic data where the true model is known. next we apply mi and control methods to reconstruct regulatory networks centered at the transcription factor myc from a published high throughput microarray dataset. the learning results are then evaluated numerically and biologically. learning myc centered transcriptional regulatory network represents an ideal test case for mi as myc is a well characterized transcriptional regulator that acts in tandem with a finite set of co-effectors and regulates the expression of a large group of genes. myc has been well investigated and online databases of myc targets are available for validation purpose. despite these efforts, many cofactors and targets remain unidentified, and corresponding regulatory mechanisms unknown. as a result, an integrated understanding of myc dependent transcriptional regulation has remained out of reach. in this study, we use mi to derive an accurate transcriptomic map surrounding myc from the same gene expression dataset used to identify myc targets. the approaches used here are general and can be directly used for any transcriptional regulator given sufficient gene expression data.  mi validation with synthetic data we validated mi against other commonly used methods listed in table including a discrete version of mi, two-way mutual information and a log conditional probability score used in bayesian network learning. each control method represents one of the three major issues described in the introduction, with details given in table learning was carried out using data sampled from a synthetic regulatory network, described in additional file 1: supplementary figure and supplementary table where the true network structure is known. we learned the best two-parent regulatory model for each dependent node by exhaustively searching through each possible parent set and scoring with each metric. note that each control method compares to and validates mi in one of the three major aspects described in introduction: data discretization; pairwise testing; emphasis on correlation over causality. all scores are calculated based on continuous nonparametric probability density estimation, except dmi based on discretization using bins of equal size. † in this paper, log conditional probability and bn are used interchangeably. # performance rank for real data experiment is based on qualitative comparison. the resulting best scoring network from a representative experiment is shown in figure using the mi score, we recovered the true models for all dependent variables with exactly two parents, including u u and u for variables with fewer or more than two parents, i.e. u u and u mi detected the best two-parent representative of the true models. continuous mi outperformed dmi as dmi identified poor models for u u and u the bn tended to select confounding yet correlative models with low or negative coordination between the two parents. for example, the bn score selected u2+u and x3+u over x1+x as the top models for u therefore, the coordinative component in mi is necessary to differentiate the true parent set from the confounding set. compared to mi mi as well as log conditional probability consistently gave more accurate models whenever there was a difference, demonstrating their advantage in capturing higher order relationships. the existence of two way edges or edges with reversed direction showed that mi could not identify direction of causality between variables. in addition, the two parents for nodes u u u and u picked by mi have highly negative coordination with each other. these results demonstrate that, among the methods tested, mi most accurately identified the underlying regulatory network for both linear and nonlinear relationships between variables. next we quantitatively compared the performance of mi to other commonly used methods in terms of both sensitivity and precision. in figure only sensitivity curves are shown because the precision curves are essentially the same but shifted. figure 3a provides the absolute performance, while 2b shows the relative performance. the relative performance is a more meaningful comparison, given that the number of parents was fixed, although both results are quite similar. the absolute sensitivity and precision mi algorithm achieved were and respectively, and the relative levels are both. in this comparison, mi consistently outperformed dmi across all different sample sizes. also mi was more robust than dmi in that the sensitivity and precision curves have smaller error bars. in addition mi always outperformed the correlative bn. mi consistently demonstrated the lowest performance by a large margin as long as the sample size was greater than all methods reached a plateau at ~ samples, indicating that the sample default used in this paper is appropriate for all methods to learn two parent regulatory models. finally, all four methods were ranked in terms of performance in table overall, mi always gave the highest true positive and the lowest false positive rate, and significantly outperformed all control methods, details of statistical tests shown given in additional file 1: supplementary table while experiment with the above-mentioned small synthetic network clearly proved the principles of mi method, we also scaled up to a synthetic network that has the size of large gene regulatory networks yet still allows exhaustive search of two-parent models. this large synthetic network with nodes and edges was created by tiling copies of the small network plus cross-tile edges. experiments with data sampled from this large synthetic network yielded similar results to those from the small synthetic study above. mi applied to high throughput microarray data we used mi and control methods to infer regulatory network centered at myc transcription factor from a human b cell microarray dataset. note that the same dataset had been generated and used for identifying myc target genes by another group. instead of doing an exhaustive search of co-regulator pairs for each target as in the synthetic data, we fixed one of the regulators to be myc and the target to be a known myc target, and searched for the second regulator. this constraint imposed by our specific biological focus made the analysis more tractable and our results more testable, because we only need to select and test the second regulator. notice that this simplified problem is a sub-case of the synthetic problem. we are still using the same scoring metrics and following the same procedure, except that one parent node is fixed by introducing extra literature data. in this sense, all methods are still comparable. experiments with synthetic data showed that such simplification does not change the final results as long as we are introducing a real parent of the target with enough marginal dependency, i.e. i >, for mi dmi and bn. for mi fixing r = myc does change the results, but it makes sense when taken as prior knowledge introduction. we pre-filtered myc targets, t, with i ≥ to prevent bias upon fixing r = myc, and to speed up analysis similar to candidate parent set selection in the sparse candidate algorithm. the verified targets were retrieved from the myc target gene database available online. after pre-filtering using the constraint i ≥, myc targets remained as shown in additional file 1: supplementary table for each filtered target of myc we selected top cofactor models using mi or control methods. because for each target gene, there are usually multiple models which score almost the same and are equally interesting biologically. for example, several coregulated myc cofactors are involved in regulation of a target gene, any one of them can be selected as the best r or multiple genes in a pathway/complex represent the same regulatory action equally well, all of them are sensible coregulators for a myc target. this is slightly different from the synthetic experiment, where only there is true or best model for each target. nonetheless, keeping only top model led to almost the same lists of most frequently selected cofactor as the list based on top models, except that the number of targets mapped to individual cofactors was smaller. all other comparisons between mi and control methods led to the same results when top models were used. top highest scoring cofactors are counted for each target. cofactors in bold font are involved in myc dependent or general transcriptional regulation, those in italics are in the list of verified myc targets with i ≥. mi and dmi selected models with significant coordination i, whereas the bn and mi selected models with high two-way dependency or i is constant because r is fixed to myc) shown by figure models inferred by all methods showed distinct patterns when plotted in three dimensional space. these patterns suggest that two parents together explain the target expression well. the difference is that bn and mi models showed distinct two dimensional patterns as well, while the mi and dmi models did not. what mi and dmi captured are 3-way interactions in that neither of the two parents alone can describe the target well enough. in contrast, the relationships bn and mi captured are essentially two-way, and as such do not require both parents. this outcome is not surprising in that the mi metric favors strong three way interactions, while the bn and mi methods have no such favor and as such would be expected to include confounding two-way models more frequently. bn and mi models had low or negative 3-way coordination, and are likely confounding models. the relationship r2~r is similar to t~r and t~r follows a nearly perfect linear pattern. such high similarity between r and t is unlikely true regulation but rather coregulation relationship when considering various other factors affecting the target gene expression that are not counted by transcription level of the regulator, such as mrna to protein translation, protein modification, and localization of the regulator. in other words, real regulators do not correlate so well with their targets. we expect that the r factors predicted by the bn and mi methods is most often another myc target tightly coregulated with t instead of a coregulator, and indeed many top r are myc targets. next we collected the top cofactors and ranked each cofactor according to its frequency of being selected. table lists the top most frequently selected cofactors using the four methods. transcriptional regulatory networks centered at myc were constructed based on the top cofactors and corresponding targets, as shown in figure literature validation was focused on these top cofactor lists. top cofactor based transcriptional regulatory networks constructed by mi and dmi were larger and covered more myc targets than the networks created by bn and mi. out of myc targets, mi places % of these targets while dmi places %, bn places %, and mi places % of the targets. in other words, more myc target genes are regulated by the top mechanisms inferred by mi or dmi which is more consistent with the current mechanistic understanding of myc dependent transcription that myc regulates a large number of targets as a global transcriptional regulator yet only interacts with a small set of cofactors. biologically, top myc cofactor list selected by mi was more consistent with the literature than the lists created by the control methods. seven out of ten mi top myc cofactors are involved in myc dependent or general transcriptional regulation. gcn5l, ash2l, mizf, cbx are chromatin structure modifiers, which change chromatin structure around target genes through chemical modification hence activate or repress their transcription. chromatin structure modification by gcn5l and similar enzymes is a well documented mechanism for myc dependent transcriptional regulation. znf and mizf are transcriptional factors. csk phosphorylates and activates gsk-3beta directly and indirectly, while gsk-3beta phosphorylates, deactivates myc and promotes its degradation. shoc complexes with ras and raf and enhances map kinase activation, which in turn positively regulates myc stability/activity by phosphorylation. in contrast, only, and top myc cofactors selected by dmi bn and mi respectively are transcriptional regulators based on gene ontology and literature. moreover, and top myc cofactors selected by bn and mi respectively are actually from the verified myc targets. other top cofactors selected by bn and mi are likely'unverified' myc targets given that they either share function annotations or have similar expression profile with these questionable cofactors. in other words, bn and mi frequently produced confounding models where target genes were mistaken as myc cofactors, while mi and dmi produced no confounding models. in figure 6d, the two-way edges between red nodes suggest that mi not only confounded coregulators with targets, but also failed to tell the causal direction of the relationships. combined with numerical comparison in figure these biological results show that unlike bn and mi scores, mi score effectively differentiates true causal models from confounding models because it takes the interaction between regulators into account. discussion in this study, we have used mi to identify mechanistically plausible relationships from gene expression data. for synthetic data, mi recovered all true two-parent models, or the best representatives of the true models, and showed superior performance over the commonly used probability based methods including bayesian networks and classical two-way mutual information. for experimental data, myc cofactors identified by mi are either true or strongly supported by literature, while cofactors identified by control methods make little sense. notably, the same microarray dataset has been used to identify myc targets based on two-way mutual information. mi uses three strategies to improve its predictions. first, mi does not require data discretization, and as such retains more of the information in the data. this continuous method enhanced the learning quality significantly, as shown by the synthetic example in figure second, we extended classical two-way mutual information to three-way, which allows mi to capture more complex relationships between regulators and targets. third, the mi score considers high order interaction or coordination and better differentiates causal relationships from confounding relationships as was shown by both the synthetic and myc problems. myc cofactors predicted by mi details agree with the established literature. notably, four of the top cofactors selected by mi are chromatin structure modifier genes, suggesting that chromatin structure modification is the primary mechanism for myc dependent transcriptional regulation. this inference is directly supported by the independent experimental results of knoepfler et al, which provides further evidence of the role of myc on chromatin structure modification via histone acetylation and methylation. among the top myc cofactors identified by mi gcn5l, csk, and shoc are known or presumed coregulators for myc transcriptional activity. all other seven myc cofactors selected by mi are novel, although their connections to myc or transcription are well documented. all these results demonstrate that mi is an accurate and powerful method to infer regulatory models from microarray data. in contrast, top myc cofactors inferred using control methods make much less sense biologically. fewer of them are known transcriptional regulators and none of them is directly connected to myc function. the fact that multiple myc targets were mistaken as top myc coregulators suggests that bn and mi methods have difficulty inferring true causal relationships from high throughput gene expression data. generally speaking, it is sensible that some myc targets can be its cofactors as seen in feedback loops. however, it is not likely that these myc targets taken as co-regulators are real co-regulators because of feedback loops, since almost all of them are not functionally related to transcriptional regulation or myc regulation activity. similar confounding regulators were selected by control methods in the synthetic example. figure show why such confounding models occurred. there are likely feedback loops in myc regulation, however these feedback relationships could only be identified with knockout data or time series data. in this work we only consider the general case where non-sequential observational gene expression data are available. learning from high throughput microarray data was different from learning from the small synthetic dataset. differences between methods were larger for the microarray data, compared to the synthetic experiment. for the microarray data, mi and dmi were closer, whereas for the synthetic data bn and mi were closer. this change in ranking suggests that the coordinative component was more significant than the difference made by using continuous versus discrete metric or 3-way versus 2-way metric for microarray data, but not for synthetic data. these differences between microarray data and synthetic data can be ascribed to the fact that large numbers of highly correlative confounding models exist for the microarray data due to the large number of variables, especially coexpressed genes, while the synthetic data contained relatively fewer possible confounding models. the high order mutual information framework presented here is generally applicable, although we have only described and used three-way mutual information. the same set of strategies can be used to model arbitrarily high order relationships. to learn a regulatory model with d dimensions or nodes by exhaustive searching through a system with v variables, we need ~10*5d data samples for nonparametric probability density estimation, and computation time is o. although 10*5d is conservative compared to sufficient sample size indicated in the performance curve, ~ for d =, undoubtedly, both the required dataset size and computational time exponentially increase with d. therefore, 4-way or 5-way relationships require more samples than currently available microarray chips. through the use of mi we have demonstrated that tailored probability based metrics can outperform more standard methods used in systems biology for identifying mechanistic regulatory relationships. we expect that future enhancements to these scoring metrics are possible to identify larger sets of regulators while making fewer assumptions during the analysis.  mi is a novel method for learning probabilistic graphical models and addresses three common issues in previous methods simultaneously: handling of continuous variables, detection of more complex three-way relationships and better differentiation of causal versus confounding relationships. simulation experiments show that mi consistently and significantly outperforms the frequently used control methods such as discrete probabilistic networks, two-way mutual information networks and bayesian networks. when applied to a human b cell microarray dataset, mi recovered cofactors for myc transcription factor and revealed the major regulatory mechanisms involved in myc dependent transcriptional regulation, which are directly verified or strongly supported by literature. overall, mi provides a powerful method for inferring mechanistic relationships underlying biological and other complex systems. 
in this work we introduce a new statistical learning strategy, mi that addresses three common issues in previous methods simultaneously: handling of continuous variables, detection of more complex three-way relationships and better differentiation of causal versus confounding relationships. unlike control methods, mi effectively differentiated true causal models from confounding models. in addition, mi significantly outperformed the control methods, including bayesian networks, classical two-way mutual information and a discrete version of mi we then used mi and control methods to infer a regulatory network centered at the myc transcription factor from a published microarray dataset. the mi network showed that limited sets of regulatory mechanisms are employed repeatedly to control the expression of large number of genes. mi recovered major myc cofactors, and revealed major mechanisms involved in myc dependent transcriptional regulation, which are strongly supported by literature. overall, our work demonstrates that mi outperforms the frequently used control methods, and provides a powerful method for inferring mechanistic relationships underlying biological and other complex systems. probability based statistical learning methods such as mutual information and bayesian networks have emerged as a major category of tools for reverse engineering mechanistic relationships from quantitative biological data. we test the mi algorithm using both synthetic and experimental data. models selected by mi were numerically and biologically distinct from those selected by control methods.
 probability based statistical learning methods such as mutual information and bayesian networks have emerged as a major category of tools for reverse engineering mechanistic relationships from quantitative biological data. in this work we introduce a new statistical learning strategy, mi that addresses three common issues in previous methods simultaneously: handling of continuous variables, detection of more complex three-way relationships and better differentiation of causal versus confounding relationships. with these improvements, we provide a more realistic representation of the underlying biological system.  we test the mi algorithm using both synthetic and experimental data. in the synthetic data experiment, mi achieved an absolute sensitivity/precision of / and a relative sensitivity/precision both of. in addition, mi significantly outperformed the control methods, including bayesian networks, classical two-way mutual information and a discrete version of mi we then used mi and control methods to infer a regulatory network centered at the myc transcription factor from a published microarray dataset. models selected by mi were numerically and biologically distinct from those selected by control methods. unlike control methods, mi effectively differentiated true causal models from confounding models. mi recovered major myc cofactors, and revealed major mechanisms involved in myc dependent transcriptional regulation, which are strongly supported by literature. the mi network showed that limited sets of regulatory mechanisms are employed repeatedly to control the expression of large number of genes.  overall, our work demonstrates that mi outperforms the frequently used control methods, and provides a powerful method for inferring mechanistic relationships underlying biological and other complex systems. the mi method is implemented in r in the "mi3" package, available under the gnu gpl from and from the r package archive cran. 