BACKGROUND
a general problem in systems biology is the integration of observational experimental data such as gene expression, with known pathways topologies.
ideally these two data sources should be complementary, however in practice there are few methods to systematically integrate these two kinds of information.
in this paper, we introduce a method called mechanistic bayesian networks  in an attempt to use knowledge about the topology of a pathway with gene expression data related to the same pathway.
as a sample case we used the topology of the sonic hedgehog  signaling pathway as a model pathway along with a targeted gene expression dataset.
using these data we used the mbn approach to identify regulatory targets of the shh pathway.
sonic hedgehog  pathway
the shh pathway plays a central role in organismal development and the progression of some cancers  <cit> .
because of its central role, the shh pathway is well studied providing us with an ideal test case to validate our mbn approach.
the details of shh are reviewed in detail elsewhere  <cit> , but here we will summarize the early steps of the pathway that we will use in this work.
shh is a secreted protein that acts as both a short-range contact-dependant factor and as a long-range diffusible morphogen.
the shh ligand binds its canonical receptor patched <dig> , which releases its inhibition of a second membrane-bound protein, smoothened .
derepression of smo in turn activates a signaling cascade inside the cell, eventually activating the gli transcription factors and regulating the expression of a variety of genes.
while some of the downstream targets of the pathway are known, many remain unknown.
the initial steps of the pathway are shown in figure  <dig> with the terminal node genex representing a putative downstream target.
note that in this cascade, the shh, ptch <dig>  and smo proteins can all directly or indirectly affect the target, although in different ways and to different degrees.
figure  <dig> and  <dig> show more abstracted representations of the pathway.
due to its dose-dependant effect and role in development, the shh signal requires strict spatiotemporal control.
the pathway's known targets include pathway components, promoters and antagonists, thus regulating the effects of the shh ligand over time.
cell surface proteins that promote shh signaling, for example, are initially expressed in shh-responsive cells, sensitizing cells to even low levels of the shh ligand.
as the level of shh signaling increases, downregulation of positive shh components such as gas <dig> and upregulation of negative components such as ptch <dig>  ptch <dig>  and hhip that sequester the ligand ensure a tight control over shh signaling  <cit> .
additionally, shh is known to cooperate with or antagonize other pathways such as bmp  <cit> , retinoic acid  <cit> , wnt  <cit> , ras  <cit> , and notch  <cit>  in a time-, dose-, and spatially-dependant manner.
these factors further complicate the identification of pathway targets.
the shh pathway raises a more general problem: given gene expression data from multiple samples, tissue and organs under different genetic knockout conditions, how can we identify downstream targets of a given pathway?
due to the interactions between different pathways and multiple cascades between ligand reception and eventual transcriptional regulation, it is not clear what constitutes a downstream target of a specific pathway.
common analysis approaches include significance testing between samples, differential expression, and clustering  <cit> .
while these approaches are widely used and helpful, they fail to incorporate knowledge of the underlying pathway.
when the pathway information is used, it is done while ignoring fundamental biological knowledge such as the central dogma.
dna, mrna, and proteins are conflated into one variable and specific interactions such as protein-protein or protein-dna are implicitly assumed to be detectable via gene expression data.
to circumvent this problem, we present a novel mechanistic bayesian network approach that more closely respects the meaning of both the pathway and the expression data.
methods
in the following sections we describe the theoretical underpinnings of the mbn framework and provide a method for evaluating the significance of a result.
next we describe how we tested the mbn method using gene expression data from gene knockout mouse models to identify targets of the shh pathway.
mechanistic bayesian networks: theory and definition
mechanistic bayesian networks represent a subset of a general class of analysis tools called bayesian networks .
a bn is a probabilistic graphical model that encodes dependencies between variables in a compact and descriptive manner.
a bn can be represented as a graph with nodes representing variables and edges representing dependencies or relationships between the variables.
mathematically, each node describes a conditional probability distribution  that quantitatively models the relationships between a node and its parent nodes.
note that an edge  between two variables indicates a directed relationship but does not specify the functional form of the relationship.
said another way, an edge in a bayesian network does not indicate activation, inhibition or any other specific function.
this interpretation of an edge differs from the usual definition of an edge in a signaling pathway, where edges are often assumed to have a defined effect.
the advantage of the more broad edge definition is that the relationships between nodes in a bayesian network can be more complex and include functions that are nonlinear, multimodal, or logical.
bns have been used successfully to model complex phenomenon in many fields and in systems biology, in particular, to model gene-regulatory networks, protein-interaction networks, signaling networks and to integrate heterogeneous biological data  <cit> .
methods for training and learning bns are well established  <cit>  and are available in a large number of software applications  <cit> .
in contrast to the more general bayesian network, a mechanistic bayesian network  adheres more closely to known biological mechanisms by differentiating between mrna transcripts and proteins, and by including pathway-based structural constraints.
because most studies do not measure both protein and gene expression, unknown quantities in an mbn are treated as unobserved, latent variables.
by creating models that more closely resemble the known mechanism, mbns effectively incorporate in additional information that is not available in the experimental data alone.
mechanistic bayesian network templates
when analyzing a biological system, users often want to generate and rank a set of entities that match certain constraints, such as "all downstream targets of the shh pathway" or "all proteins that participate in the crosstalk between the shh and wnt pathways".
while it is possible to devise specific sets of constraints and methods for each case, we propose a generic mbn based template approach.
each set of constraints is expressed as an mbn template from which specific mbns are instantiated and evaluated against an experimental dataset.
an mbn template is composed of two types of nodes: constant and template nodes.
a constant node is a regular node corresponding to either one variable in the dataset or an unobserved variable while a template node is a holding place for variables in the dataset.
each instantiation of the template is a regular mbn that can be scored using existing bn methods.
maximum entropy discretization of the data
before mbn can be applied to a dataset, the data must be discretized into a finite set of bins.
while there are many ways to bin data, we use the maximum entropy principle to derive our discretization scheme.
the principle states that the distribution that maximizes the information entropy is the true distribution given testable information  <cit> .
accordingly, we bin our data such that each bin contains the same number of data points as this maximizes the entropy of the distribution.
there is no theory prescribing the optimal number of bins for any given dataset or analytical method; increasing the number of bin decreases the information loss incurred during discretization but also increases the number of parameters and thus decreases the statistical power of the analysis given the same data samples.
most modelers working with systems biology data have arrived at three bins as a suitable compromise between information loss and statistical power  <cit> .
calculating the posterior probability of a model
once the data are discretized, a particular model topology can be scored as the posterior probability of a model given data, given below.
p is the posterior probability of the bn or mbn model m given the data d. said another way, p is our belief that the model m is correct after having observed the data d. p is the prior probability of the model, that is, our belief that the model m is correct before having observed the data d. the prior probability allows us to integrate model probabilities computed using other methods and data.
p is the likelihood that the data was generated from the model m while p is a scaling factor that is usually ignored.
p is calculated as the likelihood after marginalizing over all model parameters given a network structure.
in the mbn framework, a bayesian network is modeled using a multinomial representation with dirichlet priors for the relationships between variables.
this representation is convenient in that it is relatively agnostic to functional forms and has a convenient closed form solution, the bayesian dirichlet equivalent  metric, described below.
bde scoring metric: an approximation for the marginal likelihood
the mbn approach uses the bde scoring metric  <cit>  to calculate the marginal likelihood of a dataset given a model.
the bde metric is a closed form solution to the marginal likelihood for a multinomial bayesian network model with dirichlet priors.
the bde metric is expressed as:
where Γ is the gamma function, n is the number of nodes, ri is the arity of node i, qi is the arity of πi , nij is the number of samples where πi is in configuration j, nijk is the number of sampler where πi is in configuration j and the node has value k and αij and αijk are the prior counts corresponding to nij and nijk respectively.
a full derivation of the bde metric is described elsewhere  <cit> .
in practice, the bde metric is represented in log-space for computational convenience.
because of this log transformation, all bde scores in this paper will be negative.
calculating the marginal likelihood with missing data
if we have some missing values or latent variables in the dataset, we can use an altered method for computing the marginal likelihood.
because the likelihood is no longer fully factorizable into the product of the probabilities for each variable, we must calculate the marginal likelihood given every potential completion of the missing data.
the simplest method is to marginalize over all possible sets of values for the missing data and take the average.
however, this exact enumeration approach is impractical for non-trivial cases and so a heuristic must be used.
in this work we will used a modified gibbs sampling  <cit>  approach to approximate the marginal likelihood.
we alter the method to only sample missing data completions that result in a maximum entropy discretization for all variables.
this maximum entropy requirement ensures that p = p for all sets of variables and eliminates bias due to differential discretization for observed and hidden variables.
calculating the marginal likelihood with interventions
when calculating the marginal likelihood, data from interventions such as genetic knockouts are handled differently.
if a particular data sample is the result of an intervention on a set of variables, the values for those variables no longer depend on their parent sets.
they are, instead, arbitrarily set as the result of the intervention.
accordingly, the specific variable samples are ignored when constructing the multinomial table for any node that was intervened upon  <cit> .
assessing model significance with bootstrapping based p-values
to evaluate the significance of an mbn prediction, we calculate a p-value for each mbn result using nonparametric bootstrapping  <cit> .
in this approach, we generate a large number of mbn models with template nodes replaced by randomly generated variables with the same discretization schemes as the observed variables.
the scores for these mbn models are used as the null model against which mbn scores are compared to determine p-values.
due to the computational time required to score each bootstrap sample, we can only resolve p-values down to .
<dig> and cannot correct for multiple hypothesis testing.
we acknowledge that correction for multiple testing would lower the significance of each result.
using mbn templates to define and identify downstream targets of the shh pathway
to determine the downstream targets of the shh pathway, we created an mbn template based on the canonical pathway as described in literature with a terminal downstream template node as shown in figure  <dig>  specific mbn models are instantiated from the template by replacing the template node  with a specific gene from the dataset.
thus, each instantiated mbn model is a separate hypothesis that can be evaluated against the experimental data.
experimental data
to test our approach, we assembled data from three mouse models with gene knockouts in shh, ptch <dig>  and smo described in more detail elsewhere  <cit> .
briefly, samples from different embryonic tissues at varying developmental stages were assayed using the u74v <dig> affymetrix microarrays to determine the expression profile.
not all combinations of developmental stages and genetic backgrounds were available due to prenatal lethality for ptch <dig> and smo.
samples assayed include 6-8-somite-stage  with wildtype, smo-/-, and ptch1-/- backgrounds; 10- <dig> somite-stage  with wildtype, smo-/-, and ptch1-/- backgrounds;  <dig>  days post fertilization embryo samples from head, trunk, and limb bud with wildtype and shh-/- backgrounds.
the experimental design is summarized in table  <dig>
embryonic day  <dig>
embryonic day  <dig>
embryonic day  <dig>
embryonic day  <dig>
embryonic day  <dig>
experimental design for gathering gene expression data from a range of developmental stages, locations, and genetic backgrounds.
columns indicate the genetic background of the mouse model, rows indicate the developmental stage and location of the samples and the numbers within the table cell indicate the number of samples assayed.
data preprocessing
the raw gene expression data  were processed using rma  <cit>  as implemented in the bioconductor affy package  <cit> .
the data were annotated using updated probeset definitions  provided by the microarray lab at the university of michigan  <cit> .
data were discretized into  <dig> bins using the maximum entropy discretization scheme described above.
because the data contain samples from varying developmental stages and genes with significant expression variation between stages, genes were discretized separately for each tissue type.
this discretization scheme allowed the software to more easily identify relevant patterns within a tissue type.
genes were knocked out by altering the sequence to create inactive protein.
thus, for a knockout on a particular gene, the gene's expression was left unaltered but the protein expression was set to  <dig> and the value was marked as an intervention so it could be handled differently by the scoring procedure.
mbn template analysis
to test the mbn template system using the shh dataset, we developed a template that matched the topology shown in figure  <dig>  in this template, we search for a downstream target node, genex, that responds to the shh signaling cascade.
in this search, we assume that each gene has an equal prior probability of being the target node, and all relationships are scored using the bde score and a gibbs sampler, described above.
using this template, we generated  <dig> candidate mbn models -- one for each candidate target gene that showed significant differential expression.
due to the custom cdf annotation, smo was not present on the chip, thus we modeled the smo mrna expression as a hidden node.
no protein expression was observed in this experiment, so nodes representing protein concentrations or activities were also modeled as hidden.
rather than selecting the number of iterations for the gibbs sampler a priori, we calculated the posterior distribution p with  <dig> k iterations, then with  <dig> k iterations.
although the sampler had largely converged at this point as evidenced by the observation that the top  <dig> results did not change significantly, we resumed the gibbs's sampler for the  <dig> k run and sampled another  <dig> k iterations for a total of  <dig> k iterations per mbn model.
the computation took approximately  <dig>  minutes per mbn for a total of approximately  <dig> hours on a 48-node compute grid.
p-values were calculated with bootstrapping describe above using  <dig> mbn models to represent the null distribution.
all code and data required to replicate the analysis are included in additional file  <dig> .
the code can be easily modified to run similar analysis on different dataset using different mbn templates.
in addition, we also tested two alternative bayesian approaches shown in figure  <dig> and  <dig>  these approaches use topologies that are simpler than the mbn approach in  <dig>  but less closely adhere to the known mechanism of the biochemistry of the pathway.
RESULTS
in the following sections we show the results from the mbn analysis of the shh knockout dataset and discuss the biological plausibility for each finding.
next we discuss expected hedgehog pathway target genes that were not identified by the mbn and a comparison with other techniques.
the predicted target genes from the mbn analysis are shown in table  <dig>  the top scoring hits include many known shh targets as described below.
the top hit, gas <dig>  has been shown to be a negative target of shh signaling whereas knockout and subsequent gain-of-function studies have shown gas <dig> to be a positive component of the shh pathway that acts synergistically with ptch <dig> to bind the hedgehog ligand  <cit> .
the second target, gli <dig>  is one of the three gli transcription factors that regulate shh target genes and is often used as a canonical readout of shh activity.
shh signaling is known to both induce gli <dig> expression and also regulate its nuclear accumulation and therefore its activity  <cit> .
ptch <dig> is a homolog of ptch <dig> but has a slightly different expression pattern.
like ptch <dig>  ptch <dig> is known to be transcriptionally modulated by shh signaling  <cit> .
knockout studies have shown that shh signaling is required for expression of msx <dig>  a transcriptional repressor with a putative role in limb-pattern formation  <cit> .
foxc <dig> and foxd <dig> have been shown to be upregulated by shh signaling in null, conditional, and constitutively active smo mutant backgrounds  <cit> .
the columns specify the gene name and a short description; the bayesian score as described in methods; the resulting p-value for the bayesian score calculated using a bootstrap method as described in methods; a linear correlation to the shh transcript; a linear correlation to the ptch <dig> transcript; and a literature citation, if available, that specifies that the gene is indeed a true downstream target of the shh pathway.
beyond these known shh-related genes, some of the highly ranked results appear to be novel putative shh targets.
mig <dig> binds the opitz syndrome gene mid <dig> and the complex is thought to stabilize microtubules  <cit> .
mid <dig> is repressed by shh  <cit>  and this is partially borne out by our gene expression results which show a positive influence from shh and a negative influence from ptch <dig>  interestingly though, our analysis shows mid <dig> to be a poor target.
mig <dig> is not a known shh target but is a significant result in our analysis and appears to be upregulated by the shh pathway.
both cellular retinoic acid binding proteins, crabp <dig> and crabp <dig>  score well in our analysis but are not known shh targets.
in some developmental processes, retinoic acid  has been shown to work synergistically with shh to create the zone of polarizing activity,  <cit>  for example, and in adult pluripotent stem cells  <cit> .
crabp binds ra with high affinity and is thought to spatially modulate the effect of ra by adjusting the concentration of ra reaching the nucleus  <cit> .
one interpretation of the mbn result is that that shh modulates the effects of ra by regulating crabp.
known shh targets not identified by the mbn
our results, however, do fail to identify some known and expected targets.
for example, one might expect the genes foxa <dig>  gli <dig>  gli <dig>  and ihh to be identified by the mbn as these are often listed as targets of the shh pathway, but in our analysis they did not score near the top.
the mbn template used in this analysis  detects relationships between a candidate target and shh and ptch <dig>  constrained by the information bottleneck introduced by the hidden protein nodes.
when we examine the relationship between the expected but poorly scoring genes and shh and ptch <dig> , we find no consistent pattern of differential expression.
it is clear that the expression does not change significantly or consistently with the pattern of change for shh or ptch <dig>  genes identified as targets by the mbn such as gas <dig>  gli <dig>  ptch <dig>  crabp <dig> and foxd <dig>  show clear relationships that can be assessed by eye.
comparison of mbn analysis and non-bayesian bioinformatics techniques
the results we obtain using mbn are different than those from standard techniques for determining downstream pathway targets.
whereas clustering, differential expression and significance testing are, in general, unable to identify nonlinear and multimodal relationships, bayesian methods such as mbn do not have this limitation.
bayesian methods learn based on complex patterns of change rather than the magnitude of change or a simple linear correlation.
when we use differential expression as an indicator of significance, for example, we find none of the canonical targets near the top of the list .
gli <dig>  the commonly used indicator of shh activity, ranks 629th in terms of fold-change, but 2nd with respect to mbn analysis.
this difference is expected, however, because only small changes in gene expression are required to modulate the activities of transcription factors and pathway components that participate in feedback loops with the shh ligand.
we obtain slightly better but still poor results when identifying targets based on a linear correlation to either shh or ptch <dig> gene expression.
this difficulty with clustering can be seen in table  <dig>  clustering methods suffer from similar failures as linear correlation and fail to take into account known information about the initial steps of the pathway.
comparison of mbn and bn modeling of the shh pathway
to isolate the effect of the mbn template over more standard bn approaches, we ran a direct comparison between the approaches using the three topologies shown in figure  <dig>  the simplest translation of the shh pathway, shown in figure  <dig>  is to interpret the qualitative diagram as a bn and replace all arrows with a bayesian edge to produce a sequential model.
although the topology looks most similar to a biochemical diagram, for a bayesian network the topology indicates that only smo influences the target gene.
an alternative topology is shown in figure  <dig>  in this parallel topology, the target gene is influenced by all of the genes in the pathway.
however, this topology does not capture the unobserved proteomic effects, nor does the topology distinguish between the relative position of shh, ptch, and smo in the signaling pathway
the results for the three topologies depicted in figures  <dig>   <dig> and  <dig> are shown in tables  <dig>   <dig> and  <dig> respectively.
when we compare the results, we find that the target genes predicted by the mbn topology ) best match what is known about the pathway.
if we compare the number of known targets identified in the top  <dig> targets predicted by each method we find the mbn predicts 5/ <dig>  the parallel topology predicts 4/ <dig>  and the sequential topology predicts 2/ <dig>  in the mbn results, gas <dig>  gli <dig> and ptch <dig> rank 1st, 2nd, and 3rd respectively and are all well-known shh targets, while in the sequential model these genes rank 11th, 2nd and 3rd respectively, and in the parallel model these genes rank 24th, 20th and 21st respectively.
one reason for this discrepancy is that the bn models impose fewer constraints on the candidate gene which would produce more false positives -- a result in line with what we see.
the full lists of target genes for each topology is provided in additional file  <dig>
the columns specify the gene name and short description; the bayesian score as described in methods; and a literature citation, if available, that specifies that the gene in indeed a true downstream target of the shh pathway.
the columns are identical to those in table  <dig>
extension to larger mbn templates
while we used a small template pathway in this work as a proof of concept, the method can be applied to larger pathways in a similar way.
however, large numbers of hidden nodes increase the computational requirements and may make the analysis intractable using the gibbs sampler used in this work.
the mbn template algorithm's runtime complexity is o where n is the number of candidate genes, m in the number of data samples and h is the number of hidden or latent variables.
because the limiting step is the gibbs sampling's o runtime complexity, we could use a more complex but efficient heuristic sampling method such as variational learning  <cit> .
as currently implemented, mbn templates should be limited to models with a small number of hidden nodes.
extension of mbns to other data types
although this study focused on gene expression data, the mbn template method can be used with any type of observed data that can be represented in the pathway.
if a molecular entity can be represented in an mbn template while being faithful to the underlying biochemical mechanisms, data regarding the molecule's concentration or activity can be used.
possible other measurements that could be used with a similar mbn approach include protein expression, kinase activity, and mirna expression.
furthermore, mbn templates can be used to integrate observations for multiple data types assuming that the measurements are made on a common set of samples.
CONCLUSIONS
we have shown how mbns can be used to integrate gene expression data and known topological information to uncover mechanistic details about complex pathways.
although we have shown only one example with finding downstream targets of the shh pathway, a similar approach could be used on other pathway topologies, data types, and to identify targets at different locations in the pathway.
in this example, the mbn method provided better target predictions than other methods, due in part because the topology assumed by the mbn is closer to the known biochemical mechanism.
list of abbreviations
mbn: mechanistic bayesian network; bn: bayesian network
authors' contributions
as designed the mbn software, carried out the analysis, and wrote much of the paper.
tt gathered the gene expression data and commented on the results.
apm designed the experiment, oversaw the data collection process, and commented on the results.
pjw designed the mbn approach, guided the analysis, and edited the paper.
all authors have read and approved the final manuscript.
supplementary material
additional file 1
source code and data.
contains all the code and data required to replicate the analysis.
the scripts use pebl's distributed computing features to run analysis over a variety of distributed computing platforms.
additionally, the scripts are easy to modify to run similar analysis with different datasets and mbn templates.
the readme.txt file contains further information about using and modifying the scripts.
click here for file
additional file 2
complete list of predicted shh targets.
excel file showing the complete target lists ranked by their fit within the mbn and bn frameworks .
the file also contains results for two smaller studies using the affymetrix 430a and 430b chips.
click here for file
acknowledgements
pjw and as are supported by nih grant u54-da- <dig>  work in apms laboratory was supported by a grant from the nih  <dig> r <dig> ns <dig>

 a general problem in systems biology is the integration of observational experimental data such as gene expression, with known pathways topologies. ideally these two data sources should be complementary, however in practice there are few methods to systematically integrate these two kinds of information. in this paper, we introduce a method called mechanistic bayesian networks in an attempt to use knowledge about the topology of a pathway with gene expression data related to the same pathway. as a sample case we used the topology of the sonic hedgehog signaling pathway as a model pathway along with a targeted gene expression dataset. using these data we used the mbn approach to identify regulatory targets of the shh pathway. sonic hedgehog pathway the shh pathway plays a central role in organismal development and the progression of some cancers. because of its central role, the shh pathway is well studied providing us with an ideal test case to validate our mbn approach. the details of shh are reviewed in detail elsewhere, but here we will summarize the early steps of the pathway that we will use in this work. shh is a secreted protein that acts as both a short-range contact-dependant factor and as a long-range diffusible morphogen. the shh ligand binds its canonical receptor patched, which releases its inhibition of a second membrane-bound protein, smoothened. derepression of smo in turn activates a signaling cascade inside the cell, eventually activating the gli transcription factors and regulating the expression of a variety of genes. while some of the downstream targets of the pathway are known, many remain unknown. the initial steps of the pathway are shown in figure with the terminal node genex representing a putative downstream target. note that in this cascade, the shh, ptch and smo proteins can all directly or indirectly affect the target, although in different ways and to different degrees. figure and show more abstracted representations of the pathway. due to its dose-dependant effect and role in development, the shh signal requires strict spatiotemporal control. the pathway's known targets include pathway components, promoters and antagonists, thus regulating the effects of the shh ligand over time. cell surface proteins that promote shh signaling, for example, are initially expressed in shh-responsive cells, sensitizing cells to even low levels of the shh ligand. as the level of shh signaling increases, downregulation of positive shh components such as gas and upregulation of negative components such as ptch ptch and hhip that sequester the ligand ensure a tight control over shh signaling. additionally, shh is known to cooperate with or antagonize other pathways such as bmp, retinoic acid, wnt, ras, and notch in a time-, dose-, and spatially-dependant manner. these factors further complicate the identification of pathway targets. the shh pathway raises a more general problem: given gene expression data from multiple samples, tissue and organs under different genetic knockout conditions, how can we identify downstream targets of a given pathway? due to the interactions between different pathways and multiple cascades between ligand reception and eventual transcriptional regulation, it is not clear what constitutes a downstream target of a specific pathway. common analysis approaches include significance testing between samples, differential expression, and clustering. while these approaches are widely used and helpful, they fail to incorporate knowledge of the underlying pathway. when the pathway information is used, it is done while ignoring fundamental biological knowledge such as the central dogma. dna, mrna, and proteins are conflated into one variable and specific interactions such as protein-protein or protein-dna are implicitly assumed to be detectable via gene expression data. to circumvent this problem, we present a novel mechanistic bayesian network approach that more closely respects the meaning of both the pathway and the expression data. methods in the following sections we describe the theoretical underpinnings of the mbn framework and provide a method for evaluating the significance of a result. next we describe how we tested the mbn method using gene expression data from gene knockout mouse models to identify targets of the shh pathway. mechanistic bayesian networks: theory and definition mechanistic bayesian networks represent a subset of a general class of analysis tools called bayesian networks. a bn is a probabilistic graphical model that encodes dependencies between variables in a compact and descriptive manner. a bn can be represented as a graph with nodes representing variables and edges representing dependencies or relationships between the variables. mathematically, each node describes a conditional probability distribution that quantitatively models the relationships between a node and its parent nodes. note that an edge between two variables indicates a directed relationship but does not specify the functional form of the relationship. said another way, an edge in a bayesian network does not indicate activation, inhibition or any other specific function. this interpretation of an edge differs from the usual definition of an edge in a signaling pathway, where edges are often assumed to have a defined effect. the advantage of the more broad edge definition is that the relationships between nodes in a bayesian network can be more complex and include functions that are nonlinear, multimodal, or logical. bns have been used successfully to model complex phenomenon in many fields and in systems biology, in particular, to model gene-regulatory networks, protein-interaction networks, signaling networks and to integrate heterogeneous biological data. methods for training and learning bns are well established and are available in a large number of software applications. in contrast to the more general bayesian network, a mechanistic bayesian network adheres more closely to known biological mechanisms by differentiating between mrna transcripts and proteins, and by including pathway-based structural constraints. because most studies do not measure both protein and gene expression, unknown quantities in an mbn are treated as unobserved, latent variables. by creating models that more closely resemble the known mechanism, mbns effectively incorporate in additional information that is not available in the experimental data alone. mechanistic bayesian network templates when analyzing a biological system, users often want to generate and rank a set of entities that match certain constraints, such as "all downstream targets of the shh pathway" or "all proteins that participate in the crosstalk between the shh and wnt pathways". while it is possible to devise specific sets of constraints and methods for each case, we propose a generic mbn based template approach. each set of constraints is expressed as an mbn template from which specific mbns are instantiated and evaluated against an experimental dataset. an mbn template is composed of two types of nodes: constant and template nodes. a constant node is a regular node corresponding to either one variable in the dataset or an unobserved variable while a template node is a holding place for variables in the dataset. each instantiation of the template is a regular mbn that can be scored using existing bn methods. maximum entropy discretization of the data before mbn can be applied to a dataset, the data must be discretized into a finite set of bins. while there are many ways to bin data, we use the maximum entropy principle to derive our discretization scheme. the principle states that the distribution that maximizes the information entropy is the true distribution given testable information. accordingly, we bin our data such that each bin contains the same number of data points as this maximizes the entropy of the distribution. there is no theory prescribing the optimal number of bins for any given dataset or analytical method; increasing the number of bin decreases the information loss incurred during discretization but also increases the number of parameters and thus decreases the statistical power of the analysis given the same data samples. most modelers working with systems biology data have arrived at three bins as a suitable compromise between information loss and statistical power. calculating the posterior probability of a model once the data are discretized, a particular model topology can be scored as the posterior probability of a model given data, given below. p is the posterior probability of the bn or mbn model m given the data d. said another way, p is our belief that the model m is correct after having observed the data d. p is the prior probability of the model, that is, our belief that the model m is correct before having observed the data d. the prior probability allows us to integrate model probabilities computed using other methods and data. p is the likelihood that the data was generated from the model m while p is a scaling factor that is usually ignored. p is calculated as the likelihood after marginalizing over all model parameters given a network structure. in the mbn framework, a bayesian network is modeled using a multinomial representation with dirichlet priors for the relationships between variables. this representation is convenient in that it is relatively agnostic to functional forms and has a convenient closed form solution, the bayesian dirichlet equivalent metric, described below. bde scoring metric: an approximation for the marginal likelihood the mbn approach uses the bde scoring metric to calculate the marginal likelihood of a dataset given a model. the bde metric is a closed form solution to the marginal likelihood for a multinomial bayesian network model with dirichlet priors. the bde metric is expressed as: where Γ is the gamma function, n is the number of nodes, ri is the arity of node i, qi is the arity of πi, nij is the number of samples where πi is in configuration j, nijk is the number of sampler where πi is in configuration j and the node has value k and αij and αijk are the prior counts corresponding to nij and nijk respectively. a full derivation of the bde metric is described elsewhere. in practice, the bde metric is represented in log-space for computational convenience. because of this log transformation, all bde scores in this paper will be negative. calculating the marginal likelihood with missing data if we have some missing values or latent variables in the dataset, we can use an altered method for computing the marginal likelihood. because the likelihood is no longer fully factorizable into the product of the probabilities for each variable, we must calculate the marginal likelihood given every potential completion of the missing data. the simplest method is to marginalize over all possible sets of values for the missing data and take the average. however, this exact enumeration approach is impractical for non-trivial cases and so a heuristic must be used. in this work we will used a modified gibbs sampling approach to approximate the marginal likelihood. we alter the method to only sample missing data completions that result in a maximum entropy discretization for all variables. this maximum entropy requirement ensures that p = p for all sets of variables and eliminates bias due to differential discretization for observed and hidden variables. calculating the marginal likelihood with interventions when calculating the marginal likelihood, data from interventions such as genetic knockouts are handled differently. if a particular data sample is the result of an intervention on a set of variables, the values for those variables no longer depend on their parent sets. they are, instead, arbitrarily set as the result of the intervention. accordingly, the specific variable samples are ignored when constructing the multinomial table for any node that was intervened upon. assessing model significance with bootstrapping based p-values to evaluate the significance of an mbn prediction, we calculate a p-value for each mbn result using nonparametric bootstrapping. in this approach, we generate a large number of mbn models with template nodes replaced by randomly generated variables with the same discretization schemes as the observed variables. the scores for these mbn models are used as the null model against which mbn scores are compared to determine p-values. due to the computational time required to score each bootstrap sample, we can only resolve p-values down to. and cannot correct for multiple hypothesis testing. we acknowledge that correction for multiple testing would lower the significance of each result. using mbn templates to define and identify downstream targets of the shh pathway to determine the downstream targets of the shh pathway, we created an mbn template based on the canonical pathway as described in literature with a terminal downstream template node as shown in figure specific mbn models are instantiated from the template by replacing the template node with a specific gene from the dataset. thus, each instantiated mbn model is a separate hypothesis that can be evaluated against the experimental data. experimental data to test our approach, we assembled data from three mouse models with gene knockouts in shh, ptch and smo described in more detail elsewhere. briefly, samples from different embryonic tissues at varying developmental stages were assayed using the u74v affymetrix microarrays to determine the expression profile. not all combinations of developmental stages and genetic backgrounds were available due to prenatal lethality for ptch and smo. samples assayed include 6-8-somite-stage with wildtype, smo-/-, and ptch1-/- backgrounds; 10- somite-stage with wildtype, smo-/-, and ptch1-/- backgrounds; days post fertilization embryo samples from head, trunk, and limb bud with wildtype and shh-/- backgrounds. the experimental design is summarized in table embryonic day embryonic day embryonic day embryonic day embryonic day experimental design for gathering gene expression data from a range of developmental stages, locations, and genetic backgrounds. columns indicate the genetic background of the mouse model, rows indicate the developmental stage and location of the samples and the numbers within the table cell indicate the number of samples assayed. data preprocessing the raw gene expression data were processed using rma as implemented in the bioconductor affy package. the data were annotated using updated probeset definitions provided by the microarray lab at the university of michigan. data were discretized into bins using the maximum entropy discretization scheme described above. because the data contain samples from varying developmental stages and genes with significant expression variation between stages, genes were discretized separately for each tissue type. this discretization scheme allowed the software to more easily identify relevant patterns within a tissue type. genes were knocked out by altering the sequence to create inactive protein. thus, for a knockout on a particular gene, the gene's expression was left unaltered but the protein expression was set to and the value was marked as an intervention so it could be handled differently by the scoring procedure. mbn template analysis to test the mbn template system using the shh dataset, we developed a template that matched the topology shown in figure in this template, we search for a downstream target node, genex, that responds to the shh signaling cascade. in this search, we assume that each gene has an equal prior probability of being the target node, and all relationships are scored using the bde score and a gibbs sampler, described above. using this template, we generated candidate mbn models -- one for each candidate target gene that showed significant differential expression. due to the custom cdf annotation, smo was not present on the chip, thus we modeled the smo mrna expression as a hidden node. no protein expression was observed in this experiment, so nodes representing protein concentrations or activities were also modeled as hidden. rather than selecting the number of iterations for the gibbs sampler a priori, we calculated the posterior distribution p with k iterations, then with k iterations. although the sampler had largely converged at this point as evidenced by the observation that the top results did not change significantly, we resumed the gibbs's sampler for the k run and sampled another k iterations for a total of k iterations per mbn model. the computation took approximately minutes per mbn for a total of approximately hours on a 48-node compute grid. p-values were calculated with bootstrapping describe above using mbn models to represent the null distribution. all code and data required to replicate the analysis are included in additional file. the code can be easily modified to run similar analysis on different dataset using different mbn templates. in addition, we also tested two alternative bayesian approaches shown in figure and these approaches use topologies that are simpler than the mbn approach in but less closely adhere to the known mechanism of the biochemistry of the pathway.  in the following sections we show the results from the mbn analysis of the shh knockout dataset and discuss the biological plausibility for each finding. next we discuss expected hedgehog pathway target genes that were not identified by the mbn and a comparison with other techniques. the predicted target genes from the mbn analysis are shown in table the top scoring hits include many known shh targets as described below. the top hit, gas has been shown to be a negative target of shh signaling whereas knockout and subsequent gain-of-function studies have shown gas to be a positive component of the shh pathway that acts synergistically with ptch to bind the hedgehog ligand. the second target, gli is one of the three gli transcription factors that regulate shh target genes and is often used as a canonical readout of shh activity. shh signaling is known to both induce gli expression and also regulate its nuclear accumulation and therefore its activity. ptch is a homolog of ptch but has a slightly different expression pattern. like ptch ptch is known to be transcriptionally modulated by shh signaling. knockout studies have shown that shh signaling is required for expression of msx a transcriptional repressor with a putative role in limb-pattern formation. foxc and foxd have been shown to be upregulated by shh signaling in null, conditional, and constitutively active smo mutant backgrounds. the columns specify the gene name and a short description; the bayesian score as described in methods; the resulting p-value for the bayesian score calculated using a bootstrap method as described in methods; a linear correlation to the shh transcript; a linear correlation to the ptch transcript; and a literature citation, if available, that specifies that the gene is indeed a true downstream target of the shh pathway. beyond these known shh-related genes, some of the highly ranked results appear to be novel putative shh targets. mig binds the opitz syndrome gene mid and the complex is thought to stabilize microtubules. mid is repressed by shh and this is partially borne out by our gene expression results which show a positive influence from shh and a negative influence from ptch interestingly though, our analysis shows mid to be a poor target. mig is not a known shh target but is a significant result in our analysis and appears to be upregulated by the shh pathway. both cellular retinoic acid binding proteins, crabp and crabp score well in our analysis but are not known shh targets. in some developmental processes, retinoic acid has been shown to work synergistically with shh to create the zone of polarizing activity, for example, and in adult pluripotent stem cells. crabp binds ra with high affinity and is thought to spatially modulate the effect of ra by adjusting the concentration of ra reaching the nucleus. one interpretation of the mbn result is that that shh modulates the effects of ra by regulating crabp. known shh targets not identified by the mbn our results, however, do fail to identify some known and expected targets. for example, one might expect the genes foxa gli gli and ihh to be identified by the mbn as these are often listed as targets of the shh pathway, but in our analysis they did not score near the top. the mbn template used in this analysis detects relationships between a candidate target and shh and ptch constrained by the information bottleneck introduced by the hidden protein nodes. when we examine the relationship between the expected but poorly scoring genes and shh and ptch, we find no consistent pattern of differential expression. it is clear that the expression does not change significantly or consistently with the pattern of change for shh or ptch genes identified as targets by the mbn such as gas gli ptch crabp and foxd show clear relationships that can be assessed by eye. comparison of mbn analysis and non-bayesian bioinformatics techniques the results we obtain using mbn are different than those from standard techniques for determining downstream pathway targets. whereas clustering, differential expression and significance testing are, in general, unable to identify nonlinear and multimodal relationships, bayesian methods such as mbn do not have this limitation. bayesian methods learn based on complex patterns of change rather than the magnitude of change or a simple linear correlation. when we use differential expression as an indicator of significance, for example, we find none of the canonical targets near the top of the list. gli the commonly used indicator of shh activity, ranks 629th in terms of fold-change, but 2nd with respect to mbn analysis. this difference is expected, however, because only small changes in gene expression are required to modulate the activities of transcription factors and pathway components that participate in feedback loops with the shh ligand. we obtain slightly better but still poor results when identifying targets based on a linear correlation to either shh or ptch gene expression. this difficulty with clustering can be seen in table clustering methods suffer from similar failures as linear correlation and fail to take into account known information about the initial steps of the pathway. comparison of mbn and bn modeling of the shh pathway to isolate the effect of the mbn template over more standard bn approaches, we ran a direct comparison between the approaches using the three topologies shown in figure the simplest translation of the shh pathway, shown in figure is to interpret the qualitative diagram as a bn and replace all arrows with a bayesian edge to produce a sequential model. although the topology looks most similar to a biochemical diagram, for a bayesian network the topology indicates that only smo influences the target gene. an alternative topology is shown in figure in this parallel topology, the target gene is influenced by all of the genes in the pathway. however, this topology does not capture the unobserved proteomic effects, nor does the topology distinguish between the relative position of shh, ptch, and smo in the signaling pathway the results for the three topologies depicted in figures and are shown in tables and respectively. when we compare the results, we find that the target genes predicted by the mbn topology ) best match what is known about the pathway. if we compare the number of known targets identified in the top targets predicted by each method we find the mbn predicts 5/ the parallel topology predicts 4/ and the sequential topology predicts 2/ in the mbn results, gas gli and ptch rank 1st, 2nd, and 3rd respectively and are all well-known shh targets, while in the sequential model these genes rank 11th, 2nd and 3rd respectively, and in the parallel model these genes rank 24th, 20th and 21st respectively. one reason for this discrepancy is that the bn models impose fewer constraints on the candidate gene which would produce more false positives -- a result in line with what we see. the full lists of target genes for each topology is provided in additional file the columns specify the gene name and short description; the bayesian score as described in methods; and a literature citation, if available, that specifies that the gene in indeed a true downstream target of the shh pathway. the columns are identical to those in table extension to larger mbn templates while we used a small template pathway in this work as a proof of concept, the method can be applied to larger pathways in a similar way. however, large numbers of hidden nodes increase the computational requirements and may make the analysis intractable using the gibbs sampler used in this work. the mbn template algorithm's runtime complexity is o where n is the number of candidate genes, m in the number of data samples and h is the number of hidden or latent variables. because the limiting step is the gibbs sampling's o runtime complexity, we could use a more complex but efficient heuristic sampling method such as variational learning. as currently implemented, mbn templates should be limited to models with a small number of hidden nodes. extension of mbns to other data types although this study focused on gene expression data, the mbn template method can be used with any type of observed data that can be represented in the pathway. if a molecular entity can be represented in an mbn template while being faithful to the underlying biochemical mechanisms, data regarding the molecule's concentration or activity can be used. possible other measurements that could be used with a similar mbn approach include protein expression, kinase activity, and mirna expression. furthermore, mbn templates can be used to integrate observations for multiple data types assuming that the measurements are made on a common set of samples.  we have shown how mbns can be used to integrate gene expression data and known topological information to uncover mechanistic details about complex pathways. although we have shown only one example with finding downstream targets of the shh pathway, a similar approach could be used on other pathway topologies, data types, and to identify targets at different locations in the pathway. in this example, the mbn method provided better target predictions than other methods, due in part because the topology assumed by the mbn is closer to the known biochemical mechanism. list of abbreviations mbn: mechanistic bayesian network; bn: bayesian network authors' contributions as designed the mbn software, carried out the analysis, and wrote much of the paper. tt gathered the gene expression data and commented on the results. apm designed the experiment, oversaw the data collection process, and commented on the results. pjw designed the mbn approach, guided the analysis, and edited the paper. all authors have read and approved the final manuscript. supplementary material additional file 1 source code and data. contains all the code and data required to replicate the analysis. the scripts use pebl's distributed computing features to run analysis over a variety of distributed computing platforms. additionally, the scripts are easy to modify to run similar analysis with different datasets and mbn templates. the readme.txt file contains further information about using and modifying the scripts. click here for file additional file 2 complete list of predicted shh targets. excel file showing the complete target lists ranked by their fit within the mbn and bn frameworks. the file also contains results for two smaller studies using the affymetrix 430a and 430b chips. click here for file acknowledgements pjw and as are supported by nih grant u54-da- work in apms laboratory was supported by a grant from the nih r ns 
the topology of a biological pathway provides clues as to how a pathway operates, but rationally using this topology information with observed gene expression data remains a challenge. we introduce a new general-purpose analytic method called mechanistic bayesian networks that allows for the integration of gene expression data and known constraints within a signal or regulatory pathway to predict new downstream pathway targets. we demonstrate how mbns can be used by modeling the early steps of the sonic hedgehog pathway using gene expression data from different developmental stages and genetic backgrounds in mouse. the mbn framework is implemented in an open-source bayesian network learning package, the python environment for bayesian learning. the mbn approach shown here can easily be extended to other pathways and data types to yield a more mechanistic framework for learning genetic regulatory models. using the mbn approach we are able to automatically identify many of the known downstream targets of the hedgehog pathway such as gas and gli along with a short list of likely targets such as mig 
 the topology of a biological pathway provides clues as to how a pathway operates, but rationally using this topology information with observed gene expression data remains a challenge.  we introduce a new general-purpose analytic method called mechanistic bayesian networks that allows for the integration of gene expression data and known constraints within a signal or regulatory pathway to predict new downstream pathway targets. the mbn framework is implemented in an open-source bayesian network learning package, the python environment for bayesian learning. we demonstrate how mbns can be used by modeling the early steps of the sonic hedgehog pathway using gene expression data from different developmental stages and genetic backgrounds in mouse. using the mbn approach we are able to automatically identify many of the known downstream targets of the hedgehog pathway such as gas and gli along with a short list of likely targets such as mig  the mbn approach shown here can easily be extended to other pathways and data types to yield a more mechanistic framework for learning genetic regulatory models. 